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

// -----------------------------------------------------------
//> Effect of gravity (all harmonics) on a geostationary satellite's orbit:
//> - Longitudinal acceleration as a function of longitude 
//> (mainly caused by C22/S22 and to a lesser extent C31/S31 and C33/S33)
//> - semi-major axis offset / Keplerian value as a function of longitude
//> (main effect is due to J2)
//
// Author: A. Lamy
// -----------------------------------------------------------


// ------------------------------------
// Initializations 
// ------------------------------------

deff("[n] = NMAX()", "n = max(size(%CL_cs1nm))"); 

// max degree/order of potential
nmax = 10; 

// Plot option (see meaning below)
iplot = 1; 

desc_param = list(..
   CL_defParam("Max degree and order of potential (>=2 and <=" + msprintf("%d",NMAX()) + ")", nmax, accv=2:NMAX()),..
   CL_defParam("Plot: 1=longitudinal acceleration, 2=sma offset", iplot, accv=[1,2]) ..
);

[nmax, iplot] = CL_inputParam(desc_param);


// ------------------------------------
// Computations 
// ------------------------------------

// longitudes
lon = linspace(0, 2*%pi, 361); 

// Mean motion and semi major axis (Keplerian orbit)
// Reference frame = ECI
n = %CL_rotrBody; 
sma = (%CL_mu / n^2)^(1/3); 

// Position of satellite (frame = ECF)
pos_sph = [lon; zeros(lon); sma*ones(lon)]; 
pos = CL_co_sph2car(pos_sph); 

// Acceleration due to harmonics (central term not included)
// coords in ECF
acc = CL_fo_sphHarmAcc(pos, [nmax, nmax]); 

// Orbital (equinoctial) elements + time derivatives
// frame = inertial (=ECF frozen/ECI)
equin = [sma*ones(lon); zeros(4, length(lon)); lon]; 
equindot = CL_op_gauss("equin", equin, "inertial") * acc; 

// longitudinal acceleration (3rd Kepler's law)
lon2dot = (-3/2) * (n/sma) * equindot(1,:); 

// semi-major axis offset (3rd Kepler's law + change of sign)
dsma = (2/3) * (sma/n) * equindot(6,:); 
dsma_mean = CL_mean(dsma, meth="simp"); 

// computes stable (negative slope) and unstable longitudes
lonz_stab = CL_detectZero(lon, lon2dot, direct="decr"); 
lonz_unst = CL_detectZero(lon, lon2dot, direct="incr"); 
 

// ------------------------------------
// Prints a few results
// ------------------------------------

mprintf("-------------------------------------------\n"); 

mprintf("Keplerian semi-major axis (km) = %.3f\n", sma/1000); 
mprintf("Corrected semi-major axis - mean value (km) = %.3f\n", (sma+dsma_mean)/1000); 

mprintf("Min/max longitudinal acceleration (deg/day^2) = %.3g / %.3g\n", min(lon2dot*%CL_rad2deg*86400^2), max(lon2dot*%CL_rad2deg*86400^2)); 

str = msprintf("%.3f\n", lonz_stab' * %CL_rad2deg)'; 
mprintf("Stable equilibrium longitudes (deg) = %s\n", strcat(str, " ")); 

str = msprintf("%.3f\n", lonz_unst' * %CL_rad2deg)'; 
mprintf("Unstable equilibrium longitudes (deg) = %s\n", strcat(str, " ")); 
mprintf("\n"); 



// ------------------------------------
// Plots 
// ------------------------------------

f = scf(); 
f.visible = "on"; 
f.immediate_drawing = "off"; 

// Plot longitudinal acceleration
if (iplot == 1)

  y = lon2dot * %CL_rad2deg * 86400^2 / 0.001;  
  plot(lon*%CL_rad2deg, y, "k"); 

  // stable equilibrium longitudes
  plot(lonz_stab*%CL_rad2deg, zeros(lonz_stab), "bo"); 

  // unstable equilibrium longitudes
  plot(lonz_unst*%CL_rad2deg, zeros(lonz_unst), "ro"); 

  // adjust plot
  a = gca(); 

  // change background color for symbols, thickness for lines
  h = CL_g_select(a, "Polyline"); 
  CL_g_set(h, "mark_background", h.mark_foreground); 
  h.thickness = 2; 

  // change bounds
  dy = (max(y) - min(y))/20; // margin
  a.data_bounds = [0, min(y) - dy; 360, max(y)+dy]; 
  a.tight_limits = "on"; 

  a.title.text = "Longitudinal acceleration (1.e-3 deg/day^2) - EGM96s (" + string(nmax) + "x" + string(nmax) + ")"; 
  a.x_label.text = "Longitude (deg)"; 
  CL_g_stdaxes(a); 

end


// Plot semi major axis offset
if (iplot == 2)

  y = dsma / 1000; 
  plot(lon*%CL_rad2deg, y, "k"); 
  plot(lon*%CL_rad2deg, (dsma_mean/1000) * ones(lon), "b"); 

  // adjust plot
  a = gca(); 

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

  // change bounds
  dy = (max(y) - min(y))/20; // margin
  a.data_bounds = [0, min(y) - dy; 360, max(y)+dy]; 
  a.tight_limits = "on"; 

  a.title.text = "Semi-major axis offset / Keplerian value (km) - EGM96s (" + string(nmax) + "x" + string(nmax) + ")"; 
  a.x_label.text = "Longitude (deg)"; 
  CL_g_stdaxes(a); 

end

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