//  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_albedoAcc(pos, pos_sun, coefp, albedo, er, ersun, p0)
// Acceleration due to albedo (reflected solar pressure)
//
// Calling Sequence
// [acc] = CL_fo_albedoAcc(pos, pos_sun, coefp [, albedo, er, ersun, p0])
// 
// Description
// <itemizedlist>
// <listitem>
// <p>Acceleration due to constant albedo (reflected solar pressure).</p>
// <p></p></listitem>
// <listitem>
// <p>Notes: </p>
// <p>- The coordinates frame can be any frame.</p>
// <p>- The origin of the frame must be the central body, which is also the reflecting body.</p>
// <p>- The spacecraft is assumed spherical. Consequently, the acceleration vector lies in the plane 
// formed by the central body, the perturbing body and the Sun. </p>
// <p>- The function uses pre-computed data present in the file "acc_albedo.dat" located in the 
// CL_home()/data/force_models directory. </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 (same frame as pos) [m]. (3xN or 3x1)
// coefp: SRP coefficient (cp*area/mass) [m^2/kg]. (1x1 or 1xN)
// albedo: (optional) Albedo coefficient. Default value is %CL_albedo. (1x1 or 1xN)
// er: (optional) Equatorial radius of reflecting body. Default is %CL_eqRad. [m] (1x1)
// ersun: (optional) Sun radius. 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
//
// 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_albedoAcc(pos, pos_sun, coefp)
// 

// Declarations:

// Check optional inputs
if (~exists("albedo", "local")); albedo = CL__dataGetEnv("albedo"); 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;


// File acc_albedo.dat contains structure s (ang, q, val_x, val_z)
// ang: angle between body and spacecraft from planet center
// q : half view angle from planet center
// val_z: value in radial direction (// to planet-center->spacecraft) 
// => val_z should be positive
// val_x: value in perpendicular direction, in the Sun-Spacecraft plane, toward the Sun
// => val_x should be negative

fname = fullfile(CL_home(), "data", "force_models", "acc_albedo.dat"); 
if (~isfile(fname))
  CL__error("File " + fname + " not found");
end
load(fname);
  
// Check input size
[pos, pos_sun, coefp, albedo] = CL__checkInputs(pos, 3, pos_sun, 3, coefp, 1, albedo, 1); 

// SRP (no eclipses) - norm of acceleration 
// sat at center of central body (= [0;0;0])
acc_srp_norm = CL_norm(CL_fo_srpAcc([0;0;0], pos_sun, coefp, %f, er=er, ersun=ersun, p0=p0)); 
  
// Interpolation spline
CZ = splin2d(s.q, s.ang, s.val_z); 
CX = splin2d(s.q, s.ang, s.val_x); 

// q = Center angle for visibility considering spherical planet
q = real(acos(er ./ CL_norm(pos))); 
  
// ang = angle between the satellite direction and the Sun direction
ang = CL_vectAngle(pos, pos_sun); 
  
// Interpolated values
val_z = interp2d(q, ang, s.q, s.ang, CZ);
val_x = interp2d(q, ang, s.q, s.ang, CX);

// Check sign in case of interpolation errors: 
// val_z >= 0 and val_x <= 0
val_z = max(val_z, 0); 
val_x = min(val_x, 0); 

// New frame (Z => sat, X => Sun)
M = CL_rot_defFrameVec(pos, pos_sun, 3, 1); 

// The looked-for acceleration is the pre-computed acceleration (in the correct frame)
// multiplied by |SRP| * albedo
acc = CL_dMult(acc_srp_norm .* albedo, M' * [val_x; zeros(val_x); val_z]); 
  
endfunction