//  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 [acc] = CL_fo_srpAcc(pos, pos_sun, coefp, ecl, er, ersun, p0)
// Acceleration due to solar radiation pressure 
//
// Calling Sequence
// [acc] = CL_fo_srpAcc(pos, pos_sun, coefp [,ecl, er, ersun, p0])
//
// Description
// <itemizedlist>
// <listitem>
// <p>Acceleration due to solar radiation pressure (SRP).</p> 
// <p>The acceleration is directed from the Sun to the considered position and is 
// inversely proportional to the distance squared. </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 they are, the acceleration is multiplied by a factor 
// equal to 1 out of the eclipse region and less than 1 otherwise.</p>
// <p></p></listitem>
// <listitem>
// <p>Notes:</p>
// <p>- The coordinates frame can be any frame. The origin of the frame does not matter,  
// except if eclipses are taken into account, in which case the origin of the frame must be the
// center of the eclipsing body (usually the same as the central body). </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 an 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)
// acc: Acceleration [m/s^2]. (3xN)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_fo_srpPot
// 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_srpAcc(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); 

[u,d] = CL_unitVector(pos - pos_sun); 

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

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

// frac = 1 if no eclipse, less than 1 otherwise
frac = 1; 

if (ecl)
  pos_earth = [0;0;0]; 
  frac = 1-CL_gm_eclipseCheck(pos, pos_sun, pos_earth, ersun, er);
  K = K .* frac;
end

// acceleration
acc = CL_dMult(u, K); 

endfunction