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
|
include <math.h>
# IMG_AIRMASS -- Get or compute the image airmass from the image header.
# If the airmass cannot be determined from header then INDEF is returned.
#
# Airmass formulation from Allen "Astrophysical Quantities" 1973 p.125,133
# and John Ball's book on Algorithms for the HP-45.
real procedure img_airmass (im)
pointer im # IMIO pointer
real airmass, zd, ha, ra, dec, st, latitude, coszd, scale, x
int imaccf()
real imgetr()
errchk imgetr()
data scale/750.0/ # Atmospheric scale height approx
begin
# If the airmass is in the header return its value.
if (imaccf (im, "airmass") == YES)
return (imgetr (im, "airmass"))
# Compute zenith distance if not defined.
iferr (zd = imgetr (im, "zd")) {
# Compute hour angle if not defined.
iferr (ha = imgetr (im, "ha")) {
st = imgetr (im, "st")
ra = imgetr (im, "ra")
ha = st - ra
call imaddr (im, "ha", ha)
}
dec = imgetr (im, "dec")
latitude = imgetr (im, "latitude")
ha = DEGTORAD (ha) * 15
dec = DEGTORAD (dec)
latitude = DEGTORAD (latitude)
coszd = sin (latitude) * sin (dec) +
cos (latitude) * cos (dec) * cos (ha)
zd = RADTODEG (acos (coszd))
call imaddr (im, "zd", zd)
}
# Compute airmass from zenith distance.
zd = DEGTORAD (zd)
x = scale * cos (zd)
airmass = sqrt (x ** 2 + 2 * scale + 1) - x
call imaddr (im, "airmass", airmass)
return (airmass)
end
|
