aboutsummaryrefslogtreecommitdiff
path: root/noao/rv/rvfftcorr.x
diff options
context:
space:
mode:
authorJoe Hunkeler <jhunkeler@gmail.com>2015-08-11 16:51:37 -0400
committerJoe Hunkeler <jhunkeler@gmail.com>2015-08-11 16:51:37 -0400
commit40e5a5811c6ffce9b0974e93cdd927cbcf60c157 (patch)
tree4464880c571602d54f6ae114729bf62a89518057 /noao/rv/rvfftcorr.x
downloadiraf-osx-40e5a5811c6ffce9b0974e93cdd927cbcf60c157.tar.gz
Repatch (from linux) of OSX IRAF
Diffstat (limited to 'noao/rv/rvfftcorr.x')
-rw-r--r--noao/rv/rvfftcorr.x120
1 files changed, 120 insertions, 0 deletions
diff --git a/noao/rv/rvfftcorr.x b/noao/rv/rvfftcorr.x
new file mode 100644
index 00000000..a53a737d
--- /dev/null
+++ b/noao/rv/rvfftcorr.x
@@ -0,0 +1,120 @@
+include <math.h>
+include "rvpackage.h"
+include "rvflags.h"
+include "rvcont.h"
+
+# RV_FFTCORR - Parent routine for the FFT correlation. This procedure prepares
+# the data and processes the FFT.
+
+procedure rv_fftcorr (rv, plot_only)
+
+pointer rv #I RV struct pointer
+int plot_only #I Plot prepared spectra only?
+
+pointer sp
+pointer otempy, rtempy
+pointer wkobj, wkref, ans
+int npts, fnpts
+int i, ishift
+int fft_pow2()
+errchk realloc
+
+begin
+ RV_UPDATE(rv) = YES # set the update flag
+
+ npts = int ((RV_GLOB_W2(rv) - RV_GLOB_W1(rv)) / RV_OWPC(rv) + 1)
+ npts = max (npts, max (RV_NPTS(rv),RV_RNPTS(rv)))
+ fnpts = fft_pow2 (npts)
+ RV_CCFNPTS(rv) = fnpts
+ RV_FFTNPTS(rv) = fnpts
+
+ # (Re)allocate the answer vectors.
+ call realloc (RV_WKPIXX(rv), 2*fnpts, TY_REAL)
+ call realloc (RV_WKPIXY(rv), 2*fnpts, TY_REAL)
+
+ call smark (sp) # allocate some pointers
+ call salloc (wkobj, 2*fnpts, TY_REAL)
+ call salloc (wkref, 2*fnpts, TY_REAL)
+ call salloc (otempy, 2*fnpts, TY_REAL)
+ call salloc (rtempy, 2*fnpts, TY_REAL)
+ call salloc (ans, 2*fnpts, TY_REAL)
+
+ # Clear the work arrays.
+ call aclrr (Memr[wkobj], 2*fnpts)
+ call aclrr (Memr[wkref], 2*fnpts)
+ call aclrr (Memr[otempy], 2*fnpts)
+ call aclrr (Memr[rtempy], 2*fnpts)
+
+ # Prepare the data in the temp arrays before processing. First we'll
+ # do the object spectrum
+ if (OBJCONT(rv) == NO)
+ call amovr (OBJPIXY(rv,1), Memr[otempy], RV_NPTS(rv))
+ else
+ call amovr (OCONT_DATA(rv,1), Memr[otempy], RV_NPTS(rv))
+
+ if (RV_OW0(rv) == RV_GLOB_W1(rv) || RV_DCFLAG(rv) == -1)
+ ishift = 0
+ else
+ ishift = nint ((RV_OW0(rv) - RV_GLOB_W1(rv)) / RV_OWPC(rv))
+
+ call prep_spec (rv, RV_OSAMPLE(rv), npts, fnpts, RV_NPTS(rv),
+ otempy, wkobj, ishift, YES)
+
+
+ # Now do the template spectrum.
+ if (REFCONT(rv) == NO)
+ call amovr (REFPIXY(rv,1), Memr[rtempy], RV_RNPTS(rv))
+ else
+ call amovr (RCONT_DATA(rv,1), Memr[rtempy], RV_RNPTS(rv))
+
+ if (RV_RW0(rv) == RV_GLOB_W1(rv) || RV_DCFLAG(rv) == -1)
+ ishift = 0
+ else
+ ishift = nint ((RV_RW0(rv) - RV_GLOB_W1(rv)) / RV_RWPC(rv))
+
+ call prep_spec (rv, RV_RSAMPLE(rv), npts, fnpts, RV_RNPTS(rv),
+ rtempy, wkref, ishift, YES)
+
+
+ # Normalize the correlation.
+ call rv_normalize (Memr[wkobj], fnpts)
+ call rv_normalize (Memr[wkref], fnpts)
+
+ # Now do a plot of the prepared spectra if that what's requested.
+ if (plot_only == YES) {
+ call gclear (RV_GP(rv)) # clear the screen
+
+ if (RV_FILTER(rv) == OBJ_ONLY || RV_FILTER(rv) == BOTH) {
+ call realft (Memr[wkobj], fnpts, 1) # forward FFT
+ call rv_filter (rv, Memr[wkobj], fnpts) # Object
+ call realft (Memr[wkobj], fnpts, -1) # inverse FFT
+ }
+ if (RV_FILTER(rv) == TEMP_ONLY || RV_FILTER(rv) == BOTH) {
+ call realft (Memr[wkref], fnpts, 1) # forward FFT
+ call rv_filter (rv, Memr[wkref], fnpts) # Template
+ call realft (Memr[wkref], fnpts, -1) # inverse FFT
+ }
+
+ call split_plot (rv, RV_GP(rv), TOP, Memr[wkobj], fnpts,
+ OBJ_ONLY, PREPARED_PLOT)
+ call split_plot (rv, RV_GP(rv), BOTTOM, Memr[wkref], fnpts,
+ TEMP_ONLY, PREPARED_PLOT)
+
+ call sfree (sp)
+ return
+ }
+
+ # Call correlation routine to get answer vector.
+ call rv_correl (rv, Memr[wkobj], Memr[wkref], fnpts, Memr[ans])
+
+ # Load work arrays and fix wrap-around ordering of ans vector.
+ do i = 1, fnpts {
+ WRKPIXY(rv,i) = Memr[ans+i-1]
+ WRKPIXX(rv,i) = real ((i-fnpts/2-1))
+ }
+ call fft_fixwrap (WRKPIXY(rv,1), fnpts)
+
+ RV_Y1(rv) = INDEF # Reset some plot flags
+ RV_Y2(rv) = INDEF
+ call sfree (sp)
+end