//  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_srpPanelAcc(pos, pos_sun, coefp, normal, nsides, ecl, er, ersun, p0)
// Acceleration due to SRP (flat plate) 
//
// Calling Sequence
// [acc] =  CL_fo_srpPanelAcc(pos, pos_sun, coefp, normal, nsides [,ecl, er, ersun, p0])
// 
// Description
// <itemizedlist>
// <listitem>
// <p>Acceleration due to solar radiation pressure (SRP) for a flat plate.</p> 
// <p></p> 
// <p>The SRP coefficient <b>coefp</b> is defined by: </p>
// <p>coefp = [Ka; Kd; Ks] * a / M </p>
// <p>With: </p>
// <p>- Ka: proportion (between 0 and 1) of absorbed radiation, </p>
// <p>- Kd: proportion (between 0 and 1) of diffuse radiation. </p>
// <p>- Ks: proportion (between 0 and 1) of radiation that is reflected specularly, </p>
// <p>- a / M is the area to mass ratio (a: area of the flat plate, M: total mass of 
// the object subjected to the force). </p>
// <p>Note: Ka + Kd + Ks is equal to 1.  </p>
// <p></p>
// <p><b>nsides</b> specifies if only the side oriented by <b>normal</b> is considered or if both
// sides are. If nsides is 1, the acceleration is 0 if the normal does not point at 
// the Sun. If nsides is 2, either <b>+</b>normal or <b>-</b>normal is considered depending on which one 
// points at the Sun. </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 normal vector may not be a unit vector. </p> 
// <p>- The calculation of eclipses 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 ([]), an internal value is used. </p>
// <p>- <b>p0</b> is the solar pressure at 1 AU. If <b>p0</b> is empty ([]), an internal value is used, 
// equal to the total solar irradiance divided by the speed of light. </p>
// <p></p>
// <p>- Remark: </p>
// <p><literal>> CL_fo_srpAcc(pos, pos_sun, coefp) </literal></p>
// <p>is equivalent to: </p>
// <p><literal>> CL_fo_srpPanelAcc(pos, pos_sun, [coefp; 0; 0], pos_sun-pos, 1) </literal></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 [m]. (3xN or 3x1)
// pos_sun: Sun position vector [m]. (3xN or 3x1)
// normal: Normal vector to the flat plate (3xN or 3x1)
// coefp: SRP coefficients: [Ka; Kd; Ks] * area/mass [m^2/kg]. (3xN or 3x1)
// nsides: (integer) Number of sides to be considered: 1 or 2 (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 used). [m] (1x1)
// p0: (optional) Solar radiation pressure at 1 AU. Default is [] (internal value used). [N/m^2] (1x1)
// acc: Acceleration [m/s^2]. (3xN)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_fo_srpAcc
// CL_gm_eclipseCheck
//
// Examples
// pos = [1;0;0] * 1.e7; 
// pos_sun = CL_dataGet("au") * [1; 0; 0];  
// normal = [1; 1; 0]; 
// 
// // [Ka; Kd; Ks] *  1 m^2 / 100 kg
// // 50% absorbed, 50% specular
// coefp = [0.5; 0; 0.5] * 1.e-2; 
//
// CL_fo_srpPanelAcc(pos, pos_sun, coefp, normal, 1)
// CL_fo_srpPanelAcc(pos, pos_sun, coefp, -normal, 1)
// CL_fo_srpPanelAcc(pos, pos_sun, coefp, -normal, 2)


// 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;


// check size / resize
[pos, pos_sun, coefp, normal, nsides] = CL__checkInputs(pos, 3, pos_sun, 3, coefp, 3, normal, 3, nsides, 1); 

// check values
if (find(nsides <> 1 & nsides <> 2) <> [])
  CL__error("Invalid value for argument nsides (1 or 2 expected)"); 
end

if (find(coefp < 0) <> [])
  CL__error("Invalid value for argument coefp"); 
end

// Direction to Sun
[u,d] = CL_unitVector(pos_sun - pos); 

// Unit normal vector 
normal = CL_unitVector(normal); 

// SRP coefficients [Ka; Kd; Ks] * A/M
Ca = coefp(1,:); 
Cd = coefp(2,:); 
Cs = coefp(3,:); 

// SRP at actual distance from the Sun
K = p0 * ((AU ./ d).^2); 

// Eclipses (=> multiply by ratio between 0 and 1)
if (ecl)
  frac = 1; 
  pos_earth = [0;0;0]; 
  frac = 1 - CL_gm_eclipseCheck(pos, pos_sun, pos_earth, ersun, er);
  K = K .* frac;
end

// Normal orientation 
costheta = CL_dot(u, normal); 

// If 1 side & normal does not see the Sun => acceleration is 0
K(costheta < 0 & nsides == 1) = 0; 

// coefficients of acceleration in the "u" and "normal" directions:  
// coef_u: always negative
// coef_n: negative if costheta > 0, positive if costheta < 0 and nsides = 2
coef_u = - abs(costheta) .* K .* (Ca + Cd); 
coef_n = - costheta .* K .* ( (2. / 3.) * Cd + 2 * Cs .* abs(costheta) ); 

// acceleration
acc = CL_dMult(coef_u, u) + CL_dMult(coef_n, normal);

endfunction