;+
; Project     : SOHO - LASCO
;                   
; Name        : SUN_EPHEM()
;               
; Purpose     : To calculate the solar ephemeris parameters: ecliptic
;               longitude, P, B0 angles and the semi-diameter.
;  
; Category    : Util, Coords
;             
; Explanation : Allows for planetary and lunar perturbations in the
;               calculation of solar longitude and various other solar
;               positional parameters at date/time requested.
;               Uses semi-rigorous formulae to calculate the solar P (position
;               angle of pole) and B0 (latitude of point at disk centre) angles
;               and also the semi-diameter of the solar disk at the date/time 
;               requested.
;               
; Syntax      : IDL> ang = pb0r(date_time)
;
; Examples    :
;    
; Inputs      : date_time  -  the date/time specified in any CDS format 
;               
; Opt. Inputs : None
;               
; Outputs     : Function returns a 6-element array with
;                                    ang(0) = nu_c (degrees)
;                                    ang(1) = nu_p (degrees)
;                                    ang(2) = P  (degrees)
;                                    ang(3) = B0 (degrees)
;                                    ang(4) = R  semi-diameter (arcmin)
;                                    ang(5) = Apparent longitude of
;                                             Sun (degrees).
;               
; Opt. Outputs: None
;               
; Keywords    : None
;
; Common      : None
;               ;               
; Side effects: None
;                              
; History     : Based on Fortran programs by Hohenkerk and Emerson (RGO)
;
;     16-May-94,CDS/IDL version PB0R.PRO, C D Pike, RAL,Written 
;               
;     Update semi-diameter calculation, CDP, 20-May-94
;		Version 3, William Thompson, GSFC, 14 November 1994
;			Modified .DAY to .MJD
;               Simon Plunkett, UofB, Adapted from CDS routine PB0R to
;               include solar longitude in output parameters, 10 May 1995.
;
;     Version 3, 14 November 1994
;
; Contact     :
;-            

function  sun_ephem, date

; 
; get UT of observation
;
utc = anytim2utc(date)
if utc.mjd eq 0 then return,[0.,0.,0.]

;
;  number of Julian days since 2415020.0
;
jd = cds2jd(utc)
jd.int = jd.int - 2415020L

;
;  ignoring difference between UT and ET then ...
;
de = double(jd.int + jd.frac)

;
;  get the longitude of the sun etc.
;
sun_pos, de, longmed, ra, dec, appl, oblt

;
;   form aberrated longitude 
;
lambda =  longmed  -(20.5D0/3600.0D0)

;
; form longitude of ascending node of sun's equator on ecliptic
;
node  = 73.666666D0 + (50.25D0/3600.0D0)*( (de/365.25d0) + 50.0d0 )
arg  = lambda - node

;
;   calculate P, the position angle of the pole
;   nu_c is the angle between the normal to the ecliptic and
;   (terrestrial) north. This angle is a maximum at the equinoxes and
;   a minimum at the solstices.
;   nu_p is the angle between the normal to the ecliptic and the
;   projection of the Sun's axis on the sky.
;

nu_c = atan(tan(oblt*!Dtor) * cos(appl*!Dtor)) * !Radeg
nu_p = atan( -0.12722D0 * cos(arg*!Dtor)) *!Radeg
p = nu_p - nu_c

;
;    ... and B0 the tilt of the axis
;
b = asin( 0.12620D0 * sin(arg*!dtor) ) * !radeg

;
;  ... and the semi-diameter
;
;
;  Form the mean anomalies of Venus(MV),Earth(ME),Mars(MM),Jupiter(MJ)
;  and the mean elongation of the Moon from the Sun(D).
;
t = de/36525.0d0

mv = 212.6d0   + (  (58517.80D0   * t) mod 360.0D0 )
me = 358.476d0 + (  (35999.0498D0 * t) mod 360.0D0 )
mm = 319.5d0   + (  (19139.86D0   * t) mod 360.0D0 )
mj = 225.3d0   + (  ( 3034.69D0   * t) mod 360.0D0 )
d  = 350.7d0   + ( (445267.11D0   * t) mod 360.0D0 )

;
;  Form the geocentric distance(r) and semi-diameter(sd)
;
r = 1.000141d0 - (0.016748d0 - 0.0000418d0*t)*cos(me*!dtor) $
    - 0.000140d0 * cos(2.0d0*me*!dtor)                       $
    + 0.000016d0 * cos((58.3d0 + 2.0D0*mv - 2.0D0*me)*!dtor) $
    + 0.000005d0 * cos((209.1d0 + mv - me)*!dtor)            $
    + 0.000005d0 * cos((253.8d0 - 2.0d0*mm + 2.0d0*me)*!dtor)$
    + 0.000016d0 * cos(( 89.5d0 - mj + me)*!dtor)            $
    + 0.000009d0 * cos((357.1d0 - 2.0d0*mj + 2.0D0*me)*!dtor) $
    + 0.000031D0 * cos(d*!dtor)

sd = (0.2665685d0/r)*60.0d0

return,[nu_c, nu_p, p , b , sd , appl]

end