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
|
include <gset.h>
include "daoedit.h"
define FRACTION 0.10
# DP_ERPLOT -- Plot the radial profile.
procedure dp_erplot (gd, title, xwcs, ywcs, radius, intensity, npts,
rcentroid, pmean, integral, nbins, rmin, rmax, iannulus, oannulus,
apradius, skyval, skysigma, rscale, pnorm)
pointer gd # pointer to the graphics stream
char title[ARB] # the plot title
int xwcs # the x wcs type of the final plot
int ywcs # the y wcs type of the final plot
real radius[ARB] # the radius vector
real intensity[ARB] # the intensity vector
int npts # number of points in the profile
real rcentroid[ARB] # the radius centroid vector
real pmean[ARB] # the mean intensity vector
real integral[ARB] # the integral of the profile
int nbins # the number of bins
real rmin, rmax # min and max radius
real iannulus # the inner radius of the sky annulus
real oannulus # the outer radius of the sky annulus
real apradius # the aperture radius
real skyval # the sky value
real skysigma # the sigma of the sky value
real rscale # the image scale
real pnorm # the profile normalization factor
int i, j
pointer sp, pxlabel, pylabel, sxlabel, sylabel
real r1, r2, rp1, rp2, rs1, rs2, i1, i2, ip1, ip2, is1, is2, dr, fraction
begin
# Get space for the x and y labels.
call smark (sp)
call salloc (pxlabel, SZ_FNAME, TY_CHAR)
call salloc (pylabel, SZ_FNAME, TY_CHAR)
call salloc (sxlabel, SZ_FNAME, TY_CHAR)
call salloc (sylabel, SZ_FNAME, TY_CHAR)
# Clear the plot.
call gclear (gd)
# Determine the data range of the x and y axes.
r1 = rmin - FRACTION * (rmax - rmin)
r2 = rmax + FRACTION * (rmax - rmin)
switch (xwcs) {
case WCS_XPIX:
rp1 = r1
rp2 = r2
rs1 = r1 / rscale
rs2 = r2 / rscale
case WCS_XSCALE:
rp1 = r1 / rscale
rp2 = r2 / rscale
rs1 = r1
rs2 = r2
default:
rp1 = r1
rp2 = r2
rs1 = r1 / rscale
rs2 = r2 / rscale
}
fraction = max (FRACTION, 5.0 * skysigma / pnorm)
i1 = -pnorm * fraction
i2 = pnorm * (1.0 + fraction)
switch (ywcs) {
case WCS_YNORM:
ip1 = i1 / pnorm
ip2 = i2 / pnorm
is1 = i1 + skyval
is2 = i2 + skyval
case WCS_YCOUNT:
ip1 = i1 + skyval
ip2 = i2 + skyval
is1 = i1 / pnorm
is2 = i2 / pnorm
default:
ip1 = i1 / pnorm
ip2 = i2 / pnorm
is1 = i1 + skyval
is2 = i2 + skyval
}
# Draw the axes and axis labels.
call gsetr (gd, G_ASPECT, 0.0)
switch (xwcs) {
case WCS_XPIX:
call strcpy ("Radius (b=pixels, t=scale units)", Memc[pxlabel],
SZ_FNAME)
call strcpy ("", Memc[sxlabel], SZ_FNAME)
case WCS_XSCALE:
call strcpy ("Radius (b=scale units, t=pixels)", Memc[pxlabel],
SZ_FNAME)
call strcpy ("", Memc[sxlabel], SZ_FNAME)
default:
call strcpy ("Radius (b=pixels, t=scale units)", Memc[pxlabel],
SZ_FNAME)
call strcpy ("", Memc[sxlabel], SZ_FNAME)
}
switch (ywcs) {
case WCS_YCOUNT:
call strcpy ("Counts", Memc[pylabel], SZ_FNAME)
call strcpy ("Norm Counts", Memc[sylabel], SZ_FNAME)
case WCS_YNORM:
call strcpy ("Norm Counts", Memc[pylabel], SZ_FNAME)
call strcpy ("Counts", Memc[sylabel], SZ_FNAME)
default:
call strcpy ("Counts", Memc[pylabel], SZ_FNAME)
call strcpy ("Norm Counts", Memc[sylabel], SZ_FNAME)
}
call gseti (gd, G_XDRAWAXES, 1)
call gseti (gd, G_YDRAWAXES, 1)
call gswind (gd, rp1, rp2, ip1, ip2)
call glabax (gd, title, Memc[pxlabel], Memc[pylabel])
call gseti (gd, G_XDRAWAXES, 2)
call gseti (gd, G_YDRAWAXES, 2)
call gswind (gd, rs1, rs2, is1, is2)
call glabax (gd, title, Memc[sxlabel], Memc[sylabel])
# Plot the data points.
call gswind (gd, r1, r2, i1, i2)
call gpmark (gd, radius, intensity, npts, GM_PLUS, 1.0, 1.0)
# Plot the smoothed radial profile skipping any points with no data.
call gswind (gd, r1, r2, -fraction, 1.0 + fraction)
for (i = 1; i <= nbins && rcentroid[i] <= 0.0; i = i + 1)
;
call gamove (gd, rcentroid[i], pmean[i])
do j = i, nbins {
if (pmean[j] == 0.0)
next
call gadraw (gd, rcentroid[j], pmean[j])
}
# Plot the integral.
call gswind (gd, r1, r2, -fraction, 1.0 + fraction)
call gamove (gd, 0.0, 0.0)
dr = (rmax - rmin) / real (nbins - 1)
do j = 2, nbins {
call gadraw (gd, real (j - 1) * dr, integral[j])
}
# Plot the sky annuli and the aperture radius.
call gamove (gd, iannulus, -fraction)
call gadraw (gd, iannulus, 1.0 + fraction)
call gamove (gd, oannulus, -fraction)
call gadraw (gd, oannulus, 1.0 + fraction)
call gamove (gd, apradius, -fraction)
call gadraw (gd, apradius, 1.0 + fraction)
# Plot the zero level.
call gswind (gd, rp1, rp2, ip1, ip2)
switch (ywcs) {
case WCS_YCOUNT:
call gamove (gd, rp1, skyval)
call gadraw (gd, rp2, skyval)
case WCS_YNORM:
call gamove (gd, rp1, 0.0)
call gadraw (gd, rp2, 0.0)
default:
call gamove (gd, rp1, 0.0)
call gadraw (gd, rp2, 0.0)
}
call sfree (sp)
end
|
