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

// -----------------------------------------------------------
//> Solar system planet positions:
//> Positions of solar system planets in EOD (ecliptic of date), 
//> using JPL's DE405 ephemerides.
//> Valid dates are between years 1960 and 2200. 
//> Round symbols represent the body positions at the 
//> reference date.  
//
// Author: G. Azema / A. Lamy
// -----------------------------------------------------------

// ------------------------------
// Utility functions
// ------------------------------

// ---------------------------------------------------------
// Check if dates are in the correct range
// NB: cal is invalid if it contains %nan (then => %f)
// scal: calendar string, dur = [d_past, d_future] (days)
// ---------------------------------------------------------
function [ok] = isValidDate(scal, dur)  
  ok = %f; 
  cjd1 = CL_dat_cal2cjd(1960,1,1);
  cjd2 = CL_dat_cal2cjd(2201,1,1);
  cal = CL_dat_str2cal(scal); 
  if (cal(1) <> %nan)
    cjd = CL_dat_cal2cjd(cal); 
    ok = ( dur(1)+dur(2) > 0 & ..
           cjd - dur(1) >= cjd1 & cjd + dur(2) < cjd2 );
  end
endfunction

// ---------------------------------------------------------
// get body position from the Sun in EOD 
// cjd: TREF time scale
// ---------------------------------------------------------
function [pos] = getBodyPos(name, cjd, ephem)
  pos_ICRS = CL_eph_de405(name, cjd, "Sun", ephem=ephem);
  pos = CL_fr_convert("ICRS", "EOD", cjd, pos_ICRS);

endfunction

// ---------------------------------------------------------
// plots body trajectory
// col = color_index
// pos: in drawing units
// ---------------------------------------------------------
function plotBodyPos(pos, col)
  plot(pos(1,:)/UNIT, pos(2,:)/UNIT); 
  h = CL_g_select(gce(), "Polyline"); 
  h.foreground = col;
  h.thickness = 2;
endfunction

// ---------------------------------------------------------
// Plots body 
// col = color_index, sz = size (points)
// pos: in drawing units
// ---------------------------------------------------------
function plotBody(pos, col)
  plot(pos(1,:)/UNIT, pos(2,:)/UNIT);  
  h = CL_g_select(gce(), "Polyline");
  h.line_mode = "off"; 
  h.mark_mode = "on"; 
  h.mark_style = 0; // round mark
  h.mark_background = col;
  h.mark_foreground = col;
  h.mark_size_unit = "point"; 
  h.mark_size = 7;
endfunction


// ---------------------------------------------------------
// print position - print(pos), print("start"), print("end")
// ---------------------------------------------------------
function printPos(name, pos)
  if (name == "start")
    mprintf("-----------------------\n"); 
    mprintf("%s\n", "Positions (10^6 km)"); 
    mprintf("Date = %s (TREF)\n", CL_dat_cal2str(cal_ref)); 
    mprintf("Frame = %s\n", "EOD"); 
    mprintf("Center = %s\n", "Sun"); 
    mprintf("-----------------------\n"); 
  elseif (name == "end")
    mprintf("\n"); 
  else
    mprintf("%-10s", name); 
    for (i=1:3); mprintf("%15.6f ", pos(i,1)/1.e9); end
    printf("\n"); 
  end
endfunction


// ===========================================================
// MAIN 
// ===========================================================
// Name, color, size (points) for each body
names = [ "Sun", "Mercury", "Venus", "Earth", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune"];
colors = ["yellow2", "grey", "orange", "blue", "red", "brown", "green", "slateblue", "darkturquoise"]; 

// Initial start date = current time (TREF)
scal = CL_dat_cal2str(CL_dat_cjd2cal(floor(CL_dat_now()))); 
dur = [60,0]; // days before/after 
bodyset = 1; 

desc = list(..
  CL_defParam("Date (calendar format, TREF time scale)", scal, typ='cal', valid="isValidDate($x, $dur)"), ...
  CL_defParam("Arc length into the past/future", dur, dim=2, units=["day"], id="$dur", valid="$x >= 0"), ...
  CL_defParam("Planets: 1=inner pl.+Mars, 2=outer pl.", bodyset, accv=[1,2]) ...
);

[scal, dur, bodyset] = CL_inputParam(desc);

// converts to useable data. 
cal_ref = CL_dat_str2cal(scal);
cjd_ref = CL_dat_cal2cjd(cal_ref);
cjd_start = cjd_ref - dur(1);
cjd_end = cjd_ref + dur(2);

// select bodies
if (bodyset == 1); I = 1:5; else I = [1, (5:9)]; end
names = names(I);
colors = colors(I);


// -----------------------------------------------------------
// Results / plot
// -----------------------------------------------------------
f = scf();
f.visible = "on";
f.immediate_drawing = "off";

nb = size(names,"*");
UNIT = 1.e9; // millions of km 

// header 
printPos("start"); 

// legend handles
hleg = []; 

// load DE405 ephemeris file 
cjd_all = [cjd_start:cjd_end, cjd_ref]; 
ephem = CL_eph_de405LoadT(cjd_all); 

for k = 1 : nb
  // get all needed positions
  pos = getBodyPos(names(k), cjd_all, ephem);

  // Plot trajectories
  plotBodyPos(pos(:, 1:$-1), color(colors(k)));
  hleg($+1) = gce(); 
  
  // Plot body positions at ref date
  plotBody(pos(:,$), color(colors(k)));
  printPos(names(k), pos(:,$)); 
end

printPos("end"); 


// Plot adjustments

// legends
CL_g_legend(hleg, names, new_axes=%t);

a = gca();

str = tokens(CL_dat_cal2str(cal_ref)); 
a.title.text = "Positions of solar system bodies - " + str(1) + " - " + "EOD" + " (2D)";
a.x_label.text = "X (10^6 km)";
a.y_label.text = "Y (10^6 km)";

CL_g_stdaxes(a); 
a.grid_position = "background"; 
L = max(a.data_bounds(2,:) - a.data_bounds(1,:)) * 1.08; 
mid = (a.data_bounds(1,:) +a.data_bounds(2,:))/2;
a.data_bounds = [ mid(1) + [-L/2; L/2], mid(2) + [-L/2; L/2] ]; 

a.isoview = "on"; 
a.tight_limits = "on"; 
f.immediate_drawing="on";
f.visible="on";