//  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 [pot] = CL_fo_srpPot(pos, pos_sun, coefp, ecl, er, ersun, p0)
// Potential due to solar radiation pressure
//
// Calling Sequence
// [pot] = CL_fo_srpPot(pos, pos_sun, coefp [,ecl, er, ersun, p0])
//
// Description
// <itemizedlist>
// <listitem>
// <p>Potential due to solar radiation pressure (SRP).</p> 
// <p>By definition, the acceleration derived from the potential is +grad(potential).</p>
// <p></p>
// <p>The SRP coefficient (<b>coefp</b>) is defined by: coefp = cp * area / mass (with cp: 
// reflectivity coefficient, between 1 and 2). </p> 
// <p></p>
// <p>Eclipses can be taken into account or not. If so, the acceleration is multiplied by a factor 
// (equal to 1 if there is no eclipse and less than 1 otherwise).</p>
// <p></p></listitem>
// <listitem>
// <p>Notes:</p>
// <p>- The origins for the satellite and Sun position vectors must be the same. 
// If eclipses are considered, this origin must be the eclipsing body (i.e. central body).</p> 
// <p>- The coordinates frame can be any frame.</p> 
// <p>- The calculation of the eclipse uses the radius of the eclipsing body (<b>er</b>), 
// and the radius of the Sun (<b>ersun</b>). If <b>ersun</b> is empty ([]) then the internal value is used. </p>
// <p>- <b>p0</b> is the solar radius pressure at 1 AU. If <b>p0</b> is empty ([]) then the internal value is used (equal
// to the total solar irradiance divided by the speed of light). </p>
// <p></p></listitem>
// <listitem>
// <p>See <link linkend="Force models">Force models</link> for more details.</p>
// <p></p></listitem>
// </itemizedlist>
//
// Parameters
// pos: Position vector (from any position) [m]. (3xN or 3x1)
// pos_sun: Sun position (same origin as pos) [m]. (3xN or 3x1)
// coefp: SRP coefficient (cp*area/mass) [m^2/kg]. (1xN or 1x1)
// ecl: (optional, boolean) %t if eclipses are taken into account; %f otherwise. Default is %t. (1x1)
// er: (optional) Equatorial radius of eclipsing body. Default is %CL_eqRad. [m] (1x1)
// ersun: (optional) Equatorial radius of the Sun. Default is [] (internal value is used). [m] (1x1)
// p0: (optional) Solar radiation pressure at 1 AU. Default is [] (internal value is used). [N/m^2] (1x1)
// pot: Potential [m^2/s^2]. (1xN)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_fo_srpAcc
// CL_gm_eclipseCheck
//
// Examples
// pos = [7000.e3; 0; 0]; // ECI
// pos_sun = CL_eph_sun(CL_dat_cal2cjd(2000,3,21)); // ECI
// coefp = 1.5 * 10 / 1000; 
// CL_fo_srpPot(pos, pos_sun, coefp, %t)


// Declarations:
AU = CL__dataGetEnv("au", internal=%t);

// Code:
if (~exists('ecl','local')); ecl = %t; end
if (~exists("er", "local")); er = CL__dataGetEnv("eqRad"); end
if (~exists('ersun','local')); ersun = []; end
if (~exists('p0','local')); p0 = []; end

if (ersun == []); ersun = CL__dataGetEnv(["body", "Sun", "eqRad"]); end;
if (p0 == []); p0 = CL__dataGetEnv("totalSolarIrradiance", internal=%t) / CL__dataGetEnv("lightSpeed", internal=%t); end;

[pos, pos_sun, coefp] = CL__checkInputs(pos, 3, pos_sun, 3, coefp, 1); 

d = CL_norm(pos - pos_sun); 

// Solar Radiation Pressure at actual distance
// (no eclipse)

pot = -p0 * coefp .* (AU^2) ./ d;

if (ecl)
  // frac = 1 if no eclipse, less than 1 otherwise
  pos_earth = [0;0;0]; 
  frac = 1-CL_gm_eclipseCheck(pos, pos_sun, pos_earth, ersun, er);
  pot = pot .* frac;
end


endfunction