/
x_apohelio.pro
173 lines (153 loc) · 5.57 KB
/
x_apohelio.pro
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
;+
; NAME:
; x_keckhelio
; Version 1.1
;
; PURPOSE:
; Compute correction term to add to velocities to convert to heliocentric.
; This code also does the VLT and MMT. Note, it is the *NEGATIVE* of the
; the number that one applies to a wavelength solution.
;
; CALLING SEQUENCE:
; vcorr = x_keckhelio( ra, dec, [epoch], jd=, tai=, $
; longitude=, latitude=, altitude= )
;
; INPUTS:
; ra - Right ascension [degrees]
; dec - Declination [degrees]
; [epoch] - Epoch of observation for RA, DEC; default to 2000.
;
; RETURNS:
; vcorr - Velocity correction term, in km/s, to add to measured
; radial velocity to convert it to the heliocentric frame.
; OPTIONAL KEYWORDS:
; jd - Decimal Julian date. Note this should probably be
; type DOUBLE.
; tai - Number of seconds since Nov 17 1858; either JD or TAI
; must be specified. Note this should probably either
; be type DOUBLE or LONG64.
; longitude - Longitude of observatory;
; default to (360-155.47220) deg for APO
; latitute - Latitude of observatory; default to 32.780361 deg for APO
; altitude - Altitude of observatory; default to 4000 m for
; Keck
; OBS= - Observatory (default: 'keck')
;
; OUTPUTS:
;
; OPTIONAL OUTPUTS:
;
; COMMENTS:
;
; EXAMPLES:
; helio_shift = -1. * x_keckhelio(RA, DEC, 2000.0)
;
; BUGS:
;
; PROCEDURES CALLED:
; baryvel
; ct2lst
;
; REVISION HISTORY:
; 09-May-2000 Written by S. Burles & D. Schlegel
; 30-Aug-2002 Revised by JXP for Keck
; 22-Dec-2007 Revised to use observatory function by JFH.
;-
;------------------------------------------------------------------------------
function x_apohelio, ra, dec, epoch, jd = jd, tai = tai $
, longitude = longitude, latitude = latitude $
, altitude = altitude, OBS = OBS
if (N_params() LT 2) then begin
print,'Syntax - ' + $
'heliovel = x_keckhelio(ra, dec, [epoch], JD=) v(1.1)'
return, -1
endif
if (NOT keyword_set(epoch)) then epoch = 2000.0
IF KEYWORD_SET(LONGITUDE) AND KEYWORD_SET(LATITUDE) $
AND KEYWORD_SET(altitude) AND NOT KEYWORD_SET(OBS) THEN BEGIN
longitude = 360.0d - longitude
ENDIF ELSE BEGIN
IF NOT KEYWORD_SET(OBS) THEN OBS = 'keck'
if strmatch(OBS,'vlt',/fold) then begin
longitude = 360. - 70.40322
latitude = -24.6258
altitude = 2635. ; meters
endif else if strmatch(OBS,'apo',/fold) then begin
observatory, obs, obs_struct
longitude = 360.0d - 105.820
latitude = 32.780361
altitude = 2788 ;; meters
endif else begin
observatory, obs, obs_struct
longitude = 360.0d - obs_struct.longitude
latitude = obs_struct.latitude
altitude = obs_struct.altitude ;; meters
endelse
ENDELSE
; Default to location of Keck Observatory
; if (NOT keyword_set(longitude)) then longitude = 360. - 155.47220
; if (NOT keyword_set(latitude)) then latitude = 19.82656886
; if (NOT keyword_set(altitude)) then altitude = 4000. ; meters
; case obs of
; 'vlt': begin
; longitude = 360. - 70.40322
; latitude = -24.6258
; altitude = 2635. ; meters
; end
; 'mmt': begin
; longitude = 360. - 110.88456
; latitude = 31.688778
; altitude = 2600. ; meters
; end
; 'lick': begin
; longitude = 360. - 121.637222
; latitude = 37.343056
; altitude = 1283. ; meters
; end
; 'mko': begin
; longitude = 360. - 155.47220
; latitude = 19.82656886
; altitude = 4000. ; meters
; end
; else: stop
; endcase
if (NOT keyword_set(jd)) then begin
if (keyword_set(tai)) then begin
jd = 2400000.5D + tai / (24.D*3600.D)
endif else begin
message, 'Must specify either JD or TAI', /cont
return, 0
endelse
endif
DRADEG = 180.d0 / !DPI
;----------
; Compute baryocentric velocity
baryvel, jd, epoch, dvelh, dvelb
; Project velocity toward star
vbarycen = dvelb[0]*cos(dec/DRADEG)*cos(ra/DRADEG) + $
dvelb[1]*cos(dec/DRADEG)*sin(ra/DRADEG) + dvelb[2]*sin(dec/DRADEG)
;----------
; Compute rotational velocity of observer on the Earth
; LAT is the latitude in radians.
latrad = latitude / DRADEG
; Reduction of geodetic latitude to geocentric latitude (radians).
; DLAT is in arcseconds.
dlat = -(11.d0 * 60.d0 + 32.743000d0) * sin(2.d0 * latrad) + $
1.163300d0 * sin(4.d0 * latrad) -0.002600d0 * sin(6.d0 * latrad)
latrad = latrad + (dlat / 3600.d0) / DRADEG
; R is the radius vector from the Earth's center to the observer (meters).
; VC is the corresponding circular velocity
; (meters/sidereal day converted to km / sec).
; (sidereal day = 23.934469591229 hours (1986))
r = 6378160.0d0 * (0.998327073d0 + 0.00167643800d0 * cos(2.d0 * latrad) - $
0.00000351d0 * cos(4.d0 * latrad) + 0.000000008d0 * cos(6.d0 * latrad)) $
+ altitude
vc = 2.d0 * !DPI * (r / 1000.d0) / (23.934469591229d0 * 3600.d0)
; Compute the hour angle, HA, in degrees
ct2lst, LST, longitude, junk, jd
LST = 15. * LST ; convert from hours to degrees
HA = LST - ra
; Project the velocity onto the line of sight to the star.
vrotate = vc * cos(latrad) * cos(dec/DRADEG) * sin(HA/DRADEG)
return, (-vbarycen + vrotate)
end