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
|
# File rvsao/Util/vcombine.x
# March 27, 2015
# By Jessica Mink, Harvard-Smithsonian Center for Astrophysics
# Copyright(c) 1991-2015 Smithsonian Astrophysical Observatory
# You may do anything you like with this file except remove this copyright.
# The Smithsonian Astrophysical Observatory makes no representations about
# the suitability of this software for any purpose. It is provided "as is"
# without express or implied warranty.
# Combine cross-correlation and emission line velocities using algorithm
# derived from John Tonry's SAO Nova software by Bill Wyatt
include "../lib/rvsao.h"
procedure vcombine (velxc,vxerr,vr,velem,veerr,nlfit,vel,verr,debug)
double velxc # Radial velocity based on best template correlation
double vxerr # Error in cross-correlation radial velocity
double vr # R-value of best cross-correlation
double velem # Radial velocity based on emission lines
double veerr # Error in emission line radial velocity
int nlfit # Number of emission lines in velocity fit
double vel # Heliocentric radial velocity of object (returned)
double verr # Error in radial velocity (returned)
bool debug # True for printed information
double disperr # Dispersion error in km/sec
double fc,vdiff
double emerr,emerr2
double xcerr,xcerr2
double err,err2
double rlim, dindef,disperr2
begin
dindef = INDEFD
rlim = 15.d0
disperr2 = 0.d0
# disperr = 15.d0
# if (disperr == 0.0 || disperr == dindef)
# disperr2 = 225.d0
# else
# disperr2 = disperr * disperr
if (velem == dindef && velxc == dindef) {
vel = dindef
verr = dindef
return
}
else if (velem == dindef || veerr <= 0. || nlfit == 0) {
if (velxc != dindef) {
vel = velxc
verr = sqrt ((vxerr*vxerr) + disperr2)
}
else {
vel = 0.d0
verr = 0.d0
}
return
}
else if (velxc == dindef || vxerr <= 0.) {
if (velem != dindef) {
vel = velem
verr = sqrt ((veerr*veerr) + disperr2)
}
else {
vel = 0.d0
verr = 0.d0
}
return
}
else {
emerr = veerr
emerr2 = emerr * emerr
xcerr = vxerr
xcerr2 = xcerr * xcerr
vdiff = velxc - velem
fc = vdiff * vdiff / (xcerr2 + emerr2)
if (fc > 8) {
if (vr > rlim) {
vel = velxc
err = vxerr
}
else if (nlfit < 3 && vr > 4.) {
vel = velxc
err = vxerr
}
else {
vel = velem
err = veerr
}
}
else {
if (xcerr2 > 0 && emerr2 > 0) {
err2 = 1.d0 / ((1.d0/emerr2) + (1.d0/xcerr2))
vel = ((velxc / xcerr2) + (velem / emerr2)) * err2
err = sqrt (err2)
}
else if (xcerr2 > 0) {
err = xcerr
vel = velxc
}
else if (emerr2 > 0) {
err = emerr
vel = velem
}
}
}
err2 = err * err
err = sqrt (err2 + disperr2)
verr = err
if (debug) {
call printf ("VCOMBINE: vel = %.3f +- %.3f km/sec\n")
call pargd (vel)
call pargd (verr)
call printf ("VCOMBINE: xc vel = %.3f +- %.3f km/sec, R = %.2f\n")
call pargd (velxc)
call pargd (xcerr)
call pargd (vr)
call printf ("VCOMBINE: em vel = %.3f +- %.3f km/sec, %d lines\n")
call pargd (velem)
call pargd (emerr)
call pargi (nlfit)
}
end
# Dec 5 1991 Set vel = velxc and verr = vxerr if no emission velocity
# Aug 11 1992 Use nlfit from getim.com instead of nfit from emv common
# Nov 24 1992 Avoid divide by zero errors
# Nov 30 1992 Drop emv header file
# Jun 23 1994 Pass velocities as arguments, not in labelled common
# Aug 3 1994 Change header from fquot to rvsao
# Aug 3 1994 Remove quadrature addition of 30 km/sec error to emission error
# Aug 18 1994 Clean up code
# Jul 13 1995 Change R-value for cross-correlation-only from 10 to 15
# Jul 13 1995 Add debugging code
# Jul 19 1995 Deal with INDEF velocities
# Mar 27 1997 Back R-value for cross-correlation-only back to 10
# Mar 27 1997 Back R-value for cross-correlation-only back to 10
# Mar 27 1997 Only take emission alone if R < 4 and more than 2 lines
# Feb 12 1998 If only error is zero, return zero velocity
# Dec 9 2008 Make disperr dispersion error explicit to fix systematic error
# Sep 20 2011 Set disperr2 to zero because dispersion error is not
# always known. Instrument errors should still be added
# to error
# Mar 27 2015 Link to header and common files in lib/
|
