File: vcombine.x

package info (click to toggle)
iraf-rvsao 2.8.3-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 16,456 kB
  • sloc: ansic: 963; lisp: 651; fortran: 397; makefile: 27
file content (156 lines) | stat: -rw-r--r-- 4,295 bytes parent folder | download
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/