//  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, vel, acc] = CL_eph_de405(body,cjd, orig,tt_tref,ephem); 
// Position, velocity and acceleration of a body using JPL's DE405 ephemerides
//
// Calling Sequence
// [pos, vel, acc] = CL_eph_de405(body,cjd [,orig,tt_tref,ephem])
//
// Description
// <itemizedlist><listitem>
// <p>Computes the position, velocity and acceleration of a body using JPL's DE405 ephemerides.</p>
// <p>The results are given in the ICRS frame.</p>
// <p>Available bodies are: "Mercury", "Venus", "Earth-Moon-bary" (or "EMB"), "Mars", "Jupiter", "Saturn", 
//      "Uranus", "Neptune", "Pluto", "Moon", "Sun", "Earth", "solar-sys-bary" (or "SSB")</p>
// <p></p>
// <p>The computed data are:</p>
// <p>pos: vector from <b>orig</b> to <b>body</b></p>
// <p>vel: first time derivative of pos</p>
// <p>acc: second time derivative of pos</p>
// <p></p>
// <p>By default, the function looks for the appropriate JPL's DE405 ephemeris files located 
// in CL_home()/data/ephem/de405.</p>
// <p>One can also specify the ephemeris to be used with the optional argument <b>ephem</b>. 
// (See function <link linkend="CL_eph_de405Load">CL_eph_de405Load</link>)</p>
// <p></p></listitem>
// <listitem>
// <p>Notes:</p>
// <p> - For dates that are not covered by the ephemeris file(s), %nan values are returned.</p>
// <p> - The default value for the origin is the solar system barycenter (even for the Moon!)</p>
// <p> - Theoretically the time scale for DE405 is TDB. In practice, using a TT time scale is acceptable.
// See <link linkend="Dates and time scales">Dates and time scales</link> for more details.</p>
// <p> - The results returned for all planets except Earth actually are the position velocity and acceleration of the 
// barycenter of the planet system (i.e barycenter of planet + satellites).</p>
// <p> - Body names are case sensitive.</p>
// <p></p></listitem>
// <listitem>
// <p>See <link linkend="Reference frames">Reference frames</link> for more details on the ICRS frame,
// and <link linkend="Ephemerides">Ephemerides</link> for more details on ephemerides of celestial bodies.</p>
// <p></p></listitem>
// </itemizedlist>
//
// Parameters
// body: (string) Name of the body. (1x1)
// cjd: Modified (1950.0) julian day (Time scale: TREF) (1xN)
// orig: (optional, string) Name of the origin. Default is "solar-sys-bary" (1x1)
// tt_tref: (optional) TT-TREF [seconds]. Default is %CL_TT_TREF. (1xN or 1x1)
// ephem: (optional) Scilab structure of ephemeris (see function <link linkend="CL_eph_de405Load">CL_eph_de405Load</link>).
// pos: Position of body, relative to orig, in ICRS frame (3xN)
// vel: Velocity of body, relative to orig, in ICRS frame (3xN)
// acc: Acceleration of body, relative to orig, in ICRS frame (3xN)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_eph_de405Load
//
// Examples
// // Position of Mars in ICRS (Origin = Earth)
// cjd = 21915:21920;
// pos = CL_eph_de405("Mars", cjd, "Earth");
//
// // Position, velocity, acceleration of Jupiter in ICRS (Origin = Solar system barycenter)
// cjd = 21915:21920;
// [pos, vel, acc] = CL_eph_de405("Jupiter", cjd);
//

// Internal function to load the descriptor of ephemeris files (located in data folder)
function [desc] = de405_loadDesc(fpath)
    if (~isfile(fpath))
      CL__error("File " + fpath + " not found");
    end
    load(fpath, "desc");
endfunction

// Declarations

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

// Default origin
if (~exists("orig", "local")); orig = ""; end
if (orig == "" | orig == [])
  orig = "solar-sys-bary"; 
end

if (~exists("ephem", "local")); ephem = []; end
if (~exists("tt_tref", "local")); tt_tref = CL__dataGetEnv("TT_TREF"); end

if (typeof(body) <> "string" | typeof(orig) <> "string")
  CL__error("Invalid input arguments");
end

// Check sizes
CL__checkInputs(cjd,1, tt_tref,1);
N = size(cjd,2);

// Replace "SSB" and "EMB" by "solar-sys-bary" and "Earth-Moon-bary"
if (body == "SSB"); body = "solar-sys-bary"; end;
if (body == "EMB"); body = "Earth-Moon-bary"; end;
if (orig == "SSB"); orig = "solar-sys-bary"; end;
if (orig == "EMB"); orig = "Earth-Moon-bary"; end;

// Flag to compute velocity, acceleration only if requested
cvel = (argn(1) >= 2);
cacc = (argn(1) >= 3);

// Default initialization
pos = %nan * ones(3,N); 
vel = []; 
acc = []; 

if (cvel); vel = %nan * ones(3,N); end
if (cacc); acc = %nan * ones(3,N); end

// stops if cjd is empty (=> returns [])
if (N == 0)
  return; // <= RETURN
end

// Convert from TREF to TDB
// (Actually to TT: small approximation)
cjd = cjd + tt_tref/86400;


if (ephem == [])
  // If ephem not provided, load descriptor and load 
  // appropriate binary files

  fpath = fullfile(CL_home(), "data", "ephem", "de405", "desc.dat");
  desc = de405_loadDesc(fpath); 

  interv = [desc.cjd_begin, desc.cjd_end($)]; 

  // Search dates in intervals
  ind = dsearch(cjd, interv); 

  I = find(ind > 0); 
  indmin = min(ind(I)); 
  indmax = max(ind(I)); 

  for k = indmin : indmax

    I = find(ind == k); 
    if (I == []); continue; end

    // Load ephem from binary file
    fpath = fullfile(CL_home(), "data", "ephem", "de405", desc.fname(k));
    ephem = CL_eph_de405Load(fpath); 
    
    [pos1, vel1, acc1] = de405_getPVA(body, cjd(I), orig, ephem, cvel, cacc);
    pos(:,I) = pos1; 
    vel(:,I) = vel1; 
    acc(:,I) = acc1; 
  end

elseif (ephem.ephem_type == 1)
  // simple ephemeris structure is provided (type 1)
  [pos, vel, acc] = de405_getPVA(body, cjd, orig, ephem, cvel, cacc);

else
  // multiple ephemeris structure is provided (type 2)

  // date intervals in descriptor
  interv = [ephem.desc.cjd_begin, ephem.desc.cjd_end($)]; 

  // Search dates in intervals
  ind = dsearch(cjd, interv); 

  I = find(ind > 0); 
  index_desc = []; 
  if (I <> []); index_desc = min(ind(I)) : max(ind(I)); end
  
  for k = index_desc
    I = find(ind == k); 
    if (I == []); continue; end
    [pos1, vel1, acc1] = de405_getPVA(body, cjd(I), orig, ephem.eph(k), cvel, cacc);
    pos(:,I) = pos1; 
    vel(:,I) = vel1; 
    acc(:,I) = acc1; 
  end
  
end

endfunction


// ephem: jpl ephem structure
// numbody: index between 1 and 11  
// t: cjd (tbd time scale)
function [pos, vel, acc] = de405_interp(ephem, numbody, t, cvel, cacc)

N = size(t, 2); 

// NB: vel, acc initialized to [] if cvel or cacc = %f
pos = %nan * ones(3,N); 
vel = []; 
acc = []; 
if (cvel); vel = %nan * ones(3,N); end
if (cacc); acc = %nan * ones(3,N); end

tbeg = ephem.data(1,1); 
tend = ephem.data(2,$); 

// number and length (days) of intervals - based on data only
nbint = size(ephem.data,2); 
Tint = (tend-tbeg)/nbint; 

// number of sub-intervals per interval
nbsub = ephem.ptr(3, numbody); 

// total number and length(days) of sub-intervals
T = Tint / nbsub; 
nb = nbint * nbsub; 

// number of chebyshev coefs per sub-interval (x, y or z)
nbcoef = ephem.ptr(2, numbody);  

// check dates are in correct range
Iok = find(t >= tbeg & t <= tend); 
t = t(Iok); 

// all Chebyshev coef for body (x,y,z,x,y,z...)
// 1 column <=> 1 sub-interval 
i1 = ephem.ptr(1, numbody); // 1st index in each interval
i2 = i1 + 3 * nbsub * nbcoef - 1;  // last index in each interval
allcoefs = matrix(ephem.data(i1:i2, :), nbcoef, -1); 

// sub-interval indices for each date
It = max(1, min(nb, floor((t - tbeg) / T) + 1));

// normalized date: in (-1, 1) in each sub-interval
tn = 2 * (t - tbeg - (It-1) * T) / T - 1; 

// values of Chebyshev polynomials for each normalized date
// p = column vector of polynomials
p = CL_chebyshevPoly(0:nbcoef-1)'; 
y0 = horner(p, tn); 
y1 = zeros(tn); 
y2 = zeros(tn); 
if (cvel)
  p = derivat(p); 
  y1 = horner(p, tn);  
  if (cacc)
    p = derivat(p); 
    y2 = horner(p, tn);  
  end
end


// Chebyshev interpolation pour x,y,z (<=> j = 1,2,3)
for (j = 1:3)

   // indices of coefficients for each date
   J = 3 * (It-1) + j;  
   
   // multiply coeffs by polynomial values
   pos(j,Iok) = sum(y0 .* allcoefs(:,J), "r");
   if (cvel); vel(j,Iok) = sum(y1 .* allcoefs(:,J), "r"); end;
   if (cacc); acc(j,Iok) = sum(y2 .* allcoefs(:,J), "r"); end;

end

// change units => meters, seconds
// derivatives / time (and not normalized time)
// s = 2 (t-t0)/T - 1 => d/dt = d/ds * ds/dt = d/ds * (2/T)
pos = pos * 1000; 
vel = vel * 1000 * 2 / (T*86400); 
acc = acc * 1000 * (2 / (T*86400))^2; 

endfunction



function [pos, vel, acc] = de405_getPVA(body, t, orig, ephem, cvel, cacc); 

// "Earth" and  "solar-sys-bary" : added (don't exist in JPL files)
// "Earth": obtained from "Earth-Moon-bary" + "Moon")
// "Moon" : is in fact Earth->Moon in JPL file
// Other names: from solar system barycenter (same order as in file)
// Frame = ICRS (EME2000)

names = [ "Mercury", "Venus", "Earth-Moon-bary", "Mars", "Jupiter", "Saturn", ..
          "Uranus", "Neptune", "Pluto", "Moon", "Sun", ..
      "Earth", "solar-sys-bary" ]; 

numbody = find(names == body); 
if (numbody == []); CL__error("Unknown body"); end

numorg = find(names == orig); 
if (numorg == []); CL__error("Unknown origin"); end


// Earth / solar system barycenter
function [pos, vel, acc] = get_earth(t, cvel, cacc)
  [pos, vel, acc] = de405_interp(ephem, 3, t, cvel, cacc); // Earth-Moon barycenter
  [posm, velm, accm] = de405_interp(ephem, 10, t, cvel, cacc); // Earth->Moon
  emrat = ephem.emrat; 
  pos = pos - posm/(1+emrat); 
  vel = vel - velm/(1+emrat); 
  acc = acc - accm/(1+emrat); 
endfunction 

// Moon / solar system barycenter
function [pos, vel, acc] = get_moon(t, cvel, cacc)
  [pos, vel, acc] = de405_interp(ephem, 3, t, cvel, cacc); // Earth-Moon barycenter
  [posm, velm, accm] = de405_interp(ephem, 10, t, cvel, cacc); // Earth->Moon
  emrat = ephem.emrat; 
  pos = pos + posm * (emrat/(1+emrat)); 
  vel = vel + velm * (emrat/(1+emrat)); 
  acc = acc + accm * (emrat/(1+emrat)); 
endfunction 

// solar system barycenter (=> [0;0;0])
function [pos, vel, acc] = get_ssbary(t, cvel, cacc)
  // NB: vel, acc set to [] if cvel or cacc = %f
  pos = zeros(3, size(t,2)); 
  vel = zeros(3, size(t,2) * cvel); 
  acc = zeros(3, size(t,2) * cacc); 
endfunction 


// Special case : Moon and Earth
if (body == "Moon" & orig == "Earth")
   [pos, vel, acc] = de405_interp(ephem, 10, t, cvel, cacc); 

elseif (body == "Earth" & orig == "Moon")
   [pos, vel, acc] = de405_interp(ephem, 10, t, cvel, cacc); 
   pos = -pos;
   vel = -vel; 
   acc = -acc; 
   
else
   // process "special" bodies first

   if (body == "Earth")
     [pos, vel, acc] = get_earth(t, cvel, cacc); 
   elseif (body == "Moon");  
     [pos, vel, acc] = get_moon(t, cvel, cacc); 
   elseif (body == "solar-sys-bary");  
     [pos, vel, acc] = get_ssbary(t, cvel, cacc); 
   else
     [pos, vel, acc] = de405_interp(ephem, numbody, t, cvel, cacc); 
   end
   
   if (orig == "Earth")
     [pos0, vel0, acc0] = get_earth(t, cvel, cacc); 
   elseif (orig == "Moon");  
     [pos0, vel0, acc0] = get_moon(t, cvel, cacc); 
   elseif (orig == "solar-sys-bary");  
     [pos0, vel0, acc0] = get_ssbary(t, cvel, cacc); 
   else
     [pos0, vel0, acc0] = de405_interp(ephem, numorg, t, cvel, cacc); 
   end

   pos = pos - pos0; 
   vel = vel - vel0; 
   acc = acc - acc0; 
end

endfunction