//  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_solidTidesPot(pos, pos_b, mu_b, deg, er, lcoefs)
// Potential due to gravity (solid tides) 
//
// Calling Sequence
// [pot] = CL_fo_solidTidesPot(pos, pos_b, mu_b [, deg, er, lcoefs])
//
// Description
// <itemizedlist>
// <listitem>
// <p>Potential due to solid tides caused by a perturbing body.</p>
// <p>By definition, the acceleration derived from the potential is +grad(potential). </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)
// pot: Potential [m^2/s^2]. (1xN)
//
// 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_solidTidesPot(pos, pos_moon, mu_moon)

// Declarations:

// ----------------------------------------
// Internal function 1
// For degree 2 or 3 only - explicit formulas
// NB: pos, pos_b: same size
// ----------------------------------------
function [pot] = solidTidesPot1(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); 
    
    
  // Potential due to solid tides
  // k = love Factors (k2 <=> degree 2)

  pot = zeros(cosalpha); 

  if (deg >= 2)
    k2 = lcoefs(2);
    pot = pot + 0.5 * mu_b * k2 * (er^5) .* (3 * cosalpha.^2 - 1) ./ ((d.^3) .* (d_b.^3));
  end

  if (deg >= 3)
    k3 = lcoefs(3);
    pot = pot + 0.5 * mu_b * k3 * (er^7) .* ((5 * cosalpha.^2 - 3) .* cosalpha) ./ ((d.^4) .* (d_b.^4));
  end
  
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 [pot] = solidTidesPot2(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]; 

  // potential using spherical harmonics
  pot = CL_sphHarmVal(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)
  pot = solidTidesPot1(pos, pos_b, mu_b, deg, er, lcoefs); 
else
  pot = solidTidesPot2(pos, pos_b, mu_b, deg, er, lcoefs); 
end

endfunction