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
|
include <mwset.h>
include <smw.h>
# SMW_MERGE -- Merge split MWCS array to a single MWCS.
procedure smw_merge (smw)
pointer smw #U Input split WCS, output single WCS
int i, pdim, naps, format, beam, dtype, dtype1, nw, nw1
int ap, axes[3]
double w1, dw, z, w11, dw1, z1
real aplow[2], aphigh[2]
pointer sp, key, val, term, coeff, mw, mw1, mw_open()
data axes/1,2,3/
begin
if (SMW_NMW(smw) == 1)
return
call smark (sp)
call salloc (key, SZ_FNAME, TY_CHAR)
call salloc (val, SZ_LINE, TY_CHAR)
call salloc (term, 15, TY_DOUBLE)
coeff = NULL
pdim = SMW_PDIM(smw)
naps = SMW_NSPEC(smw)
mw1 = SMW_MW(smw,0)
# Determine output WCS format.
format = SMW_ES
do i = 1, naps {
call smw_gwattrs (smw, i, 1, ap, beam, dtype, w1, dw, nw, z,
aplow, aphigh, coeff)
if (i == 1) {
dtype1 = dtype
w11 = w1
dw1 = dw
z1 = z
nw1 = nw
}
if (dtype>1||dtype!=dtype1||w1!=w11||dw!=dw1||nw!=nw1||z!=z1) {
format = SMW_MS
break
}
}
# Setup WCS.
switch (format) {
case SMW_ES:
mw = mw_open (NULL, pdim)
call mw_newsystem (mw, "equispec", pdim)
call mw_swtype (mw, axes, pdim, "linear", "")
case SMW_MS:
mw = mw_open (NULL, pdim)
call mw_newsystem (mw, "multispec", pdim)
call mw_swtype (mw, axes, pdim, "multispec", "")
if (pdim > 2)
call mw_swtype (mw, 3, 1, "linear", "")
}
ifnoerr (call mw_gwattrs (mw1, 1, "label", Memc[val], SZ_LINE))
call mw_swattrs (mw, 1, "label", Memc[val])
ifnoerr (call mw_gwattrs (mw1, 1, "units", Memc[val], SZ_LINE))
call mw_swattrs (mw, 1, "units", Memc[val])
ifnoerr (call mw_gwattrs (mw1, 1, "units_display", Memc[val], SZ_LINE))
call mw_swattrs (mw, 1, "units_display", Memc[val])
call mw_gltermd (mw1, Memd[term+pdim], Memd[term], pdim)
call mw_sltermd (mw, Memd[term+pdim], Memd[term], pdim)
call mw_gwtermd (mw1, Memd[term], Memd[term+pdim],
Memd[term+2*pdim], pdim)
Memd[term] = 1.
Memd[term+pdim] = w1 / (1 + z)
Memd[term+2*pdim] = dw / (1 + z)
call mw_swtermd (mw, Memd[term], Memd[term+pdim],
Memd[term+2*pdim], pdim)
# Set the SMW structure.
call smw_open (mw, smw, NULL)
if (format == SMW_MS) {
do i = 1, SMW_NMW(mw) - 1
call mw_close (SMW_MW(mw,i))
SMW_NMW(mw) = 1
}
do i = 1, naps {
call smw_gwattrs (smw, i, 1, ap, beam, dtype, w1, dw, nw, z,
aplow, aphigh, coeff)
call smw_swattrs (mw, i, 1, ap, beam, dtype, w1, dw, nw, z,
aplow, aphigh, Memc[coeff])
call smw_gapid (smw, i, 1, Memc[val], SZ_LINE)
call smw_sapid (mw, i, 1, Memc[val])
}
call smw_close (smw)
smw = mw
call mfree (coeff, TY_CHAR)
call sfree (sp)
end
|
