diff options
Diffstat (limited to 'noao/rv/rvrebin.x')
| -rw-r--r-- | noao/rv/rvrebin.x | 155 |
1 files changed, 155 insertions, 0 deletions
diff --git a/noao/rv/rvrebin.x b/noao/rv/rvrebin.x new file mode 100644 index 00000000..4c0710f7 --- /dev/null +++ b/noao/rv/rvrebin.x @@ -0,0 +1,155 @@ +include <error.h> +include "rvpackage.h" +include "rvflags.h" +include "rvcont.h" + + +# FORCE_REBIN - Rebin the input spectrum in log space from a specified +# uniform log binning. This is to force both object and templates to the +# same dispersion if at all possible. + +int procedure force_rebin (rv) + +pointer rv #I RV struct pointer + +int npts, dcflag, which, status +int force_which(), rv_getim() +real w0, wpc, w2 +errchk realloc + +begin + # Skip this if not dispersion corrected. + if (RV_DCFLAG(rv) == -1) + return (OK) + + # First get the intermediate values to rebin the spectrum + dcflag = 1 # Rebinning to log dispersion. + if (force_which (rv, which, w0, wpc, npts) == NO) + return (OK) + + # Move the data before rebinning and increase the size of the data cache + if (which == OBJECT_SPECTRUM) { + call realloc (RV_OPIXX(rv), npts, TY_REAL) + call realloc (RV_OPIXY(rv), npts, TY_REAL) + call aclrr (OBJPIXX(rv,1), npts) + call aclrr (OBJPIXY(rv,1), npts) + } else { + call realloc (RV_RPIXX(rv), npts, TY_REAL) + call realloc (RV_RPIXY(rv), npts, TY_REAL) + call aclrr (REFPIXX(rv,1), npts) + call aclrr (REFPIXY(rv,1), npts) + } + + w2 = w0 + (npts - 1) * wpc + if (DEBUG(rv)) { + call d_printf (DBG_FD(rv), "force_rebin:\n") + call d_printf (DBG_FD(rv), "\tw0,wpc,w2,npts = %g,%g,%g,%d - %d\n") + call pargr(w0);call pargr(wpc);call pargr(w2) + call pargi(npts);call pargi(which) + call d_printf (DBG_FD(rv), "\to_w0=%g o_wpc=%g o_w2=%g\n") + call pargr(RV_OW0(rv)) ; call pargr(RV_OWPC(rv)) + call pargr(RV_OW2(rv)) + call d_printf (DBG_FD(rv), "\tr_w0=%g r_wpc=%g r_w2=%g\n") + call pargr(RV_RW0(rv)) ; call pargr(RV_RWPC(rv)) + call pargr(RV_RW2(rv)) + } + + # Now rebin to the desired dispersion. + #wpc = (10.0 ** (w0 + (npts-1)*wpc) - 10.0 ** w0) / (npts - 1) + w2 = 10.0 ** w2 + w0 = 10.0 ** w0 + if (which == OBJECT_SPECTRUM) + status = rv_getim (rv, IMAGE(rv), OBJECT_SPECTRUM, w0, w2, npts) + else + status = rv_getim (rv, RIMAGE(rv), REFER_SPECTRUM, w0, w2, npts) + + if (DEBUG(rv)) { + call d_printf (DBG_FD(rv), "\tend: w0,wpc,w2,npts = %g,%g,%g,%d\n") + call pargr(w0);call pargr(wpc);call pargr(w2);call pargi(npts) + } + + # Re-calculate the velocity dispersion and global endpoints + RV_DELTAV(rv) = wpc * CLN10 + RV_GLOB_W1(rv) = min (RV_OW0(rv), RV_RW0(rv)) + RV_GLOB_W2(rv) = max (RV_OW2(rv), RV_RW2(rv)) + + return (OK) +end + + +# FORCE_WHICH - Determine if spectra have to be rebinned to have the same +# dispersion and same starting wavelengths to within an integer shift. + +int procedure force_which (rv, which, w0, wpc, npts) + +pointer rv #I RV struct pointer +int which #O Which spectrum? +real w0 #O Rebinned W0 +real wpc #O Rebinned WPC +int npts #O Rebinned number of pixels + +int i +int stat +bool fp_equalr() + +begin + # Check if the spectra have the same dispersion. + which = OBJECT_SPECTRUM + w0 = RV_RW0(rv) + wpc = RV_RWPC(rv) + npts = RV_RNPTS(rv) + stat = NO + if (fp_equalr (w0, RV_OW0(rv)) && fp_equalr (wpc, RV_OWPC(rv))) + return (stat) + + # Determine spectrum to rebin. + switch (RV_REBIN(rv)) { + case RB_OBJ: + which = REFER_SPECTRUM + case RB_TEMP: + which = OBJECT_SPECTRUM + case RB_SMALL: + if (abs (RV_OWPC(rv)) < abs (RV_RWPC(rv))) + which = REFER_SPECTRUM + else + which = OBJECT_SPECTRUM + case RB_BIG: + if (abs (RV_OWPC(rv)) > abs (RV_RWPC(rv))) + which = REFER_SPECTRUM + else + which = OBJECT_SPECTRUM + } + + # Set the new dispersion parameters. The dispersion is the same as + # the target spectrum and the starting wavelength is adjusted by a + # fractional pixel amount into the data so that the starting + # wavelengths of the two spectra are an integer number of pixels + # apart in the common dispersion. + + switch (which) { + case REFER_SPECTRUM: + w0 = RV_OW0(rv) + wpc = RV_OWPC(rv) + if (RV_RWPC(rv) / RV_OWPC(rv) > 0) + i = (RV_RW0(rv) - w0) / wpc + 0.999 + else + i = (RV_RW0(rv) - w0) / wpc + w0 = w0 + i * wpc + npts = (RV_RW2(rv) - w0) / wpc + 1 + if (!fp_equalr (w0, RV_RW0(rv)) || !fp_equalr (wpc, RV_RWPC(rv))) + stat = YES + case OBJECT_SPECTRUM: + w0 = RV_RW0(rv) + wpc = RV_RWPC(rv) + if (RV_RWPC(rv) / RV_OWPC(rv) > 0) + i = (RV_OW0(rv) - w0) / wpc + 0.999 + else + i = (RV_OW0(rv) - w0) / wpc + w0 = w0 + i * wpc + npts = (RV_OW2(rv) - w0) / wpc + 1 + if (!fp_equalr (w0, RV_OW0(rv)) || !fp_equalr (wpc, RV_OWPC(rv))) + stat = YES + } + + return (stat) +end |