;+
; computes the projected velocity of the telescope wrt
; six coordinate systems: geo, helio, bary, lsrk, lsrd, gal
; negative velocities mean approach
;
; The standard LSR is defined as follows: the sun moves at 20.0 km/s
; toward ra=18.0h, dec=30.0 deg in 1900 epoch coords
;
; <p><b>CHDOPPLER IS NOT VECTORIZED; ONLY SINGLE NUMBERS</b>
;
; <p>
; This code came via e-mail from Carl Heiles via Tom Bania on 11/04/2004.
; Local changes:
; <UL>
; <LI> modify this documentation for use by idldoc.
; <LI> removed path argument and related code, replaced by obspos
; argument.
; <LI> default position is the GBT.
; <LI> Observatory longitude was not being passed to juldaytolmst.
; <LI> LSRD added 
; <LI> Galactocentric added
; <LI> Checked against aips++ Measures.  Differences were < 20 m/s in
; one test case (<10m/s for geo, bary, and lsrk).
; <LI> Double precision throughout.
; <LI> Relativistic addition of velocities.
; </UL>
;
; Previous revision history: carlh 29oct04
; <UL>
; <LI> from idpppler_chl; changed calculation epoch to 2000.
; </UL>
;
; @param ra {in}{required} The source ra in decimal hours, equinox 2000
; @param dec {in}{required} The source dec in decimal hours, equinox 2000
; @param julday {in}{required} The julian day
;
; @keyword obspos {in}{optional}{double [2]} observatory position
; [East longitude, latitude] in degrees.  
; Uses the GBT position if not specified.
; 
; @keyword light {in}{optional}{boolean} When set, returns the
; velocity as a fraction of c
;
;
; @returns The velocity in km/s, or as a faction of c if
; the keyword /light is specified. the result is a 6-element
; vector whose elements are [geo, helio, bary, lsrk, lsrd, gal].
;
; @uses <a href="http://idlastro.gsfc.nasa.gov/ftp/pro/astro/baryvel.pro">baryvel</a>
; @uses <a href="http://idlastro.gsfc.nasa.gov/ftp/pro/astro/precess.pro">precess</a>
; @uses <a href="juldaytolmst.html">juldaytolmst</a>
;
; @version $Id: chdoppler.pro,v 1.3 2005/04/08 20:36:06 bgarwood Exp $
;-
function chdoppler, ra, dec, julday, $
        obspos=obspos, light=light
;
; Default to GBT if obspos not provided.
if (not keyword_set(obspos)) then begin
    obspos=dblarr(2);
    obspos[0]=[-(79.D + 50.D/60.D + 23.42D/3600.D)]
    obspos[1]=[(38.D + 25.D/60.D + 59.26D/3600.D)]
endif

;------------------ORBITAL SECTION-------------------------
;GET THE COMPONENTS OF RA AND DEC, 2000u EPOCH
rasource=ra*15.d*!dtor
decsource=dec*!dtor
xxsource = dblarr(3)
xxsource[0] = cos(decsource) * cos(rasource)
xxsource[1] = cos(decsource) * sin(rasource)
xxsource[2] = sin(decsource)

;GET THE EARTH VELOCITY WRT THE SUN CENTER
baryvel, julday, 2000.,vvorbit,velb

;PROJECTED VELOCITY OF EARTH CENTER WRT SUN TO THE SOURCE
pvorbit_helio = total(vvorbit* xxsource)
pvorbit_bary = total(velb* xxsource)



;-----------------------LSRK SECTION-------------------------
;THE STANDARD LSRK IS DEFINED AS FOLLOWS: THE SUN MOVES AT 20.0 KM/S
;TOWARD RA=18.0H, DEC=30.0 DEG IN 1900 EPOCH COORDS
;using PRECESS, this works out to ra=18.063955 dec=30.004661 in 2000
;coords.
ralsrk_rad= 2.d*!pi*18.d/24.d
declsrk_rad= !dtor*30.d
precess, ralsrk_rad, declsrk_rad, 1900.d, 2000.d,/radian

;FIND THE COMPONENTS OF THE VELOCITY OF THE SUN WRT THE LSRK FRAME
xxlsrk = dblarr(3)
xxlsrk[0] = cos(declsrk_rad) * cos(ralsrk_rad)
xxlsrk[1] = cos(declsrk_rad) * sin(ralsrk_rad)
xxlsrk[2] = sin(declsrk_rad)
vvlsrk = 20.d*xxlsrk

;PROJECTED VELOCITY OF THE SUN WRT LSRK TO THE SOURCE
pvlsrk=total(vvlsrk*xxsource)


;-----------------------LSRD SECTION-------------------------
;THE LSRD IS DEFINED AS FOLLOWS: THE SUN MOVES AT 16.6 KM/S
;TOWARD RA=17:49:58.7 hours, DEC=28.07.04 DEG IN 2000 EPOCH COORDS

ralsrd_rad= 2.d*!dpi*(17.D + 49.D/60.D + 58.7D/3600.D)/24.d
declsrd_rad= !dtor*(28.D + 07.D/60.D + 04.0D/3600.D)

;FIND THE COMPONENTS OF THE VELOCITY OF THE SUN WRT THE LSRD FRAME
xxlsrd = dblarr(3)
xxlsrd[0] = cos(declsrd_rad) * cos(ralsrd_rad)
xxlsrd[1] = cos(declsrd_rad) * sin(ralsrd_rad)
xxlsrd[2] = sin(declsrd_rad)
vvlsrd = 16.6D*xxlsrd

;PROJECTED VELOCITY OF THE SUN WRT LSRD TO THE SOURCE
pvlsrd=total(vvlsrd*xxsource)

;-----------------------GALACTOCENTRIC SECTION------------------
; LSRD + 220 km/s towards RA=21:12:01.1 DEC=48.19.47 in J2000 Epoch

ragal_rad = 2.d*!dpi*(21.D + 12.D/60.D + 01.1D/3600.D)/24.d
decgal_rad = !dtor*(48.D + 19.D/60.D + 47.D/3600.D)

; Find the components of the velocity of the sun wrt this frame
xxgal = dblarr(3)
xxgal[0] = cos(decgal_rad) * cos(ragal_rad)
xxgal[1] = cos(decgal_rad) * sin(ragal_rad)
xxgal[2] = sin(decgal_rad)
vvgal = 220.D*xxgal

;PROJECTED VELOCITY OF THE SUN WRT GAL TO THE SOURCE
pvgal=total(vvgal*xxsource)

;---------------------EARTH SPIN SECTION------------------------
lst_mean= 24.d/(2.d*!pi)*juldaytolmst( julday,obslong=obspos[0])
raspin=(lst_mean- ra- 6.d)* 15.d* !dtor
decspin=0.0d*!dtor

;PROJECTED DISTANCE FROM CENTER OF EARTH
lat=obspos[1]
gclat = lat - 0.1924d * sin(2.d*lat*!dtor) ; true angle lat
rearth =( 0.99883d + 0.00167d * cos(2.d*lat*!dtor))* 6378.1d ;dist from center, km
rho=rearth*cos(gclat*!dtor)

;SPIN VELOCITY, KM/S
vspin=2.d*!pi*rho/86164.090d

xxspin = dblarr(3)
xxspin[0] = cos(decspin) * cos(raspin)
xxspin[1] = cos(decspin) * sin(raspin)
xxspin[2] = sin(decspin)
vvspin = vspin*xxspin

;PROJECTED VELOCITY OF STATION WRT EARTH CENTER TO SOURCE
pvspin=total(vvspin*xxsource)


;---------------------NOW PUT IT ALL TOGETHER------------------

vtotal= dblarr( 6)
vtotal[0]= -pvspin
vtotal[1]= -shiftvel(pvspin,pvorbit_helio)
vtotal[2]= -shiftvel(pvspin,pvorbit_bary)
vtotal[3]= -shiftvel(-vtotal[2],pvlsrk)
vtotal[4]= -shiftvel(-vtotal[2],pvlsrd)
vtotal[5]= -shiftvel(-vtotal[4],pvgal)

if keyword_set(light) then vtotal=vtotal/(!gc.light_speed*1.d3)

return,vtotal
end