//  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 and Moon directions in ECI frame. 
//> - Sun: at the specified time (round mark) + trajectory over 350 days. 
//> - Moon: at the specified time (round mark) + trajectory over 26 days. 
//> - Earth meridians: at the specified time (every 15 degrees) so that
//>   the longitude of the Sun and Moon can be read on the graph. 'G' is the 
//>   Greenwich meridian (longitude = 0 deg).  
//
// Author: A. Lamy
// -----------------------------------------------------------


// -----------------------------------------------------------
// Computation
// -----------------------------------------------------------
function print_info(t)
  s = CL_eph_sun(t);  // Frame = ECI
  s_sph = CL_co_car2sph(s);
  m = CL_eph_moon(t);   // Frame = ECI
  m_sph = CL_co_car2sph(m);
  tsid = CL_rMod(CL_mod_siderealTime(t),0,2*%pi);

  printf("Time (TREF): %s\n", CL_dat_cal2str(CL_dat_cjd2cal(t)));
  printf("Coordinates relative to ECI (=%s) frame:\n", CL_configGet("ECI_FRAME"));
  printf("Sun  ra/decl (deg):  %7.2f  %7.2f \n", s_sph(1)*%CL_rad2deg, s_sph(2)*%CL_rad2deg);
  printf("Moon ra/decl (deg):  %7.2f  %7.2f \n", m_sph(1)*%CL_rad2deg, m_sph(2)*%CL_rad2deg);
  printf("Sidereal time (deg): %7.2f\n", tsid*%CL_rad2deg);
  printf("\n");
endfunction


// Calendar date - TREF 
scal = []; 

desc_param = list(..
   CL_defParam("Date (calendar format, TREF)", typ= "cal", scal) ..
);

[scal] = CL_inputParam(desc_param);

t0 = CL_dat_cal2cjd(CL_dat_str2cal(scal)); 
d1 = 350; // days for Sun propagation
d2 = 26; // days for Moon propagation
t1 = linspace(t0, t0+d1, round(d1/2)+2); // considered as UT
t2 = linspace(t0, t0+d2, round(d2*5)+2); // considered as UT

print_info(t0);


// -----------------------------------------------------------
// plot
// -----------------------------------------------------------
// name of body (capitalized)
function [name] = body_name(body)
  if (part(body,1:1) == 'm')
    name = "Moon"; 
  else
    name = "Sun"; 
  end
endfunction

// xtring interface
function String(x,y,text,col,sz,style)
  xstring(x,y,text);
  e=gce();
  e.font_foreground = col;
  e.font_size = sz;
  e.font_style = style;
endfunction

// thickness after plot
function Thickness(th)
  e=gce();
  e.children.thickness = th;
endfunction

// marks for sun and moon 
function plot_mark(t,body,col,fg,sz)
  if (body == "sun")
    d = CL_eph_sun(t);  
  elseif (body == "moon")
    d = CL_eph_moon(t); 
  end
  d_sph = CL_co_car2sph(d);
  plot(d_sph(1)*%CL_rad2deg, d_sph(2)*%CL_rad2deg, "o");
  e=gce();
  e.children.mark_background=col;
  e.children.mark_foreground=fg;
  e.children.mark_size=sz;
  name = body_name(body);
  String(d_sph(1)*%CL_rad2deg, d_sph(2)*%CL_rad2deg, " "+name, col,3,5);
endfunction


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

col1 = color("orange"); // Sun
col2 = color("steelblue"); // Moon
col3 = color("darkseagreen"); // meridians
col4 = color("forestgreen"); // meridians

db = [0,-30;360,30]; // data bounds

// ---------------------
// Greenwich meridian
// ---------------------

tsid = CL_rMod(CL_mod_siderealTime(t0),0,2*%pi);

for dl=0:15:359
  lon_deg = CL_rMod(dl + tsid*%CL_rad2deg, 0, 360); 
  plot2d([lon_deg,lon_deg], [-30, 30], style=col3);
  if modulo(dl,60) == 0
     Thickness(2);
     String(lon_deg, -29.9, string(dl), col4, 1, 6);
  end
end

plot2d([tsid, tsid]*%CL_rad2deg, [-30, 30], style=col4);
Thickness(2);
String(tsid*%CL_rad2deg, 0, ' G', col4, 3, 4);

// ---------------------
// Sun
// ---------------------

d = CL_eph_sun(t1);  
CL_plot_ephem(d, win_id=f, color_id=col1, data_bounds=db, thickness=2);
plot_mark(t0,"sun",col1,1,10)

// ---------------------
// Moon
// ---------------------

d = CL_eph_moon(t2);  
CL_plot_ephem(d, win_id=f, color_id=col2, data_bounds=db, thickness=2);
plot_mark(t0,"moon",col2,1,10)


// various settings

a.x_ticks = tlist("ticks", 0:60:360, string(0:60:360));
a.y_ticks = tlist("ticks", -30:10:30, string(-30:10:30));
a.x_label.text = "Right ascension (deg)";
a.y_label.text = "Declination (deg)";
months = ["Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"];  
cal = CL_dat_cjd2cal(t0); 

str_day = sprintf("%d %s %d", cal(3), months(cal(2)), cal(1));
str_hour = sprintf("%02d:%02d", cal(4), cal(5)); 
 
ECI_frame = CL_configGet("ECI_FRAME");
a.title.text = "Sun and Moon in ECI (=" + ECI_frame + ") - " + str_day + " " + str_hour + " (TREF)"; 

a.data_bounds = db;
CL_g_stdaxes(a);
a.grid = [-1,0];


f.immediate_drawing="on";
f.visible="on";