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

// -----------------------------------------------------------
//> Trajectory of objects from TLEs
//> 
//> The TLE are read from a file specified by a path or an URL. By default (empty name), 
//> "tle_examples.txt" in CelestLab data directory is used. 
//> 
//> For up-to-date data, go for instance to http://www.celestrak.com/NORAD/elements. 
//
// Author: A. Lamy
// -----------------------------------------------------------

fname=""; 
objname = "SIRIUS"; 
dur = 1; 
step = 60; 
iframe = 1; 

// validate the file name (does not check URL is valid)
function [ok] = validFname(fname)
  ok = ( stripblanks(fname,%t) == ""  | ..
         part(fname, 1:5) == "http:" | ..
         isfile(fname) ); 
endfunction


desc_param = list(..
  CL_defParam("TLE file path or URL (empty = default file)", fname, typ="s", valid='validFname($x)'),..
  CL_defParam("Part of object(s) name(s)", objname, typ="s", valid='stripblanks($x,%t) <> """"'),..
  CL_defParam("Propagation duration", dur, units = ["day"], valid ='$x>0 & $x < 100'),..
  CL_defParam("Propagation step", step, units = ["s"], valid='$x >0'),..
  CL_defParam("Frame: 1=ECI, 2=ECF, 3=TEME", iframe, accv = [1,2,3])..
);

[fname, objname, dur, step, iframe] = CL_inputParam(desc_param);

FRAMES = ["ECI", "ECF", "TEME"]; 
frame = FRAMES(iframe); 

// -----------------------------
// Computation 
// -----------------------------
// Internal function: 
// Load a TLE file
// may return one or more TLE
function [tle] = loadTLE(fname)
  // default output
  tle = CL_tle_new(0); 
  
  if (stripblanks(fname,%t) == "")
    fpath = CL_path("tle_examples.txt", fullfile(CL_home(), "data", "misc"));
  elseif(part(fname, 1:5) == "http:")
    fpath = getURL(fname, fullfile(TMPDIR, "CL_demos_TLE.txt")); 
  elseif (isfile(fname))
    fpath = fname; 
  else
    fpath = []; 
  end

  if (fpath == [])
    error("Invalid TLE file path or URL"); 
  end
  tle = CL_tle_load(fpath);
endfunction


// Initialize TLE
tle = loadTLE(fname); 
I = grep(tle.desc, objname); 
if (I == []); error("Specified TLE not found"); end  
tle = tle(I); 


// -----------------------------------------------------------
// compute trajectories
// -----------------------------------------------------------
list_traj = list(); 

// time relative to TLE epoch
delta_t = 0 : step : dur * 86400; 

for k = 1:size(tle)
  cjd = tle(k).epoch_cjd + delta_t/86400; // considered as TREF
  [pos,vel] = CL_tle_genEphem(tle(k), cjd, frame); 
  res = struct(); 
  res.cjd = cjd; 
  res.pos = pos; 
  list_traj(k) = []; // to avoid possible Scilab bug
  list_traj(k) = res; 
end


// -----------------------------------------------------------
// display info
// -----------------------------------------------------------
// print desc and 
mprintf("Number of elements: %d\n", size(tle)); 
for k = 1:size(tle)
  str = stripblanks(tle(k).desc, %t); 
  if (str == ""); str = "?"; end
  mprintf("%d: %s [%s]\n", k, tle(k).intldesg, str); 
end  
mprintf("\n"); 
  
  
// -----------------------------------------------------------
// plot
// -----------------------------------------------------------

function plot_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 plot_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 = color("grey90");
  e.foreground = color("grey80");

  // origin meridian / equator
  plot_traj(R*[cos(u);zeros(u);sin(u)], F=1, col=color("grey70"), th=1);  
  plot_traj(R*[cos(v);sin(v);zeros(v)], F=1, col=color("grey70"), th=1);  

endfunction


f = scf();
f.visible = "on"; 
f.immediate_drawing = "off"; 
ech = 1000; // unit = 1 km
	
// sphere representing the Earth, 
plot_sphere(%CL_eqRad/ech);

// orbit
vmin = +%inf * ones(3,1); // min x, y, z
vmax = -%inf * ones(3,1); // max x, y, z
cols = [2, 5, 13, 24, 32, 28]; // color indices
nbcols = length(cols); 
for k = 1:size(tle)
  pos = list_traj(k).pos; 
  vmin = min(vmin, min(pos, "c")); 
  vmax = max(vmax, max(pos, "c")); 
  param3d(pos(1,:)/ech, pos(2,:)/ech, pos(3,:)/ech);
  h = CL_g_select(gce(), "Polyline"); 
  h.foreground = cols(modulo(k-1,nbcols)+1); 
  h.thickness = 2; 
end

a = gca();
L = max(vmax-vmin); 
cen = (vmax+vmin)/2; 
a.data_bounds = [cen - (L/2) * 1.05, cen + (L/2) * 1.05]' / ech; 
a.tight_limits = "on"; 
a.isoview = 'on'; 
a.rotation_angles = [80, 20]; 
  
a.title.text = "Object(s): " + objname + " - Frame = " + frame; 
CL_g_stdaxes(); 

// CL_g_legend(a, tle.desc');

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