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_polyfit.x | 244 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 244 insertions(+) create mode 100644 pkg/utilities/t_polyfit.x (limited to 'pkg/utilities/t_polyfit.x') diff --git a/pkg/utilities/t_polyfit.x b/pkg/utilities/t_polyfit.x new file mode 100644 index 00000000..003935ee --- /dev/null +++ b/pkg/utilities/t_polyfit.x @@ -0,0 +1,244 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include + +define POLY_INSTRUMENTAL 1 +define POLY_UNIFORM 2 +define POLY_STATISTICAL 3 + +define MAX_ITEMS 1000 # maximum number of data elements + +# POLYFIT -- Fit a polynomial to the list of input pairs (x,y, sigmay) +# A polynomial fit of user specifiable order is made to the data. +# +# y = a0 + a1*x + a2*x**2 + ... +# +# The values for the coeficients a0 - aN are printed on the +# first line of the standard output. The uncertainties in the +# coeficients are then printed on the next line. +# +# Optionally (verbose = yes) the values for chi-square, ftest, +# and correlation coefficient are printed along with the calculated +# values from the fit for the dependent variable. +# +# If listdata = yes, then the only output will be pairs of +# X,Yc values. Yc is the value of dependent variable as Calculated +# from the fit. This option allows piping to the GRAPH task. +# +# The data are taken from STDIN, a file, or a list of files. +# In the latter case, each data file results in an independent +# set of results. + +# The routines REGRES and MATINV from Bevington are used to perfrom +# the fit. + +procedure t_polyfit() + +char fname[SZ_FNAME], weights[SZ_FNAME] +bool verbose, listdata +int filelist, order, weighting +int clpopni(), clgfil(), clgeti(), clgwrd() +bool clgetb() +define exit_ 91 + +begin + # Input can come from the standard input, a file, or a list of files. + # The following procedure makes both cases look like a list of files. + + filelist = clpopni ("input") + order = clgeti ("order") + weighting = clgwrd ("weighting", weights, SZ_FNAME, + ",instrumental,uniform,statistical,") + verbose = clgetb ("verbose") + listdata = clgetb ("listdata") + + while (clgfil (filelist, fname, SZ_FNAME) != EOF) + call pf_fitdatalist (fname, order, weighting, verbose, listdata) + + call clpcls (filelist) +end + + +# PF_FITDATALIST -- Perform polynomial fit to data from named file. + +procedure pf_fitdatalist (listfile, order, weighting, verbose, listdata) + +char listfile[SZ_FNAME] # input file +int order # polynomial order +int weighting # mode of weighting fit +bool verbose, listdata # printout options + +int i, in, item, m[50], mode, nterms, line_number +real x[MAX_ITEMS], y[MAX_ITEMS], yfit[MAX_ITEMS], sigy[MAX_ITEMS] +real a0, a[50], siga0, siga[50], r[50], rmul, chisqr, ftest +real stdev + +extern pf_fctn() +bool fp_equalr() +int fscan(), nscan(), open() +errchk open, fscan, printf + +begin + # Set term selection array for REGRES. + for (i=1; i <= 50; i=i+1) + m[i] = i + + in = open (listfile, READ_ONLY, TEXT_FILE) + + # Read successive X,Y, SIGMAY triples from the standard input, + # accumulating the values in the arrays X, Y, weight. Skip list + # elements containing less than two numbers. The maximum number of + # elements that can be read is fixed. + + item = 1 + line_number = 0 + + while ((fscan (in) != EOF) && (item < MAX_ITEMS+1)) { + line_number = line_number + 1 + + call gargr (x[item]) + call gargr (y[item]) + + # There must be two items per entry for the x,y pair. + if (nscan() < 2) { + call eprintf ("Bad entry in list on line:%d - item ignored\n") + call pargi (line_number) + next + } + + # Set undefined errors to 0.0. + if (weighting == POLY_INSTRUMENTAL) { + call gargr (sigy[item]) + if (nscan() < 3) { + call eprintf ("Undefined sigmay on line:%d - item ignored\n") + call pargi (line_number) + next + } else if (fp_equalr (sigy[item], 0.0)) { + call eprintf ("Zero-valued sigmay on line:%d - item ignored\n") + call pargi (line_number) + next + } + + } else + sigy[item] = 0 + + item = item + 1 + } + + item = item - 1 + if (item > MAX_ITEMS) { + call printf ("Number of data elements exceeded - max=%d\n") + call pargi (MAX_ITEMS) + } + + if (item <= order) { + call eprintf ("Not enough data for fit: order=%d, items=%d\n") + call pargi (order) + call pargi (item) + goto exit_ + } + + # It is necessary to scale the dependent variable values + # to 1.0 on average to minimize the dynamic range during + # matrix inversion. + + nterms = order + switch (weighting) { + case POLY_INSTRUMENTAL: + mode = 1 + case POLY_UNIFORM: + mode = 0 + case POLY_STATISTICAL: + mode = -1 + default: + mode = 0 + } + + call pf_regres (x, y, sigy, item, nterms, m, mode, yfit, + a0, a, siga0, siga, r, rmul, chisqr, ftest, pf_fctn) + + # Compute standard deviation of residuals from reduced chi-sqr. + switch (weighting) { + case POLY_STATISTICAL, POLY_INSTRUMENTAL: + stdev = 0.0 + do i = 1, item + stdev = stdev + (y[i] - yfit[i]) ** 2 + stdev = sqrt (stdev / (item - 1)) + case POLY_UNIFORM: + stdev = sqrt ((item - nterms - 1) * chisqr / (item - 1)) + default: + stdev = sqrt ((item - nterms - 1) * chisqr / (item - 1)) + } + + # Print coefficients scaled back to input Y levels + if (!listdata) { + call printf ("%12.7g") + call pargr (a0) + + for (i=1; i <= order; i=i+1) { + call printf (" %12.7g") + call pargr (a[i]) + } + call printf ("\n") + + # Print sigmas. + call printf ("%12.7g") + call pargr (siga0) + + for (i=1; i <= order; i=i+1) { + call printf (" %12.7g") + call pargr (siga[i]) + } + call printf ("\n") + + # If verbose option specified, also print chi-square, ftest + # correlation coefficient, standard deviation of residuals, + # number of input pairs, and calculated y-values. + # ***25Nov85,SeH - ftest is undefined when correlation = 1. + + if (verbose) { + + if (fp_equalr (ftest, -99999.)) { + call printf ("\nchi sqr: %7g ftest: UNDEF ") + call pargr (chisqr) + if (fp_equalr (rmul, -99999.)) { + call printf (" correlation: UNDEF\n") + call pargr (rmul) + } else { + call printf (" correlation: %7g\n") + call pargr (rmul) + } + } else { + call printf ("\nchi sqr: %7g ftest: %7g ") + call pargr (chisqr) + call pargr (ftest) + if (fp_equalr (rmul, -99999.)) { + call printf (" correlation: UNDEF\n") + call pargr (rmul) + } else { + call printf ("correlation: %7g\n") + call pargr (rmul) + } + } + + call printf (" nr pts: %7g std dev res: %0.6g\n") + call pargi (item) + call pargr (stdev) + + call printf ("\nx(data) y(calc) y(data) sigy(data)\n") + } + } + + if (verbose || listdata) { + for (i=1; i <= item; i=i+1) { + call printf ("%7g %7g %7g %7g\n") + call pargr (x[i]) + call pargr (yfit[i]) + call pargr (y[i]) + call pargr (sigy[i]) + } + } + +exit_ + call close (in) +end -- cgit