pro calhrs_hel,ih,wave,log
;
;+
;			calhrs_hel
; Convert to heliocentric wavelengths by correcting for the earth's
; motion around the sun.
;
; CALLING SEQUENCE:
;	calhrs_help,ih,wave,log
;
; INPUT/OUTPUTS:
;
;	ih - 128 x N integer*2 header vectors
;	wave - wavelength vectors (npts x N)
;	log - processing log (with HST DG2 data included)
;
;
; DESCRIPTION
;	Use the low-precision formulae for the Sun's coordiantes described
;	in the Astronimical Almanac (1984), page C24.  The velocities
;	are obtained by takeing the derivatives of the coordinates.  The
;	velocity vector is in the equatorial coordiante system of epoch
;	J2000.  THis algorithm does not include the Earth-Moon Motion, Sun-
;	barycenter motion, light time correction from the earth to the Sun.
;	It should be accurate to approx. 0.025 km/sec. within the epoch range
;	of 1900 to 2100 AD.
;
; HISTORY:
;	version 1  D. Lindler  Nov 1989
;       Nov 21 1991  (VER 1.1) JHB/CSC  - Modified to handle either old or 
;			new style PODPS keywords.
;	11-jun-1992	JKF/ACC		- version 1.2 - removed warning
;					messages for old RA/DEC. Added 
;					error message if DEC was not found.
;	30-jun-1992	JKF/ACC		- fixed bug in DECLNTRG condition.
;
;-             
;----------------------------------------------------------------------------
	VERSION = 1.2
	if !prelaunch then return
 
;
; extract RA and DEC from DG2 area of the LOG -------------------------------
;
	;
	; RTASNTRG changed to RA_TARG in new style PODPS keywords (JHB 11-21-91)
	;
        ra=sxpar(log,'RA_TARG')
        if !err lt 0 then begin
          RA = sxpar(log,'RTASNTRG')
  	  if !err lt 0 then begin
	     message,' Keyword (RA_TARG/RTASNTRG) not found in log',/cont
	     message,' WARNING: No heliocentric wavelength correction done',/cont
	     !err = -1
	     return
	  endif
	endif
	;
	; DECLNTRG changed to DEC_TARG in new style PODPS keywords
	;	(JHB 11-21-91)
	;
        dec = sxpar(log,'DEC_TARG')
        if !err lt 0 then begin
          DEC = sxpar(log,'DECLNTRG')
  	  if !err lt 0 then begin
	     message,' Keyword (DEC_TARG/DECLNTRG) not found in log',/cont
	     message,' WARNING: No heliocentric wavelength correction done',/cont
	     !err = -1
	     return
	  endif
	end
;
; get the time ---------------------------------------------------------------
;
	time = string(byte(ih,108,24))	;DD-MMM-YYYY HH:MM:SS.SS
	day = long(strmid(time,0,2))
	mon = strmid(time,3,3)
	year = long(strmid(time,7,4))
	hour = long(strmid(time,12,2))
	min = long(strmid(time,15,2))
	sec = double(strmid(time,18,5))
;
; convert month to number
;
	months = ['   ','JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG', $
		  'SEP','OCT','NOV','DEC']
	month = 0
	for i=1,12 do if mon eq months(i) then month=i
	if month eq 0 then begin
		print,'CALHRS_HEL-- Invalid packet time ='+time
		retall
	endif
;
; convert to fraction of a day
;
	day = sec/(24.0D0*3600.0D0) + min/(24.0D0*60.0D0) + hour/24.0D0 + $
		day
;
; Convert to modified Julian Date
;
	intexp = 367*year - 7*(year+(month+9)/12)/4 - $
		 	3*((year+(month-9)/7)/100+1)/4 + $
			275*month/9 + 1721029 - 2400001
	MJD = intexp + day
;
; compute velocity vector of the earth  ----------------------------------------
;
	AU = 1.4959787D8
	RADN = 57.29577951308232088D0
	MJ2000 = 51544.5D0
;
; coefficients for mean ecliptic longitude of the Earth, corrected for
; abberation:
;
	L0 = (280.460 - 180.0)/RADN
	L1 = (0.985647D0 - 50.29/(3600.0 *365.25))/RADN ;2nd term is precision
	L2 = 1.915/RADN
	l3 = 0.02/RADN
;
; coefficients for mean anomaly
;
	G0 = 357.528D0/RADN
	G1 = 0.9856003D0/RADN
;
; coefficients for distance from the Sun in km.
;
	R0 = 1.00014D0 * AU
	R1 = -0.01671D0 * AU
	R2 = -0.00014D0 * AU
;
; days since J2000
;
	TDAYS = MJD - MJ2000
;
; obliquity of ecliptic at J2000
;
	OBL = 23.439291 / RADN
;
; mean anomaly (radians) and its derivative (radians/day)
;
	G = G0 + G1*TDAYS
	DG = G1
 
	SING = sin(G)
	COSG = cos(G)
	SIN2G = 2.0*SING*COSG
	COS2G = COSG*COSG - SING*SING
;
; ecliptic longitude of the Earth (radians): distance from the sun (km)
;
	LAMBDA = L0 + L1*TDAYS + L2*SING + L3*SIN2G
	R = R0 + R1*COSG + R2*COS2G
;
; derivativs in radians
;
	DLAM = (L1 + L2*DG*COSG + 2.*L3*DG*COS2G)/86400.
	DR = (-R1*DG*SING - 2.0*R2*DG*SIN2G)/86400.
 
	SINLAM = sin(LAMBDA)
	COSLAM = cos(LAMBDA)
	VX = DR*COSLAM - R*SINLAM*DLAM
	VY = (DR*SINLAM + R*COSLAM*DLAM) * COS(OBL)
	VZ = (DR*SINLAM + R*COSLAM*DLAM) * SIN(OBL)
;
; compute velocity (km/sec) towards the target
;
	DEC = DEC/RADN
	RA = RA/RADN
	V = VX*COS(DEC)*COS(RA) + VY*COS(DEC)*SIN(RA) + VZ*SIN(DEC)
;
; correct wavelengths -------------------------------------------------------
;
	WAVE = WAVE*(1.0D0 + V/2.99792458D5)
;
; update history
;
	hist = strarr(2)
	hist(0) = 'CALHRS_HEL version'+string(version,'(f5.2)')+ $
		': Conversion to Heliocentric wavelengths'
	hist(1) = "    Earth's velocity towards target = " + $
		string(v,'(f7.2)') + ' km/sec.'
	sxaddhist,hist,log
	if !dump gt 0 then printf,!textunit,hist
	ih(67,0) = ih(67,*) or 16
	!err = 1
	return
end