//  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 [pos_planet_EOD,vel_planet_EOD] = CL__eph_planet_VSOP87(planet,jd, nbcoef,cvel)
// Planets positions in EOD frame (VSOP87 theory)
//
// Calling Sequence
// [pos_planet_EOD,vel_planet_EOD] = CL__eph_planet_VSOP87(planet,jd [,nbcoef, cvel])
//
// Description
// <itemizedlist><listitem>
// <p>Computes various planets positions and velocity relative to the ecliptic of date frame (EOD) 
// using VSOP87 theory.</p>
// <p>The center of the frame is the Sun.</p>
// <p>Available planets are Mercury,Venus,Earth,Mars,Jupiter,Saturn,Uranus,Neptune and Pluto</p>
// <p></p>
// <p>Each variable (longitude, latitude and distance) are computed using :
// L0 + L1*T + L2*T^2 + L3*T^3 + L4*T^4 + L5*T^5 and each Ln is the sum of harmonics terms</p>
// <p>The number of harmonics terms to be used for each series can be tweaked using nbcoef. (default is maximum)</p>
// <p>It is a (3x6) matrix with the row being the variable (longitude,latitude,distance) and the 
// column being the power of T (from 0 to 5)</p>
// <p>For example :</p>
// <p>- nbcoef(1,1) is the number of harmonics to be used in the computation of longitude, for T^0</p>
// <p>- nbcoef(3,6) is the number of harmonics to be used in the computation of distance, for T^5</p>
// </listitem>
// </itemizedlist>
//
// Parameters
// planet: (string) Name of planet ("Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus" or "Neptune") (1x1)
// jd : Julian date (TT) (1xN)
// nbcoef : (integer, optional) Number of harmonics terms to be used. Default is the maximum. (3x6)
// cvel: (optional, boolean) %t if velocity should be computed. Default is %t. (1x1)
// pos_planet_EOD : Position relative to EOD frame [m] (3xN)
// vel_planet_EOD : Velocity relative to EOD frame [m/s] (3xN)
//
// Bibliography
// Planetary Theories in rectangular and spherical variables: VSOP87 solution. Bretagnon P., Francou G.
//
// Authors
// CNES - DCT/SB
//
// Examples
// // Julian date, in TT time scale
// jd = 21512210.3;
//
// // Default = highest precision
// pos_mars_EOD = CL__eph_planet_VSOP87("Mars", jd)
//
// // Lower precision : only 10 harmonics for each serie
// nbcoef = 10 * ones(3,6);
// pos_mars_EOD = CL__eph_planet_VSOP87("Mars", jd, nbcoef);


  function series = load_vsop87(planet)
  // Loads the series of a planet from a file
  // 
  // planet: (string) Name of planet: Mercury,Venus,Earth,Mars,Jupiter,Saturn,Uranus or Neptune
  // series: structure of series:
  // contains 3 sub-structures ("L", "B", "R") each containing a list of coefficients for powers of T (T^0 up to T^5)
  // Each sub-structure is a Px3 matrix. First column = A, Second column = B, Third column = C
  //
  // Note on units :
  // tjy: time expressed in Thousands of Julian Years (tjy)    
  // alpha : degree of time in-between 0 and 5.
  // A is expressed in au/(tjy^alpha) for the distances and in rad/(tjy^alpha) for angular variables.
  // The phase B is expressed in radians.
  // The frequency C is expressed in rad/tjy.
  //
  // See VSOP87 theory for more details on the meaning of these terms.

    fname = planet + ".dat";
    fpath = fullfile(CL_home(), "data", "ephem", "VSOP87_series", fname);
    
    if (~isfile(fpath)); CL__error("File " + fpath + "does not exist"); end;
    
    // Loads a variable named "series"
    load(fpath, "series");
    
  endfunction


  function [pos_sum, vel_sum] = sum_series(series,T,TT,nbcoef)
  // Computes L0 + L1*T + L2*T^2 + L3*T^3 + L4*T^4 + L5*T^5
  // With each Ln being the sum of the series : An*cos(Bn+CnT)*
  // Also computes its derivative with respect to time
  //
  // series: structure containing matrixes ABC (Px3) for powers of T: T^0, T^1 ... T^5
  // T: Time in millenia from epoch 2000 12h (TT time scale)
  // TT: T^0; T^1; ... T^5 (5xN)
  // nbcoef: number of terms to be used in each series (1x6)
  // pos_sum: L0 + L1*T + L2*T^2 + L3*T^3 + L4*T^4 + L5*T^5 (1xN)
  // vel_sum: derivative of pos_sum with respect to T (1xN)
  //          = L0p + L1+L1p*T + 2*L2*T+L2pT^2 + 3*L3*T^2+L3pT^3 + 4*L4*T^3+L4pT^4 + 5*L5*T^4+L5pT^5
  
  pos_sum = zeros(T);
  vel_sum = zeros(T);
  
  // Compute series L and Lp
  for k = 1 : 6
    ABC = series(k);
    nbcoefs = min( nbcoef(k) , size(ABC,1) ); // nbcoef should not be superior to total number of series
    
    if (nbcoefs > 0)
      A = ABC(1:nbcoefs,1) ;
      B = ABC(1:nbcoefs,2) ;
      C = ABC(1:nbcoefs,3) ;
      Lk = sum (A*ones(T) .* cos(B*ones(T) + C * T ) , "r");
      pos_sum = pos_sum + Lk .* TT(k,:);
      
      // If velocity required
      if (argn(1) == 2)
        Lpk = sum ((A.*C)*ones(T) .* -sin(B*ones(T) + C * T ) , "r"); 
        if (k==1)
          vel_sum = Lpk;
        else
          vel_sum = vel_sum + (k-1) * Lk .* TT(k-1,:)  +  Lpk .* TT(k,:);
        end
      end
    end
  end
  
  endfunction


// Declarations:
AU = CL__dataGetEnv("au", internal=%t); 

// Code:
if (argn(2) < 2 | argn(2) > 4) 
  CL__error("Invalid number of input arguments");
end

if (~exists('nbcoef','local')); nbcoef = %inf * ones(3,6); end // Default = maximum number of harmonics
if (~exists("cvel","local")); cvel = %t; end;

// Compute velocity only if needed
if (argn(1) == 1)
  cvel = %f;
end

// Time in millenia from epoch 2000 12h (TT time scale)
T   = (jd - 2451545.0) / 365250.0;
TT = [ ones(T) ; T ; T.^2 ; T.^3 ; T.^4 ; T.^5 ];

// Loads series of the planet
series = load_vsop87(planet);

// Compute series for longitude, latitude, distance, with various number of harmonics
if (~cvel)
  lon = sum_series(series("L"),T,TT,nbcoef(1,:)); 
  lat = sum_series(series("B"),T,TT,nbcoef(2,:));
  dist = sum_series(series("R"),T,TT,nbcoef(3,:));
  
  pos_planet_EOD = CL_co_sph2car([lon;lat;dist]);
  pos_planet_EOD = pos_planet_EOD * AU; // convert AU to m
  vel_planet_EOD = [];
  
// Compute with velocity
else
  [lon,lonp] = sum_series(series("L"),T,TT,nbcoef(1,:)); 
  [lat,latp] = sum_series(series("B"),T,TT,nbcoef(2,:));
  [dist,distp] = sum_series(series("R"),T,TT,nbcoef(3,:));
  
  [pos_planet_EOD,vel_planet_EOD] = CL_co_sph2car([lon;lat;dist] , [lonp;latp;distp]);
  pos_planet_EOD = pos_planet_EOD * AU; // convert AU to m
  vel_planet_EOD = vel_planet_EOD * AU / (365250 * 86400); // convert AU/millenia to m/s
end

endfunction