diff options
Diffstat (limited to 'math/curfit/cvacpts.gx')
| -rw-r--r-- | math/curfit/cvacpts.gx | 186 |
1 files changed, 186 insertions, 0 deletions
diff --git a/math/curfit/cvacpts.gx b/math/curfit/cvacpts.gx new file mode 100644 index 00000000..56a36cb2 --- /dev/null +++ b/math/curfit/cvacpts.gx @@ -0,0 +1,186 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <math/curfit.h> + +$if (datatype == r) +include "curfitdef.h" +$else +include "dcurfitdef.h" +$endif + +# CVACPTS -- Procedure to add a set of points to the normal equations. +# The inner products of the basis functions are calculated and +# accumulated into the CV_ORDER(cv) by CV_NCOEFF(cv) matrix MATRIX. +# The main diagonal of the matrix is stored in the first row of +# MATRIX followed by the remaining non-zero diagonals. This method +# of storage is particularly efficient for the large symmetric +# banded matrices produced during spline fits. The inner product +# of the basis functions and the data ordinates are stored in the +# CV_NCOEFF(cv)-vector VECTOR. The array LEFT stores the +# indices which show which elements of MATRIX and VECTOR are +# to receive the inner products. + +$if (datatype == r) +procedure cvacpts (cv, x, y, w, npts, wtflag) +$else +procedure dcvacpts (cv, x, y, w, npts, wtflag) +$endif + +pointer cv # curve descriptor +PIXEL x[npts] # array of abcissa +PIXEL y[npts] # array of ordinates +PIXEL w[npts] # array of weights +int npts # number of data points +int wtflag # type of weighting + +int i, ii, j, k +pointer sp +pointer vzptr, vindex, mzptr, mindex, bptr, bbptr +pointer bw, rows + +begin + + # increment the number of points + CV_NPTS(cv) = CV_NPTS(cv) + npts + + # remove basis functions calculated by any previous cvrefit call + if (CV_BASIS(cv) != NULL) { + + call mfree (CV_BASIS(cv), TY_PIXEL) + call mfree (CV_WY(cv), TY_PIXEL) + + CV_BASIS(cv) = NULL + CV_WY(cv) = NULL + + if (CV_LEFT(cv) != NULL) { + call mfree (CV_LEFT(cv), TY_INT) + CV_LEFT(cv) = NULL + } + } + + # calculate weights + switch (wtflag) { + case WTS_UNIFORM: + call amovk$t (PIXEL(1.0), w, npts) + case WTS_SPACING: + if (npts == 1) + w[1] = 1. + else + w[1] = abs (x[2] - x[1]) + do i = 2, npts - 1 + w[i] = abs (x[i+1] - x[i-1]) + if (npts == 1) + w[npts] = 1. + else + w[npts] = abs (x[npts] - x[npts-1]) + case WTS_USER: + # user supplied weights + case WTS_CHISQ: + # data assumed to be scaled to photons with Poisson statistics + do i = 1, npts { + if (y[i] > PIXEL(0.0)) + w[i] = PIXEL(1.0) / y[i] + else if (y[i] < PIXEL(0.0)) + w[i] = -PIXEL(1.0) / y[i] + else + w[i] = PIXEL(0.0) + } + default: + call amovk$t (PIXEL(1.0), w, npts) + } + + + # allocate space for the basis functions + call smark (sp) + call salloc (CV_BASIS(cv), npts * CV_ORDER(cv), TY_PIXEL) + + # calculate the non-zero basis functions + switch (CV_TYPE(cv)) { + case LEGENDRE: + call $tcv_bleg (x, npts, CV_ORDER(cv), CV_MAXMIN(cv), + CV_RANGE(cv), BASIS(CV_BASIS(cv))) + case CHEBYSHEV: + call $tcv_bcheb (x, npts, CV_ORDER(cv), CV_MAXMIN(cv), + CV_RANGE(cv), BASIS(CV_BASIS(cv))) + case SPLINE3: + call salloc (CV_LEFT(cv), npts, TY_INT) + call $tcv_bspline3 (x, npts, CV_NPIECES(cv), -CV_XMIN(cv), + CV_SPACING(cv), BASIS(CV_BASIS(cv)), + LEFT(CV_LEFT(cv))) + case SPLINE1: + call salloc (CV_LEFT(cv), npts, TY_INT) + call $tcv_bspline1 (x, npts, CV_NPIECES(cv), -CV_XMIN(cv), + CV_SPACING(cv), BASIS(CV_BASIS(cv)), + LEFT(CV_LEFT(cv))) + case USERFNC: + call $tcv_buser (cv, x, npts) + } + + + # allocate temporary storage space for matrix accumulation + call salloc (bw, npts, TY_PIXEL) + call salloc (rows, npts, TY_INT) + + # one index the pointers + vzptr = CV_VECTOR(cv) - 1 + mzptr = CV_MATRIX(cv) + bptr = CV_BASIS(cv) + + switch (CV_TYPE(cv)) { + + case LEGENDRE, CHEBYSHEV, USERFNC: + + # accumulate the new right side of the matrix equation + do k = 1, CV_ORDER(cv) { + call amul$t (w, BASIS(bptr), Mem$t[bw], npts) + vindex = vzptr + k + do i = 1, npts + VECTOR(vindex) = VECTOR(vindex) + Mem$t[bw+i-1] * y[i] + bbptr = bptr + ii = 0 + do j = k, CV_ORDER(cv) { + mindex = mzptr + ii + do i = 1, npts + MATRIX(mindex) = MATRIX(mindex) + Mem$t[bw+i-1] * + BASIS(bbptr+i-1) + ii = ii + 1 + bbptr = bbptr + npts + } + bptr = bptr + npts + mzptr = mzptr + CV_ORDER(cv) + } + + case SPLINE1,SPLINE3: + + call amulki (LEFT(CV_LEFT(cv)), CV_ORDER(cv), Memi[rows], npts) + call aaddki (Memi[rows], CV_MATRIX(cv), Memi[rows], npts) + call aaddki (LEFT(CV_LEFT(cv)), vzptr, LEFT(CV_LEFT(cv)), npts) + + # accumulate the new right side of the matrix equation + do k = 1, CV_ORDER(cv) { + call amul$t (w, BASIS(bptr), Mem$t[bw], npts) + do i = 1, npts { + vindex = LEFT(CV_LEFT(cv)+i-1) + k + VECTOR(vindex) = VECTOR(vindex)+ Mem$t[bw+i-1] * y[i] + } + bbptr = bptr + ii = 0 + do j = k, CV_ORDER(cv) { + do i = 1, npts { + mindex = Memi[rows+i-1] + ii + MATRIX(mindex) = MATRIX(mindex) + Mem$t[bw+i-1] * + BASIS(bbptr+i-1) + } + ii = ii + 1 + bbptr = bbptr + npts + } + bptr = bptr + npts + call aaddki (Memi[rows], CV_ORDER(cv), Memi[rows], npts) + } + } + + # release the space + call sfree (sp) + CV_BASIS(cv) = NULL + CV_LEFT(cv) = NULL +end |