From 40e5a5811c6ffce9b0974e93cdd927cbcf60c157 Mon Sep 17 00:00:00 2001 From: Joe Hunkeler Date: Tue, 11 Aug 2015 16:51:37 -0400 Subject: Repatch (from linux) of OSX IRAF --- pkg/utilities/curfit.gx | 216 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 216 insertions(+) create mode 100644 pkg/utilities/curfit.gx (limited to 'pkg/utilities/curfit.gx') diff --git a/pkg/utilities/curfit.gx b/pkg/utilities/curfit.gx new file mode 100644 index 00000000..588bc1d1 --- /dev/null +++ b/pkg/utilities/curfit.gx @@ -0,0 +1,216 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include +include +include +include +include "curfit.h" + +define VERBOSE_OUTPUT 1 +define LIST_OUTPUT 2 +define DEFAULT_OUTPUT 3 + +define CF_UNIFORM 1 +define CF_USER 2 +define CF_STATISTICAL 3 +define CF_INSTRUMENTAL 4 + +# CF_FIT -- Called once for each curve to be fit. + +$for (rd) +procedure cf_fit$t (ic, gt, x, y, wts, nvalues, nmax, device, interactive, ofmt, + power) + +pointer ic # ICFIT pointer +pointer gt # Graphics tools pointer +PIXEL x[nmax] # X data values +PIXEL y[nmax] # Y data values +PIXEL wts[nmax] # Weights +int nvalues # Number of data points +int nmax # Maximum number of data points +char device[SZ_FNAME] # Output graphics device +int interactive # Fit curve interactively? +int ofmt # Type of output listing +bool power # Convert coeff to power series? + +int ncoeff, i +PIXEL xmin, xmax +pointer sp, gp, cv, coeff, tty +pointer gopen(), ttyodes() +int fstati(), $tcvstati() + +begin + # Determine data range and set up curve fitting limits. + call alim$t (x, nvalues, xmin, xmax) + call ic_putr (ic, "xmin", real (xmin)) + call ic_putr (ic, "xmax", real (xmax)) + + if (interactive == YES) { + gp = gopen (device, NEW_FILE, STDGRAPH) + call icg_fit$t (ic, gp, "cursor", gt, cv, x, y, wts, nvalues) + call gclose (gp) + } else + # Do fit non-interactively + call ic_fit$t (ic, cv, x, y, wts, nvalues, YES, YES, YES, YES) + + # Output answers to STDOUT + if (ofmt != LIST_OUTPUT) { + if (fstati (STDOUT, F_REDIR) == NO) { + tty = ttyodes ("terminal") + call ttyclear (STDOUT, tty) + call ttycdes (tty) + } + + #call ic_show (ic, "STDOUT", gt) + call ic_vshow$t (ic, "STDOUT", cv, x, y, wts, nvalues, gt) + + if (ofmt == VERBOSE_OUTPUT) { + call printf ( + "\n# \t X \t Yc \t Y \t W\n") + call cf_listxy$t (cv, x, y, wts, nvalues) + } + } else + call cf_listxy$t (cv, x, y, wts, nvalues) + + # Convert coefficients if requested for legendre or chebyshev + if (power && ofmt != LIST_OUTPUT) { + # Calculate and print coefficients + ncoeff = $tcvstati (cv, CVNCOEFF) + call smark (sp) + call salloc (coeff, ncoeff, TY_PIXEL) + call $tcvpower (cv, Mem$t[coeff], ncoeff) + call printf ("# Power series coefficients would be:\n") + call printf ("# \t\tcoefficient\n") + do i = 1, ncoeff { + call printf ("# \t%d \t%14.7e\n") + call pargi (i) + call parg$t (Mem$t[coeff+i-1]) + } + call sfree (sp) + } + +$if (datatype == r) + call cvfree (cv) +$else + call $tcvfree (cv) +$endif + #call ic_close$t (ic) +end + + +# CF_LISTXY -- Print answers to STDOUT as x,y pairs. + +procedure cf_listxy$t (cv, xvals, yvals, wts, nvalues) + +pointer cv # Pointer to curfit structure +int nvalues # Number of data values +PIXEL xvals[nvalues] # Array of x data values +PIXEL yvals[nvalues] # Array of y data values +PIXEL wts[nvalues] # Array of weights + +int i +PIXEL $tcveval() + +begin + do i = 1, nvalues { + call printf ("\t%14.7e \t%14.7e \t%14.7e \t%14.7e\n") + call parg$t (xvals[i]) + call parg$t ($tcveval (cv, xvals[i])) + call parg$t (yvals[i]) + call parg$t (wts[i]) + } +end + +# IM_PROJECTION -- Given an image section of arbitrary dimension, compute +# the projection along a single axis by taking the average over the other +# axes. We do not know about bad pixels. + +procedure im_projection$t (im, x, y, w, npix, weighting, axis) + +pointer im # Pointer to image header structure +PIXEL x[npix] # Index of projection vector +PIXEL y[npix] # Receives the projection vector +PIXEL w[npix] # Receives the weight vector +int weighting # Weighting of the individual points +int npix # Length of projection vector +int axis # The axis to be projected to (x=1) + +int i, lastv +long v[IM_MAXDIM], nsum, totpix +pointer pix +PIXEL asum$t() +pointer imgnl$t() +errchk imgnl$t + +begin + if (im == NULL) + call error (1, "Image projection operator called with null im") + if (axis < 1 || axis > IM_NDIM(im)) + call error (2, "Attempt to take projection over nonexistent axis") + + + # Set the y projection vector + call aclr$t (y, npix) + call amovkl (long(1), v, IM_MAXDIM) + + switch (axis) { + case 1: + # Since the image is read line by line, it is easy to compute the + # projection along the x-axis (axis 1). We merely sum all of the + # image lines. + + while (imgnl$t (im, pix, v) != EOF) + call aadd$t (Mem$t[pix], y, y, npix) + + default: + # Projecting along any other axis when reading the image line + # by line is a bit difficult to understand. Basically, the + # element 'axis' of the V vector (position of the line in the + # image) gives us the index into the appropriate element of + # y. When computing the projection over multiple dimensions, + # the same output element will be referenced repeatedly. All + # of the elmenents of the input line are summed and added into + # this output element. + + for (lastv=v[axis]; imgnl$t (im, pix, v) != EOF; lastv=v[axis]) { + i = lastv + if (i <= npix) + y[i] = y[i] + asum$t (Mem$t[pix], IM_LEN(im,1)) + } + } + + # Now compute the number of pixels contributing to each element + # of the output vector. This is the number of pixels in the image + # divided by the length of the projection. + + totpix = 1 + do i = 1, IM_NDIM(im) + if (i == axis) + totpix = totpix * min (npix, IM_LEN(im,i)) + else + totpix = totpix * IM_LEN(im,i) + nsum = totpix / min (npix, IM_LEN(im,axis)) + + # Compute the average by dividing by the number if pixels summed at + # each point. + call adivk$t (y, PIXEL (nsum), y, npix) + + # Set the x and weight vectors + do i = 1, npix { + x[i] = i + switch (weighting) { + case CF_STATISTICAL: + if (y[i] > 0.0) + w[i] = 1.0 / y[i] + else if (y[i] < 0.0) + w[i] = abs (1.0 / y[i]) + else + w[i] = 1.0 + case CF_UNIFORM: + w[i] = 1. + default: + w[i] = 1. + } + } +end +$endfor -- cgit