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
|
# PAIRMASS -- plot the airmass for a given RA & Dec on a given date.
procedure t_pairmass()
pointer sp, observat, timesys, ut, air, title
pointer obs
double ra, dec, epoch, ra0, dec0, epoch0
double longitude, latitude, zone, st, ha, utd, resolution
int day, month, year, nsteps, i, tsys
real amin, amax
pointer obsopen()
double clgetd(), obsgetd(), ast_mst(), airmass()
bool clgetb()
int clgeti(), clgwrd()
begin
call smark (sp)
call salloc (observat, SZ_FNAME, TY_CHAR)
call salloc (timesys, SZ_FNAME, TY_CHAR)
ra0 = clgetd ("ra")
dec0 = clgetd ("dec")
epoch0 = clgetd ("epoch")
year = clgeti ("year")
month = clgeti ("month")
day = clgeti ("day")
# Set time to plot.
tsys = clgwrd ("timesys", Memc[timesys], SZ_FNAME,
"|Universal|Standard|Siderial|")
# Get observatory information.
call clgstr ("observatory", Memc[observat], SZ_FNAME)
obs = obsopen (Memc[observat])
#call obslog (obs, "PAIRMASS", "latitude longitude timezone", STDOUT)
call obsgstr (obs, "name", Memc[observat], SZ_FNAME)
longitude = obsgetd (obs, "longitude")
latitude = obsgetd (obs, "latitude")
zone = obsgetd (obs, "timezone")
call obsclose (obs)
resolution = clgetd ("resolution")
nsteps = nint (24 * resolution) + 1
call salloc (ut, 3*nsteps, TY_REAL)
call salloc (air, 3*nsteps, TY_REAL)
call salloc (title, SZ_LINE, TY_CHAR)
call ast_date_to_epoch (year, month, day, 12.d0, epoch)
call ast_precess (ra0, dec0, epoch0, ra, dec, epoch)
do i = 1, nsteps {
utd = (i-1) / resolution
call ast_date_to_epoch (year, month, day, utd, epoch)
st = ast_mst (epoch, longitude)
ha = st - ra
switch (tsys) {
case 1:
Memr[ut+i-1] = utd
case 2:
if (utd < zone)
Memr[ut+i-1] = utd - zone + 24
else
Memr[ut+i-1] = utd - zone
case 3:
Memr[ut+i-1] = st
}
Memr[air+i-1] = real (airmass (ha, dec, latitude))
Memr[ut+i-1+nsteps] = Memr[ut+i-1] - 24
Memr[ut+i-1+2*nsteps] = Memr[ut+i-1] + 24
Memr[air+i-1+nsteps] = Memr[air+i-1]
Memr[air+i-1+2*nsteps] = Memr[air+i-1]
}
call xt_sort2 (Memr[ut], Memr[air], 3*nsteps)
call alimr (Memr[air], 3*nsteps, amin, amax)
call sprintf (Memc[title], SZ_LINE,
"Airmass for %d/%d/%d\n%s\nRA=%h, Dec=%h (%g)")
call pargi (month)
call pargi (day)
call pargi (year)
call pargstr (Memc[observat])
call pargd (ra0)
call pargd (dec0)
call pargd (epoch0)
if (clgetb ("listout")) {
call printf ("%s\nTime System=%s\n\n")
call pargstr (Memc[title])
call pargstr (Memc[timesys])
amax = clgetd ("wy2")
if (amax < 1)
amax = 5
do i = 1, 3*nsteps {
if (Memr[ut+i-1] < 0. || Memr[ut+i-1] >= 24.)
next
if (Memr[air+i-1] > amax)
next
call printf ("%6.0m\t%8.4f\n")
call pargr (Memr[ut+i-1])
call pargr (Memr[air+i-1])
}
} else
call draw_vector (Memc[title], Memc[timesys], Memr[ut], Memr[air],
3*nsteps, 0., 24., amin, amax)
call sfree (sp)
end
|
