//  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 [b, bdot] = CL_mod_geomagField(cjd, pos)
// Earth's magnetic field (IGRF model)
//
// Calling Sequence
// [b, bdot] = CL_mod_geomagField(cjd, pos)
//
// Description
// <itemizedlist>
// <listitem>
// <p>Computes the 3 components of the Earth's magnetic field according to IGRF model. Also returns their time derivatives. </p> 
// <p></p></listitem>
// <listitem>
// <p>Notes:</p>
// <p>- The model coefficients are given at the beginning of year (January 1st). A linear interpolation 
// is performed in order to determine the coefficients between 2 reference dates. </p>
// <p>- The implicit time scale is assumed to be UTC but TREF can be considered as acceptable. </p>
// <p>- The implicit reference frame is assumed to be ECF. </p>
// <p>- %nan is returned for dates outside the range: [1900/1/1, 2020/1/1]. The results may be
// approximate for years beyond 2015. </p>
// <p>- The model is inaccurate for radii greater than 6.6 Earth radii (around 42000 km). </p>
// <p></p></listitem>
// <listitem>
// <p>Go to http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html for more details on the model.</p>
// <p></p></listitem>
// </itemizedlist>
//
// Parameters
// cjd: Date (modified julian day from 1950.0) at which the magnetic field is computed. (1xN or 1x1)
// pos: Position at which the magnetic field is computed, in cartesian coordinates [m]. (3xN or 3x1)
// b: Magnetic field in cartesian coordinates [tesla]. (3xN)
// bdot: Magnetic field time derivative in cartesian coordinates [tesla/sec]. (3xN)
//
// Authors
// CNES - DCT/SB
//
// Examples
// cjd = CL_dat_cal2cjd(2010:2013,1,1); 
// // position - geodetic coordinates (lon,lat,alt)
// pos_ell = [0*%pi/180; 45*%pi/180; 0];
// // magnetic field in topocentric North frame - nanotesla
// b = CL_mod_geomagField(cjd, CL_co_ell2car(pos_ell)); 
// b_topo = CL_fr_topoNMat(pos_ell) * b * 1.e9


// Declarations:

// ------------------------------------------------------
// Computation for a given reference date
// cjd, pos: same as in main function (1xN) and (3xN)
// cjd_ref: date at which znm and znmd are given (1x1)
// znm: spherical harmonics  
// znmd: spherical harmonics time derivatives (per day)
// er: value of Earth radius to be used
// coef: multiplying coef for potential function to be used
// b: magnetic field (tesla)
// bdot: derivative (tesla/s)
function [b,bdot] = geomagField(cjd, pos, cjd_ref, znm, znmd, er, coef)

  // magnetic field (tesla) and derivative  
  // NB : - grad()
  gr = -CL_sphHarmGrad(pos, er, coef, znm); 
  grd = -CL_sphHarmGrad(pos, er, coef, znmd); 

  // value at cjd (linear interpolation from cjd_ref)
  b = gr + CL_dMult(grd, cjd - cjd_ref); 
  bdot = grd / 86400; // => per second

endfunction


// Code

// check argument size / resize
[cjd, pos] = CL__checkInputs(cjd,1,pos,3); 

if (cjd == [] | pos == [])
  CL__error("Invalid input arguments"); 
end

// default initialization
// (%nan is case dates out of range)
b = %nan * ones(pos); 
bdot = b; 

// load model coefficients.
// variable "mod" contains: 
// - er: Earth radius specific value (m) (1x1)
// - coef: Coef for spherical harmonics (m^2) (1x1)
// - cjd_ref: ref dates (1xN)
// - cjd_max: max date considered for validity (1x1)
// - l: list (size N):
// l(k).znm: harmonic coef (tesla) at date cjd_ref(k)
// l(k).znmd: harmonic coef derivative (tesla/day) at date cjd_ref(k)
// Note: l(k).znm and l(k).znmd may not have the same size

fpath = fullfile(CL_home(), "data", "environment", "igrf11.dat"); 
if (~isfile(fpath))
  CL__error("File not found:" + fpath); 
end
load(fpath, "mod"); 

// date intervals (mod.cjd_max = greatest valid date)
cjd_ref = [mod.cjd, mod.cjd_max]; 

// find appropriate interval for each input date
// NB1: dseach returns 0 if out of range
// NB2: dsearch(b, [a,b,c]) => returns 1 and not 2!
// "unique" is not used as supposed less efficient and dates
// should be close to one another

K = dsearch(cjd, cjd_ref);
kmin = max(min(K),1); 
kmax = min(max(K), length(cjd_ref) - 1); 
 
for (k = kmin : kmax)
  I = find(K == k);
  // computes results for date interval "k"
  if (I <> [])
    [b(:,I), bdot(:,I)] = geomagField(cjd(I), pos(:,I), ..
        mod.cjd(k), mod.l(k).znm, mod.l(k).znmd, mod.er, mod.coef); 
  end
end


endfunction