//  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_zonHarmAcc(pos, nz, er, mu, j1jn)
// Acceleration due to gravity (zonal harmonics)
//
// Calling Sequence
// [acc] = CL_fo_zonHarmAcc(pos, [nz, er, mu, j1jn])
//
// Description
// <itemizedlist>
// <listitem>
// <p>Acceleration due to zonal harmonics.</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. (1x1)
// mu: (optional) Gravitational constant [m^3/s^2]. Default value is %CL_mu. (1x1)
// j1jn: (optional) Zonal harmonics. Default value is %CL_j1jn.
// acc: Acceleration [m/s^2]. (3xN)
//
// Authors
// CNES - DCT/SB
//
// Bibliography
// D. Vallado - Fundamentals of Astrodynamics and Applications, 2nd edition
//
// See also
// CL_fo_zonHarmPot
// CL_sphHarmGrad
//
// Examples
// pos = [[7000.e3; 0; 0], [0; 7000.e3; 0], [0; 0; 7000.e3]];
// CL_fo_zonHarmAcc(pos, 2)
// CL_fo_zonHarmAcc(pos, 2:6)
// CL_fo_zonHarmAcc(pos, [2,4,6])

// Declarations:

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

  r = CL_norm(pos); 
  zr = pos(3,:) ./ r; 
  zr2 = zr .* zr; 
  zr4 = zr2 .* zr2; 
  zr6 = zr4 .* zr2; 
  
  // Jn: only those considered are <> 0
  Jn = zeros(6,1);
  Jn(nz) = j1jn(nz); 

  acc = zeros(pos); 

  // effect of J2
  if (Jn(2) <> 0)
    K = (-3.0 / 2) * mu * Jn(2) * er^2 ./ r.^5; 
    Kxy = K .* (1 - 5 * zr2); 
    Kz = K .* (3 - 5 * zr2); 
    acc = acc + [ Kxy .* pos(1,:); Kxy .* pos(2,:); Kz .* pos(3,:) ]; 
  end

  // effect of J3
  if (Jn(3) <> 0)
    K = (-1.0 / 2) * mu * Jn(3) * er^3 ./ r.^5; 
    Kxy = K .* zr .* (15 - 35 * zr2) ./ r; 
    Kz = K .* (-3 + 30 * zr2 - 35 * zr4); 
    acc = acc + [ Kxy .* pos(1,:); Kxy .* pos(2,:); Kz ]; 
  end

  // effect of J4
  if (Jn(4) <> 0)
    K = (-5.0 / 8) * mu * Jn(4) * er^4 ./ r.^7;
    Kxy =  K .* (-3 + 42 * zr2 - 63 * zr4);  
    Kz = K .* (-15 + 70 * zr2 - 63 * zr4);
    acc = acc + [ Kxy .* pos(1,:); Kxy .* pos(2,:); Kz .* pos(3,:) ]; 
  end

  // effect of J5
  if (Jn(5) <> 0)
    K = (3.0 / 8) * mu * Jn(5) * er^5 ./ r.^7;
    Kxy =  K .* zr .* (35 - 210 * zr2 + 231 * zr4) ./ r;  
    Kz = K .* (-5 + 105 * zr2 - 315 * zr4 + 231 * zr6);
    acc = acc + [ Kxy .* pos(1,:); Kxy .* pos(2,:); Kz ]; 
  end

  // effect of J6
  if (Jn(6) <> 0)
    K = (-1.0 / 16) * mu * Jn(6) * er^6 ./ r.^9;
    Kxy =  K .* (35 - 945 * zr2 + 3465 * zr4 - 3003 * zr6);  
    Kz = K .* (245 - 2205 * zr2 + 4851 * zr4 - 3003 * zr6);
    acc = acc + [ Kxy .* pos(1,:); Kxy .* pos(2,:); Kz .* pos(3,:)]; 
  end

endfunction


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

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

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

  acc = CL_sphHarmGrad(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)
  acc = zonHarmAcc1(pos, nz, er, mu, j1jn); 
else
  acc = zonHarmAcc2(pos, nz, er, mu, j1jn); 
end

endfunction