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

// -----------------------------------------------------------
//> Short periods due to J2, J3... either on: 
//> - orbital elements (adapted to circular orbits)
//> - components in "qsw" local frame
//>  
//> x-axis: argument of latitude (degrees)
//> Blue line: mean value
//> 
//> NB: 
//> - LydanneLp model is used to compute the short periods. 
//> - The mean orbit parameters are all constant except the mean 
//>   argument of latitude (which varies linearly with time).  
//
// Author: A. Lamy
// -----------------------------------------------------------


sma = %CL_eqRad + 700.e3; // sma minus equatorial radius
pom = %pi/2; 
inc = 98*%CL_deg2rad; 
ecc = 1.e-3; 
iplot = 1; // 1: orbital elements   2: components in "qsw" 

desc_param = list(..
   CL_defParam("Mean semi-major axis", sma, units=['m', 'km'], valid='$x>0'),..
   CL_defParam("Mean eccentricity", ecc,  valid="$x>=0 & $x<1"  ),..
   CL_defParam("Mean inclination", inc, units=['rad', 'deg'], valid="$x>=0 & $x<=180"  ),..
   CL_defParam("Mean argument of perigee", pom, units=['rad', 'deg']),..
   CL_defParam("Show short periods on: (1=orbital elements, 2=QSW components)", iplot, accv=1:2)..
);

[sma, ecc, inc, posm, iplot] = CL_inputParam(desc_param);

ex = ecc*cos(pom); 
ey = ecc*sin(pom); 
circ0 = [sma; ex; ey; inc; 0; 0]; 

T = CL_kp_params('per', sma)/86400;  // orbital period in seconds
t0 = 0; 
t = linspace(0,T,200); // final times (days, origin = 1950.0)

circ = CL_ex_kepler(t0, circ0, t);
pso = (t/T) * 2*%pi; 
circ(6,:) = pso; 

// NB: no propagation !
[mean_cir,osc_cir] = CL_ex_propagate("lydlp", "cir", t, circ, t, "mo"); 

dcirc = osc_cir - mean_cir; 
dcirc(5:6,:) = CL_rMod(dcirc(5:6,:),-%pi, %pi);

// QSW computed from mean elements
[p_nom,v_nom] = CL_oe_cir2car(mean_cir);
[p,v] = CL_oe_cir2car(osc_cir);
[dp_qsw] = CL_fr_inertial2qsw(p_nom, v_nom, p-p_nom);


// -----------------------------------------------------------
// plot 1 (Orbital elements)
// -----------------------------------------------------------
function Plot1(pso, dcirc)

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

  unit_coef = [1.e-3,1.e3,1.e3,%CL_rad2deg,%CL_rad2deg,%CL_rad2deg];
  unit = ["(km)", "(1.e-3)", "(1.e-3)", "(deg)", "(deg)", "(deg)"]; 
  param = ["Semi-major axis", "ex", "ey", "Inclination", "RAAN", "Arg. of latitude"]; 

  [xmi, xma, nb1] = CL_graduate(min(pso*%CL_rad2deg), max(pso*%CL_rad2deg));
  xgrad = linspace(xmi, xma, nb1+1);
  ordre = [1,4,2,3,5,6]; // plot order

  for i = 1:6
    subplot(3,2,i);
    k = ordre(i); 
    a=gca();
    y = dcirc(k,:)*unit_coef(k); 
    plot2d(pso*%CL_rad2deg, y);
    plot2d(pso*%CL_rad2deg, mean(y)*ones(pso), style=2);

    a.title.text = param(k) + " " + unit(k);
    //a.x_label.text = "Argument of latitude (deg) ";
    a.x_ticks = tlist("ticks", xgrad, string(xgrad));

    // tight limits for x-axis only. 
    db = a.data_bounds; 
    db(:,2) = [min(y) - (max(y)-min(y)) * 0.1; max(y) + (max(y)-min(y)) * 0.1]; 
    a.data_bounds = db; 
    a.tight_limits(1) = "on"; 
    
    a.margins = [0.15,0.1,0.18,0.18];
    CL_g_stdaxes(a, fg=1, ft=2, fl=1);
    a.sub_ticks = [0,0]; 
  end

  h = CL_g_select(f, "Polyline");
  h.thickness = 2;

  f.immediate_drawing="on";

endfunction

// -----------------------------------------------------------
// plot 2 (QSW)
// -----------------------------------------------------------

function Plot2(pso, dp_qsw)

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

  param = ["Radial axis (m)", "Tangential axis (m)", "Normal axis (m)"]; 
  [xmi, xma, nb1] = CL_graduate(min(pso*%CL_rad2deg), max(pso*%CL_rad2deg));
  xgrad = linspace(xmi, xma, nb1+1);

  for i=1:3
    k = i;

    subplot(3,1,k);
    a=gca();
    y = dp_qsw(k,:); 

    plot2d(pso*%CL_rad2deg, y);
    plot2d(pso*%CL_rad2deg, mean(y)*ones(pso), style=2);

    a.title.text = param(k);
    a.x_ticks = tlist("ticks", xgrad, string(xgrad));
   
    // tight limits for x-axis only. 
    db = a.data_bounds; 
    db(:,2) = [min(y) - (max(y)-min(y)) * 0.1; max(y) + (max(y)-min(y)) * 0.1]; 
    a.data_bounds = db; 
    a.tight_limits(1) = "on"; 

    a.margins = [0.15,0.1,0.18,0.18];
    CL_g_stdaxes(a, fg=1, ft=2, fl=1);
    a.sub_ticks = [0,0]; 
  end

  h = CL_g_select(f, "Polyline");
  h.thickness = 2;

  f.immediate_drawing="on";

endfunction

// change 'visible' status when everything is drawn 
// (which activates pointer location)
if (iplot == 1); Plot1(pso, dcirc); end
if (iplot == 2); Plot2(pso, dp_qsw); end