Initial commit

1 files changed, 188 insertions, 0 deletions

diff --git a/noao/rv/rvidlines/idvelocity.x b/noao/rv/rvidlines/idvelocity.x
new file mode 100644
index 00000000..c62e5c89
--- /dev/null
+++ b/noao/rv/rvidlines/idvelocity.x
@@ -0,0 +1,188 @@
+include	<smw.h>
+include	<units.h>
+include	"identify.h"
+
+
+# ID_VELOCITY -- Compute velocity.
+
+procedure id_velocity (id, interactive)
+
+pointer	id			# ID pointer
+int	interactive		# Called interactively?
+
+int	i, n
+double	z, sumz, sumz2, sumw, zerr, zhelio, v, verr, id_zval()
+
+begin
+	sumz = 0
+	sumw = 0
+	n = 0
+	for (i=1; i <= ID_NFEATURES(id); i = i + 1) {
+	    if (IS_INDEFD (USER(id,i)) || WTS(id,i) == 0.)
+		next
+	    z = id_zval (id, FIT(id,i), USER(id,i))
+	    sumz = sumz + WTS(id,i) * z
+	    sumw = sumw + WTS(id,i)
+	    n = n + 1
+	}
+
+	if (sumw > 0.) {
+	    zhelio = ID_ZHELIO(id)
+	    sumz = sumz / sumw
+
+	    sumz2 = 0.
+	    for (i=1; i <= ID_NFEATURES(id); i = i + 1) {
+		if (IS_INDEFD (USER(id,i)) || WTS(id,i) == 0.)
+		    next
+		z = id_zval (id, FIT(id,i), USER(id,i))
+		sumz2 = sumz2 + WTS(id,i) * (z - sumz) ** 2
+	    }
+	    if (sumz2 > 0.)
+		sumz2 = sqrt (sumz2 / sumw)
+	    else
+		sumz2 = 0.
+	    zerr = sumz2
+	    if (n > 1)
+		zerr = zerr / sqrt (n - 1.)
+
+	    if (interactive == YES) {
+		v = (sumz + zhelio) * VLIGHT
+		verr = zerr * VLIGHT
+
+		if (zhelio == 0D0)
+		    call printf (
+			"%s%s: Lines=%3d, Vobs=%.5g (%.5g), Zobs=%.5g (%.5g)\n")
+		else
+		    call printf (
+		    "%s%s: Lines=%3d, Vhelio=%.5g (%.5g), Zhelio=%.5g (%.5g)\n")
+
+		    call pargstr (Memc[ID_IMAGE(id)])
+		    call pargstr (Memc[ID_SECTION(id)])
+		    call pargi (n)
+		    call pargd (v)
+		    call pargd (verr)
+		    call pargd (sumz + zhelio)
+		    call pargd (zerr)
+	    }
+	    ID_REDSHIFT(id) = sumz
+	    ID_RMSRED(id) = sumz2
+	    ID_ZHELIO(id) = zhelio
+	}
+end
+
+
+# ID_ZVAL -- Compute Z value.
+
+double procedure id_zval (id, x, xref)
+
+pointer	id		#I Identify pointer
+double	x		#I Coordinate
+double	xref		#I Reference coordinate
+double	z		#O Z value
+
+double	y, yref
+pointer	un
+
+begin
+	y = x
+	yref = xref
+
+	un = UN(ID_SH(id))
+	if (UN_LOG(un) == YES) {
+	    y = 10D0 ** y
+	    yref = 10D0 ** yref
+	}
+	if (UN_INV(un) == YES) {
+	    y = 1D0 / y
+	    yref = 1D0 / yref
+	}
+
+	switch (UN_CLASS(un)) {
+	case UN_WAVE:
+	    z = (y - yref) / yref
+	case UN_FREQ, UN_ENERGY:
+	    z = (yref - y) / y
+	case UN_VEL:
+	    y = sqrt ((1 + y) / (1 - y))
+	    yref = sqrt ((1 + yref) / (1 - yref))
+	    z = (y - yref) / yref
+	case UN_DOP:
+	    y = y + 1
+	    yref = yref + 1
+	    z = (y - yref) / yref
+	}
+
+	return (z)
+end
+
+
+
+# ID_ZSHIFTD -- Shift coordinate by redshift.
+
+double procedure id_zshiftd (id, x, dir)
+
+pointer	id		#I Identify pointer
+double	x		#I Coordinate
+int	dir		#I Direction (0=to rest, 1=from rest)
+double	y		#O Shifted coordinate
+
+pointer	un
+
+begin
+	y = x
+
+	un = UN(ID_SH(id))
+	if (UN_LOG(un) == YES)
+	    y = 10D0 ** y
+	if (UN_INV(un) == YES)
+	    y = 1D0 / y
+
+	switch (UN_CLASS(un)) {
+	case UN_WAVE:
+	    if (dir == 0)
+		y = y / (1 + ID_REDSHIFT(id))
+	    else
+		y = y * (1 + ID_REDSHIFT(id))
+	case UN_FREQ, UN_ENERGY:
+	    if (dir == 0)
+		y = y * (1 + ID_REDSHIFT(id))
+	    else
+		y = y / (1 + ID_REDSHIFT(id))
+	case UN_VEL:
+	    y = sqrt ((1 + y) / (1 - y))
+	    if (dir == 0)
+		y = y / (1 + ID_REDSHIFT(id))
+	    else
+		y = y * (1 + ID_REDSHIFT(id))
+	    y = y ** 2
+	    y = (y - 1) / (y + 1)
+	case UN_DOP:
+	    y = (y + 1)
+	    if (dir == 0)
+		y = (y + 1) / (1 + ID_REDSHIFT(id)) - 1
+	    else
+		y = (y + 1) * (1 + ID_REDSHIFT(id)) - 1
+	}
+
+	if (UN_INV(un) == YES)
+	    y = 1D0 / y
+	if (UN_LOG(un) == YES)
+	    y = log10 (y)
+
+	return (y)
+end
+
+
+# ID_ZSHIFTR -- Shift coordinate by redshift.
+
+real procedure id_zshiftr (id, x, dir)
+
+pointer	id		#I Identify pointer
+real	x		#I Coordinate
+int	dir		#I Direction (0=to rest, 1=from rest)
+
+double	id_zshiftd()
+
+begin
+	return (real (id_zshiftd (id, double(x), dir)))
+end