PRO co_aberration, jd, ra, dec, d_ra, d_dec, eps=eps
;+
; NAME:
; CO_ABERRATION
; PURPOSE:
; Calculate changes to Ra and Dec due to the effect of annual aberration
; EXPLANATION:
; as described in Meeus, Chap 23.
; CALLING SEQUENCE:
; co_aberration, jd, ra, dec, d_ra, d_dec, [EPS = ]
; INPUTS
; jd : Julian Date [scalar or vector]
; ra, dec : Arrays (or scalars) of the ra and dec's in degrees
; Note: if jd is a vector, then ra and dec must either be scalars, or
; vectors of the same length.
;
; OUTPUTS
; d_ra, d_dec: the corrections to ra and dec due to aberration in
; arcseconds. (These values can be added to the true RA
; and dec to get the apparent position). Note that d_ra
; is *not* multiplied by cos(dec), so that
; apparent_ra = ra + d_ra/3600.
; OPTIONAL INPUT KEYWORD:
; eps : set this to the true obliquity of the ecliptic (in radians), or
; it will be set for you if you don't know it (in that case, set it to
; an empty variable).
; EXAMPLE:
; Compute the change in RA and Dec of Theta Persei (RA = 2h46m,11.331s, Dec =
; 49d20',54.54") due to aberration on 2028 Nov 13.19 TD
;
; IDL> jdcnv,2028,11,13,.19*24,jd ;Get Julian date
; IDL> co_aberration,jd,ten(2,46,11.331)*15,ten(49,20,54.54),d_ra,d_dec
;
; ==> d_ra = 30.045" (=2.003s) d_dec = 6.697"
; NOTES:
; These formula are from Meeus, Chapters 23. Accuracy is much better than 1
; arcsecond.
;
; The maximum deviation due to annual aberration is 20.49" and occurs when the
; Earth velocity is perpendicular to the direction of the star.
;
; REVISION HISTORY:
; Written, June 2002, Chris O'Dell, U. of Wisconsin
; Fix error with vector input W. Landsman June 2009
; June 2009 update fixed case where JD was scalar but RA,Dec were vectors, but
; broke the case when both JD and RA,Dec were vectors Aug 2012 W. Landsman
; Further fix when JD is 1 element vector W. Landsman
;-
compile_opt idl2
d2r = !dpi/180.
if N_elements(jd) EQ 1 then jd = jd[0]
T = (jd -2451545.0)/36525.0 ; julian centuries from J2000 of jd.
if n_elements(eps) eq 0 then begin ; must calculate obliquity of ecliptic
njd = n_elements(jd)
d_psi = dblarr(njd)
d_epsilon = d_psi
for i=0L,njd-1 do begin
nutate, jd[i], dp, de ; d_psi and d_epsilon in degrees
d_psi[i] = dp
d_epsilon[i] = de
endfor
eps0 = ten(23,26,21.448)*3600.d - 46.8150*T - 0.00059*T^2 + $
0.001813*T^3
eps = (eps0 + d_epsilon)/3600.*d2r ; true obliquity of the ecliptic
; in radians
endif
sunpos, jd, sunra, sundec, sunlon
; Earth's orbital eccentricity
e = 0.016708634d - 0.000042037d*T - 0.0000001267d*T^2
; longitude of perihelion, in degrees
pi = 102.93735 + 1.71946*T + 0.00046*T^2
k = 20.49552 ;constant of aberration, in arcseconds
;Useful Trig Functions
cd = cos(dec*d2r) & sd = sin(dec*d2r)
if N_elements(eps) EQ 1 then eps = eps[0] ;Special scalar case
ce = cos(eps) & te = tan(eps)
cp = cos(pi*d2r) & sp = sin(pi*d2r)
cs = cos(sunlon*d2r) & ss = sin(sunlon*d2r)
ca = cos(ra*d2r) & sa = sin(ra*d2r)
term1 = (ca*cs*ce+sa*ss)/cd
term2 = (ca*cp*ce+sa*sp)/cd
term3 = (cs*ce*(te*cd-sa*sd)+ca*sd*ss)
term4 = (cp*ce*(te*cd-sa*sd)+ca*sd*sp)
d_ra = -k * term1 + e*k * term2
d_dec = -k * term3 + e*k * term4
END