//  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_sun_EOD, vel_sun_EOD] = CL__eph_sun_Meeus(jd, cvel)
// Sun position and velocity in EOD frame (Meeus algorithm)
//
// Calling Sequence
// [pos_sun_EOD, vel_sun_EOD] = CL__eph_sun_Meeus(jd [,cvel])
//
// Description
// <itemizedlist><listitem>
// <p>Computes Moon position and velocity relative to the ecliptic of date frame (EOD) using Meeus algorithm.</p>
// <p>The frame is centered on Earth. (pos_sun_EOD is the Earth-Sun vector).</p>
// <p></p>
// <p>Note: if cvel is set to %t but only one output argument is required, the velocity is not computed.</p>
// </listitem>
// </itemizedlist>
//
// Parameters
// jd : Julian date (TT) (1xN)
// cvel: (optional, boolean) %t if velocity should be computed. Default is %t. (1x1)
// pos_sun_EOD : Position relative to EOD frame (m) (3xN)
// vel_sun_EOD : Velocity relative to EOD frame (m/s) (3xN)
//
// Bibliography
// Astronomical Algorithms J.Meeus - Second edition 1998, Chapter 25
//
// Authors
// CNES - DCT/SB
//
// Examples
// // Julian date, in TT time scale
// jd = 21512210.3;
//
// pos_sun_EOD = CL__eph_sun_Meeus(jd)
//


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

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

if (~exists("cvel","local")); cvel = %t; end;

// Handle [] case
if (jd == [])
  pos_sun_EOD = [];
  vel_sun_EOD = [];
  return; // <-- RETURN
end

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

// --------------------------
// ephemeris computation
// --------------------------

// Time in centuries from epoch 2000 12h (TT time scale)
TT   = (jd - 2451545.0) / 36525.0;
TT2   = TT .* TT ;
TT3   = TT2 .* TT ;
TT4   = TT2 .* TT2 ;

// Sun mean longitude
L0 = 280.46646 + 36000.76983 * TT + 0.0003032 * TT2; // deg
L0 = pmodulo(L0 * DEG2RAD, 2*%pi);

// Sun Mean Anomaly
M = 357.52911 + 35999.05029 * TT - 0.0001537 * TT2; // deg
M = pmodulo(M * DEG2RAD, 2*%pi);

// Eccentricity of Earth orbit
e = 0.016708634 - 4.2037e-05 * TT - 1.267e-07 * TT2;

// Equation of Sun center
C = (1.914602 - 0.004817 * TT - 0.000014 * TT2) .* sin(M) + ..
       (0.019993 - 0.000101 * TT) .* sin(2*M) + 0.000289 .* sin(3*M); // deg
C = pmodulo(C * DEG2RAD, 2*%pi);

// Sun longitude    
lon = L0 + C;    
 
// True anomaly
v = M + C;

// Sun radius vector (m)
R = AU * 1.000001018 * (1-e.^2) ./ (1 + e .* cos(v)) ;

// ---------------------  
// Velocity computation
// ---------------------
if (cvel)
  // Derivative of Sun mean longitude (rad/century)
  L0_d = (36000.76983 + 0.0003032 * 2*TT) * DEG2RAD;
  
  // Derivative of Sun Mean Anomaly (rad/century)
  M_d = (35999.05029 - 0.0001537 * 2*TT) * DEG2RAD ;
  
  // Derivative of eccentricity of Earth orbit
  e_d = -4.2037e-05  - 1.267e-07 * 2*TT;
  
  // Derivative of equation of Sun center (rad/century)
  C_d = (( -0.004817 - 0.000014 * 2*TT) .* sin(M) + ..
        (1.914602 - 0.004817 * TT - 0.000014 * TT2) .* M_d .* cos(M) + ..
        -0.000101 .* sin(2*M) + (0.019993 - 0.000101 * TT) .* (2*M_d) .* cos(2*M) + ..
        0.000289 .* (3*M_d) .* cos(3*M)) * DEG2RAD;
         
  // Derivative of Sun longitude (rad/century) 
  lon_d = L0_d + C_d;
  
  // Derivative of true anomaly (rad/century)
  v_d = M_d + C_d;
  
  // Derivative of Sun radius vector (m/century)
  // (u/v)'=(vu'-uv'/v2)
  R_d = AU * 1.000001018 * ((1 + e .* cos(v)) .* (-2 * e_d .* e) -               ...   // vu'
                               ((1-e.^2) .* (e_d.*cos(v) - e .* v_d .* sin(v)))) ./  ...   // uv'                                                    // uv'
                               ((1 + e .* cos(v)).^2);                                     // v2
end

// ------------------   
// Output
// ------------------ 

// Position (in EOD frame)
if (~cvel)
  pos_sun_EOD = CL_co_sph2car([lon;zeros(lon);R]);
  vel_sun_EOD = [];
  
// Position and velocity (in EOD frame)
else
  [pos_sun_EOD,vel_sun_EOD] = CL_co_sph2car([lon;zeros(lon);R],[lon_d;zeros(lon_d);R_d]/(36525*86400));
end
           
endfunction