//  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_sphHarmPot(pos, maxnm, er, mu, j1jn, cs1nm, inc00)
// Potential due to gravity (spherical harmonics)   
//
// Calling Sequence
// [pot] = CL_fo_sphHarmPot(pos, maxnm [, er, mu, j1jn, cs1nm, inc00])
//
// Description
// <itemizedlist>
// <listitem>
// <p>Potential defined by spherical harmonics.</p> 
// <p>By definition, the acceleration derived from the potential is +grad(potential).</p> 
// <p></p></listitem>
// <listitem>
// <p>Notes:</p>
// <p>- The origin for the position vector must be the central body.</p>
// <p>- The coordinates frame must be the frame in which the potential is defined.</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)
// maxnm : (integer) Maximum degree/order to consider: [nmax,mmax] (1x2) or "all"
// 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. (Nzx1)
// cs1nm : (complex, optional) Tesseral and sectorial harmonics: cs1nm(n,m) = Cnm + %i * Snm. Default value is %CL_cs1nm (NnxNm)
// inc00: (optional) %t if the effect of the central term is included. Default is %f. 
// pot: Potential [m^2/s^2]. (1xN)  
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_fo_sphHarmAcc
// CL_sphHarmVal
//
// Examples
// pos = [[7000.e3; 0; 0], [0; 7000.e3; 0], [0; 0; 7000.e3]];
// CL_fo_sphHarmPot(pos)
// CL_fo_sphHarmPot(pos, [6,6])

// Code: 

if (~exists("maxnm", "local")); maxnm = "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 (~exists("cs1nm", "local")); cs1nm = CL__dataGetEnv("cs1nm"); end
if (~exists("inc00", "local")); inc00 = %f; end

// j1jn: column vector
if (size(j1jn,2) <> 1)
  CL__error("Invalid size for j1jn");
end

// maxnm == "all" => maximum sizes
if (maxnm == "all")
  maxnm = max(size(j1jn), size(cs1nm))
end

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

I = find(maxnm <> round(maxnm) | maxnm < 0);
if (I <> [])
  CL__error("Invalid values for maxnm");
end

nmax = maxnm(1); 
mmax = maxnm(2); 
n = min(size(j1jn,1), nmax);
ni = min(size(cs1nm,1), nmax); 
nj = min(size(cs1nm,2), mmax); 

znm = zeros(nmax+1, mmax+1); 
znm(1,1) = 1; 
znm(2:n+1,1) = -j1jn(1:n); 
znm(2:ni+1,2:nj+1) = cs1nm(1:ni,1:nj); 

pot = CL_sphHarmVal(pos, er, mu, znm, maxnm, inc00); 

endfunction