//  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_zonHarmPot(pos, nz, er, mu, j1jn)
// Potential due to gravity (zonal harmonics)   
//
// Calling Sequence
// [pot] = CL_fo_zonHarmPot(pos, [nz, er, mu, j1jn])
//
// Description
// <itemizedlist>
// <listitem>
// <p>Potential due to zonal harmonics.</p> 
// <p>By definition, the acceleration derived from the potential is +grad(potential).</p> 
// <p></p></listitem>
// <listitem>
// <p>Notes:</p>
// <p>- nz is an array of the zonal term numbers to consider.</p> 
// <p>- The origin for the position vector must be the central body.</p>
// <p>- The coordinates frame can be any frame whose Z axis is the revolution axis of the perturbation field.</p> 
// <p>- If nz contains numbers in the range [2,6], explicit formulas are used (faster). </p>
// <p>- If nz is equal to -1 or "all", all the harmonics are used. </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 the central body [m]. (3xN)
// nz: (optional, integer) Zonal term numbers to considerer (1xP) or "all" or -1
// er: (optional) Equatorial radius [m]. Default value is %CL_eqRad.
// mu: (optional) Gravitational constant [m^3/s^2]. Default value is %CL_mu.
// j1jn: (optional) Zonal harmonics. Default value is %CL_j1jn.
// pot: Potential [m^2/s^2]. (1xN)  
//
// Authors
// CNES - DCT/SB
//
// Bibliography
// D. Vallado - Fundamentals of Astrodynamics and Applications, 2nd edition
//
// See also
// CL_fo_zonHarmAcc
// CL_sphHarmVal
//
// Examples
// pos = [[7000.e3; 0; 0], [0; 7000.e3; 0], [0; 0; 7000.e3]];
// CL_fo_zonHarmPot(pos) // all harmonics
// CL_fo_zonHarmPot(pos, 2:6)
// CL_fo_zonHarmPot(pos, [2,4,6])

// Declarations:

// ----------------------------------------
// Internal function 1 (J2-J6 - analytical formulas)
// ----------------------------------------
function [pot] = zonHarmPot1(pos, nz, er, mu, j1jn)

  // Jn: only those considered are <> 0
  Jn = zeros(6,1);
  Jn(nz) = j1jn(nz); 

  r = CL_norm(pos); 
  z = pos(3,:); // "z" coordinate 
  zr = z./r; 
  
  // Note: 
  // pot = grad(acceleration)
  // for all terms: 
  // pot = - Jk * (mu/r) * (R/r)^k * Pk(z/r)
  // Pk: Legendre polynomial of degree k

  pot = zeros(size(pos,2)); 

  // effect of J2
  if (Jn(2) <> 0)
    pot = pot - (1.0 / 2) * Jn(2) * (mu./r) .* (er./r).^2 .* (3*zr.^2 - 1);
  end

  // effect of J3
  if (Jn(3) <> 0)
    pot = pot - (1.0 / 2) * Jn(3) * (mu./r) .* (er./r).^3 .* (5*zr.^3 - 3*zr); 
  end

  // effect of J4
  if (Jn(4) <> 0)
    pot = pot - (1.0 / 8) * Jn(4) * (mu./r) .* (er./r).^4 .* (35*zr.^4 - 30*zr.^2 + 3); 
  end

  // effect of J5
  if (Jn(5) <> 0)
    pot = pot - (1.0 / 8) * Jn(5) * (mu./r) .* (er./r).^5 .* (63*zr.^5 - 70*zr.^3 + 15*zr); 
  end

  // effect of J6
  if (Jn(6) <> 0)
    pot = pot - (1.0 / 16) * Jn(6) * (mu./r) .* (er./r).^6 .* (231*zr.^6 - 315*zr.^4 + 105*zr.^2 - 5); 
  end

endfunction


// ----------------------------------------
// Internal function 2 (all cases)
// ----------------------------------------
function [pot] = zonHarmPot2(pos, nz, er, mu, j1jn)

  // Cn0 = -Jn: only those considered are <> 0
  // Warning: Cn0 = znm(n+1)
  nmax = max(nz); 

  znm = zeros(nmax+1,1);
  znm(nz+1) = -j1jn(nz); 

  pot = CL_sphHarmVal(pos, er, mu, znm, [nmax,0]); 

endfunction


// Code:
if (~exists('nz','local')); nz = "all"; end
if (~exists("er", "local")); er = CL__dataGetEnv("eqRad"); end
if (~exists("mu", "local")); mu = CL__dataGetEnv("mu"); end
if (~exists("j1jn", "local")); j1jn = CL__dataGetEnv("j1jn"); end

// if nz == "all" or -1, all coefs are used
if (nz == "all" | nz == -1); nz = 1:length(j1jn); end

I = find(nz <> round(nz) | nz < 1 | nz > length(j1jn));
if (I <> [])
  CL__error("Invalid zonal terms numbers");
end

// check sizes as algorithms could not work
// nz: row vector
if (size(nz,1) <> 1)
  CL__error("Wrong size for the number of coefficients");
end

// j1jn: column vector
if (size(j1jn,2) <> 1)
  CL__error("Wrong size for the zonal coefficients");
end


nzmin = min(nz);
nzmax = max(nz);

if (nzmin >= 2 & nzmax <= 6)
  pot = zonHarmPot1(pos, nz, er, mu, j1jn); 
else
  pot = zonHarmPot2(pos, nz, er, mu, j1jn); 
end

endfunction