aboutsummaryrefslogtreecommitdiff
path: root/noao/rv/deblend.x
diff options
context:
space:
mode:
Diffstat (limited to 'noao/rv/deblend.x')
-rw-r--r--noao/rv/deblend.x776
1 files changed, 776 insertions, 0 deletions
diff --git a/noao/rv/deblend.x b/noao/rv/deblend.x
new file mode 100644
index 00000000..ee8c0e8f
--- /dev/null
+++ b/noao/rv/deblend.x
@@ -0,0 +1,776 @@
+include <error.h>
+include <mach.h>
+include <gset.h>
+include "rvpackage.h"
+include "rvflags.h"
+
+# DEBLEND -- Deblend up to 4 lines in a spectral region.
+
+procedure deblend (rv, gp, x1, x2, dx, wx1, wy1, pix, ans, nans)
+
+pointer rv #I RV struct pointer
+pointer gp #I GIO file descriptor
+real x1, x2, dx #I Coordinate scale
+real wx1, wy1 #I Cursor position
+real pix[ARB] #I Spectrum data
+char ans[2*SZ_LINE,4] #O Answer strings
+int nans #O Number of answer strings
+
+int i, j, i1, npts, nlines, maxlines, wc, key, op, stat
+
+double vobs, vhelio, verr
+real w, wxc, wyc, wx, wy, wx2, wy2, a[14], waves[4]
+real slope, height, flux, cont, sigma, eqw, scale, chisq
+real serr, shift, fwhm
+bool fit
+pointer sp, cmd, x, y, anti
+
+int scan(), clgcur(), clgkey(), rv_rvcorrect()
+errchk dofit
+
+include "fitcom.com"
+define done_ 99
+define HELP "noao$lib/scr/deblend.key"
+define OP "Option (a=0p1s, b=1p1s, c=np1s, d=0pns, e=1pns, f=npns, q=quit):"
+define SQ2PI 2.5066283
+
+begin
+ call smark (sp)
+ call salloc (cmd, SZ_FNAME, TY_CHAR)
+
+ # Input cursor is first continuum point now get second continuum point.
+ call printf ("d again:")
+ if (clgcur ("cursor", wx2, wy2, wc, key, Memc[cmd], SZ_FNAME) == EOF) {
+ call sfree (sp)
+ return
+ }
+ call gctran (gp, wx2, wy2, wx2, wy2, wc, 2)
+ if (RV_FITDONE(rv) == YES) {
+ call rv_erase_fit (rv, false)
+ RV_FITDONE(rv) = NO
+ IS_DBLSTAR(rv) = NO
+ }
+
+ # Set pixel indices and determine number of points to fit.
+ call fixx (wx1, wx2, wy1, wy2, x1, x2)
+ call pixind (x1, dx, wx1, i1)
+ call pixind (x1, dx, wx2, j)
+ npts = j - i1 + 1
+ RV_IEND(rv) = j
+ RV_ISTART(rv) = i1
+ if (npts < 3) {
+ call rv_errmsg ("At least 3 points are required\n")
+ call sfree (sp)
+ return
+ }
+
+ # Allocate space for the points to be fit.
+ call salloc (x, npts, TY_REAL)
+ call salloc (y, npts, TY_REAL)
+
+ # Subtract the continuum and scale the data.
+ wxc = wx1
+ wyc = wy1
+ slope = (wy2-wy1) / (wx2-wx1)
+ scale = 0.
+ do i = 1, npts {
+ w = x1 + (i1+i-2) * dx
+ Memr[y+i-1] = pix[i1+i-1] - (wyc + slope * (w-wxc))
+ scale = max (scale, abs (Memr[y+i-1]))
+ Memr[x+i-1] = w
+ }
+ call adivkr (Memr[y], scale, Memr[y], npts)
+
+ # Select the lines to be fit. If no lines return.
+ maxlines = 4
+ nlines = 0
+ call printf ("Lines ('m' to mark, 't' to type, 'q' to quit):")
+ while (clgcur ("cursor", wx, wy, wc, key, Memc[cmd], SZ_FNAME) != EOF) {
+ switch (key) {
+ case 'm':
+ call gctran (gp, wx, wy, wx, wy, wc, 2)
+ case 't':
+ if (RV_DCFLAG(rv) == -1) {
+ call printf ("shift: ")
+ call flush (STDOUT)
+ if (scan() != EOF)
+ call gargr (wx)
+ } else {
+ call printf ("velocity: ")
+ call flush (STDOUT)
+ if (scan() != EOF)
+ call gargr (wx)
+ wx = wx / RV_DELTAV(rv)
+ }
+ call printf ("Lines ('m' to mark, 't' to type, 'q' to quit):")
+ case 'q':
+ call printf ("\n")
+ break
+ case 'I':
+ call fatal (0, "Interrupt")
+ default:
+ call printf (
+ "Lines ('m' to mark, 't' to type, 'q' to quit):")
+ next
+ }
+ for (i = 1; i <= nlines && wx != waves[i]; i = i + 1)
+ ;
+ if (i > nlines) {
+ nlines = nlines + 1
+ waves[nlines] = wx
+ call gmark (gp, wx, wy, GM_VLINE, 4., 4.)
+ call gflush (gp)
+ }
+ if (nlines == maxlines) {
+ call printf ("\n")
+ break
+ }
+ }
+ if (nlines == 0)
+ goto done_
+
+ # Do fits.
+ fit = false
+ call printf (OP)
+ while (clgcur ("cursor", wx, wy, wc, op, Memc[cmd], SZ_FNAME) != EOF) {
+ switch (op) {
+ case '?':
+ call gpagefile (gp, HELP, "Rvxcor Deblending Options")
+ call printf (OP)
+ next
+ case 'a', 'b', 'c', 'd', 'e', 'f':
+ case 'q':
+ call printf ("\n")
+ break
+ case 'I':
+ call fatal (0, "Interrupt")
+ default:
+ call printf ("%s")
+ call pargstr (OP)
+ next
+ }
+
+ # Erase old deblended fit in case we've been here before. Fit is
+ # erased above from when we first entered.
+ if (IS_DBLSTAR(rv) == YES) {
+ call gseti (gp, G_PLTYPE, GL_CLEAR)
+ call rv_plt_deblend (rv, gp, NO)
+ call gseti (gp, G_PLTYPE, GL_SOLID)
+ }
+
+ # Save some variables for later plotting.
+ DBL_X1(rv) = wx1
+ DBL_X2(rv) = wx2
+ DBL_Y1(rv) = wy1
+ DBL_Y2(rv) = wy2
+ DBL_I1(rv) = i1
+ DBL_NFITP(rv) = npts
+ DBL_SCALE(rv) = scale
+ DBL_SLOPE(rv) = slope
+
+ # Convert line postions to relative to first line.
+ a[1] = waves[1]
+ a[2] = 0.25 * abs (Memr[x+npts-1] - Memr[x]) / nlines
+ do i = 1, nlines {
+ call pixind (x1, dx, waves[i], j)
+ a[3*i] = (pix[j] - (wyc + slope * (waves[i]-wxc))) / scale
+ a[3*i+1] = waves[i] - waves[1]
+ a[3*i+2] = 1.
+ }
+
+ switch (op) {
+ case 'a':
+ iferr {
+ call dofit ('a', Memr[x], Memr[y], npts, a, nlines, chisq)
+ call dofit ('a', Memr[x], Memr[y], npts, a, nlines, chisq)
+ } then {
+ call erract (EA_WARN)
+ next
+ }
+ case 'b':
+ iferr {
+ call dofit ('a', Memr[x], Memr[y], npts, a, nlines, chisq)
+ call dofit ('b', Memr[x], Memr[y], npts, a, nlines, chisq)
+ } then {
+ call erract (EA_WARN)
+ next
+ }
+ case 'c':
+ iferr {
+ call dofit ('a', Memr[x], Memr[y], npts, a, nlines, chisq)
+ call dofit ('b', Memr[x], Memr[y], npts, a, nlines, chisq)
+ call dofit ('c', Memr[x], Memr[y], npts, a, nlines, chisq)
+ } then {
+ call erract (EA_WARN)
+ next
+ }
+ case 'd':
+ iferr {
+ call dofit ('a', Memr[x], Memr[y], npts, a, nlines, chisq)
+ call dofit ('d', Memr[x], Memr[y], npts, a, nlines, chisq)
+ } then {
+ call erract (EA_WARN)
+ next
+ }
+ case 'e':
+ iferr {
+ call dofit ('a', Memr[x], Memr[y], npts, a, nlines, chisq)
+ call dofit ('b', Memr[x], Memr[y], npts, a, nlines, chisq)
+ call dofit ('e', Memr[x], Memr[y], npts, a, nlines, chisq)
+ } then {
+ call erract (EA_WARN)
+ next
+ }
+ case 'f':
+ iferr {
+ call dofit ('a', Memr[x], Memr[y], npts, a, nlines, chisq)
+ call dofit ('b', Memr[x], Memr[y], npts, a, nlines, chisq)
+ call dofit ('c', Memr[x], Memr[y], npts, a, nlines, chisq)
+ call dofit ('f', Memr[x], Memr[y], npts, a, nlines, chisq)
+ } then {
+ call erract (EA_WARN)
+ next
+ }
+ }
+ fit = true
+ RV_FITDONE(rv) = YES
+ DBL_NSHIFTS(rv) = nlines
+ call amovr (a, DBL_COEFFS(rv,1), 3*nlines+2)
+
+ # Update parameters in the fitting common for the output log
+ nfit = npts
+ nfitpars = 3*nlines+2
+ binshift = INDEFI
+ niter = 3
+ chisqr = INDEF
+ ccfvar = INDEF
+ mresid = INDEF
+ sresid = INDEF
+
+ # Compute model spectrum with continuum and plot.
+ IS_DBLSTAR(rv) = YES
+ call rv_plt_deblend (rv, gp, NO)
+
+ # Print computed values on status line.
+ i = 1
+ key = ''
+ repeat {
+ call flush (STDOUT)
+ switch (key) {
+ case '-':
+ i = i - 1
+ if (i < 1)
+ i = nlines
+ case '+':
+ i = i + 1
+ if (i > nlines)
+ i = 1
+ case 'q':
+ call printf ("\n")
+ break
+ }
+
+ height = scale * a[3*i]
+ w = a[1] + a[3*i+1]
+ sigma = abs (a[2]*a[3*i+2])
+ flux = sigma * height * SQ2PI
+ cont = wyc + slope * (w - wxc)
+ if (cont > 0.)
+ eqw = abs (flux) / cont
+ else
+ eqw = INDEF
+
+ if (key == 'r') {
+ call printf ("\nrms = %8.4g")
+ call pargr (scale * sqrt (chisq / npts))
+ } else if (key == 'I') {
+ call fatal (0, "Interrupt")
+ } else if (key == 'v') {
+ serr = 0.0
+ shift = w
+ stat = rv_rvcorrect (rv, shift, serr, vobs, vhelio, verr)
+ call printf (
+ "\n%d: shift = %8.4f Vo = %8.3f Vh = %8.3f fwhm = %6.4f")
+ call pargi (i)
+ call pargr (shift)
+ call pargd (vobs)
+ call pargd (vhelio)
+ call pargr (2.35482 * sigma * RV_DELTAV(rv))
+ } else {
+ call printf (
+ "\n%d: center = %8.6g, flux = %8.4g, eqw = %6.4g, fwhm = %6.4g")
+ call pargi (i)
+ call pargr (w)
+ call pargr (flux)
+ call pargr (eqw)
+ call pargr (2.35482 * sigma)
+ }
+
+ call printf (" (+,-,v,r,q):")
+ call flush (STDOUT)
+ } until (clgkey ("ukey", key, Memc[cmd], SZ_FNAME) == EOF)
+
+ # Log computed values
+ nans = nlines
+ do i = 1, nlines {
+ w = a[1] + a[3*i+1]
+ cont = wyc + slope * (w - wxc)
+ height = scale * a[3*i]
+ sigma = abs (a[2]*a[3*i+2])
+ flux = sigma * height * SQ2PI
+ if (cont > 0.)
+ eqw = abs (flux) / cont
+ else
+ eqw = INDEF
+
+ call sprintf (ans[1,i], 2*SZ_LINE,
+ " %9.7g %9.7g %9.6g %9.4g %9.6g %9.4g %9.4g\n")
+ call pargr (w)
+ call pargr (cont)
+ call pargr (flux)
+ call pargr (eqw)
+ call pargr (height)
+ call pargr (sigma)
+ call pargr (2.35482 * sigma)
+
+ # Now calculate and save the velocity information
+ serr = 0.0
+ if (RV_DCFLAG(rv) != -1) {
+ stat = rv_rvcorrect (rv, w, serr, vobs, vhelio, verr)
+ call salloc (anti, RV_CCFNPTS(rv), TY_REAL)
+ fwhm = 2.35482 * sigma
+ call rv_antisym (rv, w, height, fwhm, WRKPIXY(rv,1),
+ RV_CCFNPTS(rv), Memr[anti], ccfvar, verr, DBL_R(rv,i))
+ if (IS_INDEFD(vobs))
+ DBL_VOBS(rv,i) = INDEFR
+ else
+ DBL_VOBS(rv,i) = real (vobs)
+ if (IS_INDEFD(vhelio))
+ DBL_VHELIO(rv,i) = INDEFR
+ else
+ DBL_VHELIO(rv,i) = real (vhelio)
+ if (IS_INDEFD(verr))
+ DBL_VERR(rv,i) = INDEFR
+ else
+ DBL_VERR(rv,i) = real (verr)
+ DBL_FWHM(rv,i) = 2.35482 * sigma * RV_DELTAV(rv)
+ } else {
+ DBL_VOBS(rv,i) = INDEFR
+ DBL_VHELIO(rv,i) = INDEFR
+ DBL_VERR(rv,i) = INDEFR
+ DBL_FWHM(rv,i) = 2.35482 * sigma
+ }
+ DBL_HEIGHT(rv,i) = height
+ DBL_SHIFT(rv,i) = w
+ }
+ call printf (OP)
+ }
+
+
+done_ call sfree (sp)
+end
+
+
+# SUBBLEND -- Subtract last fit.
+
+procedure subblend (rv, gp, pix, x1, x2, dx, wx1, wy1)
+
+pointer rv #I RV struct pointer
+pointer gp #I Graphics descriptor
+real pix[ARB] #I CCF array
+real x1, x2, dx #I Coordinate scale
+real wx1, wy1 #I Cursor position
+
+int i, j, i1, wc, npts, key
+real w, wx2, wy2
+pointer sp, cmd
+
+int clgcur()
+real model()
+
+begin
+ call smark (sp)
+ call salloc (cmd, SZ_FNAME, TY_CHAR)
+
+ # Subtract continuum subtracted curve from spectrum
+ if (RV_FITDONE(rv) == NO) {
+ call sfree (sp)
+ return
+ }
+
+ # Determine fit range
+ call printf ("- again:")
+ call flush (STDOUT)
+ if (clgcur ("cursor", wx2, wy2, wc, key, Memc[cmd], SZ_FNAME) == EOF) {
+ call sfree (sp)
+ return
+ }
+
+ call fixx (wx1, wx2, wy1, wy2, x1, x2)
+ call pixind (x1, dx, wx1, i1)
+ call pixind (x1, dx, wx2, j)
+ npts = j - i1 + 1
+
+ do i = 1, npts {
+ w = x1 + (i1+i-2) * dx
+ pix[i1+i-1] = pix[i1+i-1] - DBL_SCALE(rv) * model (w,
+ DBL_COEFFS(rv,1), 3*DBL_NSHIFTS(rv)+2)
+ }
+
+ # Plot subtracted curve
+ call gvline (gp, pix[i1], npts, wx1, wx2)
+ call gflush (gp)
+
+ RV_FITDONE(rv) = NO
+ call sfree (sp)
+end
+
+
+# DOFIT -- Perform nonlinear iterative fit for the specified parameters.
+# This uses the Levenberg-Marquardt method from NUMERICAL RECIPES.
+
+procedure dofit (key, x, y, npts, a, nlines, chisq)
+
+int key #I Fitting option
+real x[npts] #I X data
+real y[npts] #I Y data
+int npts #I Number of points
+real a[ARB] #I Fitting parameters
+int nlines #I Number of lines
+real chisq #O Chi squared
+
+int i, np, nfit
+real mr, chi2
+pointer sp, flags
+errchk mr_solve
+
+begin
+ # Number of terms is 3 for each line plus common center and sigma.
+ np = 3 * nlines + 2
+
+ call smark (sp)
+ call salloc (flags, np, TY_INT)
+
+ # Peaks are always fit.
+ switch (key) {
+ case 'a': # Solve one sigma.
+ nfit = 1 + nlines
+ Memi[flags] = 2
+ do i = 1, nlines
+ Memi[flags+i] = 3 * i
+ case 'b': # Solve one position and one sigma.
+ nfit = 2 + nlines
+ Memi[flags] = 1
+ Memi[flags+1] = 2
+ do i = 1, nlines
+ Memi[flags+1+i] = 3 * i
+ case 'c': # Solve independent positions and one sigma.
+ nfit = 1 + 2 * nlines
+ Memi[flags] = 2
+ do i = 1, nlines {
+ Memi[flags+2*i-1] = 3 * i
+ Memi[flags+2*i] = 3 * i + 1
+ }
+ case 'd': # Solve for sigmas.
+ nfit = 2 * nlines
+ do i = 1, nlines {
+ Memi[flags+2*i-2] = 3 * i
+ Memi[flags+2*i-1] = 3 * i + 2
+ }
+ case 'e': # Solve for one position and sigmas.
+ nfit = 1 + 2 * nlines
+ Memi[flags] = 1
+ do i = 1, nlines {
+ Memi[flags+2*i-1] = 3 * i
+ Memi[flags+2*i] = 3 * i + 2
+ }
+ case 'f': # Solve for positions and sigmas.
+ nfit = 3 * nlines
+ do i = 1, nfit
+ Memi[flags+i-1] = i + 2
+ }
+
+
+ mr = -1.
+ i = 0
+ chi2 = MAX_REAL
+ repeat {
+ call mr_solve (x, y, npts, a, Memi[flags], np, nfit, mr, chisq)
+ if (chi2 - chisq > 1.)
+ i = 0
+ else
+ i = i + 1
+ chi2 = chisq
+ } until (i == 3)
+
+ mr = 0.
+ call mr_solve (x, y, npts, a, Memi[flags], np, nfit, mr, chisq)
+
+ call sfree (sp)
+end
+
+
+# MODEL -- Compute model from fitted parameters.
+#
+# I(x) = I(i) exp {[(x - xc - dx(i)) / (sig sig(i))] ** 2 / 2.}
+#
+# where the parameters are xc, sig, I(i), dx(i), and sig(i) (i=1,nlines).
+
+real procedure model (x, a, na)
+
+real x #I X value to be evaluated
+real a[na] #I Parameters
+int na #I Number of parameters
+
+int i
+real y, arg
+
+begin
+ y = 0.
+ do i = 3, na, 3 {
+ arg = (x - a[1] - a[i+1]) / (a[2] * a[i+2])
+ if (abs (arg) < 7.)
+ y = y + a[i] * exp (-arg**2 / 2.)
+ }
+ return (y)
+end
+
+
+# DERIVS -- Compute model and derivatives for MR_SOLVE procedure.
+#
+# I(x) = I(i) exp {[(x - xc - dx(i)) / (sig sig(i))] ** 2 / 2.}
+#
+# where the parameters are xc, sig, I(i), dx(i), and sig(i) (i=1,nlines).
+
+procedure derivs (x, a, y, dyda, na)
+
+real x #I X value to be evaluated
+real a[na] #I Parameters
+real y #O Function value
+real dyda[na] #O Derivatives
+int na #I Number of parameters
+
+int i
+real sig, arg, ex, fac
+
+begin
+ y = 0.
+ dyda[1] = 0.
+ dyda[2] = 0.
+ do i = 3, na, 3 {
+ sig = a[2] * a[i+2]
+ arg = (x - a[1] - a[i+1]) / sig
+ if (abs (arg) < 7.)
+ ex = exp (-arg**2 / 2.)
+ else
+ ex = 0.
+ fac = a[i] * ex * arg
+
+ y = y + a[i] * ex
+ dyda[1] = dyda[1] + fac / sig
+ dyda[2] = dyda[2] + fac * arg / a[2]
+ dyda[i] = ex
+ dyda[i+1] = fac / sig
+ dyda[i+2] = fac * arg / a[i+2]
+ }
+end
+
+
+# MR_SOLVE -- Levenberg-Marquardt nonlinear chi square minimization.
+#
+# Use the Levenberg-Marquardt method to minimize the chi squared of a set
+# of paraemters. The parameters being fit are indexed by the flag array.
+# To initialize the Marquardt parameter, MR, is less than zero. After that
+# the parameter is adjusted as needed. To finish set the parameter to zero
+# to free memory. This procedure requires a subroutine, DERIVS, which
+# takes the derivatives of the function being fit with respect to the
+# parameters. There is no limitation on the number of parameters or
+# data points. For a description of the method see NUMERICAL RECIPES
+# by Press, Flannery, Teukolsky, and Vetterling, p523.
+
+procedure mr_solve (x, y, npts, params, flags, np, nfit, mr, chisq)
+
+real x[npts] #I X data array
+real y[npts] #I Y data array
+int npts #I Number of data points
+real params[np] #U Parameter array
+int flags[np] #I Flag array indexing parameters to fit
+int np #I Number of parameters
+int nfit #I Number of parameters to fit
+real mr #O MR parameter
+real chisq #O Chi square of fit
+
+int i
+real chisq1
+pointer new, a1, a2, delta1, delta2
+
+errchk mr_invert
+
+begin
+ # Allocate memory and initialize.
+ if (mr < 0.) {
+ call mfree (new, TY_REAL)
+ call mfree (a1, TY_REAL)
+ call mfree (a2, TY_REAL)
+ call mfree (delta1, TY_REAL)
+ call mfree (delta2, TY_REAL)
+
+ call malloc (new, np, TY_REAL)
+ call malloc (a1, nfit*nfit, TY_REAL)
+ call malloc (a2, nfit*nfit, TY_REAL)
+ call malloc (delta1, nfit, TY_REAL)
+ call malloc (delta2, nfit, TY_REAL)
+
+ call amovr (params, Memr[new], np)
+ call mr_eval (x, y, npts, Memr[new], flags, np, Memr[a2],
+ Memr[delta2], nfit, chisq)
+ mr = 0.001
+ }
+
+ # Restore last good fit and apply the Marquardt parameter.
+ call amovr (Memr[a2], Memr[a1], nfit * nfit)
+ call amovr (Memr[delta2], Memr[delta1], nfit)
+ do i = 1, nfit
+ Memr[a1+(i-1)*(nfit+1)] = Memr[a2+(i-1)*(nfit+1)] * (1. + mr)
+
+ # Matrix solution.
+ call mr_invert (Memr[a1], Memr[delta1], nfit)
+
+ # Compute the new values and curvature matrix.
+ do i = 1, nfit
+ Memr[new+flags[i]-1] = params[flags[i]] + Memr[delta1+i-1]
+ call mr_eval (x, y, npts, Memr[new], flags, np, Memr[a1],
+ Memr[delta1], nfit, chisq1)
+
+ # Check if chisq has improved.
+ if (chisq1 < chisq) {
+ mr = 0.1 * mr
+ chisq = chisq1
+ call amovr (Memr[a1], Memr[a2], nfit * nfit)
+ call amovr (Memr[delta1], Memr[delta2], nfit)
+ call amovr (Memr[new], params, np)
+ } else
+ mr = 10. * mr
+
+ if (mr == 0.) {
+ call mfree (new, TY_REAL)
+ call mfree (a1, TY_REAL)
+ call mfree (a2, TY_REAL)
+ call mfree (delta1, TY_REAL)
+ call mfree (delta2, TY_REAL)
+ }
+end
+
+
+# MR_EVAL -- Evaluate curvature matrix. This calls procedure DERIVS.
+
+procedure mr_eval (x, y, npts, params, flags, np, a, delta, nfit, chisq)
+
+real x[npts] #I X data array
+real y[npts] #I Y data array
+int npts #I Number of data points
+real params[np] #I Parameter array
+int flags[np] #I Flag array indexing parameters to fit
+int np #I Number of parameters
+real a[nfit,nfit] #U Curvature matrix
+real delta[nfit] #U Delta array
+int nfit #I Number of parameters to fit
+real chisq #U Chi square of fit
+
+int i, j, k
+real ymod, dy, dydpj, dydpk
+pointer sp, dydp
+
+begin
+ call smark (sp)
+ call salloc (dydp, np, TY_REAL)
+
+ do j = 1, nfit {
+ do k = 1, j
+ a[j,k] = 0.
+ delta[j] = 0.
+ }
+
+ chisq = 0.
+ do i = 1, npts {
+ call derivs (x[i], params, ymod, Memr[dydp], np)
+ dy = y[i] - ymod
+ do j = 1, nfit {
+ dydpj = Memr[dydp+flags[j]-1]
+ delta[j] = delta[j] + dy * dydpj
+ do k = 1, j {
+ dydpk = Memr[dydp+flags[k]-1]
+ a[j,k] = a[j,k] + dydpj * dydpk
+ }
+ }
+ chisq = chisq + dy * dy
+ }
+
+ do j = 2, nfit
+ do k = 1, j-1
+ a[k,j] = a[j,k]
+
+ call sfree (sp)
+end
+
+
+# MR_INVERT -- Solve a set of linear equations using Householder transforms.
+
+procedure mr_invert (a, b, n)
+
+real a[n,n] #I Input matrix and returned inverse
+real b[n] #U Input RHS vector and returned solution
+int n #I Dimension of input matrices
+
+int krank
+real rnorm
+pointer sp, h, g, ip
+
+begin
+ call smark (sp)
+ call salloc (h, n, TY_REAL)
+ call salloc (g, n, TY_REAL)
+ call salloc (ip, n, TY_INT)
+
+ call hfti (a, n, n, n, b, n, 1, 0.001, krank, rnorm,
+ Memr[h], Memr[g], Memi[ip])
+
+ call sfree (sp)
+end
+
+
+# FIXX - Check for bounds on x's.
+
+procedure fixx (eqx1, eqx2, eqy1, eqy2, x1, x2)
+
+real eqx1, eqx2, eqy1, eqy2, x1, x2
+
+real temp
+
+begin
+ if ((x1 - x2) * (eqx1 - eqx2) < 0.) {
+ temp = eqx2
+ eqx2 = eqx1
+ eqx1 = temp
+
+ temp = eqy2
+ eqy2 = eqy1
+ eqy1 = temp
+ }
+
+ eqx1 = max (min (x1, x2), min (max (x1, x2), eqx1))
+ eqx2 = max (min (x1, x2), min (max (x1, x2), eqx2))
+end
+
+
+# PIXIND -- Compute pixel index.
+
+procedure pixind (x1, dx, valx, i1)
+
+real x1, dx, valx
+int i1
+
+begin
+# i1 = aint ((valx-x1)/dx +0.5) + 1
+ i1 = (valx - x1) / dx + 1.5
+end
generated by cgit (git 2.34.1) at 2025-09-20 07:24:39 -0400