From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- pkg/utilities/t_curfit.x | 446 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 446 insertions(+) create mode 100644 pkg/utilities/t_curfit.x (limited to 'pkg/utilities/t_curfit.x') diff --git a/pkg/utilities/t_curfit.x b/pkg/utilities/t_curfit.x new file mode 100644 index 00000000..af4dea06 --- /dev/null +++ b/pkg/utilities/t_curfit.x @@ -0,0 +1,446 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include +include +include +include +include +include + +define VERBOSE_OUTPUT 1 +define LIST_OUTPUT 2 +define DEFAULT_OUTPUT 3 +define IMAGE_OP 1 +define LIST_OP 2 + +define CF_UNIFORM 1 +define CF_USER 2 +define CF_STATISTICAL 3 +define CF_INSTRUMENTAL 4 + +define NADD 20 # Number of points that can be added by ICFIT + +# T_CURFIT -- cl interface to the curfit package. Task CURFIT provides +# four fitting options: legendre, chebyshev, cubic spline or linear spline. +# The output can be printed in default, verbose or tabular formats. The +# user can also choose to interactively fit the curve. + +procedure t_curfit () + +pointer x, y, w, gt, fcn, fname, flist, dev, str, sp, ic +bool listdata, verbose, power, redir +int fd, ofmt, interactive, datatype +int axis, nvalues, nmax, weighting +pointer gt_init() +bool clgetb() +int imtopen(), clgeti(), cf_operand(), cf_rimage(), cf_rlist() +int imtgetim(), clgwrd() +int fstati() + +begin + # Allocate space for string buffers + call smark (sp) + call salloc (fcn, SZ_FNAME, TY_CHAR) + call salloc (fname, SZ_FNAME, TY_CHAR) + call salloc (flist, SZ_LINE, TY_CHAR) + call salloc (dev, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_FNAME, TY_CHAR) + + # First get cl parameters. Check to see if input has been redirected. + redir = false + if (fstati (STDIN, F_REDIR) == YES) { + redir = true + call strcpy ("STDIN", Memc[fname], SZ_FNAME) + } else { + call clgstr ("input", Memc[flist], SZ_LINE) + fd = imtopen (Memc[flist]) + } + + listdata = clgetb ("listdata") + verbose = clgetb ("verbose") + ofmt = DEFAULT_OUTPUT + if (listdata) + ofmt = LIST_OUTPUT + else if (verbose) + ofmt = VERBOSE_OUTPUT + + # Determine the calculation datatype. + switch (clgwrd ("calctype", Memc[dev], SZ_FNAME, "|real|double|")) { + case 1: + datatype = TY_REAL + case 2: + datatype = TY_DOUBLE + } + + if (clgetb ("interactive")) { + interactive = YES + call clgstr ("device", Memc[dev], SZ_FNAME) + } else { + interactive = ALWAYSNO + call strcpy ("", Memc[dev], SZ_FNAME) + } + + power = clgetb ("power") + + call ic_open (ic) + call clgstr ("function", Memc[fcn], SZ_FNAME) + call ic_pstr (ic, "function", Memc[fcn]) + call ic_puti (ic, "order", clgeti ("order")) + weighting = clgwrd ("weighting", Memc[str], SZ_FNAME, + "|uniform|user|statistical|instrumental|") + + gt = gt_init () + + repeat { + if (!redir) { + if (imtgetim (fd, Memc[fname], SZ_FNAME) == EOF) + break + } + + if (cf_operand (Memc[fname]) == IMAGE_OP) { + axis = clgeti ("axis") + nvalues = cf_rimage (Memc[fname], axis, x, y, w, weighting, + datatype) + call gt_sets (gt, GTTYPE, "line") + } else { + nvalues = cf_rlist (Memc[fname], x, y, w, weighting, datatype) + call gt_sets (gt, GTTYPE, "mark") + + # For list input only, order the input array. The + # rg_ranges package requires an x ordered input array, or else + # points will be excluded from the fit. This test can be + # removed when/if the ordering restriction is removed from + # rg_xranges. Sorted data is required even when no sampling + # is done, as in the default case of sample=*. (ShJ 6-24-88) + + switch (datatype) { + case TY_REAL: + call xt_sort3 (Memr[x], Memr[y], Memr[w], nvalues) + case TY_DOUBLE: + call xt_sort3d (Memd[x], Memd[y], Memd[w], nvalues) + } + } + + # Allow for adding points. + nmax = nvalues + NADD + call realloc (x, nmax, datatype) + call realloc (y, nmax, datatype) + call realloc (w, nmax, datatype) + + call gt_sets (gt, GTTITLE, Memc[fname]) + + switch (datatype) { + case TY_REAL: + call cf_fitr (ic, gt, Memr[x], Memr[y], Memr[w], nvalues, + nmax, Memc[dev], interactive, ofmt, power) + case TY_DOUBLE: + call cf_fitd (ic, gt, Memd[x], Memd[y], Memd[w], nvalues, + nmax, Memc[dev], interactive, ofmt, power) + } + + call flush (STDOUT) + call mfree (x, datatype) + call mfree (y, datatype) + call mfree (w, datatype) + + if (redir) + break + } + + switch (datatype) { + case TY_REAL: + call ic_closer (ic) + case TY_DOUBLE: + call ic_closed (ic) + } + + if (!redir) + call imtclose (fd) + call gt_free (gt) + call sfree (sp) +end + + +define IMAGE_OP 1 +define LIST_OP 2 + +# CF_OPERAND -- Determine whether the operand argument is an image section +# or a list. If the string is STDIN, it is a list; if a subscript is +# present, it is an image; otherwise we must test whether or not it is a +# binary file and make the decision based on that. + +int procedure cf_operand (operand) + +char operand[ARB] # Input list + +int first, last, ip +int access(), strncmp() + +begin + # Strip off any whitespace at the beginning or end of the string. + for (ip=1; IS_WHITE(operand[ip]); ip=ip+1) + ; + first = ip + for (last=ip; operand[ip] != EOS; ip=ip+1) + if (!IS_WHITE(operand[ip])) + last = ip + + if (first == last) + return (LIST_OP) + else if (strncmp (operand[first], "STDIN", 5) == 0) + return (LIST_OP) + else if (operand[last] == ']') + return (IMAGE_OP) + else if (access (operand, 0, TEXT_FILE) == YES) + return (LIST_OP) + else + return (IMAGE_OP) +end + +define SZ_BUF 1000 + +# CF_RLIST -- Read a list of two dimensional data pairs into two type +# datatype arrays in memory. Return pointers to the arrays and a count of the +# number of pixels. + +int procedure cf_rlist (fname, x, y, w, weighting, datatype) + +char fname[ARB] # Name of list file +pointer x # Pointer to x data values (returned) +pointer y # Pointer to y data values (returned) +pointer w # Pointer to weight values (returned) +int weighting # Type of weighting +int datatype # Datatype of x and Y values + +int buflen, n, fd, ncols, lineno +pointer sp, lbuf, ip + +int getline(), nscan(), open() +real cf_divzr(), cfz_divzr() +double cf_divzd(), cfz_divzd() +extern cf_divzr(), cf_divzd() +extern cfz_divzr(), cfz_divzd() +errchk open, sscan, getline, malloc + +begin + call smark (sp) + call salloc (lbuf, SZ_LINE, TY_CHAR) + + fd = open (fname, READ_ONLY, TEXT_FILE) + + n = 0 + ncols = 0 + lineno = 0 + + while (getline (fd, Memc[lbuf]) != EOF) { + # Skip comment lines and blank lines. + lineno = lineno + 1 + if (Memc[lbuf] == '#') + next + for (ip=lbuf; IS_WHITE(Memc[ip]); ip=ip+1) + ; + if (Memc[ip] == '\n' || Memc[ip] == EOS) + next + + if (n == 0) { + buflen = SZ_BUF + iferr { + call malloc (x, buflen, datatype) + call malloc (y, buflen, datatype) + call malloc (w, buflen, datatype) + } then + call erract (EA_FATAL) + } else if (n + 1 > buflen) { + buflen = buflen + SZ_BUF + call realloc (x, buflen, datatype) + call realloc (y, buflen, datatype) + call realloc (w, buflen, datatype) + } + + # Decode the points to be plotted. + call sscan (Memc[ip]) + switch (datatype) { + case TY_REAL: + call gargr (Memr[x+n]) + call gargr (Memr[y+n]) + call gargr (Memr[w+n]) + case TY_DOUBLE: + call gargd (Memd[x+n]) + call gargd (Memd[y+n]) + call gargd (Memd[w+n]) + } + + # The first line determines whether we have an x,y list or a + # y-list. It is an error if only one value can be decoded when + # processing a two column list. + + if (ncols == 0 && nscan() > 0) + ncols = nscan() + + switch (nscan()) { + case 0: + call eprintf ("no args; %s, line %d: %s\n") + call pargstr (fname) + call pargi (lineno) + call pargstr (Memc[lbuf]) + next + case 1: + if (ncols >= 2) { + call eprintf ("only 1 arg; %s, line %d: %s\n") + call pargstr (fname) + call pargi (lineno) + call pargstr (Memc[lbuf]) + next + } else { + switch (datatype) { + case TY_REAL: + Memr[y+n] = Memr[x+n] + Memr[x+n] = n + 1.0 + Memr[w+n] = 1.0 + case TY_DOUBLE: + Memd[y+n] = Memd[x+n] + Memd[x+n] = n + 1.0 + Memd[w+n] = 1.0d0 + } + } + case 2: + if (ncols == 3) { + call eprintf ("only 2 args; %s, line %d: %s\n") + call pargstr (fname) + call pargi (lineno) + call pargstr (Memc[lbuf]) + next + } else { + switch (datatype) { + case TY_REAL: + Memr[w+n] = 1.0 + case TY_DOUBLE: + Memd[w+n] = 1.0d0 + } + } + } + + n = n + 1 + } + + call realloc (x, n, datatype) + call realloc (y, n, datatype) + call realloc (w, n, datatype) + + switch (weighting) { + case CF_UNIFORM: + if (datatype == TY_REAL) + call amovkr (1.0, Memr[w], n) + else + call amovkd (1.0d0, Memd[w], n) + case CF_USER: + ; + case CF_STATISTICAL: + if (datatype == TY_REAL) { + call aabsr (Memr[y], Memr[w], n) + call arczr (1.0, Memr[w], Memr[w], n, cf_divzr) + } else { + call aabsd (Memd[y], Memd[w], n) + call arczd (1.0d0, Memd[w], Memd[w], n, cf_divzd) + } + case CF_INSTRUMENTAL: + if (datatype == TY_REAL) { + call apowkr (Memr[w], 2, Memr[w], n) + call arczr (1.0, Memr[w], Memr[w], n, cfz_divzr) + } else { + call apowkd (Memd[w], 2, Memd[w], n) + call arczd (1.0d0, Memd[w], Memd[w], n, cfz_divzd) + } + } + + call close (fd) + call sfree (sp) + return (n) +end + + +# CF_RIMAGE -- Read an image section and compute the projection about +# one dimension, producing x and y vectors as output. + +int procedure cf_rimage (imsect, axis, x, y, w, weighting, datatype) + +char imsect[ARB] # Name of image section +pointer x # Pointer to x data values +pointer y # Pointer to y data values +pointer w # Pointer to weight values +int weighting # Type of weighting +int axis # Axis about which projection is taken +int datatype # Datatype of data values + +int npix +pointer im +pointer immap() +errchk immap, im_projectionr, im_projectiond, malloc + +begin + im = immap (imsect, READ_ONLY, 0) + + if (axis < 1 || axis > IM_NDIM(im)) + call error (2, "Attempt to take projection over nonexistent axis") + npix = IM_LEN(im,axis) + + call malloc (x, npix, datatype) + call malloc (y, npix, datatype) + call malloc (w, npix, datatype) + + switch (datatype) { + case TY_REAL: + call im_projectionr (im, Memr[x], Memr[y], Memr[w], npix, weighting, + axis) + case TY_DOUBLE: + call im_projectiond (im, Memd[x], Memd[y], Memd[w], npix, weighting, + axis) + } + + call imunmap (im) + return (npix) +end + + +# CF_DIVZR -- Procedure to return a real number in case of a divide by zero. + +real procedure cf_divzr (a) + +real a # real number + +begin + return (a) +end + + +# CF_DIVZD -- Procedure to return a double number in case of a divide by zero. + +double procedure cf_divzd (a) + +double a # double precision number number + +begin + return (a) +end + + +# CFZ_DIVZR -- Procedure to return a real number in case of a divide by zero. + +real procedure cfz_divzr (a) + +real a # real number + +begin + return (0.0) +end + + +# CFZ_DIVZD -- Procedure to return a double number in case of a divide by zero. + +double procedure cfz_divzd (a) + +double a # double precision number number + +begin + return (0.0d0) +end -- cgit