aboutsummaryrefslogtreecommitdiff
path: root/math/gsurfit/gs_deval.gx
diff options
context:
space:
mode:
Diffstat (limited to 'math/gsurfit/gs_deval.gx')
-rw-r--r--math/gsurfit/gs_deval.gx241
1 files changed, 241 insertions, 0 deletions
diff --git a/math/gsurfit/gs_deval.gx b/math/gsurfit/gs_deval.gx
new file mode 100644
index 00000000..38b90dac
--- /dev/null
+++ b/math/gsurfit/gs_deval.gx
@@ -0,0 +1,241 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+# GS_DPOL -- Procedure to evaluate the polynomial derivative basis functions.
+
+procedure $tgs_dpol (x, npts, order, nder, k1, k2, basis)
+
+PIXEL x[npts] # array of data points
+int npts # number of points
+int order # order of new polynomial, order = 1, constant
+int nder # order of derivative, order = 0, no derivative
+PIXEL k1, k2 # normalizing constants
+PIXEL basis[ARB] # basis functions
+
+int bptr, k, kk
+PIXEL fac
+
+begin
+ # Optimize for oth and first derivatives.
+ if (nder == 0) {
+ call $tgs_bpol (x, npts, order, k1, k2, basis)
+ return
+ } else if (nder == 1) {
+ call $tgs_bpol (x, npts, order, k1, k2, basis)
+ do k = 1, order {
+ call amulk$t(basis[1+(k-1)*npts], PIXEL (k),
+ basis[1+(k-1)*npts], npts)
+ }
+ return
+ }
+
+ # Compute the polynomials.
+ bptr = 1
+ do k = 1, order {
+ if (k == 1)
+ call amovk$t (PIXEL(1.0), basis, npts)
+ else if (k == 2)
+ call amov$t (x, basis[bptr], npts)
+ else
+ call amul$t (basis[bptr-npts], x, basis[bptr], npts)
+ bptr = bptr + npts
+ }
+
+ # Apply the derivative factor.
+ bptr = 1
+ do k = 1, order {
+ if (k == 1) {
+ fac = PIXEL(1.0)
+ do kk = 2, nder
+ fac = fac * PIXEL (kk)
+ } else {
+ fac = PIXEL(1.0)
+ do kk = k + nder - 1, k, -1
+ fac = fac * PIXEL(kk)
+ }
+ call amulk$t (basis[bptr], fac, basis[bptr], npts)
+ bptr = bptr + npts
+ }
+end
+
+
+# GS_DCHEB -- Procedure to evaluate the chebyshev polynomial derivative
+# basis functions using the usual recursion relation.
+
+procedure $tgs_dcheb (x, npts, order, nder, k1, k2, basis)
+
+PIXEL x[npts] # array of data points
+int npts # number of points
+int order # order of polynomial, order = 1, constant
+int nder # order of derivative, order = 0, no derivative
+PIXEL k1, k2 # normalizing constants
+PIXEL basis[ARB] # basis functions
+
+int i, k
+pointer fn, dfn, xnorm, bptr, fptr
+PIXEL fac
+
+begin
+ # Optimze the no derivatives case.
+ if (nder == 0) {
+ call $tgs_bcheb (x, npts, order, k1, k2, basis)
+ return
+ }
+
+ # Allocate working space for the basis functions and derivatives.
+ call calloc (fn, npts * (order + nder), TY_PIXEL)
+ call calloc (dfn, npts * (order + nder), TY_PIXEL)
+
+ # Compute the normalized x values.
+ call malloc (xnorm, npts, TY_PIXEL)
+ call alta$t (x, Mem$t[xnorm], npts, k1, k2)
+
+ # Compute the current solution.
+ bptr = fn
+ do k = 1, order + nder {
+ if (k == 1)
+ call amovk$t (PIXEL(1.0), Mem$t[bptr], npts)
+ else if (k == 2)
+ call amov$t (Mem$t[xnorm], Mem$t[bptr], npts)
+ else {
+ call amul$t (Mem$t[xnorm], Mem$t[bptr-npts], Mem$t[bptr], npts)
+ call amulk$t (Mem$t[bptr], PIXEL(2.0), Mem$t[bptr], npts)
+ call asub$t (Mem$t[bptr], Mem$t[bptr-2*npts], Mem$t[bptr], npts)
+ }
+ bptr = bptr + npts
+ }
+
+ # Compute the derivative basis functions.
+ do i = 1, nder {
+
+ # Compute the derivatives.
+ bptr = fn
+ fptr = dfn
+ do k = 1, order + nder {
+ if (k == 1)
+ call amovk$t (PIXEL(0.0), Mem$t[fptr], npts)
+ else if (k == 2) {
+ if (i == 1)
+ call amovk$t (PIXEL(1.0), Mem$t[fptr], npts)
+ else
+ call amovk$t (PIXEL(0.0), Mem$t[fptr], npts)
+ } else {
+ call amul$t (Mem$t[xnorm], Mem$t[fptr-npts], Mem$t[fptr],
+ npts)
+ call amulk$t (Mem$t[fptr], PIXEL(2.0), Mem$t[fptr], npts)
+ call asub$t (Mem$t[fptr], Mem$t[fptr-2*npts], Mem$t[fptr],
+ npts)
+ fac = PIXEL (2.0) * PIXEL (i)
+ call awsu$t (Mem$t[bptr-npts], Mem$t[fptr], Mem$t[fptr],
+ npts, fac, PIXEL(1.0))
+
+ }
+ bptr = bptr + npts
+ fptr = fptr + npts
+ }
+
+ # Make the derivatives the old solution
+ if (i < nder)
+ call amov$t (Mem$t[dfn], Mem$t[fn], npts * (order + nder))
+ }
+
+ # Copy the solution into the basis functions.
+ call amov$t (Mem$t[dfn+nder*npts], basis[1], order * npts)
+
+ call mfree (xnorm, TY_PIXEL)
+ call mfree (fn, TY_PIXEL)
+ call mfree (dfn, TY_PIXEL)
+end
+
+
+# GS_DLEG -- Procedure to evaluate the Legendre polynomial derivative basis
+# functions using the usual recursion relation.
+
+procedure $tgs_dleg (x, npts, order, nder, k1, k2, basis)
+
+PIXEL x[npts] # number of data points
+int npts # number of points
+int order # order of new polynomial, 1 is a constant
+int nder # order of derivate, 0 is no derivative
+PIXEL k1, k2 # normalizing constants
+PIXEL basis[ARB] # array of basis functions
+
+int i, k
+pointer fn, dfn, xnorm, bptr, fptr
+PIXEL ri, ri1, ri2, fac
+
+begin
+ # Optimze the no derivatives case.
+ if (nder == 0) {
+ call $tgs_bleg (x, npts, order, k1, k2, basis)
+ return
+ }
+
+ # Allocate working space for the basis functions and derivatives.
+ call calloc (fn, npts * (order + nder), TY_PIXEL)
+ call calloc (dfn, npts * (order + nder), TY_PIXEL)
+
+ # Compute the normalized x values.
+ call malloc (xnorm, npts, TY_PIXEL)
+ call alta$t (x, Mem$t[xnorm], npts, k1, k2)
+
+ # Compute the basis functions.
+ bptr = fn
+ do k = 1, order + nder {
+ if (k == 1)
+ call amovk$t (PIXEL(1.0), Mem$t[bptr], npts)
+ else if (k == 2)
+ call amov$t (Mem$t[xnorm], Mem$t[bptr], npts)
+ else {
+ ri = k
+ ri1 = (PIXEL(2.0) * ri - PIXEL(3.0)) / (ri - PIXEL(1.0))
+ ri2 = - (ri - PIXEL(2.0)) / (ri - PIXEL(1.0))
+ call amul$t (Mem$t[xnorm], Mem$t[bptr-npts], Mem$t[bptr], npts)
+ call awsu$t (Mem$t[bptr], Mem$t[bptr-2*npts], Mem$t[bptr],
+ npts, ri1, ri2)
+ }
+ bptr = bptr + npts
+ }
+
+ # Compute the derivative basis functions.
+ do i = 1, nder {
+
+ # Compute the derivatives.
+ bptr = fn
+ fptr = dfn
+ do k = 1, order + nder {
+ if (k == 1)
+ call amovk$t (PIXEL(0.0), Mem$t[fptr], npts)
+ else if (k == 2) {
+ if (i == 1)
+ call amovk$t (PIXEL(1.0), Mem$t[fptr], npts)
+ else
+ call amovk$t (PIXEL(0.0), Mem$t[fptr], npts)
+ } else {
+ ri = k
+ ri1 = (PIXEL(2.0) * ri - PIXEL(3.0)) / (ri - PIXEL(1.0))
+ ri2 = - (ri - PIXEL(2.0)) / (ri - PIXEL(1.0))
+ call amul$t (Mem$t[xnorm], Mem$t[fptr-npts], Mem$t[fptr],
+ npts)
+ call awsu$t (Mem$t[fptr], Mem$t[fptr-2*npts], Mem$t[fptr],
+ npts, ri1, ri2)
+ fac = ri1 * PIXEL (i)
+ call awsu$t (Mem$t[bptr-npts], Mem$t[fptr], Mem$t[fptr],
+ npts, fac, PIXEL(1.0))
+
+ }
+ bptr = bptr + npts
+ fptr = fptr + npts
+ }
+
+ # Make the derivatives the old solution
+ if (i < nder)
+ call amov$t (Mem$t[dfn], Mem$t[fn], npts * (order + nder))
+ }
+
+ # Copy the solution into the basis functions.
+ call amov$t (Mem$t[dfn+nder*npts], basis[1], order * npts)
+
+ call mfree (xnorm, TY_PIXEL)
+ call mfree (fn, TY_PIXEL)
+ call mfree (dfn, TY_PIXEL)
+end
generated by cgit (git 2.34.1) at 2025-09-20 07:09:17 -0400