//  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'.

// -----------------------------------------------------------
//> The moon is shown each day as it appears as seen from some
//> location on Earth. 
//> The time each day is the time (between reference time and 
//> reference time + 1 day) when the elevation of the moon is 
//> at its maximum.  
//> Cases where the moon is partly visible (partly below the
//> horizon, eclipsed by the Sun...) are not considered.  
//> The background color indicates whether it's daylight or not. 
//>
//> Note: The time scale is TREF (see help for details). 

// Author: A. Lamy
// -----------------------------------------------------------

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

// elevation of body from location loc at times t
function [elev] = Elevation(t, loc, body)
  // Compute position of body in ECI frame
  if (body == "sun")
    pos_body = CL_eph_sun(t);  
  elseif (body == "moon")
    pos_body = CL_eph_moon(t); 
  end  
  pos_ECF = CL_fr_convert("ECI", "ECF", t,pos_body);
  loc_car = CL_co_ell2car(loc) * ones(1,size(pos_ECF,2));
  pos_topoN = CL_fr_topoNMat(loc) * (pos_ECF-loc_car);
  pos_topoN_sph = CL_co_car2sph(pos_topoN);  
  elev = pos_topoN_sph(2,:); 
endfunction

// maximum elevation of body between t0 and t0+1
function [tmax] = Date_max_elev(t0, loc, body)
  t = linspace(t0, t0+1-1.e-5, 1440); 
  [elevmav, k] = max(Elevation(t, loc, body)); 
  tmax = t(k);
endfunction

// defines the color for each facet, considering dot product 
// of each direction the the light (cosang)
function [col] = Color(cosang, icol)
  c1 = -0.03; 
  c2 = 0.03; 
  c3 = 1; 
  col = zeros(cosang);
  I = find(cosang <= c1);
  col(I) = icol(1); 
  I = find(cosang > c1 & cosang < c2);
  col(I) = icol(1) + round(((cosang(I)-c1)/(c2-c1)) * (icol(2)-icol(1))); 
  I = find(cosang >= c2);
  col(I) = icol(2) + round(((cosang(I)-c2)/(c3-c2)) * (icol(3)-icol(2))); 
endfunction

// draw the Moon: 
// P: transformation matrix defining the new frame
// icol: array of color indices for color interpolation
function Draw_moon(P, icol)

  u = linspace(-%pi/2, %pi/2, 100);
  v = linspace(0, 2*%pi, 100);
  x = cos(u)'*cos(v);
  y = cos(u)'*sin(v);
  z = sin(u)'*ones(v);
  
  plot3d2(x, y, z); 
  e = gce();

  d = e.data;

  // The light comes from +z
  // d.z is the scalar product of the vector with the direction 
  // of the light => color
  col = Color(d.z, icol); 

  // change frame of the points defining the sphere
  dim = size(d.x); 
  coord = P * [d.x(:)'; d.y(:)'; d.z(:)']; 
  x = matrix(coord(1,:), dim);
  y = matrix(coord(2,:), dim);
  z = matrix(coord(3,:), dim);

  tl = tlist(["3d" "x" "y" "z" "color"],x,y,z,col);
  e.data = tl;
  e.color_flag = 3; // interpolated
  e.color_mode = -1;

endfunction

// P:transf matrix from frame with sun towards z to local frame
// local frame: 
// X: horizontal (towards the right when looking at the moon
// Y: perpendicular to the Moon direction (towards the North 
// in the northern hemisphere, and South in the southern)
function [P] = Obs_matrix(t, loc)

// Sun in ECI (t: TREF)
pos_sun = CL_eph_sun(t);

// Moon in ECI (t: TREF)
pos_moon = CL_eph_moon(t);

// Observer in ECI
pos_ECF = CL_co_ell2car(loc);
vertical_ECF = CL_co_sph2car([loc(1);loc(2);1]);
pos_ECI = CL_fr_convert("TIRS","ECI",t,pos_ECF);
vertical_ECI = CL_fr_convert("TIRS","ECI",t, vertical_ECF);

// observer local frame
// NB: pos ~ local vertical
Mobs = CL_rot_defFrameVec(pos_ECI-pos_moon,vertical_ECI,3,2);

// frame : Z -> sun
M = CL_rot_defFrameVec(pos_sun-pos_moon,pos_sun-pos_moon,3,3);

P = Mobs*M'; 

endfunction 


// -----------------------------------------
// main
// -----------------------------------------

// longitude/latitude of location
lon = 0 * %CL_deg2rad;       
lat = 42 * %CL_deg2rad; 

// scal = reference date (string, calendar)
scal = [];  
ndays = 0:0;  // numbers of relative days / t0

desc_param = list(..
   CL_defParam("Longitude of location", lon, units=["rad", "deg"], valid='$x>=0 & $x <=360'), ..
   CL_defParam("Latitude of location", lat, units=["rad", "deg"], valid='$x>=-90 & $x <=90'), ..
   CL_defParam("Reference time (calendar format, TREF)", scal, typ="cal"), ..
   CL_defParam("Day numbers relative to reference (empty => ref time)", ndays, dim=-1, valid='$x>=-100 & $x<=100 & round($x)==$x')..
);

[lon,lat,scal,ndays] = CL_inputParam(desc_param);

loc = [lon;lat;0]; 
t0 = CL_dat_cal2cjd(CL_dat_str2cal(scal)); 

// visualization instants
if (ndays == [])
  tab_t = t0; 
else
  tab_t = []; 
  tprec = -%inf; 
  for idays = ndays
    t = Date_max_elev(t0 + idays, loc, "moon"); 
    if (abs(t-tprec) < 0.5); continue; end
    tab_t = [tab_t, t];
    tprec = t; 
  end
end


// -----------------------------------------
// plot
// -----------------------------------------

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

c0 = 0.3*[1,1,1]; // grey
c1 = [1,0.85,0.1]; // dark yellow
c2 = [1,1,0.5]; // light yellow

p = (0:30)'/30;
Cmap1 = CL_dMult(1-p, c0) + CL_dMult(p, c1);
Cmap2 = CL_dMult(1-p, c1) + CL_dMult(p, c2);

f.color_map = [[1,1,1]; Cmap1; Cmap2]; 
icol = 1 + [1, 31, 62];  

months = ["Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"]; 


// loop
for t = tab_t

tic();

f.immediate_drawing = "off"; 
clf(f);

if (Elevation(t, loc, "moon") >= 0); 
   P = Obs_matrix(t, loc); 
   Draw_moon(P, icol);
end

// graphic settings
a = gca();

a.axes_visible = ["off", "off", "off"];
a.box = "off"; 
a.x_label.text = ""; 
a.y_label.text = ""; 
a.z_label.text = ""; 
a.isoview = "on";

// rotation angles: very important! 
// screen = (X,Y) frame
// X axis: positive towards the right
// Y axis: positive towards the top 
a.rotation_angles = [0, -90]; 

cal = CL_dat_cjd2cal(t); 
str_day = sprintf("%d %s %d", cal(3), months(cal(2)), cal(1));
str_hour = sprintf("%02d:%02d", cal(4), cal(5)); 
 
a.title.text = str_day + " " + str_hour + " (TREF)"; 
a.title.font_size = 4;
a.title.font_style = 4;

col_bg = "darkblue"; 
col_fg = "lightblue"; 
if (Elevation(t, loc, "sun") >= 0); 
  col_bg = "lightblue"; 
  col_fg = "darkblue"; 
end

f.background = color(col_bg); 
a.background = color(col_bg); 
a.title.font_foreground = color(col_fg);

f.immediate_drawing = "on"; 

dt = (0.8-toc()) * 1000;
sleep(max(dt,1)); // sleeps at most 0.8 s

end