diff options
Diffstat (limited to 'pkg/xtools/icfit/icfvshow.gx')
| -rw-r--r-- | pkg/xtools/icfit/icfvshow.gx | 164 |
1 files changed, 164 insertions, 0 deletions
diff --git a/pkg/xtools/icfit/icfvshow.gx b/pkg/xtools/icfit/icfvshow.gx new file mode 100644 index 00000000..458c0664 --- /dev/null +++ b/pkg/xtools/icfit/icfvshow.gx @@ -0,0 +1,164 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <math/curfit.h> +include "icfit.h" +include "names.h" + +# IC_FVSHOW -- Show fit parameters in verbose mode. + +procedure ic_fvshow$t (ic, cv, x, y, wts, npts, fd) + +pointer ic # ICFIT pointer +pointer cv # Curfit pointer +PIXEL x[ARB] # Ordinates +PIXEL y[ARB] # Abscissas +PIXEL wts[ARB] # Weights +int npts # Number of data points +int fd # Output descriptor + +int i, n, deleted, ncoeffs +PIXEL chisqr, rms +pointer sp, fit, wts1, coeffs, errors + +int $tcvstati() +PIXEL ic_rms$t() + +begin + # Do the standard ic_show option, then add on the verbose part. + call ic_fshow (ic, fd) + + if (npts == 0) { + call eprintf ("# Incomplete output - no data points for fit\n") + return + } + + # Determine the number of coefficients and allocate memory. + + ncoeffs = $tcvstati (cv, CVNCOEFF) + call smark (sp) + call salloc (coeffs, ncoeffs, TY_PIXEL) + call salloc (errors, ncoeffs, TY_PIXEL) + + if (npts == IC_NFIT(ic)) { + # Allocate memory for the fit. + + n = npts + call salloc (fit, n, TY_PIXEL) + call salloc (wts1, n, TY_PIXEL) + + # Eliminate rejected points and count deleted points. + + call amov$t (wts, Mem$t[wts1], n) + if (IC_NREJECT(ic) > 0) { + do i = 1, npts { + if (Memi[IC_REJPTS(ic)+i-1] == YES) + Mem$t[wts1+i-1] = 0. + } + } + deleted = 0 + do i = 1, n { + if (wts[i] == 0.) + deleted = deleted + 1 + } + + # Get the coefficients and compute the errors. + + call $tcvvector (cv, x, Mem$t[fit], n) + call $tcvcoeff (cv, Mem$t[coeffs], ncoeffs) + call $tcverrors (cv, y, Mem$t[wts1], Mem$t[fit], n, chisqr, + Mem$t[errors]) + rms = ic_rms$t (x, y, Mem$t[fit], Mem$t[wts1], n) + } else { + # Allocate memory for the fit. + + n = IC_NFIT(ic) + call salloc (fit, n, TY_PIXEL) + call salloc (wts1, n, TY_PIXEL) + + # Eliminate rejected points and count deleted points. + + call amov$t (Mem$t[IC_WTSFIT(ic)], Mem$t[wts1], n) + if (IC_NREJECT(ic) > 0) { + do i = 1, npts { + if (Memi[IC_REJPTS(ic)+i-1] == YES) + Mem$t[wts1+i-1] = 0. + } + } + deleted = 0 + do i = 1, n { + if (wts[i] == 0.) + deleted = deleted + 1 + } + + # Get the coefficients and compute the errors. + + call $tcvvector (cv, Mem$t[IC_XFIT(ic)], Mem$t[fit], n) + rms = ic_rms$t (Mem$t[IC_XFIT(ic)], Mem$t[IC_YFIT(ic)], + Mem$t[fit], Mem$t[wts1], n) + call $tcvcoeff (cv, Mem$t[coeffs], ncoeffs) + call $tcverrors (cv, Mem$t[IC_YFIT(ic)], Mem$t[wts1], Mem$t[fit], + n, chisqr, Mem$t[errors]) + } + + # Print the error analysis. + + call fprintf (fd, "# total points = %d\n# sample points = %d\n") + call pargi (npts) + call pargi (n) + call fprintf (fd, "# nrejected = %d\n# deleted = %d\n") + call pargi (IC_NREJECT(ic)) + call pargi (deleted) + call fprintf (fd, "# RMS = %10.7g\n") + call parg$t (rms) + call fprintf (fd, "# square root of reduced chi square = %10.7g\n") + call parg$t (sqrt (chisqr)) + + call fprintf (fd, "# \t coefficent\t error\n") + do i = 1, ncoeffs { + call fprintf (fd, "# \t%14.7e\t%14.7e\n") + call parg$t (Mem$t[coeffs+i-1]) + call parg$t (Mem$t[errors+i-1]) + } + + # Free allocated memory. + + call sfree (sp) +end + + +# IC_FXYSHOW -- List data as x, y, fit, weight lines on output. + +procedure ic_fxyshow$t (ic, cv, x, y, w, npts, fd) + +pointer ic # ICFIT pointer +pointer cv # Pointer to curfit structure +PIXEL x[npts] # Array of x data values +PIXEL y[npts] # Array of y data values +PIXEL w[npts] # Array of weight data values +int npts # Number of data values +int fd # Output file descriptor + +int i +PIXEL $tcveval() + +begin + # List the data being fit (not necessarily the input data). + call fprintf (fd, "# X Y Y FIT WEIGHT\n") + if (npts == IC_NFIT(ic)) { + do i = 1, npts { + call fprintf (fd, "%8g %8g %8g %8g\n") + call parg$t (x[i]) + call parg$t (y[i]) + call parg$t ($tcveval (cv, x[i])) + call parg$t (w[i]) + } + } else { + do i = 1, IC_NFIT(ic) { + call fprintf (fd, "%8g %8g %8g %8g\n") + call parg$t (Mem$t[IC_XFIT(ic)+i-1]) + call parg$t (Mem$t[IC_YFIT(ic)+i-1]) + call parg$t ($tcveval (cv, Mem$t[IC_XFIT(ic)+i-1])) + call parg$t (Mem$t[IC_WTSFIT(ic)+i-1]) + } + } +end |