Viewing contents of file '../idllib/ssw/allpro/baryvel.pro'
pro baryvel, dje, deq, dvelh, dvelb
;+
; NAME:
; BARYVEL
; PURPOSE:
; Calculates heliocentric and barycentric velocity components of Earth.
;
; EXPLANATION:
; BARYVEL takes into account the Earth-Moon motion, and is useful for
; radial velocity work to an accuracy of ~1 m/s.
;
; CALLING SEQUENCE:
; BARYVEL, dje, deq, dvelh, dvelb
;
; INPUTS:
; DJE - (scalar) julian ephemeris date.
; DEQ - (scalar) epoch of mean equinox of dvelh and dvelb. If deq=0
; then deq is assumed to be equal to dje.
; OUTPUTS:
; DVELH: (vector(3)) heliocentric velocity component. in km/s
; DVELB: (vector(3)) barycentric velocity component. in km/s
;
; The 3-vectors DVELH and DVELB are given in a right-handed coordinate
; system with the +X axis toward the Vernal Equinox, and +Z axis
; toward the celestial pole.
;
; PROCEDURE CALLED:
; Function PREMAT() -- computes precession matrix
;
; NOTES:
; Algorithm taken from FORTRAN program of Stumpff (1980, A&A Suppl, 41,1)
; Stumpf claimed an accuracy of 42 cm/s for the velocity. A
; comparison with the JPL FORTRAN planetary ephemeris program PLEPH
; found agreement to within about 65 cm/s between 1986 and 1994
;
; EXAMPLE:
; Compute the radial velocity of the Earth toward Altair on 15-Feb-1994
;
; IDL> jdcnv, 1994, 2, 15, 0, jd ;==> JD = 2449398.5
; IDL> baryvel, jd, 2000, vh, vb
; ==> vh = [-17.07809, -22.80063, -9.885281] ;Heliocentric km/s
; ==> vb = [-17.08083, -22.80471, -9.886582] ;Barycentric km/s
;
; IDL> ra = ten(19,50,46.77)*15/!RADEG ;RA in radians
; IDL> dec = ten(08,52,3.5)/!RADEG ;Dec in radians
; IDL> v = vb(0)*cos(dec)*cos(ra) + $ ;Project velocity toward star
; vb(1)*cos(dec)*sin(ra) + vb(2)*sin(dec)
;
; REVISION HISTORY:
; Jeff Valenti, U.C. Berkeley Translated BARVEL.FOR to IDL.
; W. Landsman, Cleaned up program sent by Chris McCarthy (SfSU) June 1994
;-
On_Error,2
if N_params() LT 4 then begin
print,'Syntax: BARYVEL, dje, deq, dvelh, dvelb'
print,' dje - input Julian ephemeris date
print,' deq - input epoch of mean equinox of dvelh and dvelb
print,' dvelh - output vector(3) heliocentric velocity comp in km/s
print,' dvelb - output vector(3) barycentric velocity comp in km/s
return
endif
;Define constants
dc2pi = 2*!DPI
cc2pi = 2*!PI
dc1 = 1.0D0
dcto = 2415020.0D0
dcjul = 36525.0D0 ;days in julian year
dcbes = 0.313D0
dctrop = 365.24219572D0 ;days in tropical year (...572 insig)
dc1900 = 1900.0D0
AU = 1.4959787D8
;Constants dcfel(i,k) of fast changing elements.
dcfel = [1.7400353D00, 6.2833195099091D02, 5.2796D-6 $
,6.2565836D00, 6.2830194572674D02, -2.6180D-6 $
,4.7199666D00, 8.3997091449254D03, -1.9780D-5 $
,1.9636505D-1, 8.4334662911720D03, -5.6044D-5 $
,4.1547339D00, 5.2993466764997D01, 5.8845D-6 $
,4.6524223D00, 2.1354275911213D01, 5.6797D-6 $
,4.2620486D00, 7.5025342197656D00, 5.5317D-6 $
,1.4740694D00, 3.8377331909193D00, 5.6093D-6 ]
dcfel = reform(dcfel,3,8)
;constants dceps and ccsel(i,k) of slowly changing elements.
dceps = [4.093198D-1, -2.271110D-4, -2.860401D-8 ]
ccsel = [1.675104E-2, -4.179579E-5, -1.260516E-7 $
,2.220221E-1, 2.809917E-2, 1.852532E-5 $
,1.589963E00, 3.418075E-2, 1.430200E-5 $
,2.994089E00, 2.590824E-2, 4.155840E-6 $
,8.155457E-1, 2.486352E-2, 6.836840E-6 $
,1.735614E00, 1.763719E-2, 6.370440E-6 $
,1.968564E00, 1.524020E-2, -2.517152E-6 $
,1.282417E00, 8.703393E-3, 2.289292E-5 $
,2.280820E00, 1.918010E-2, 4.484520E-6 $
,4.833473E-2, 1.641773E-4, -4.654200E-7 $
,5.589232E-2, -3.455092E-4, -7.388560E-7 $
,4.634443E-2, -2.658234E-5, 7.757000E-8 $
,8.997041E-3, 6.329728E-6, -1.939256E-9 $
,2.284178E-2, -9.941590E-5, 6.787400E-8 $
,4.350267E-2, -6.839749E-5, -2.714956E-7 $
,1.348204E-2, 1.091504E-5, 6.903760E-7 $
,3.106570E-2, -1.665665E-4, -1.590188E-7 ]
ccsel = reform(ccsel,3,17)
;Constants of the arguments of the short-period perturbations.
dcargs = [5.0974222D0, -7.8604195454652D2 $
,3.9584962D0, -5.7533848094674D2 $
,1.6338070D0, -1.1506769618935D3 $
,2.5487111D0, -3.9302097727326D2 $
,4.9255514D0, -5.8849265665348D2 $
,1.3363463D0, -5.5076098609303D2 $
,1.6072053D0, -5.2237501616674D2 $
,1.3629480D0, -1.1790629318198D3 $
,5.5657014D0, -1.0977134971135D3 $
,5.0708205D0, -1.5774000881978D2 $
,3.9318944D0, 5.2963464780000D1 $
,4.8989497D0, 3.9809289073258D1 $
,1.3097446D0, 7.7540959633708D1 $
,3.5147141D0, 7.9618578146517D1 $
,3.5413158D0, -5.4868336758022D2 ]
dcargs = reform(dcargs,2,15)
;Amplitudes ccamps(n,k) of the short-period perturbations.
ccamps = $
[-2.279594E-5, 1.407414E-5, 8.273188E-6, 1.340565E-5, -2.490817E-7 $
,-3.494537E-5, 2.860401E-7, 1.289448E-7, 1.627237E-5, -1.823138E-7 $
, 6.593466E-7, 1.322572E-5, 9.258695E-6, -4.674248E-7, -3.646275E-7 $
, 1.140767E-5, -2.049792E-5, -4.747930E-6, -2.638763E-6, -1.245408E-7 $
, 9.516893E-6, -2.748894E-6, -1.319381E-6, -4.549908E-6, -1.864821E-7 $
, 7.310990E-6, -1.924710E-6, -8.772849E-7, -3.334143E-6, -1.745256E-7 $
,-2.603449E-6, 7.359472E-6, 3.168357E-6, 1.119056E-6, -1.655307E-7 $
,-3.228859E-6, 1.308997E-7, 1.013137E-7, 2.403899E-6, -3.736225E-7 $
, 3.442177E-7, 2.671323E-6, 1.832858E-6, -2.394688E-7, -3.478444E-7 $
, 8.702406E-6, -8.421214E-6, -1.372341E-6, -1.455234E-6, -4.998479E-8 $
,-1.488378E-6, -1.251789E-5, 5.226868E-7, -2.049301E-7, 0.E0 $
,-8.043059E-6, -2.991300E-6, 1.473654E-7, -3.154542E-7, 0.E0 $
, 3.699128E-6, -3.316126E-6, 2.901257E-7, 3.407826E-7, 0.E0 $
, 2.550120E-6, -1.241123E-6, 9.901116E-8, 2.210482E-7, 0.E0 $
,-6.351059E-7, 2.341650E-6, 1.061492E-6, 2.878231E-7, 0.E0 ]
ccamps = reform(ccamps,5,15)
;Constants csec3 and ccsec(n,k) of the secular perturbations in longitude.
ccsec3 = -7.757020E-8
ccsec = [1.289600E-6, 5.550147E-1, 2.076942E00 $
,3.102810E-5, 4.035027E00, 3.525565E-1 $
,9.124190E-6, 9.990265E-1, 2.622706E00 $
,9.793240E-7, 5.508259E00, 1.559103E01 ]
ccsec = reform(ccsec,3,4)
;Sidereal rates.
dcsld = 1.990987D-7 ;sidereal rate in longitude
ccsgd = 1.990969E-7 ;sidereal rate in mean anomaly
;Constants used in the calculation of the lunar contribution.
cckm = 3.122140E-5
ccmld = 2.661699E-6
ccfdi = 2.399485E-7
;Constants dcargm(i,k) of the arguments of the perturbations of the motion
; of the moon.
dcargm = [5.1679830D0, 8.3286911095275D3 $
,5.4913150D0, -7.2140632838100D3 $
,5.9598530D0, 1.5542754389685D4 ]
dcargm = reform(dcargm,2,3)
;Amplitudes ccampm(n,k) of the perturbations of the moon.
ccampm = [ 1.097594E-1, 2.896773E-7, 5.450474E-2, 1.438491E-7 $
,-2.223581E-2, 5.083103E-8, 1.002548E-2, -2.291823E-8 $
, 1.148966E-2, 5.658888E-8, 8.249439E-3, 4.063015E-8 ]
ccampm = reform(ccampm,4,3)
;ccpamv(k)=a*m*dl,dt (planets), dc1mme=1-mass(earth+moon)
ccpamv = [8.326827E-11, 1.843484E-11, 1.988712E-12, 1.881276E-12]
dc1mme = 0.99999696D0
;Time arguments.
dt = (dje - dcto) / dcjul
tvec = [1d0, dt, dt*dt]
;Values of all elements for the instant(aneous?) dje.
temp = (tvec # dcfel) mod dc2pi
dml = temp(0)
forbel = temp(1:7)
g = forbel(0) ;old fortran equivalence
deps = total(tvec*dceps) mod dc2pi
sorbel = (tvec # ccsel) mod dc2pi
e = sorbel(0) ;old fortran equivalence
;Secular perturbations in longitude.
dummy=cos(2.0)
sn = sin((tvec(0:1) # ccsec(1:2,*)) mod cc2pi)
;Periodic perturbations of the emb (earth-moon barycenter).
pertl = total(ccsec(0,*) * sn) + dt*ccsec3*sn(2)
pertld = 0.0
pertr = 0.0
pertrd = 0.0
for k=0,14 do begin
a = (dcargs(0,k)+dt*dcargs(1,k)) mod dc2pi
cosa = cos(a)
sina = sin(a)
pertl = pertl + ccamps(0,k)*cosa + ccamps(1,k)*sina
pertr = pertr + ccamps(2,k)*cosa + ccamps(3,k)*sina
if k lt 11 then begin
pertld = pertld + (ccamps(1,k)*cosa-ccamps(0,k)*sina)*ccamps(4,k)
pertrd = pertrd + (ccamps(3,k)*cosa-ccamps(2,k)*sina)*ccamps(4,k)
endif
endfor
;Elliptic part of the motion of the emb.
phi = (e*e/4d0)*(((8d0/e)-e)*sin(g) +5*sin(2*g) +(13/3d0)*e*sin(3*g))
f = g + phi
sinf = sin(f)
cosf = cos(f)
dpsi = (dc1 - e*e) / (dc1 + e*cosf)
phid = 2*e*ccsgd*((1 + 1.5*e*e)*cosf + e*(1.25 - 0.5*sinf*sinf))
psid = ccsgd*e*sinf / sqrt(dc1 - e*e)
;Perturbed heliocentric motion of the emb.
d1pdro = dc1+pertr
drd = d1pdro * (psid + dpsi*pertrd)
drld = d1pdro*dpsi * (dcsld+phid+pertld)
dtl = (dml + phi + pertl) mod dc2pi
dsinls = sin(dtl)
dcosls = cos(dtl)
dxhd = drd*dcosls - drld*dsinls
dyhd = drd*dsinls + drld*dcosls
;Influence of eccentricity, evection and variation on the geocentric
; motion of the moon.
pertl = 0.0
pertld = 0.0
pertp = 0.0
pertpd = 0.0
for k = 0,2 do begin
a = (dcargm(0,k) + dt*dcargm(1,k)) mod dc2pi
sina = sin(a)
cosa = cos(a)
pertl = pertl + ccampm(0,k)*sina
pertld = pertld + ccampm(1,k)*cosa
pertp = pertp + ccampm(2,k)*cosa
pertpd = pertpd - ccampm(3,k)*sina
endfor
;Heliocentric motion of the earth.
tl = forbel(1) + pertl
sinlm = sin(tl)
coslm = cos(tl)
sigma = cckm / (1.0 + pertp)
a = sigma*(ccmld + pertld)
b = sigma*pertpd
dxhd = dxhd + a*sinlm + b*coslm
dyhd = dyhd - a*coslm + b*sinlm
dzhd= -sigma*ccfdi*cos(forbel(2))
;Barycentric motion of the earth.
dxbd = dxhd*dc1mme
dybd = dyhd*dc1mme
dzbd = dzhd*dc1mme
for k=0,3 do begin
plon = forbel(k+3)
pomg = sorbel(k+1)
pecc = sorbel(k+9)
tl = (plon + 2.0*pecc*sin(plon-pomg)) mod cc2pi
dxbd = dxbd + ccpamv(k)*(sin(tl) + pecc*sin(pomg))
dybd = dybd - ccpamv(k)*(cos(tl) + pecc*cos(pomg))
dzbd = dzbd - ccpamv(k)*sorbel(k+13)*cos(plon - sorbel(k+5))
endfor
;Transition to mean equator of date.
dcosep = cos(deps)
dsinep = sin(deps)
dyahd = dcosep*dyhd - dsinep*dzhd
dzahd = dsinep*dyhd + dcosep*dzhd
dyabd = dcosep*dybd - dsinep*dzbd
dzabd = dsinep*dybd + dcosep*dzbd
;Epoch of mean equinox (deq) of zero implies that we should use
; Julian ephemeris date (dje) as epoch of mean equinox.
if deq eq 0 then begin
dvelh = AU * ([dxhd, dyahd, dzahd])
dvelb = AU * ([dxbd, dyabd, dzabd])
return
endif
;General precession from epoch dje to deq.
deqdat = (dje-dcto-dcbes) / dctrop + dc1900
prema = premat(deqdat,deq, /FK4)
dvelh = AU * ( prema # [dxhd, dyahd, dzahd] )
dvelb = AU * ( prema # [dxbd, dyabd, dzabd] )
return
end