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
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
|
include <gset.h>
include "rvpackage.h"
include "rvflags.h"
# RV_FIT - Fit the CCF in the specified region. Return the exact pixel
# shift and sigma of the fit.
procedure rv_fit (rv, xcf, ycf, ledge, redge, npts, ishift, shift, sigma)
pointer rv #I RV struct pointer
real xcf[ARB], ycf[ARB] #I Array of correlation peaks
int ledge #I Left edge to fit
int redge #I Right edge to fit
int npts #I Npts between edges
int ishift #I Initial index of shift
real shift #O Computed shift
real sigma #O Sigma of fit
pointer gp
real hght, init, peak, c[4]
real a, b, thresh, fwhm
int i, tnum
real rv_width(), center1d()
int rv_getshift()
include "fitcom.com"
include "rvsinc.com"
define NPARS 4
begin
# Erase the old fit first.
call rv_erase_fit (rv, false)
RV_FITDONE(rv) = NO
RV_ERRCODE(rv) = OK
IS_DBLSTAR(rv) = NO
# Do some window bounds checking.
if (ledge < (RV_WINCENTER(rv)-RV_WINDOW(rv)) ||
redge > (RV_WINCENTER(rv)+RV_WINDOW(rv))) {
call rv_err_comment (rv,
"WARNING: Some points in fit are outside window bounds.", "")
if (RV_INTERACTIVE(rv) == YES) {
call rv_errmsg (
"Warning: Some points in fit are outside window bounds.")
call tsleep (1)
}
}
# Save some info
RV_ISHIFT(rv) = ishift
RV_ISTART(rv) = ledge
RV_IEND(rv) = redge
gp = RV_GP(rv)
# Initialize some variables
tnum = RV_TEMPNUM(rv)
# Do the fitting
switch (RV_FITFUNC(rv)) {
case GAUSSIAN: # call gaussian fitting
call rv_fgauss (rv, xcf, ycf, ledge, redge, npts, ishift, c, sigma)
if (RV_ERRCODE(rv) == ERR_FIT) {
if (c[3] <= 0.0) {
return
} else {
call rv_draw_fit (rv, gp, NO)
return
}
}
shift = c[2]
if (IS_INDEF(RV_BACKGROUND(rv)))
hght = c[1] + c[4]
else
hght = c[1] + RV_BACKGROUND(rv)
case LORENTZIAN: # call lorentzian fitting
call rv_fgauss (rv, xcf, ycf, ledge, redge, npts, ishift, c, sigma)
#c[3] = abs (c[3])
if (RV_ERRCODE(rv) == ERR_FIT) {
call rv_draw_fit (rv, gp, NO)
return
}
shift = c[2]
call lorentz (shift, 3, c, 4, hght)
hght = 2.0 * c[1] / c[3]
if (!IS_INDEF(RV_BACKGROUND(rv)))
hght = hght + RV_BACKGROUND(rv)
case PARABOLA: # call parabola fitting routine
call rv_fparab (rv, xcf, ycf, ledge, redge, npts, ishift, c, sigma)
shift = -c[2] / (2. * c[3])
hght = c[1] + (c[2] * shift) + (c[3] * shift * shift)
peak = c[1] - c[2]**2 / (4. * c[3])
RV_FWHM(rv) = sqrt ( abs(-peak / (2. * c[3])))
case CENTER1D:
init = real (rv_getshift (ycf[ledge], npts, MAXIMUM))
call alimr (ycf[ledge], npts, a, b)
#thresh = (b - a) - 0.01
thresh = 0.0
peak = center1d (init, ycf[ledge], npts, npts/5., 1, 2., thresh)
RV_HEIGHT(rv) = INDEF
RV_FWHM(rv) = INDEF
RV_ERROR(rv) = INDEFD
RV_R(rv) = INDEF
RV_DISP(rv) = INDEF
if (RV_GP(rv) != NULL && RV_INTERACTIVE(rv) == YES) {
call gpmark (gp, xcf[ledge], ycf[ledge], npts, 4, 2., 2.)
call gflush (gp)
}
if (IS_INDEF(peak)) { # check for an error
RV_ERRCODE(rv) = ERR_FIT
return
} else {
shift = peak + (xcf[ledge] - 1.0)
if (RV_GP(rv) != NULL && RV_INTERACTIVE(rv) == YES) {
call gline (gp, shift, ycf[ishift], shift, ycf[ishift]+0.1)
call gflush (gp)
}
}
RV_FITDONE(rv) = YES
return
case SINC:
call rv_sinc (rv, shift, fwhm, hght)
RV_FWHM(rv) = fwhm
RV_SHIFT(rv) = shift
RV_ERROR(rv) = INDEFD
RV_R(rv) = INDEF
RV_DISP(rv) = INDEF
if (RV_GP(rv) != NULL && RV_INTERACTIVE(rv) == YES)
call gpmark (gp, xcf[ledge], ycf[ledge], npts, 4, 2., 2.)
RV_FITDONE(rv) = YES
if (IS_INDEF(fwhm)) { # no fwhm computed, so leave here
call rv_draw_fit (rv, gp, NO)
return
}
}
call amovr (c, COEFF(rv,1), NPARS)
RV_FWHM(rv) = rv_width (rv)
RV_HEIGHT(rv) = hght
# Redraw the new points fit if they were changed. Also save new fit
# window parameters
if (ledge != RV_ISTART(rv) && RV_FITFUNC(rv) != SINC) {
if (gp != NULL && RV_INTERACTIVE(rv) == YES) {
call gseti (gp, G_PMLTYPE, GL_CLEAR)
call gpmark (gp, xcf[RV_ISTART(rv)], ycf[RV_ISTART(rv)],
(RV_IEND(rv)-RV_ISTART(rv)+1), 4, 2., 2.)
call gseti (gp, G_PMLTYPE, GL_SOLID)
call gpmark (gp, xcf[ledge], ycf[ledge], npts, 4, 2., 2.)
call gpline (gp, xcf[RV_ISTART(rv)], ycf[RV_ISTART(rv)],
(RV_IEND(rv)-RV_ISTART(rv)+1))
call gflush (gp)
}
RV_ISTART(rv) = ledge
RV_IEND(rv) = ledge + npts
}
nfit = npts
# Compute the antisymmetric part of correlation and velocity error
call realloc (RV_ANTISYM(rv), RV_CCFNPTS(rv), TY_REAL)
if (!IS_INDEF(RV_FWHM(rv))) {
call rv_antisym (rv, shift, hght, RV_FWHM(rv), ycf, RV_CCFNPTS(rv),
ANTISYM(rv,1), ccfvar, RV_ERROR(rv), RV_R(rv))
} else {
RV_R(rv) = INDEF
if (RV_DCFLAG(rv) != -1)
RV_ERROR(rv) = sigma * RV_DELTAV(rv)
else
RV_ERROR(rv) = sigma
}
# Now get the dispersion of the peak
if (RV_DCFLAG(rv) != -1 && !IS_INDEF(RV_FWHM(rv)))
RV_DISP(rv) = RV_FWHM(rv) * RV_DELTAV(rv)
else
RV_DISP(rv) = INDEF
# Debugging info
if (DBG_DEBUG(rv) == YES && DBG_LEVEL(rv)>=2) {
call d_printf(DBG_FD(rv), "rvfitfunc:\n")
call d_printf(DBG_FD(rv), "\tledge=%d redge=%d npts=%d ishift=%d\n")
call pargi(ledge); call pargi(redge)
call pargi(npts); call pargi(ishift)
call d_printf(DBG_FD(rv),
"\tshift=%.4g sigma=%.4g fwhm=%.4g disp=%.4g hght=%.4g peak=%.4g\n")
call pargr(shift); call pargr(sigma); call pargr(RV_FWHM(rv))
call pargr(RV_DISP(rv));call pargr(hght); call pargr(peak)
do i = 1, NPARS {
call d_printf (DBG_FD(rv), "\t c[%d]=%g +/- %g\n")
call pargi(i); call pargr(c[i]); call pargr (ECOEFF(rv,i))
}
call flush (DBG_FD(rv))
}
# Put stuff in the common for the log
binshift = xcf[ishift]
if (RV_ERRCODE(rv) == OK) {
RV_FITDONE(rv) = YES
RV_UPDATE(rv) = YES
}
# Plot the computed fit.
call rv_draw_fit (rv, gp, NO)
# Mark the background level.
if (RV_FITFUNC(rv) == GAUSSIAN || RV_FITFUNC(rv) == LORENTZIAN) {
if (RV_INTERACTIVE(rv) == YES)
call rv_draw_background (rv, gp)
}
end
# RV_WIDTH - Procedure to compute the width of the CCF.
real procedure rv_width (rv)
pointer rv #I RV struct pointer
real fwhm, h, l, r, peak, shift
int gstati()
include "fitcom.com"
begin
# Now correct it for a fixed baseline
switch (RV_FITFUNC(rv)) {
case GAUSSIAN:
fwhm = sqrt (COEFF(rv,3)) * 2.35482
l = COEFF(rv,2) - fwhm / 2.
r = COEFF(rv,2) + fwhm / 2.
call cgauss1d (l, 1, COEFF(rv,1), nfitpars, h)
case LORENTZIAN:
fwhm = 2. * abs (COEFF(rv,3))
fwhm = abs (COEFF(rv,3)) # for new lorentzian
l = COEFF(rv,2) - fwhm / 2.
r = COEFF(rv,2) + fwhm / 2.
call lorentz (l, 1, COEFF(rv,1), nfitpars, h)
case PARABOLA:
peak = COEFF(rv,1) - COEFF(rv,2)**2 / (4. * COEFF(rv,3))
fwhm = sqrt ( abs(- peak / (2. * COEFF(rv,3))))
#return (fwhm)
shift = -COEFF(rv,2) / (2.*COEFF(rv,3))
l = shift - fwhm / 2.
r = shift + fwhm / 2.
h = COEFF(rv,1) + COEFF(rv,2) * l + COEFF(rv,3) * l**2
case CENTER1D:
return (INDEF)
case SINC:
# The structure parameters were computed before, we just need
# this to draw the marker.
l = RV_SHIFT(rv) - RV_FWHM(rv) / 2.
r = RV_SHIFT(rv) + RV_FWHM(rv) / 2.
h = RV_FWHM_Y(rv)
fwhm = RV_FWHM(rv)
}
RV_FWHM_Y(rv) = h
# Now draw the line showing the width
if (RV_GP(rv) != NULL && RV_INTERACTIVE(rv) == YES) {
if (gstati(RV_GP(rv),G_PLCOLOR) != C_BACKGROUND) {
call gseti (RV_GP(rv), G_PLTYPE, GL_DASHED)
call gseti (RV_GP(rv), G_PLCOLOR, RV_LINECOLOR(rv))
}
call gline (RV_GP(rv), l, h, r, h)
if (gstati(RV_GP(rv),G_PLCOLOR) != C_BACKGROUND) {
call gseti (RV_GP(rv), G_PLTYPE, GL_SOLID)
call gseti (RV_GP(rv), G_PLCOLOR, C_FOREGROUND)
}
call gflush (RV_GP(rv))
}
return (fwhm)
end
# RV_XFIT - Set the fitting endpoints as described for the 'g' keystroke
# command.
procedure rv_xfit (rv, x, do_correction)
pointer rv #I RV struct pointer
real x #I Current cursor x position
int do_correction #I Do heliocentric correction?
real sregion, eregion, y
real shift, sigma
int istart, iend, ishift, npts, stat
int rv_getshift(), rv_rvcorrect()
include "fitcom.com"
begin
sregion = x # get endpoints
call rv_getpts (rv, eregion, y, 2)
npts = RV_CCFNPTS(rv) # Fit the region
call rv_fixx (sregion, eregion, WRKPIXX(rv,1), WRKPIXX(rv,npts))
istart = int (npts/2 + 1 + sregion)
iend = int (npts/2 + 1 + eregion)
npts = int (iend - istart + 1)
nfit = npts
ishift = rv_getshift (WRKPIXY(rv,istart), npts, MAXIMUM)
# now jump into the fitting routines
call rv_fit (rv, WRKPIXX(rv,1), WRKPIXY(rv,1), istart, iend, npts,
ishift+istart-1, shift, sigma)
if (RV_ERRCODE(rv) == ERR_FIT) {
if (RV_INTERACTIVE(rv) == YES)
call rv_errmsg ("Fit did not converge")
else
call rv_err_comment (rv, "Fit did not converge", "")
return
}
RV_SHIFT(rv) = shift
RV_SIGMA(rv) = sigma
if (do_correction == YES) {
stat = rv_rvcorrect (rv, shift, sigma, RV_VOBS(rv), RV_VCOR(rv),
RV_ERROR(rv))
if (stat != OK) {
call rv_err_comment (rv,
"WARNING: Heliocentric correction not done properly.", "")
}
if (RV_INTERACTIVE(rv) == YES)
call rv_writeln (rv, STDOUT)
}
end
# RV_YFIT - Fit the CCF based on the Y value of the cursor, as described
# for the 'y' keystroke command or the HEIGHT parameter.
procedure rv_yfit (rv, y, do_correction)
pointer rv #I RV struct pointer
real y #I Current Y cursor value
int do_correction #I Do heliocentric correction?
real sregion, eregion
real shift, sigma, center
int istart, iend, ishift, npts, stat, i
int rv_getshift(), rv_rvcorrect()
include "fitcom.com"
begin
# Search the array for the closest points in y
npts = RV_WINR(rv) - RV_WINL(rv) + 1
center = RV_CCFNPTS(rv)/2 + 1 + WRKPIXX(rv,RV_WINCENTER(rv))
i = int (center - RV_WINDOW(rv))
ishift = rv_getshift (WRKPIXY(rv,i), npts, MAXIMUM) + i - 1
i = 0
while (WRKPIXY(rv,ishift-i) > y && i <= npts) {
sregion = WRKPIXX(rv, ishift-i)
i = i + 1
}
i = 0
while (WRKPIXY(rv, ishift+i) > y && i <= npts) {
eregion = WRKPIXX(rv, ishift+i)
i = i + 1
}
# Pick up at fitting routines
npts = RV_CCFNPTS(rv)
istart = int (npts/2 + 1 + sregion)
iend = int (npts/2 + 1 + eregion)
npts = (iend - istart + 1)
# Do the minwidth/maxwidth/window constraints
call rv_fix_window (rv, 1., real(RV_CCFNPTS(rv)), y, WRKPIXX(rv,1),
WRKPIXY(rv,1), istart, iend, ishift, npts)
# Go ahead and fit this puppy
call rv_fit (rv, WRKPIXX(rv,1), WRKPIXY(rv,1), istart, iend,
npts, ishift, shift, sigma)
if (RV_ERRCODE(rv) == ERR_FIT) {
if (RV_INTERACTIVE(rv) == YES)
call rv_errmsg ("Fit did not converge")
else
call rv_err_comment (rv, "Fit did not converge", "")
return
}
RV_SHIFT(rv) = shift
RV_SIGMA(rv) = sigma
if (do_correction == YES) {
stat = rv_rvcorrect (rv, shift, sigma, RV_VOBS(rv), RV_VCOR(rv),
RV_ERROR(rv))
if (RV_INTERACTIVE(rv) == YES)
call rv_writeln (rv, STDOUT)
}
end
# RV_FIX_WINDOW - Resize the fit window according to the minwidth/maxwidth
# constraint parameters. This routine also recenters the window on the
# initial guess at the shift so points are evenly spaced about the peak.
# Does a bounds check to avoid segmentation violations.
procedure rv_fix_window (rv, x1, x2, y, xcf, ycf, istart, iend, ishift, npts)
pointer rv #I RV struct pointer
real x1, x2 #I Bounds check
real y #I Threshold level
real xcf[ARB], ycf[ARB] #I CCF array
int istart #U Start pixel of fit
int iend #U End pixel of fit
int ishift #U Peak pixel of fit
int npts #U Npts in between
int i, np1, np2
begin
if (npts < RV_MINWIDTH(rv)) {
np1 = RV_MINWIDTH(rv) - npts # Pad some points
istart = istart - (np1 / 2)
iend = iend + (np1 / 2)
if (mod(np1,2) == 1)
iend = iend + 1
} else if (npts > RV_MAXWIDTH(rv)) {
np1 = npts - RV_MAXWIDTH(rv) # Delete some points
istart = istart + (np1 / 2)
iend = iend - (np1 / 2)
if (mod(np1,2) == 1)
iend = iend - 1
}
npts = int (iend - istart + 1)
# Next, we have to make sure we honor the original constraint that
# all points are above a certain level
if (npts > RV_MINWIDTH(rv)) {
i = istart
while (ycf[i] < y && i <= ishift) # Fix left side
i = i + 1
if (i != istart)
istart = i #+ 1
i = iend
while (ycf[i] < y && i >= ishift) # Fix right side
i = i - 1
if (i != iend)
iend = i #- 1
}
npts = int (iend - istart + 1)
# Now recenter the window on the peak
np1 = ishift - istart
np2 = iend - ishift
if ((np1 - np2) < -1) { # Peak is left of center
istart = istart - abs(np1 - np2) / 2
iend = iend - abs(np1 - np2) / 2
} else if ((np1 - np2) > 1) { # Peak is right of center
istart = istart + abs(np1 - np2) / 2
iend = iend + abs(np1 - np2) / 2
}
npts = int (iend - istart + 1)
# Lastly, make sure we aren't out of bounds after all this work
if (istart < x1) {
np1 = (x1 - istart)
istart = int (x1)
iend = iend + np1
}
if (iend > x2) {
np1 = (iend - x2)
iend = int (x2)
istart = istart - np1
}
npts = int (iend - istart + 1) # Update npts and return
end
|
