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

function [output] = CL_op_ssoJ2(type_output,par1,par2,er,mu,j2,rotr_pla_sun)
// Semi major axis, eccentricity or inclination for a sun synchronous orbit (J2)
//
// Calling Sequence
// [output] = CL_op_ssoJ2(type_output,par1,par2 [,er,mu,j2,rotr_pla_sun])
// [sma] = CL_op_ssoJ2("a",ecc,inc [,er,mu,j2,rotr_pla_sun])
// [ecc] = CL_op_ssoJ2("e",sma,inc [,er,mu,j2,rotr_pla_sun])
// [inc] = CL_op_ssoJ2("i",sma,ecc [,er,mu,j2,rotr_pla_sun])
//
// Description
// <itemizedlist><listitem>
// <p>Computes orbital parameters (semi major axis, eccentricity or inclination) ensuring sun-synchronicity so that 
// the orbit's ascending node precession rate is equal to the Sun apparent rotation rate <b>rotr_pla_sun</b>.(See equation (1))</p>
// <p> - If <b>type_output</b>="a", <b>par1</b> is eccentricity, <b>par2</b> is inclination and semi major axis is computed (<b>output</b>).</p>
// <p> - If <b>type_output</b>="e", <b>par1</b> is semi major axis, <b>par2</b> is inclination and eccentricity is computed (<b>output</b>).</p>
// <p> - If <b>type_output</b>="i", <b>par1</b> is semi major axis, <b>par2</b> is eccentricity and inclination is computed (<b>output</b>)</p>
// <para/><inlinemediaobject><imageobject><imagedata fileref="helio.gif"/></imageobject></inlinemediaobject></listitem>
// <p><b>Note:</b> %nan is returned when sun-synchronicity is not possible</p>
// </itemizedlist>
//
// Parameters
// type_output: String defining the parameter to be computed (output). It can be "a" for semi major axis, "e" for eccentricity or "i" for inclination
// par1: Semi major axis if type_output is "e" or "i" ; eccentricity if type_output is "a" [m] (1xN or 1x1)
// par2: Inclination if type_output is "a" or "e" ; eccentricity if type_output is "i" [rad] (1xN or 1x1)
// er: (optional) Equatorial radius [m] (default is %CL_eqRad)
// mu: (optional) Gravitational constant [m^3/s^2] (default value is %CL_mu)
// j2: (optional) Second zonal harmonic term (default is %CL_j1jn(2))
// rotr_pla_sun : (optional) Mean apparent rotation rate of the Sun around the planet (default is %CL_rotrBodySun) (1 x 1)
// output: Semi major axis [m] if type_output is "a" ; eccentricity if type_output is "e" ; inclination [rad] if type_output is "i" (1xN). Output value is %nan if no valid result exists. 
//
// Authors
// CNES - DCT/SB
//
// Examples
// // Compute inclination for sun-synchronicity
// sma = 7078.e3;
// ecc = 0.01;
// inc = CL_op_ssoJ2 ("i", sma, ecc)
//
// // Compute eccentricity for sun-synchronicity
// sma = 7078.e3;
// inc = CL_deg2rad(98.15);
// ecc = CL_op_ssoJ2 ("e", sma, inc)
//
// // Compute semi major axis for sun-synchronicity
// ecc = 0.01;
// inc = CL_deg2rad(97);
// sma = CL_op_ssoJ2 ("a", ecc, inc)

// Declarations:

// Code:
if (~exists("er", "local")); er = CL__dataGetEnv("eqRad"); end
if (~exists("mu", "local")); mu = CL__dataGetEnv("mu"); end
if (~exists("j2", "local")); j2 = CL__dataGetEnv("j1jn", 2); end
if (~exists("rotr_pla_sun", "local")); rotr_pla_sun = CL__dataGetEnv("rotrBodySun"); end

// Check input arguments and resize if necessary
[par1, par2, N] = CL__checkInputs(par1,1, par2,1);

K = -3.0 * j2 * (er.^2) * sqrt(mu) / (2.0 * rotr_pla_sun) // auxiliary (scalar) variable

// Calculation of semi major axis
if (type_output == "a")
  ecc = par1;
  inc = par2;

  I = find (ecc < 0 | ecc >= 1); 
  if (I <> []); CL__error("Eccentricity out of range"); end
  I = find (inc < 0 | inc > %pi+%eps); 
  if (I <> []); CL__error("Inclination out of range"); end

  // output = semi major axis
  tmp = K * cos(inc) ./ ((1-ecc.^2).^2);
  I = find (tmp > 0); 
  output = %nan * ones(1:N);  
  output(I) = tmp(I) .^ (1/3.5); 

// Calculation of eccentricity
elseif (type_output == "e")
  sma = par1;
  inc = par2;

  I = find (sma <= 0); 
  if (I <> []); CL__error("Semi major axis out of range"); end
  I = find (inc < 0 | inc > %pi+%eps); 
  if (I <> []); CL__error("Inclination out of range"); end

  // output = eccentricity
  tmp = K * cos(inc) ./ (sma .^ (3.5));
  I = find (tmp >= 0 & tmp < 1); 
  output = %nan * ones(1:N);  
  output(I) = sqrt(1-sqrt(tmp(I))); 

// Calculation of inclination
elseif (type_output == "i")
  sma = par1;
  ecc = par2;

  I = find (sma <= 0); 
  if (I <> []); CL__error("Semi major axis out of range"); end
  I = find (ecc < 0 | ecc >= 1); 
  if (I <> []); CL__error("Eccentricity out of range"); end

  // output = inclination
  tmp = (sma.^ (3.5) .* (1-ecc.^2).^2) / K;
  I = find (tmp >= -1 & tmp <= 1); 
  output = %nan * ones(1:N);
  output(I) = acos(tmp(I)); 

  
else
  CL__error("Unknown value for type_output");
end

endfunction