From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- math/surfit/sf_feval.x | 280 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 280 insertions(+) create mode 100644 math/surfit/sf_feval.x (limited to 'math/surfit/sf_feval.x') diff --git a/math/surfit/sf_feval.x b/math/surfit/sf_feval.x new file mode 100644 index 00000000..58b2f765 --- /dev/null +++ b/math/surfit/sf_feval.x @@ -0,0 +1,280 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +# SF_EVCHEB -- Procedure to evaluate a Chebyshev polynomial assuming that +# the coefficients have been calculated. + +procedure sf_evcheb (coeff, x, y, zfit, npts, xterms, xorder, yorder, k1x, k2x, + k1y, k2y) + +real coeff[ARB] # 1D array of coefficients +real x[npts] # x values of points to be evaluated +real y[npts] +real zfit[npts] # the fitted points +int npts # number of points to be evaluated +int xterms # cross terms ? +int xorder,yorder # order of the polynomials in x and y +real k1x, k2x # normalizing constants +real k1y, k2y + +int i, k, j +int ytorder, cptr +pointer sp +pointer xb, yb, accum +pointer ybzptr, ybptr, xbzptr + +begin + # fit a constant + if (xorder == 1 && yorder == 1) { + call amovkr (coeff[1], zfit, npts) + return + } + + # allocate temporary space for the basis functions + call smark (sp) + call salloc (xb, xorder * npts, TY_REAL) + call salloc (yb, yorder * npts, TY_REAL) + call salloc (accum, npts, TY_REAL) + + # calculate basis functions + call sf_bcheb (x, npts, xorder, k1x, k2x, Memr[xb]) + call sf_bcheb (y, npts, yorder, k1y, k2y, Memr[yb]) + + # clear the accumulator + call aclrr (zfit, npts) + + # accumulate the values + cptr = 0 + ybzptr = yb - 1 + xbzptr = xb - 1 + ytorder = yorder + + do i = 1, xorder { + call aclrr (Memr[accum], npts) + ybptr = ybzptr + do k = 1, ytorder { + do j = 1, npts + Memr[accum+j-1] = Memr[accum+j-1] + coeff[cptr+k] * + Memr[ybptr+j] + ybptr = ybptr + npts + } + do j = 1, npts + zfit[j] = zfit[j] + Memr[accum+j-1] * Memr[xbzptr+j] + + if (xterms == NO) + ytorder = 1 + + cptr = cptr + yorder + xbzptr = xbzptr + npts + } + + # free temporary space + call sfree (sp) +end + + +# SF_EVLEG -- Procedure to evaluate a Chebyshev polynomial assuming that +# the coefficients have been calculated. + +procedure sf_evleg (coeff, x, y, zfit, npts, xterms, xorder, yorder, k1x, k2x, + k1y, k2y) + +real coeff[ARB] # 1D array of coefficients +real x[npts] # x values of points to be evaluated +real y[npts] +real zfit[npts] # the fitted points +int npts # number of points to be evaluated +int xterms # cross terms ? +int xorder,yorder # order of the polynomials in x and y +real k1x, k2x # normalizing constants +real k1y, k2y + +int i, k, j +int ytorder, cptr +pointer sp +pointer xb, yb, accum +pointer ybzptr, ybptr, xbzptr + +begin + # fit a constant + if (xorder == 1 && yorder == 1) { + call amovkr (coeff[1], zfit, npts) + return + } + + # allocate temporary space for the basis functions + call smark (sp) + call salloc (xb, xorder * npts, TY_REAL) + call salloc (yb, yorder * npts, TY_REAL) + call salloc (accum, npts, TY_REAL) + + # calculate basis functions + call sf_bleg (x, npts, xorder, k1x, k2x, Memr[xb]) + call sf_bleg (y, npts, yorder, k1y, k2y, Memr[yb]) + + # clear the accumulator + call aclrr (zfit, npts) + + # accumulate the values + cptr = 0 + ybzptr = yb - 1 + xbzptr = xb - 1 + ytorder = yorder + + do i = 1, xorder { + call aclrr (Memr[accum], npts) + ybptr = ybzptr + do k = 1, ytorder { + do j = 1, npts + Memr[accum+j-1] = Memr[accum+j-1] + coeff[cptr+k] * + Memr[ybptr+j] + ybptr = ybptr + npts + } + do j = 1, npts + zfit[j] = zfit[j] + Memr[accum+j-1] * Memr[xbzptr+j] + + if (xterms == NO) + ytorder = 1 + + cptr = cptr + yorder + xbzptr = xbzptr + npts + } + + # free temporary space + call sfree (sp) +end + + +# SF_EVSPLINE3 -- Procedure to evaluate a piecewise linear spline function +# assuming that the coefficients have been calculated. + +procedure sf_evspline3 (coeff, x, y, zfit, npts, nxpieces, nypieces, k1x, k2x, + k1y, k2y) + +real coeff[ARB] # array of coefficients +real x[npts] # array of x values +real y[npts] # array of y values +real zfit[npts] # array of fitted values +int npts # number of data points +int nxpieces, nypieces # number of fitted points minus 1 +real k1x, k2x # normalizing constants +real k1y, k2y + +int i, j, k, cindex +pointer xb, xbzptr, yb, ybzptr, ybptr +pointer accum, leftx, lefty +pointer sp + +begin + # allocate temporary space for the basis functions + call smark (sp) + call salloc (xb, 4 * npts, TY_REAL) + call salloc (yb, 4 * npts, TY_REAL) + call salloc (accum, npts, TY_REAL) + call salloc (leftx, npts, TY_INT) + call salloc (lefty, npts, TY_INT) + + # calculate basis functions + call sf_bspline3 (x, npts, nxpieces, k1x, k2x, Memr[xb], Memi[leftx]) + call sf_bspline3 (y, npts, nypieces, k1y, k2y, Memr[yb], Memi[lefty]) + + # set up the indexing + call amulki (Memi[leftx], (nypieces+4), Memi[leftx], npts) + call aaddi (Memi[leftx], Memi[lefty], Memi[lefty], npts) + + # clear the accumulator + call aclrr (zfit, npts) + + # accumulate the values + + ybzptr = yb - 1 + xbzptr = xb - 1 + + do i = 1, 4 { + call aclrr (Memr[accum], npts) + ybptr = ybzptr + do k = 1, 4 { + do j = 1, npts { + cindex = k + Memi[lefty+j-1] + Memr[accum+j-1] = Memr[accum+j-1] + coeff[cindex] * + Memr[ybptr+j] + } + ybptr = ybptr + npts + } + do j = 1, npts + zfit[j] = zfit[j] + Memr[accum+j-1] * Memr[xbzptr+j] + + xbzptr = xbzptr + npts + call aaddki (Memi[lefty], (nypieces+4), Memi[lefty], npts) + } + + # free temporary space + call sfree (sp) +end + + +# SF_EVSPLINE1 -- Procedure to evaluate a piecewise linear spline function +# assuming that the coefficients have been calculated. + +procedure sf_evspline1 (coeff, x, y, zfit, npts, nxpieces, nypieces, k1x, k2x, + k1y, k2y) + +real coeff[ARB] # array of coefficients +real x[npts] # array of x values +real y[npts] # array of y values +real zfit[npts] # array of fitted values +int npts # number of data points +int nxpieces, nypieces # number of fitted points minus 1 +real k1x, k2x # normalizing constants +real k1y, k2y + +int i, j, k, cindex +pointer xb, xbzptr, yb, ybzptr, ybptr +pointer accum, leftx, lefty +pointer sp + +begin + # allocate temporary space for the basis functions + call smark (sp) + call salloc (xb, 2 * npts, TY_REAL) + call salloc (yb, 2 * npts, TY_REAL) + call salloc (accum, npts, TY_REAL) + call salloc (leftx, npts, TY_INT) + call salloc (lefty, npts, TY_INT) + + # calculate basis functions + call sf_bspline1 (x, npts, nxpieces, k1x, k2x, Memr[xb], Memi[leftx]) + call sf_bspline1 (y, npts, nypieces, k1y, k2y, Memr[yb], Memi[lefty]) + + # set up the indexing + call amulki (Memi[leftx], (nypieces+2), Memi[leftx], npts) + call aaddi (Memi[leftx], Memi[lefty], Memi[lefty], npts) + + # clear the accumulator + call aclrr (zfit, npts) + + # accumulate the values + + ybzptr = yb - 1 + xbzptr = xb - 1 + + do i = 1, 2 { + call aclrr (Memr[accum], npts) + ybptr = ybzptr + do k = 1, 2 { + do j = 1, npts { + cindex = k + Memi[lefty+j-1] + Memr[accum+j-1] = Memr[accum+j-1] + coeff[cindex] * + Memr[ybptr+j] + } + ybptr = ybptr + npts + } + do j = 1, npts + zfit[j] = zfit[j] + Memr[accum+j-1] * Memr[xbzptr+j] + + xbzptr = xbzptr + npts + call aaddki (Memi[lefty], (nypieces+2), Memi[lefty], npts) + } + + # free temporary space + call sfree (sp) +end -- cgit