//  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 [gravmod] = CL_fo_gravityLoad(fname, incel)
// Loads a gravity model
//
// Calling Sequence
// [gravmod] = CL_fo_gravityLoad(fname, [, incel])
//
// Description
// <itemizedlist>
// <listitem>
// <p>Loads a gravity model from a file.</p> 
// <p>The function creates and returns a structure containing the data that describe the gravity 
// model (spherical harmonics). The file must conform to the syntax used for gravity models. </p> 
// <p></p>
// <p>If the <b>incel</b> ("in CelestLab") argument is true, the files are looked for in 
// the directory: CL_home()/data/potential. </p>
// <p></p></listitem>
// </itemizedlist>
//
// Parameters
// fname: (string) Name (or path) of gravity file (1x1)
// incel: (boolean, optional) %t if the file is looked for in CelestLab database. Default is %f.
// gravmod: (structure) Gravity model (spherical harmonics).
//
// Authors
// CNES - DCT/SB
//
// Examples
// gravmod = CL_fo_gravityLoad("egm96s.scd", %t)

// read potential
// creates potential structure
// unnormalize coefficients
// Cnm, Snm => gravmod.csnm(n+1,m+1)

if (~exists("incel", "local")); incel = %f; end

data = struct(); 

if (incel)
  fname = fullfile(CL_home(), "data", "potential", fname); 
end

err = exec(fname, "errcatch", -1); 
if (err <> 0)
  CL__error("Error reading potential file"); 
end

// initialize structure fields
gravmod = struct();
gravmod.name = basename(fname); 
 
gravmod.mu = data.mu; 
gravmod.er = data.er; 

maxdeg = max(max(data.coef(:,1)), max(data.coef(:,1))); 
gravmod.maxdeg = maxdeg; 

gravmod.csnm = zeros(maxdeg+1,maxdeg+1); 
  

n = data.coef(:,1); 
m = data.coef(:,2); 

// check that m <= n and n >= 0 and m >= 0
I = find(m > n | n < 0 | m < 0)
if (I <> []) 
  CL__error("Invalid values for n or m in file"); 
end

// check that Sn0 is 0
I = find(m == 0 & data.coef(:,4) <> 0);
if (I <> [])
  CL__warning(gravmod.name + " gravity model: Values Sn0 enforced to 0"); 
  data.coef(I,4) = 0; 
end  

// unnormalization factor
k = 2 * ones(m); 
I = find(m == 0); 
k(I) = 1; 
// fact = ( (n-m)! * k * (2n+1) ) / (n+m)! 
// k=1 if m == 0, k=2 otherwise

fact = sqrt(factorial(n-m) .* (2*n+1) .* k ./ factorial(n+m)); 

// unnormalized coefficients
cnm = data.coef(:,3) .* fact; 
snm = data.coef(:,4) .* fact; 

// fills matrix (complex numbers)
// Cnm,Snm => stored in gravmod.csnm(n+1,m+1)
K = sub2ind(size(gravmod.csnm), n+1, m+1); 
gravmod.csnm(K) = cnm + imult(snm); // cnm + %i*snm

endfunction