diff options
Diffstat (limited to 'noao/rv/rvidlines/idvelocity.x')
| -rw-r--r-- | noao/rv/rvidlines/idvelocity.x | 188 |
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 |