//  Copyright (c) CNES  2008
//
//  This software is part of CelestLab, a CNES toolbox for Scilab
//
//  This software is governed by the CeCILL  license under French law and
//  abiding by the rules of distribution of free software.  You can  use,
//  modify and/ or redistribute the software under the terms of the CeCILL
//  license as circulated by CEA, CNRS and INRIA at the following URL
//  'http://www.cecill.info'.

// -----------------------------------------------------------
//> Sun reflection point (glint): 
//> The orbit is circular and the Sun is fixed. 
//>
//> Graph: 
//> - Sun zenith angle: angle between zenith and Sun direction 
//>   (from reflection point)
//> - Tilt to vertical: angle between nadir and reflection point
//>   (from satellite position) 
//> - Dev. from GT: angle between reflection point and 
//>   satellite trajectory (from Earth centre)
//> 3D view: 
//> - Yellow lines: Sun rays (truncated arbitrarily)
//> - Pink lines: Reflection point's vertical (truncated arbitrarily)
//> - Green lines: From satellite to reflection point 
//> - Red curve: Trajectory of the reflection point 
//> - Black curve: Trajectory of the satellite 
//
// Author: A. Lamy
// -----------------------------------------------------------

betaa = 50 * %CL_deg2rad; // beta angle
alt = 1000.e3; // altitude  
showopt = [1,1,1]; // 3D options

desc_param = list(..
   CL_defParam("Satellite altitude", alt, units=['m', 'km'], valid='$x>0'),..
   CL_defParam("Sun beta angle", betaa, units=['rad', 'deg'], valid='$x>=0 & $x<=90'),..
   CL_defParam("Show lines from glint to: vertical?, satellite?, Sun? (1=yes, 0=no)", showopt, dim=3, accv =[0,1])..
);

[alt, betaa, showopt] = CL_inputParam(desc_param);

// -----------------------------------------------
// Computation
// -----------------------------------------------

// satellite trajectory
rsat = alt + %CL_eqRad; 
alpha = linspace(-%pi, %pi, 91); // every 4 degrees
pos_sat = rsat * [cos(alpha); sin(alpha); zeros(alpha)]; 

// Sun direction
dir_sun = [cos(betaa); 0; sin(betaa)] * ones(alpha); 

// reflection point 
[pos_glint, incid] = CL_gm_reflectionPtSph(pos_sat,dir_sun,%CL_eqRad,pos2_inf=%t);

// spherical coordinates
pos_glint_sph = CL_co_car2sph(pos_glint); 

// tilt angle to vertical
tilt = CL_gm_visiParams(rsat, %CL_eqRad, 'incid', incid, ['sat']);

I = find(incid >= %pi/2); 
pos_satg = pos_sat; 
pos_glintg = pos_glint; 
pos_glintg(:,I) = %nan; 
pos_satg(:,I) = %nan; 


// -----------------------------------------------
// internal functions
// -----------------------------------------------

// segments: A -> B
function trace_segs(A,B,col)
  N = max(size(A,2), size(B,2));
  if (size(A,2) < N); A = A * ones(1:N); end  
  if (size(B,2) < N); B = B * ones(1:N); end  
  xsegs([A(1,:); B(1,:)], [A(2,:); B(2,:)], [A(3,:); B(3,:)]); 
  e=gce();
  e.segs_color = col * ones(1:N); 
endfunction

// trajectory
function trace_traj(traj,F,col,th)
  param3d(F*traj(1,:), F*traj(2,:), F*traj(3,:)); 
  e=gce();
  e.foreground=col;
  e.thickness=th;
endfunction

// sphere
function trace_sphere(R)
  F = 0.999;
  u = linspace(-%pi/2,%pi/2,37); // every 5 deg
  v = linspace(0,2*%pi,73); // every 5 deg
  x = F*R*cos(u)'*cos(v);
  y = F*R*cos(u)'*sin(v);
  z = F*R*sin(u)'*ones(v);
  plot3d2(x,y,z);
  e=gce();
  e.color_flag = 0; 
  e.color_mode = 12;
  e.foreground = 18;

  // origin meridian / equator
  trace_traj(R*[cos(u);zeros(u);sin(u)], F=1, col=16, th=1);  
  trace_traj(R*[cos(v);sin(v);zeros(v)], F=1, col=16, th=1);  

endfunction



// -----------------------------------------------
// PLOT 1
// -----------------------------------------------

f1=scf();
f1.visible = "on"; 
f1.immediate_drawing = "off"; 

// Earth center -> glint 
if (showopt(1) == 1) 
  trace_segs([0;0;0], 1.5*pos_glintg, col=31); 
  trace_traj(1.5*pos_glintg, F=1, col=31, th=1); 
end

// sat -> glint 
if (showopt(2) == 1) 
  trace_segs(pos_satg, pos_glintg, col=13);
end

// glint -> sun 
if (showopt(3) == 1) 
  trace_segs(pos_glintg, pos_glintg+%CL_eqRad*dir_sun, col=32);
  trace_traj(pos_glintg+%CL_eqRad*dir_sun, F=1, col=32, th=1); 
end
 
// sat trajectory
trace_traj(pos_sat, F=1, col=1, th=3); 

// glint trajectory 
trace_traj(pos_glintg, F=1.001, col=20, th=3); 

// Earth
trace_sphere(%CL_eqRad);


a=gca(); 
CL_g_stdaxes(a);
a.isoview = "on"; 
a.x_label.text = "";
a.y_label.text = "";
a.z_label.text = "";
a.rotation_angles = [95, 75]; 
a.title.text = "Sun glint - 3D view"; 
a.grid = [-1,-1,-1];
a.axes_visible = ["off", "off", "off"]; 
a.box = "off"; 

f1.immediate_drawing = "on"; 


// -----------------------------------------------
// PLOT 2
// -----------------------------------------------
f2=scf();
f2.figure_position = f1.figure_position + [200,200]; 
f2.visible = "on"; 
f2.immediate_drawing = "off"; 

a=gca();
I = find(incid <= %pi/2); 

plot2d(alpha(I)*%CL_rad2deg, incid(I) * %CL_rad2deg, style=20); 
plot2d(alpha(I)*%CL_rad2deg, tilt(I) * %CL_rad2deg, style=13); 
plot2d(alpha(I)*%CL_rad2deg, pos_glint_sph(2,I) * %CL_rad2deg, style=1); 

h = CL_g_select(a, "Polyline"); 
h.thickness = 2; 

xtitle("Sun glint - characteristic angles (deg)", "Satellite position / local orbital noon (deg)"); 

CL_g_legend(a, ["Sun zenith angle", "Tilt to vertical", "Dev. from GT"]); 

CL_g_stdaxes();
a.data_bounds = [-180,0;180,90];
a.tight_limits = "on"; 

f2.immediate_drawing = "on"; 


// makes everything visible
show_window(f1);
f1.visible = "on"; 
f2.visible = "on";