aboutsummaryrefslogtreecommitdiff
path: root/math/surfit/sf_feval.x
diff options
context:
space:
mode:
Diffstat (limited to 'math/surfit/sf_feval.x')
-rw-r--r--math/surfit/sf_feval.x280
1 files changed, 280 insertions, 0 deletions
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
generated by cgit (git 2.34.1) at 2025-09-20 07:11:45 -0400