//  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_solidTidesAcc(pos, pos_b, mu_b, deg, er, lcoefs)
// Acceleration due to gravity (solid tides)  
//
// Calling Sequence
// [acc] = CL_fo_solidTidesAcc(pos, pos_b, mu_b [, deg, er, lcoefs])
//
// Description
// <itemizedlist>
// <listitem>
// <p>Acceleration due to solid tides caused by a perturbing body.</p>
// <p></p>
// <p>The computation is based on the Love model. 
// The selection of the degree is done via the input parameter <b>deg</b>.</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. </p>
// <p>- Explicit formulas are used for degrees 2 and 3. </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 central body) [m]. (3xN or 3x1)
// pos_b: Position vector of perturbing body (from central body) [m]. (3xN or 3x1)
// mu_b: Gravitational constant of perturbing body [m^3/s^2]. (1x1)
// deg: (integer, optional) Expansion degree of the Love model. Default is -1 which means "biggest possible"
// er: (optional) Equatorial radius of central body [m]. Default value is %CL_eqRad. (1x1)
// lcoefs: (optional) Tidal Love coefficients: [k1(=0); k2; ...]. Default value is %CL_tidalLoveCoefs. (Px1)
// acc: Acceleration [m/s^2]. (3xN) 
//
// Authors
// CNES - DCT/SB
//
// Examples
// pos = [7000.e3; 0; 0]; // ECI
// pos_moon = CL_eph_moon(CL_dat_cal2cjd(2000,3,21)); // ECI
// mu_moon = CL_dataGet("body.Moon.mu");
// CL_fo_solidTidesAcc(pos, pos_moon, mu_moon)

// Declarations:

// ----------------------------------------
// Internal function 1
// For degree 2 or 3 only - explicit formulas
// NB: pos, pos_b: same size
// ----------------------------------------
function [acc] = solidTidesAcc1(pos, pos_b, mu_b, deg, er, lcoefs)

  // unit vectors + norms
  [upos, d] = CL_unitVector(pos); 
  [upos_b, d_b] = CL_unitVector(pos_b); 

  // cos(angle between satellite and perturbing body)
  cosalpha = CL_dot(upos, upos_b); 
    
  // Acceleration due to solid tides
  // k = love Factor
  // acc_tide1 = acc magnitude in the direction of the satellite due to Perturbing body
  // acc_tide2 = acc magnitude in the direction of the Perturbing body

  acc_tide1 = zeros(cosalpha); 
  acc_tide2 = zeros(cosalpha); 

  if (deg >= 2)
    k2 = lcoefs(2);
    K = mu_b * 1.5 * k2 * (er^5) ./ ((d.^4) .* (d_b.^3));
    acc_tide1 = acc_tide1 + K .* (1 - 5 * cosalpha.^2);
    acc_tide2 = acc_tide2 + K .* (2 * cosalpha);
  end
  
  if (deg >= 3)
    k3 = lcoefs(3);
    K = mu_b * 0.5 * k3 * (er^7) ./ ((d.^5) .* (d_b.^4));
    acc_tide1 = acc_tide1 + K .* (15 - 35 * cosalpha.^2) .* cosalpha;
    acc_tide2 = acc_tide2 + K .* (15 * cosalpha.^2 - 3);
  end

  // acceleration vector 
  acc = CL_dMult(acc_tide1, upos) + CL_dMult(acc_tide2, upos_b);

endfunction

// ----------------------------------------
// Internal function 2
// For any degree - use of spherical harmonics. 
// Potential can be written: 
// pot = (mub*R/rb) * (1/r) * sum(kn * ((R^2/rb)/r)^n Pn(cos(ang))
// -> if frame is changed, cos(ang) can become sin(lat)
// -> kn = Love coef of degree n
// NB: pos, pos_b: same size
// ----------------------------------------
function [acc] = solidTidesAcc2(pos, pos_b, mu_b, deg, er, lcoefs)

  rb = CL_norm(pos_b); 
  a = (er * er) ./ rb; 
  f = mu_b * er ./ rb; 

  // new frame: body -> z axis, other axes: doesn't matter
  // => Pn(cos(ang)) in potential expansion = Pn(sin(lat)); 
  M = CL_rot_defFrameVec(pos_b, [0;0;1], 3, 1); 
  
  // coefficients: index 1 <=> order 0
  znm = [0; lcoefs]; 

  // acceleration using spherical harmonics
  // computation in new frame. Result converted to initial frame
  acc = (M') * CL_sphHarmGrad(M*pos, a, f, znm, [deg, 0], inc00=%f); 

endfunction
  

// Code

// Check inputs
if (~exists("deg", "local")); deg = -1; end
if (~exists("er", "local")); er = CL__dataGetEnv("eqRad"); end
if (~exists("lcoefs", "local")); lcoefs = CL__dataGetEnv("tidalLoveCoefs"); end

if (deg <> round(deg) | (deg <= 0 & deg <> -1))
  CL__error("Invalid value for argument deg");
end

if (size(lcoefs, 2) <> 1 | size(lcoefs, 1) < deg)
  CL__error("Invalid size for argument lcoefs");
end

// check size / resize
[pos, pos_b] = CL__checkInputs(pos, 3, pos_b, 3); 

// special case for deg: "all"
if (deg == -1)
  deg = size(lcoefs, 1); 
end


// computation (using explict formulas whenever possible)
if (deg >= 2 & deg <= 3)
  acc = solidTidesAcc1(pos, pos_b, mu_b, deg, er, lcoefs); 
else
  acc = solidTidesAcc2(pos, pos_b, mu_b, deg, er, lcoefs); 
end


endfunction