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
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
|
include <mach.h>
include <math.h>
include "rvpackage.h"
include "rvflags.h"
# RV_SINC - Do the Fourier (sinc) interpolation to determine the peak center
# and FWHM values. Height of the ccf at the peak is also returned.
procedure rv_sinc (rv, shift, fwhm, height)
pointer rv #I RV struct pointer
real shift #O Shift of the peak
real fwhm #O FWHM of the peak
real height #O Height of peak at center
int i, j, k, il, ir
real x, y, back, lhp, rhp, hpower, ipeak
real brent(), sinc_interp(), rv_maxpix()
errchk realloc, mfree
include "rvsinc.com"
begin
# Initialize.
il = RV_ISTART(rv)
ir = RV_IEND(rv)
ipeak = WRKPIXX(rv,RV_ISHIFT(rv))
snfit = ir - il + 1
# Allocate the pointers in the common
call realloc (sx, snfit, TY_REAL)
call realloc (sy, snfit, TY_REAL)
call realloc (splx, snfit*10, TY_REAL)
call realloc (sply, snfit*10, TY_REAL)
# Now move the part of the ccf we're fitting into the arrays, but
# change the sign of the ccf because we're finding a minimum with
# the algorithm used
call amovr (WRKPIXX(rv,il), Memr[sx], snfit)
call amulkr (WRKPIXY(rv,il), -1.0, Memr[sy], snfit)
# Now find the peak center and height.
height = brent (ipeak-1., ipeak, ipeak+1., RV_TOLERANCE(rv),
RV_MAXITERS(rv), shift)
height = - height
# Compute the sinc interpolant over the ccf to see how close we came.
# It will be plotted later so we don't free the pointers right away.
do i = 1, snfit-1 {
do j = 1, 10 {
k = ((i-1)*10+j)
Memr[splx+k-1] = Memr[sx+i-1] + ((j-1)/10.)
Memr[sply+k-1] = - sinc_interp (Memr[splx+k-1])
}
}
height = rv_maxpix (Memr[sply], snfit*10)
# Compute the FWHM through some rather brute force methods.
back = RV_BACKGROUND(rv)
hpower = back + (height - back) / 2.
if (IS_INDEF(RV_BACKGROUND(rv))) {
fwhm = INDEF # don't compute a fwhm
RV_FWHM_Y(rv) = INDEF
} else if (WRKPIXY(rv,RV_ISTART(rv)) > hpower ||
WRKPIXY(rv,RV_IEND(rv)) > hpower) {
fwhm = INDEF # don't compute a fwhm
RV_FWHM_Y(rv) = INDEF
} else {
RV_FWHM_Y(rv) = hpower
for (i=RV_ISHIFT(rv); WRKPIXY(rv,i)>hpower&&i>=1; i=i-1)
x = WRKPIXX(rv,i)
for (y=-sinc_interp(x); abs(y-hpower)>0.005; x=x-0.005)
y = - sinc_interp (x-0.005);
lhp = x - 0.005
for (i=RV_ISHIFT(rv); WRKPIXY(rv,i)>hpower && i<=RV_CCFNPTS(rv);
i=i+1)
x = WRKPIXX(rv,i)
for (y=-sinc_interp(x); abs(y-hpower)>0.005; x=x+0.005)
y = - sinc_interp (x+0.005);
rhp = x + 0.005
fwhm = abs (rhp - lhp)
}
# Clean up a little
call mfree (sx, TY_REAL)
call mfree (sy, TY_REAL)
end
# SINC_INTERP - Function subroutine to do the sine (fourier) interpolation and
# return the value of the correlation function for any x value.
#
# h(t) = Sum(all n) [ h_n * sin(pi*(t-n)) /(pi*(t-n))]
#
# The interval between samples is assumed to be unity.
real procedure sinc_interp (x)
real x #I Point to be evaluated
real sval, tmp
int i
include "rvsinc.com"
begin
# Check for an integer x. If present, return the ccf value and don't
# interpolate.
tmp = abs (float(nint(x)) - x)
if (tmp < EPSILON && x >= Memr[sx] && x <= Memr[sx+snfit-1]) {
do i = 1, snfit {
if (abs(Memr[sx+i-1]-x) < EPSILON) # find the y point
sval = Memr[sy+i-1]
}
return (sval)
}
# Do the evaluation.
tmp = sin (PI*(x-Memr[sx]))
sval = Memr[sy] * tmp / (PI * (x - Memr[sx]))
do i = 1, snfit-1 {
tmp = -tmp
sval = sval + Memr[sy+i] * tmp / (PI * (x - Memr[sx+i]))
}
return (sval)
end
# BRENT - Given a function F(), and given a bracketing triplet of abscissas
# AX, BX, CX (such that BX is between AX and CX, and F(bx) is less than both
# F(AX) and F(BX)), this routine isolates the minimum to a fractional precision
# of about TOL using Brent's Method. The abscissa of the minimum is returned
# as XMIN, and the minimum function value is the function return value.
real procedure brent (ax, bx, cx, tol, itmax, xmin)
real ax, bx, cx #I Interp. bracketing points
real tol #I Tolerance
int itmax #I Max no. of iterations
real xmin #O Minimum point
real a, b, d, etemp, fu, fv, fw, fx # local variables
real p, q, r, tol1, tol2, u, v, w, x, xm
real e
int iter
real sinc_interp()
define CGOLD .3819660 # golden ration
define ZEPS 1.0e-10 # a small number
define SHIFT {$1=$2;$2=$3;$3=$4}
begin
a = min (ax, cx) # initialize
b = max (ax, cx)
v = bx
w = v
x = v
e = 0.
fx = sinc_interp (x)
fv = fx
fw = fx
e = 0.0
do iter = 1, itmax {
xm = 0.5 * (a + b)
tol1 = tol * abs (x) + ZEPS
tol2 = 2. * tol1
if (abs(x-xm) <= (tol2-0.5*(b-a))) { # test for convergence
xmin = x
return (fx)
}
if (abs(e) > tol1) { # construct a trial
r = (x-w) * (fx-fv) # parabolic fit
q = (x-v) * (fx-fw)
p = (x-v) * q - (x-w) * r
q = 2. * (q-r)
if (q > 0.)
p = -p
q = abs(q)
etemp = e
e = d
# Determine the acceptability of the parabolic fit. Here we
# take the golden section step into the larger of the two
# segments.
if (abs(p) >= abs(.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x)) {
if (x >= xm)
e = a - x
else
e = b - x
d = CGOLD * e
} else {
d = p/q # take parabolic step
u = x+d
if (u-a < tol2 || b-u < tol2)
d = sign (tol1, xm-x)
}
} else {
if (x >= xm)
e = a - x
else
e = b - x
d = CGOLD * e
}
if (abs(d) >= tol1)
u = x + d
else
u = x + sign (tol1,d)
fu = sinc_interp (u)
if (fu <= fx) {
if (u >= x)
a = x
else
b = x
SHIFT(v,w,x,u)
SHIFT(fv,fw,fx,fu)
} else {
if (u < x)
a = u
else
b = u
if (fu <= fw || w == x) {
v = w
fv = fw
w = u
fw = fu
} else if (fu <= fv || v == x || v == w) {
v = u
fv = fu
}
}
}
call rv_errmsg ("BRENT: Exceeded maximum number of iterations.\n")
xmin = x
return (fx)
end
|
