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
|
# MKCWWCS -- MaKe Celestial, Wavelength 3D World Coordinate System
procedure mkcwwcs (wcsname)
file wcsname {prompt="WCS to create"}
file wcsref = "" {prompt="WCS reference\n"}
real equinox = INDEF {prompt="Equinox (years)"}
real ra = INDEF {prompt="RA (hours)"}
real dec = INDEF {prompt="DEC (degrees)"}
real scale = INDEF {prompt="Celestial pixel scale (arcsec/pix)"}
real pa = 0. {prompt="Position angle (deg)"}
bool lefthanded = yes {prompt="Left-handed system?"}
string projection = "tan" {prompt="Celestial projection\n",
enum="linear|tan|sin"}
real wave = INDEF {prompt="Wavelength"}
real wscale = INDEF {prompt="Wavelength scale\n"}
real rapix = INDEF {prompt="RA reference pixel"}
real decpix = INDEF {prompt="DEC reference pixel"}
real wpix = INDEF {prompt="Wavelength reference pixel"}
begin
int wcsdim = 3
real c, s, lh
file name, ref, wcs
# Determine the input and reference images.
name = wcsname
if (fscan (wcsref, ref) > 0)
wcscopy (name, ref)
# Set the axes.
wcsedit (name, "axtype", "ra", "1", wcsdim=wcsdim,
wcs="world", interactive-, verbose-, update+)
wcsedit (name, "axtype", "dec", "2", wcsdim=wcsdim,
wcs="world", interactive-, verbose-, update+)
wcsedit (name, "wtype", projection, "1,2", wcsdim=wcsdim,
wcs="world", interactive-, verbose-, update+)
# Set the celestial equinox if desired. Note this is not WCS.
if (equinox != INDEF)
hedit (name, "equinox", equinox,
add+, addonly-, verify-, show-, update+)
# Set the reference point if desired.
if (ra != INDEF)
wcsedit (name, "crval", ra*15, "1", wcsdim=wcsdim,
wcs="world", interactive-, verbose-, update+)
if (dec != INDEF)
wcsedit (name, "crval", dec, "2", wcsdim=wcsdim,
wcs="world", interactive-, verbose-, update+)
if (wave != INDEF)
wcsedit (name, "crval", wave, "3", wcsdim=wcsdim,
wcs="world", interactive-, verbose-, update+)
# Set the scales and celestial position angle.
if (scale != INDEF) {
if (pa != INDEF) {
c = cos (pa * 3.14159 / 180.) / 3600.
s = sin (pa * 3.14159 / 180.) / 3600.
} else {
c = 1.
s = 0.
}
if (lefthanded) {
wcsedit (name, "cd", -scale*c, "1", "1", wcsdim=wcsdim,
wcs="world", interactive-, verbose-, update+)
wcsedit (name, "cd", -scale*s, "1", "2", wcsdim=wcsdim,
wcs="world", interactive-, verbose-, update+)
wcsedit (name, "cd", -scale*s, "2", "1", wcsdim=wcsdim,
wcs="world", interactive-, verbose-, update+)
wcsedit (name, "cd", scale*c, "2", "2", wcsdim=wcsdim,
wcs="world", interactive-, verbose-, update+)
} else {
wcsedit (name, "cd", scale*c, "1", "1", wcsdim=wcsdim,
wcs="world", interactive-, verbose-, update+)
wcsedit (name, "cd", -scale*s, "1", "2", wcsdim=wcsdim,
wcs="world", interactive-, verbose-, update+)
wcsedit (name, "cd", scale*s, "2", "1", wcsdim=wcsdim,
wcs="world", interactive-, verbose-, update+)
wcsedit (name, "cd", scale*c, "2", "2", wcsdim=wcsdim,
wcs="world", interactive-, verbose-, update+)
}
}
if (wscale != INDEF)
wcsedit (name, "cd", wscale, "3", "3", wcsdim=wcsdim,
wcs="world", interactive-, verbose-, update+)
# Set reference pixel if desired.
if (rapix != INDEF)
wcsedit (name, "crpix", rapix, "1", wcsdim=wcsdim,
wcs="world", interactive-, verbose-, update+)
if (decpix != INDEF)
wcsedit (name, "crpix", decpix, "2", wcsdim=wcsdim,
wcs="world", interactive-, verbose-, update+)
if (wpix != INDEF)
wcsedit (name, "crpix", wpix, "3", wcsdim=wcsdim,
wcs="world", interactive-, verbose-, update+)
end
|
