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

// -----------------------------------------------------------
//> Prediction of visible passes using 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.   
//> Only one TLE should be selected using the selection criteria (name). 
//> For up-to-date data, go to http://www.celestrak.com/NORAD/elements, and in particular: 
//> http://www.celestrak.com/NORAD/elements/stations.txt for the ISS. 
//> 
//> The conditions that should be met for the object to be visible are: 
//> - object elevation above minimum value as seen from the observer at specified location, 
//> - object is totally or partly lit by the Sun, 
//> - observer is in the dark (Sun elevation from his location is less than 
//>   the minimum elevation). 
//> 
//> The default setting (minimum visibility elevation and Sun elevation at beginning 
//> of night) seem to be consistent with http://www.heavens-above.com. 
//> 
//> The Sun elevation threshold (-6 deg) defines the beginning/end of civil twilight. 
//> See http://en.wikipedia.org/wiki/Twilight#Civil_twilight for more details. 
//> 
//> The looked-for visibility periods of time are plotted in yellow. 
//> 
//> Note: UT1 and UTC are supposed identical. 
//> 
// Author: A. Lamy
// -----------------------------------------------------------

// -----------------------------
// Inputs 
// -----------------------------

// TLE file name
fname = ""; 

// name of object
objname = "ISS"; 

// location (geodetic coordinates)
lon = 1.499*%CL_deg2rad; 
lat = 43.43*%CL_deg2rad;
alt = 150; // m

// min elevation from location for visibility
elevmin = 10 * %CL_deg2rad; 

// Sun elevation maximum value 
elevmax_sun = -6 * %CL_deg2rad; // civil twilight

// Prediction time span (days) from TLE epoch
pred_start = 0; 
pred_len = 7;  

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

desc = list(..
  CL_defParam("TLE file path or URL (empty = default file)", fname, typ="s", valid='validFname($x)'),..
  CL_defParam("Part of object''s name", objname, typ="s", valid='stripblanks($x,%t) <> """"'),..
  CL_defParam("Location geodetic longitude", lon, units=['rad', 'deg'], valid='$x >= 0 & $x <= 360'),..
  CL_defParam("Location geodetic latitude", lat, units=['rad', 'deg'], valid='$x >= -90 & $x <= 90'),..
  CL_defParam("Location geodetic altitude", alt, units=['m'], valid='$x >= -10.e3 & $x <= 10.e3'),..
  CL_defParam("Minimum elevation for visibility", elevmin, units=['rad', 'deg'], valid='$x >= -90 & $x < 90'),..
  CL_defParam("Sun elevation at beginning/end of night", elevmax_sun, units=['rad', 'deg'], valid='$x >= -90 & $x < 90'),..
  CL_defParam("Prediction start / TLE epoch", pred_start, units=["day"], valid='$x >= 0 & $x <= 100'),..
  CL_defParam("Prediction time length", pred_len, units=["day"], valid='$x >= 1 & $x <= 100')..
);

[fname, objname, lon, lat, alt, elevmin, elevmax_sun, pred_start, pred_len] = CL_inputParam(desc);


// -----------------------------
// 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:7) == "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(1)); 


// Time instants for computation (TREF)
cjd0 = floor(tle.epoch_cjd + pred_start); // considered here as TREF 
cjdf = ceil(cjd0 + pred_len); // considered here as TREF 
cjd = [(cjd0 : 60 / 86400 : cjdf - 0.1), cjdf];
 
// Object position (ECF)
[pos_ecf] = CL_tle_genEphem(tle, cjd, "ECF"); 

// Sun position (ECF)
pos_sun_eci = CL_eph_sun(cjd, frame="ECI"); 
pos_sun_ecf = CL_fr_convert("ECI", "ECF", cjd, pos_sun_eci); 

// location in geodetic coordinates (ECF)
location = [lon; lat; alt]; 

// Prediction interval 
interv_predict = [cjd(1); cjd($)]; 

// Geometrical visibility intervals of object from location
interv_obj_visib = CL_ev_stationVisibility(cjd, pos_ecf, location, elevmin); 

// Visibility intervals of Sun from location (<=> day or night)
interv_daylight = CL_ev_stationVisibility(cjd, pos_sun_ecf, location, elevmax_sun); 
interv_night = CL_intervInv(interv_predict, interv_daylight); 

// Intervals when object is (partly or totally) lit by the Sun 
interv_obj_ecl = CL_ev_eclipse(cjd, pos_ecf, pos_sun_ecf, typ="umb"); 
interv_obj_lit = CL_intervInters(interv_obj_visib, CL_intervInv(interv_predict, interv_obj_ecl)); 

// Visibility periods
// => intersection of all intervals
interv = CL_intervInters(interv_obj_visib, interv_night, interv_obj_lit); 


// -----------------------------
// Function for processing of results 
// -----------------------------

// computes elevation, azimuth(+ towards east), aspect angle
function [elev, azim, aspect] = direction(t)
  pos_ecf_t = CL_interpLagrange(cjd, pos_ecf, t); 
  pos_sun_ecf_t = CL_interpLagrange(cjd, pos_sun_ecf, t); 
  [elev, azim] = CL_gm_stationPointing(location, pos_ecf_t, res= ["elev", "azim"]); 
  azim = CL_rMod(-azim, 0, 2*%pi); // => positive towards East
  pos_loc = CL_co_ell2car(location) * ones(t); 
  aspect = CL_vectAngle(pos_loc - pos_ecf_t, pos_sun_ecf_t - pos_ecf_t); 
endfunction
  

// determines cardinal direction (S, W ...) 
function [s] = str_azim(azim)
  AZIM = [-%inf, linspace(%pi/16, 2*%pi - %pi/16, 16), %inf]; 
  STR = ["N", "NNE", "NE", "ENE", "E", "ESE", "SE", "SSE", ..
         "S", "SSW" "SW", "WSW", "W", "WNW", "NW", "NNW", "N", ""]; 
  k = find(azim >= AZIM(1:$-1) & azim < AZIM(2:$)); 
  s = STR(k); 
endfunction


// plots an interval (elevation plot)
function plot_interv(cjd0, interv, hminmax, colorname, thickness)
  if (interv == []); return; end

  hminmax = hminmax * ones(interv(1,:));  
  rects = [interv(1,:) - cjd0; hminmax(2,:); interv(2,:)-interv(1,:); hminmax(2,:)-hminmax(1,:)]; 
  xrects(rects, color(colorname) * ones(rects(1,:))); 
endfunction


// Plot elevation
function plot_elev(cjd0, interv, delevdtmax, colorname, thickness)
  if (interv == []); return; end
   
  steps = min(1 * %CL_deg2rad / delevdtmax, 60/86400); 
  d = interv(2,:) - interv(1,:); 
  nb = max(2,max(round(d ./ steps))); 
  t = CL_intervLinspace(interv, nb); 
  t = matrix([t; %nan * ones(t(1,:))], 1, -1); 
  elev = direction(t);  
  
  plot(t-cjd0, elev * %CL_rad2deg); 
  e = CL_g_select(gce(), "Polyline"); 
  e.foreground = color(colorname); 
  e.thickness = thickness; 
endfunction


// -----------------------------
// Print results 
// -----------------------------

str = [ ..
  "----------------------------------------------------------------"; ..
  msprintf("Visible passes for object %s [%s]\n", tle.intldesg, stripblanks(tle.desc, %t)); ..
  "Date (TREF), El (deg), Dir  [Az (deg), Aspect angle (deg)]"; ..
  "----------------------------------------------------------------" ..
]; 
       
t = matrix(interv(:), 1, -1); 
[elev, azim, aspect] = direction(t); 
cal_str = CL_dat_cal2str(CL_dat_cjd2cal(t));

cal_str = matrix(cal_str, 2, -1); 
elev = matrix(elev, 2, -1);
azim = matrix(azim, 2, -1);
aspect = matrix(aspect, 2, -1);

for k = 1 : size(interv,2)
  for (i = 1:2)
    str($+1) = msprintf("%s %5.1f %-3s  [%5.1f %5.1f]\n",  part(cal_str(i,k), 1:19), elev(i,k)*%CL_rad2deg, .. 
                       str_azim(azim(i,k)), azim(i,k)*%CL_rad2deg, aspect(i,k)*%CL_rad2deg); 
  end
  str($+1) = msprintf("%.1f min\n", (interv(2,k) - interv(1,k)) * 1440);                             
  str($+1) = ""; 
  
end


mprintf("%s\n", str); 


// max derivative of elevation (rad/day) for computation of step 
elev = direction(cjd); 
delevdtmax = max(abs((elev(2:$) - elev(1:$-1)) ./ (cjd(2:$)-cjd(1:$-1)))); 


// -----------------------------
// Plot 
// -----------------------------

f = scf(); 
f.figure_size = [900, 500]; 
f.visible = "on"; 
f.immediate_drawing = "off";

a=gca();
a.clip_state = "clipgrf";

plot(0,0); 

plot_interv(cjd0, interv_daylight, [elevmin; %pi/2] * %CL_rad2deg, "slategray1"); 
plot_interv(cjd0, interv_night, [elevmin; %pi/2] * %CL_rad2deg, "slategray3"); 

plot_elev(cjd0, interv_obj_visib, delevdtmax, "darkblue", 1); 
plot_elev(cjd0, interv, delevdtmax, "yellow", 3); 

a.data_bounds = [cjd(1)-cjd0, elevmin*%CL_rad2deg; cjd($)-cjd0, 90]; 
a.tight_limits = "on"; 

step = max(round((cjd($) - cjd(1))/7), 1); 
x_ticks = 0 : step : cjd($) - cjd(1); 
cal = CL_dat_cjd2cal(cjd0 + x_ticks); 
a.x_ticks = tlist(["ticks", "locations", "labels"], x_ticks, ..
                    string(cal(1,:))+"/"+string(cal(2,:))+"/"+string(cal(3,:))); 

a.title.text = "Visible passes - object = " + msprintf("%s (%s)", objname, tle.intldesg) + " - " + ..
               msprintf("lon=%.2f deg, lat=%.2f deg", lon*%CL_rad2deg, lat*%CL_rad2deg); 
cal0 = CL_dat_cjd2cal(cjd0); 
a.x_label.text = "Year / Month / Day (TREF)";  
a.y_label.text = "Elevation (deg)"; 

CL_g_stdaxes(a); 
a.sub_ticks(1) = 0; 
a.font_size = 2; 

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