//  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, ICRS or ITRF, using
//> JPL's DE405 ephemerides.
//> Valid dates are between years 1960 and 2200. 
//> Bodies are referred to by names (beginning of name with
//> first letter capitalized so that it is not ambiguous). 
//> The availables bodies are: 
//> Sun, Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus,
//> Neptune, Pluto, Moon.
//> The center of the frame can be any body. 
//> Round symbols represent the body positions at the 
//> reference date.  
//
// Authors: 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

// ---------------------------------------------------------
// find (full) body name from beginning of name 
// => [] if does not exist
// NB: Use global variables "Body_names"
// ---------------------------------------------------------
function [name] = getBodyName(str)
  name = []; 
  str = stripblanks(str, %t);
  k = find(str == part(Body_names, 1:length(str))); 
  if (length(k) == 1); 
    name = Body_names(k); 
  end
endfunction

// ---------------------------------------------------------
// find (full) body names from beginning of names
// => [] if does not exist
// ---------------------------------------------------------
function [names] = getBodyNames(str)
  names = []; 
  s = stripblanks(tokens(str), %t); 
  for (i = 1 : size(s, "*"))
    name = getBodyName(s(i)); 
    if (name == []); names = []; return; end
    names = [names, name]; 
  end
endfunction

// ---------------------------------------------------------
// indices from (full) body names 
// ---------------------------------------------------------
function [ind] = getBodyIndices(names)
  ind = zeros(names); 
  for (i = 1 : size(names, "*"))
    ind(i) = find(Body_names == getBodyName(names(i))); 
  end
endfunction

// ---------------------------------------------------------
// body name => color index
// NB: Use global variables "Body_names", "Body_colors"
// ---------------------------------------------------------
function [icol] = getBodyColor(name)
  k = find(Body_names == name); 
  icol = color(Body_colors(k)); 
endfunction

// ---------------------------------------------------------
// body name => size (points)
// NB: Use global variables "Body_names", "Body_sizes" 
// ---------------------------------------------------------
function [sz] = getBodySize(name)
  k = find(Body_names == name); 
  sz = Body_sizes(k); 
endfunction

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

// ---------------------------------------------------------
// plots body trajectory
// col = color_index
// pos: in drawing units
// ---------------------------------------------------------
function plotBodyPos(pos, col)
  param3d(pos(1,:)/UNIT, pos(2,:)/UNIT, pos(3,:)/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, sz)
  param3d(pos(1,:)/UNIT, pos(2,:)/UNIT, pos(3,:)/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 = sz;
endfunction

// ---------------------------------------------------------
// ajust plot
// vmin, vmax: in drawing units (3x1)
// ---------------------------------------------------------
function adjustPlot(vmin, vmax)
  // db = data_bounds = [xmin, ymin, zmin; xmax, ymax, zmax]
  // 8% margin taken (longest side)
  L = max(max((vmax-vmin)*1.08, "r")); 
  db = [(vmin+vmax)/2 - L/2, (vmin+vmax)/2 + L/2]' / UNIT; 
  a = gca(); 
  a.data_bounds = db; 
  a.tight_limits = "on"; 
  a.isoview = "on"; 
  a.view = "2d"; // XY by default 
  CL_g_stdaxes(a);
  a.grid_position = "background"; 
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", frame); 
    mprintf("Center = %s\n", cen); 
    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
Body_names = [ "Sun", "Mercury", "Venus", "Earth", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune", "Pluto", "Moon"];
Body_colors = ["yellow2", "grey", "orange", "blue", "red", "brown", "green", "slateblue", "darkturquoise", "black", "lightblue"]; 

Frames = ["EOD", "ICRS", "ITRS"];

// Initial start date = current time (TREF)
scal = CL_dat_cal2str(CL_dat_cjd2cal(floor(CL_dat_now()))); 
dur = [60,0]; // days before/after 
names = "Su Me V E Ma J Sa"; // body names
cen = "Sun"; // center (body name)
iframe = 1; 
tsteps = [1, 0]; // steps (trajectory/marker) - days

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("Body names", names, typ='s', valid="getBodyNames($x)<>[]"), ...
  CL_defParam("Frame (1=EOD, 2=ICRS, 3=ITRF)", iframe, accv = 1:3), ...
  CL_defParam("Frame center (body name)", cen, typ='s', valid="getBodyName($x)<>[]"), ...
  CL_defParam("Time steps for trajectory/markers (0=no plot)", tsteps, dim=2, units=["day"], valid="$x>=0") ...
);

[scal, dur, names, iframe, cen, tsteps] = 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);
names = getBodyNames(names); 
cen = getBodyName(cen); 
frame = Frames(iframe); 
traj_step = tsteps(1);
mark_step = tsteps(2);
ind = getBodyIndices(names); 

// dates for trajectory or markers
cjd_traj = cjd_ref; 
if (traj_step > 0)
  cjd_traj = cjd_start : traj_step : cjd_end;
end
cjd_mark = []; 
if (mark_step > 0)
  cjd_mark = [(cjd_ref-mark_step : -mark_step : cjd_start), (cjd_ref+mark_step : mark_step : cjd_end)];
end

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

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

// min/max position values for scale adjustments
vmin =  %inf * ones(3,1); 
vmax = -%inf * ones(3,1); 

// header 
printPos("start"); 

// legend handles
hleg = []; 

// load DE405 ephemeris file 
cjd_all = [cjd_traj, cjd_mark, cjd_ref]; 
ephem = CL_eph_de405LoadT(cjd_all); 

// Plot trajectories
for k = 1 : nb
  // get all needed positions
  pos = getBodyPos(names(k), cjd_all, cen, frame, ephem);
  
  // Plot trajectories
  plotBodyPos(pos(:, 1:length(cjd_traj)), getBodyColor(names(k)));
  hleg($+1) = gce(); 
  vmin = min([vmin, pos], "c"); 
  vmax = max([vmax, pos], "c"); 
  
  // Plot time markers
  if (cjd_mark <> [])
    plotBody(pos(:, length(cjd_traj)+1:$-1), color("black"), 3); 
  end
  
  // Plot body positions at ref date
  // size of biggest body = 10 points
  plotBody(pos(:,$), getBodyColor(names(k)), 5);
  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) + " - " + frame + " (3D)";
a.x_label.text = "X (10^6 km)";
a.y_label.text = "Y (10^6 km)";

adjustPlot(vmin, vmax); 

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