diff options
Diffstat (limited to 'math/gsurfit/gsrejectr.x')
| -rw-r--r-- | math/gsurfit/gsrejectr.x | 153 |
1 files changed, 153 insertions, 0 deletions
diff --git a/math/gsurfit/gsrejectr.x b/math/gsurfit/gsrejectr.x new file mode 100644 index 00000000..fea86cef --- /dev/null +++ b/math/gsurfit/gsrejectr.x @@ -0,0 +1,153 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <math/gsurfit.h> +include "gsurfitdef.h" + +# GSREJ-- Procedure to reject a point from the normal equations. +# The inner products of the basis functions are calculated and +# accumulated into the GS_NCOEFF(sf) ** 2 matrix MATRIX. +# The main diagonal of the matrix is stored in the first row of +# MATRIX followed by the remaining non-zero diagonals. +# The inner product +# of the basis functions and the data ordinates are stored in the +# NCOEFF(sf)-vector VECTOR. + +procedure gsrej (sf, x, y, z, w, wtflag) + +pointer sf # surface descriptor +real x # x value +real y # y value +real z # z value +real w # weight +int wtflag # type of weighting + +int ii, j, k, l +int maxorder, xorder, xxorder, xindex, yindex, ntimes +pointer sp, vzptr, mzptr, xbptr, ybptr +real byw, bw + +begin + # increment the number of points + GS_NPTS(sf) = GS_NPTS(sf) - 1 + + # remove basis functions calculated by any previous gsrefit call + if (GS_XBASIS(sf) != NULL || GS_YBASIS(sf) != NULL) { + + if (GS_XBASIS(sf) != NULL) + call mfree (GS_XBASIS(sf), TY_REAL) + GS_XBASIS(sf) = NULL + if (GS_YBASIS(sf) != NULL) + call mfree (GS_YBASIS(sf), TY_REAL) + GS_YBASIS(sf) = NULL + if (GS_WZ(sf) != NULL) + call mfree (GS_WZ(sf), TY_REAL) + GS_WZ(sf) = NULL + } + + # calculate weight + switch (wtflag) { + case WTS_UNIFORM: + w = 1. + case WTS_USER: + # user supplied weights + default: + w = 1. + } + + # allocate space for the basis functions + call smark (sp) + + # calculate the non-zero basis functions + switch (GS_TYPE(sf)) { + case GS_LEGENDRE: + call salloc (GS_XBASIS(sf), GS_XORDER(sf), TY_REAL) + call salloc (GS_YBASIS(sf), GS_YORDER(sf), TY_REAL) + call rgs_b1leg (x, GS_XORDER(sf), GS_XMAXMIN(sf), + GS_XRANGE(sf), XBASIS(GS_XBASIS(sf))) + call rgs_b1leg (y, GS_YORDER(sf), GS_YMAXMIN(sf), + GS_YRANGE(sf), YBASIS(GS_YBASIS(sf))) + case GS_CHEBYSHEV: + call salloc (GS_XBASIS(sf), GS_XORDER(sf), TY_REAL) + call salloc (GS_YBASIS(sf), GS_YORDER(sf), TY_REAL) + call rgs_b1cheb (x, GS_XORDER(sf), GS_XMAXMIN(sf), + GS_XRANGE(sf), XBASIS(GS_XBASIS(sf))) + call rgs_b1cheb (y, GS_YORDER(sf), GS_YMAXMIN(sf), + GS_YRANGE(sf), YBASIS(GS_YBASIS(sf))) + case GS_POLYNOMIAL: + call salloc (GS_XBASIS(sf), GS_XORDER(sf), TY_REAL) + call salloc (GS_YBASIS(sf), GS_YORDER(sf), TY_REAL) + call rgs_b1pol (x, GS_XORDER(sf), GS_XMAXMIN(sf), + GS_XRANGE(sf), XBASIS(GS_XBASIS(sf))) + call rgs_b1pol (y, GS_YORDER(sf), GS_YMAXMIN(sf), + GS_YRANGE(sf), YBASIS(GS_YBASIS(sf))) + default: + call error (0, "GSACCUM: Illegal curve type.") + } + + # one index the pointers + vzptr = GS_VECTOR(sf) - 1 + mzptr = GS_MATRIX(sf) - 1 + xbptr = GS_XBASIS(sf) - 1 + ybptr = GS_YBASIS(sf) - 1 + + switch (GS_TYPE(sf)) { + + case GS_LEGENDRE, GS_CHEBYSHEV, GS_POLYNOMIAL: + + maxorder = max (GS_XORDER(sf) + 1, GS_YORDER(sf) + 1) + ntimes = 0 + xorder = GS_XORDER(sf) + do l = 1, GS_YORDER(sf) { + + byw = w * YBASIS(ybptr+l) + do k = 1, xorder { + bw = byw * XBASIS(xbptr+k) + VECTOR(vzptr+k) = VECTOR(vzptr+k) - bw * z + ii = 1 + xindex = k + yindex = l + xxorder = xorder + do j = k + ntimes, GS_NCOEFF(sf) { + MATRIX(mzptr+ii) = MATRIX(mzptr+ii) - bw * + XBASIS(xbptr+xindex) * YBASIS(ybptr+yindex) + if (mod (xindex, xxorder) == 0) { + xindex = 1 + yindex = yindex + 1 + switch (GS_XTERMS(sf)) { + case GS_XNONE: + xxorder = 1 + case GS_XHALF: + if ((yindex + GS_XORDER(sf)) > maxorder) + xxorder = xxorder - 1 + default: + ; + } + } else + xindex = xindex + 1 + ii = ii + 1 + } + mzptr = mzptr + GS_NCOEFF(sf) + } + + vzptr = vzptr + xorder + ntimes = ntimes + xorder + switch (GS_XTERMS(sf)) { + case GS_XNONE: + xorder = 1 + case GS_XHALF: + if ((l + GS_XORDER(sf) + 1) > maxorder) + xorder = xorder - 1 + default: + ; + } + } + + default: + call error (0, "GSACCUM: Unknown curve type.") + } + + # release the space + call sfree (sp) + GS_XBASIS(sf) = NULL + GS_YBASIS(sf) = NULL +end |