File: t_bcvcorr1.x

package info (click to toggle)
iraf-rvsao 2.8.3-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 16,456 kB
  • sloc: ansic: 963; lisp: 651; fortran: 397; makefile: 27
file content (340 lines) | stat: -rw-r--r-- 9,246 bytes parent folder | download
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
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
# File rvsao/Util/t_bcvcorr.x
# July 5, 1995
# By Doug Mink, Harvard-Smithsonian Center for Astrophysics

# Copyright(c) 1995 Smithsonian Astrophysical Observatory
# You may do anything you like with this file except remove this copyright.
# The Smithsonian Astrophysical Observatory makes no representations about
# the suitability of this software for any purpose.  It is provided "as is"
# without express or implied warranty.

# BCVCORR is an IRAF task which computes heliocentric and solar system
# barycentric corrections for spectra


include	<imhdr.h>
include	<fset.h>
include <imhdr.h>
include <imio.h>

define	LEN_USER_AREA	100000

procedure t_bcvcorr ()

pointer	specim		# image structure for spectrum
char	specfile[SZ_PATHNAME]	# Object spectrum file name
char	specdir[SZ_PATHNAME]	# Directory for object spectra
char	specpath[SZ_PATHNAME]	# Object spectrum path name
char	obsname[SZ_LINE]	# Observatory name
char	keyword[SZ_LINE]	# Observatory name
pointer	speclist	# List of spectrum files
bool	savebcv,savebcv0	# Save velocity correction in data file header
int	mode		# Mode (READ_ONLY, READ_WRITE)
pointer	obs		# pointer to observatory structure
bool	gottime		# True if observation time obtained
bool	gotpos		# True if position of observed object obtained
#bool	debug
bool	verbose, clgetb()
int	mm, dd, yyyy, ldir, ndim, k
double	ra, dec, eq, exp, dindef
double	ut1,ut2,ut, dlong, dlat, dalt, gjd
double	dc[3],dcc[3],dprema[3,3],dvelh[3],dvelb[3]
double	djd,dpi,daukm,dctrop,dcbes,dc1900,dct0,dtr,dst
double	dlongs, dlats, deqt, dras, decs, dha, dra2, dec2
double	dgcvel, dbcvel, dhcvel
real	rbcv, rhcv
int	strlen(), imtgetim(), imaccess(), strncmp()
double	julday(), clgetd()
pointer	imtopenp()
pointer	immap()
errchk	immap()
double	obsgetd()
pointer	obsopen()
errchk	obsopen

define	newspec_	10
define	endvc_	90

begin
	dpi = 3.1415926535897932d0
	daukm  =1.4959787d08
	dctrop = 365.24219572d0
	dcbes = 0.313d0
	dc1900 = 1900.0d0
	dctrop = 365.24219572d0
	dcbes = 0.313d0
	dc1900 = 1900.0d0
	dtr = dpi / 180.0d0

# Spectra for which to compute solar system center velocity corrections
	speclist = imtopenp ("spectra")
	call clgstr ("specdir",specdir,SZ_PATHNAME)
	verbose = clgetb ("verbose")
	savebcv0 = clgetb ("savebcv")

# Get next object spectrum file name from the list
newspec_
	savebcv = savebcv0
	if (imtgetim (speclist, specfile, SZ_PATHNAME) == EOF)
	   go to endvc_

# Check for readability of object spectrum
	ldir = strlen (specdir)
	if (ldir > 0) {
	    if (specdir[ldir] != '/') {
		specdir[ldir+1] = '/'
		specdir[ldir+2] = EOS
		}
	    call strcpy (specdir,specpath,SZ_PATHNAME)
	    call strcat (specfile,specpath,SZ_PATHNAME)
	    }
	else
	    call strcpy (specfile,specpath,SZ_PATHNAME)
	if (imaccess (specpath, READ_ONLY) == NO) {
	    call eprintf ("BCVCORR: cannot read spectrum file %s \n")
		call pargstr (specpath)
	    go to newspec_
	    }

# Load spectrum header
	mode = READ_WRITE
        iferr (specim = immap (specpath, mode, LEN_USER_AREA)) {
            mode = READ_ONLY
	    iferr (specim = immap (specpath, mode, LEN_USER_AREA)) {
		if (specim != NULL) call imunmap (specim)
		call eprintf ("BCVCORR: Cannot read image %s\n")
		call pargstr (specpath)
		specim = ERR
		go to endvc_
		}
	    call eprintf ("BCVCORR: Cannot write to image %s\n")
	    call pargstr (specpath)
	    savebcv = FALSE
	    }
        ndim = IM_NDIM(specim)
	if (verbose) {
	    call printf ("%s: %d x %d x %d %d-D image\n")
		call pargstr (specpath)
		call pargi (IM_LEN(specim,1))
		call pargi (IM_LEN(specim,2))
		call pargi (IM_LEN(specim,3))
		call pargi (ndim)
	    }

	dindef = INDEFD

# Direction for velocity correction
	call clgstr ("keyra",keyword,SZ_LINE)
	ra = dindef
 	call imgdpar (specim, keyword, ra)
	call clgstr ("keydec",keyword,SZ_LINE)
	dec = dindef
	call imgdpar (specim,keyword, dec)
	call clgstr ("keyeqnx",keyword,SZ_LINE)
	eq = 1950.0d0
	call imgdpar (specim,keyword,eq)
	if (ra != dindef && dec != dindef)
	    gotpos = TRUE
	else
	    gotpos = FALSE

# Midtime of observation
	call clgstr ("keyjd",keyword,SZ_LINE)
	gjd = dindef
	call imgdpar (specim, keyword, gjd)
	if (gjd == dindef) {
	    call clgstr ("keymid",keyword,SZ_LINE)
	    ut = dindef
	    call imgdpar (specim,keyword, ut)
	    if (ut == dindef) {
		call clgstr ("keyend",keyword,SZ_LINE)
		ut2 = dindef
		call imgdpar (specim, keyword, ut2)
		ut1 = dindef
		call clgstr ("keystart",keyword,SZ_LINE)
		call imgdpar (specim,keyword, ut1)
		exp = 0.
		call clgstr ("keyexp",keyword,SZ_LINE)
		call imgdpar (specim,keyword, exp)
		if (ut1 != dindef && ut2 != dindef) {
		    ut = (ut1 + ut2) * 0.5
		    gottime = TRUE
		    }
		else if (ut2 != dindef && exp > 0.d0) {
		    ut = ut2 - (exp * 0.5 / 3600.)
		    gottime = TRUE
		    }
		else if (ut1 != dindef && exp > 0.d0) {
		    ut = ut1 + (exp * 0.5 / 3600.)
		    gottime = TRUE
		    }
		else {
		    ut = dindef
		    gottime = FALSE
		    }
		}
	    else {
		ut = dindef
		gottime = FALSE
		}

# Date of observation, if midtime is not Julian Date
	    mm = INDEFI
	    dd = INDEFI
	    yyyy = INDEFI
	    call clgstr ("keydate",keyword,SZ_LINE)
	    call imgdate (specim, keyword, mm, dd, yyyy)
	    if (verbose ) {
		call printf("%d/%d/%d %h UT  ra: %h, dec: %h\n")
		    call pargi (dd)
		    call pargi (mm)
		    call pargi (yyyy)
		    call pargd (ut)
		    call pargd (ra)
		    call pargd (dec)
		}
	    if (mm != dindef && dd != dindef && yyyy != INDEFI)
		gottime = TRUE
	    else
		gottime = FALSE
	    }
	else {
	    if (verbose ) {
		call printf("Julian date: %.4f\n")
		    call pargd (gdj)
		}
	    gottime = TRUE
	    }

# Location of observatory
	obsname[1] = EOS
	dlong = dindef
	dlat = dindef
	dalt = dindef
	call clgstr ("obsname",obsname,SZ_LINE)
	if (strncmp (obsname,"file",4) == 0) {
	    call clgstr ("keyobs",keyword,SZ_LINE)
 	    call imgspar (specim,keyword,obsname,SZ_LINE)
	    call printf ("BCVCORR: %s = %s\n")
		call pargstr (keyword)
		call pargstr (obsname)
	    if (obsname[1] != EOS) {
		call clgstr ("keylong",keyword,SZ_LINE)
 		call imgdpar (specim,keyword,dlong)
		call clgstr ("keylong",keyword,SZ_LINE)
 		call imgdpar (specim,keyword,dlat)
		call clgstr ("keyalt",keyword,SZ_LINE)
 		call imgdpar (specim,keyword,dalt)
		}
	    }
	else {
            obs = obsopen (obsname)
            call obslog (obs, "BCVCORR", "latitude longitude altitude", STDOUT)
            dlat = obsgetd (obs, "latitude")
            dlong = obsgetd (obs, "longitude")
            dalt = obsgetd (obs, "altitude")
            call obsclose (obs)
	    }
	if (dlong == dindef && dlat == dindef &&dalt == dindef) {
	    dlong = clgetd ("obslong")
	    dlat = clgetd ("obslat")
	    dalt = clgetd ("obsalt")
	    }
	if (verbose) {
	    call printf("%s lat %h , long %h, alt %.1fm\n")
		call pargstr (obsname)
		call pargd (dlat)
		call pargd (dlong)
		call pargd (dalt)
	    }
         
# reformat so juldat and heliocvel and use the info
	rhcv = 0.
	rbcv = 0.
	if (gottime && gotpos) {

	# Calculate julian date
            djd = julday(mm,dd,yyyy) - 0.5d0
            djd = djd + ut / 24.0d0

	# Calculate local sidereal time
	    dlongs = dlong*dtr
	    call sidtim (djd,dlongs,dst)

	# Precess r.a. and dec. to mean equator and equinox of date (deqt)
	    deqt = (djd - dct0 - dcbes) / dctrop + dc1900
	    dras = ra * 15.0d0 * dtr
	    decs = dec * dtr
	    dc(1) = dcos(dras) * dcos(decs)
	    dc(2) = dsin(dras) * dcos(decs)
	    dc(3) =              dsin(decs)
	    call pre (eq,deqt,dprema)
	    do k = 1, 3 {
		dcc[k]=dc[1]*dprema[k,1]+dc[2]*dprema[k,2]+dc[3]*dprema[k,3]
		}

	    dra2 = datan2 (dcc[1],dcc[2])
	    dec2 = dasin(dcc(3])

	# Calculate hour angle = local sidereal time - r.a.
	    dha = dst - dra2
	    dha = dmod (dha + 2.0d0*dpi , 2.0d0*dpi)

	# Calculate observer's geocentric velocity
	# (altitude assumed to be zero)
	    dlats = dlat*dtr
	    call geovel (dlats,dalt,dec2,-dha,dgcvel)

	# Calculate components of earth's barycentric veolcity,
	# dvelb(i), i=1,2,3  in units of a.u./s
	    call barvel (djd,deqt,dvelh,dvelb)

	# Project barycentric velocity to the direction of the star, and
	# convert to km/s
	    dbcvel = 0.0d0
	    dhcvel = 0.0d0
	    do k = 1, 3 {
		dbcvel = dbcvel + dvelb[k]*dcc[k]*daukm
		dhcvel = dhcvel + dvelh[k]*dcc[k]*daukm
		}
	    if (verbose) {
		call printf ("bcvel = %.4f  hcvel = %.4f  geovel = %.4f\n")
		    call pargd (dbcvel)
		    call pargd (dhcvel)
		    call pargd (dgcvel)
		}

	# Add up both corrections
	    rbcv = dbcvel + dgcvel
	    rhcv = dhcvel + dgcvel
	    if (verbose) {
		call printf ("bcv = %.4f  hcv = %.4f\n")
		    call pargr (rbcv)
		    call pargr (rhcv)
		}

# Check for writability if results are being saved in image header
	    if (savebcv) {
		if (imaccess (specpath, READ_WRITE) == NO) {
		call eprintf ("BCVCORR: cannot write to %s; not saving results\n")
		    call pargstr (specpath)
		    IM_UPDATE(specim) = NO
		    }
		else {
		    call imaddr (specim,"BCV     ",rbcv)
		    IM_UPDATE(specim) = YES
		    }
		}
	    else
		IM_UPDATE(specim) = NO
	    }

# Close current image and move on to next image
	if (specim != ERR && specim != NULL)
	    call imunmap (specim)
	go to newspec_

#  Close spectrum list
endvc_	call imtclose (speclist)

end