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

// -----------------------------------------------------------
//> Evaluation of orbit lifetime (based on STELA propagation). 
//> 
//> Note that lifetime may be very sensitive to initial conditions
//> and that a statistical analysis might be required. 
//> 
//> Input orbital elements are mean orbital elements (in the "STELA" sense). 
//> The reference frame in which the orbit is defined is ECI. 
//> 
//> "Altitude" here means: norm of radius vector minus equatorial radius. 
//> It is computed using the mean orbital elements:  
//> - perigee altitude = sma * (1-ecc)   
//> - apogee altitude = sma * (1+ecc)  
//> 
//> Among the parameters that are implicitly defined: 
//> - The initial mean anomaly is 0 (<=> at the perigee)
//> - Thirdbody perturbation is included as well as Earth gravity 
//> up to degree and order 7. 
//> - If solar activity is constant, solar flux (F10.7) and geomagnetic
//> index (Ap) values are 140 and 9 respectively, which corresponds to a medium 
//> solar activity level. 
//> - A "-1" value for the areodynamic coefficient means that it is 
//> automatically computed as a function of altitude. 
//> 
//> An simulation file for STELA (tool) is optionally created.  
//> Note that the physical constants used in CelestLab depend on the user 
//> configuration and may not be the same as the ones used by STELA (tool).    
//
// Author: A. Lamy
// -----------------------------------------------------------

// -----------------------------------------------------------
// Input of parameters
// -----------------------------------------------------------

function [ok] = isValidDate(scal)  
  ok = %f; 
  cal = CL_dat_str2cal(scal); 
  if (cal(1) <> %nan)
    ok = ( cal(1) >= 2000 & cal(2) <= 2100);
  end
endfunction

// Input parameters default values
cal = CL_dat_now("cal");
scal = CL_dat_cal2str([cal(1);1;1]);  // january 1st
hp = 350.e3; 
ha = 900.e3; 
inc = 78 * %CL_deg2rad; 
pom = 90 * %CL_deg2rad; 
mltan = 12; 
area2mass = 1.e-2; 
cd = 2.2; 
cr = 1.5; 
solact = 1; 
maxdur = 10; // years
save_stela = 0; 

desc = list(..
   CL_defParam("Date (calendar format, TREF time scale)", scal, typ='cal', valid="isValidDate($x)"),..
   CL_defParam("Perigee altitude", hp, units=['m', 'km'], id='$hp', valid='$hp>=100' ),..
   CL_defParam("Apogee altitude",  ha, units=['m', 'km'], id='$ha', valid='$ha>=$hp' ),..
   CL_defParam("Inclination", inc, units = ["rad", "deg"], valid='$x >= 0 & $x <= 180'),..
   CL_defParam("Argument of perigee", pom, units = ["rad", "deg"], valid='$x >= 0 & $x <= 360'),..
   CL_defParam("Mean local time of the ascending node", mltan, units = ["h"], valid='$x >= 0 & $x <= 24'),..
   CL_defParam("Area to mass ratio (for drag and SRP)", area2mass, units = ["m^2/kg"], valid='$x >= 0'),..
   CL_defParam("Aerodynamic coefficient (for drag), -1 -> automatic", cd, valid='$x >= 0 | $x == -1'),..
   CL_defParam("Reflectivity coefficient (for SRP)", cr, valid='$x >= 0'),..
   CL_defParam("Solar activity: 1=constant (medium), 2=variable", solact, accv=1:2),..
   CL_defParam("Maximum length of simulation period", maxdur, units = ["yr"], valid='$x >= 1 & $x <= 100'),..
   CL_defParam("Create simulation file for STELA tool: 1=yes, 0=no", save_stela, accv=0:1)..
);

[scal, hp, ha, inc, pom, mltan, area2mass, cd, cr, solact, maxdur, save_stela] = CL_inputParam(desc);

SA_TYPES = ["constant", "variable"]; 
sa_type = SA_TYPES(solact); 


// -----------------------------------------------------------
// Computation
// -----------------------------------------------------------
// Initial date (cjd, TREF)
cjd0 = CL_dat_cal2cjd(CL_dat_str2cal(scal));

// Keplerian mean orbital elements
[sma, ecc] = CL_op_rarp2ae(%CL_eqRad+ha, %CL_eqRad+hp); 
gom = CL_op_locTime(cjd0, "mlh", mltan, "ra"); 
mean_kep0 = [sma; ecc; inc; pom; gom; 0];

// Computation dates (step = integration step = 1 day)
cjd = cjd0 + (0 : maxdur * 365.25);

// STELA model parameters 
params = CL_stela_params();
params.drag_area = area2mass * params.mass; 
params.drag_solarActivityType = sa_type; 
params.drag_solarActivityFlux = 140; 
params.drag_solarActivityAp = 9; 
if (cd == -1)
  params.drag_coefType = "variable"; 
else
  params.drag_coefType = "constant"; 
  params.drag_coef = cd; 
end
params.srp_area = area2mass * params.mass; 
params.srp_coef = cr; 

// Propagation
[mean_kep, info] = CL_stela_extrap("kep", cjd0, mean_kep0, cjd, params, ["m", "i"]);

// perigee/apogee altitudes
hp = mean_kep(1,:) .* (1 - mean_kep(2,:)) - %CL_eqRad; 
ha = mean_kep(1,:) .* (1 + mean_kep(2,:)) - %CL_eqRad; 


// -----------------------------------------------------------
// Save input file for STELA tool
// -----------------------------------------------------------
function [fname] = createStelaTmpName()
  fname = ""; 
  // 1000 tries
  for (n = 1 : 1000)
    fname = fullfile(TMPDIR, "CLdemos_lifetime_" + msprintf("%d", n) + "_sim.xml");
    if (~isfile(fname)); return; end // OK
  end
endfunction

if (save_stela)
  fname = createStelaTmpName(); 
  if (fname <> "")
    mputl(info.sim_xml, fname); 
    mprintf("\n=> Name of simulation file for STELA tool: \n%s\n\n", getlongpathname(fname)); 
  end
end


// -----------------------------------------------------------
// Plot
// -----------------------------------------------------------
// Plot ha/hp 
f = scf();
f.visible = "on"; 
f.immediate_drawing = "off"; 

time_unit = 365.25; // time unit for plot (1 year by default)
I = find(~isnan(ha)); 
if (I <> [])
  if (cjd(I($))-cjd(I(1)) < time_unit); time_unit = 1; end // unit = 1 day
  plot((cjd(I) - cjd0)/time_unit, hp(I) / 1000, "b", "thickness", 2); 
  plot((cjd(I) - cjd0)/time_unit, ha(I) / 1000, "r", "thickness", 2); 
end

if (time_unit == 1); text = "days"; else text = "years"; end

a = gca(); 
a.title.text = "Orbit lifetime"; 
a.x_label.text = "Time since beginning of simulation (" + text + ")"; 
a.y_label.text = "Perigee / Apogee altitude (km)"; 

CL_g_legend(a, ["Perigee", "Apogee"]); 
CL_g_stdaxes(a); 

f.immediate_drawing = "on";