aboutsummaryrefslogtreecommitdiff
path: root/math/curfit/curfit.sem
diff options
context:
space:
mode:
Diffstat (limited to 'math/curfit/curfit.sem')
-rw-r--r--math/curfit/curfit.sem708
1 files changed, 708 insertions, 0 deletions
diff --git a/math/curfit/curfit.sem b/math/curfit/curfit.sem
new file mode 100644
index 00000000..28058626
--- /dev/null
+++ b/math/curfit/curfit.sem
@@ -0,0 +1,708 @@
+# Semi-code for curfit.h
+
+# define the permitted types of curves
+
+define CHEBYSHEV 1
+define LEGENDRE 2
+define L2SPLINE4 3
+
+# define the weighting flags
+
+define NORMAL 1 # user enters weights
+define UNIFORM 2 # equal weights, weight 1.0
+define SPACING 3 # weigth proportional to spacing of data points
+
+define SPLINE_ORDER 4
+
+# set up the curve fitting structure
+
+define LEN_CVSTRUCT
+
+struct curfit {
+
+define CV_TYPE Memi[] # Type of curve to be fitted
+define CV_ORDER Memi[] # Order of the fit
+define CV_NPIECES Memi[] # Number of polynomial pieces, spline
+define CV_NCOEFF Memi[] # Number of coefficients
+define CV_XMAX Memr[] # Maximum x value
+define CV_XMIN Memr[] # Minimum x value
+define CV_RANGE Memr[] # Xmax minus xmin
+define CV_MAXMIN Memr[] # Xmax plus xmin
+define CV_SPACING Memr[] # Knot spacing for splines
+define CV_YNORM Memr[] # Norm of the Y vector
+define CV_NPTS Memi[] # Number of data points
+
+define CV_MATRIX Memi[] # Pointer to original matrix
+define CV_CHOFAC Memi[] # Pointer to Cholesky factorization
+define CV_BASIS Memi[] # Pointer to basis functions
+define CV_VECTOR Memi[] # Pointer to vector
+define CV_COEFF Memi[] # Pointer to coefficient vector
+define CV_LEFT Memi[] #
+
+}
+
+# matrix and vector element definitions
+
+define MATRIX Memr[$1+($2-1)*NCOEFF(cv)] # Matrix element
+define CHOFAC Memr[$1+($2-1)*NCOEFF(cv)] # Triangular matrix
+define VECTOR Memr[$1+$2] # Right side
+define COEFF Memr[$1+$2] # Coefficient vector
+define LEFT Memi[$1+$2]
+
+# matrix and vector definitions
+
+define MAT Memr[$1]
+define CHO Memr[$1]
+define VECT Memr[$1]
+define COF Memr[$1]
+
+# semi-code for the initialization procedure
+
+include "curfit.h"
+
+# CVINIT -- Procedure to set up the curve descriptor.
+
+procedure cvinit (cv, curve_type, order, xmin, xmax)
+
+pointer cv # pointer to curve descriptor structure
+int curve_type # type of curve to be fitted
+int order # order of curve to be fitted, or in the case of the
+ # spline the number of polynomial pieces to be fit
+real xmin # minimum value of x
+real xmax # maximum value of x
+
+begin
+ # allocate space for the curve descriptor
+ call smark (sp)
+ call salloc (cv, LEN_CVSTRUCT, TY_STRUCT)
+
+ if (order < 1)
+ call error (0, "CVINIT: Illegal order.")
+
+ if (xmax <= xmin)
+ call error (0, "CVINIT: xmax <= xmin.")
+
+ switch (curve_type) {
+ case CHEBYSHEV, LEGENDRE:
+ CV_ORDER(cv) = order
+ CV_NCOEFF(CV) = order
+ CV_RANGE(cv) = xmax - xmin
+ CV_MAXMIN(cv) = xmax + xmin
+ case L2SPLINE4:
+ CV_ORDER(cv) = SPLINE_ORDER
+ CV_NCOEFF(cv) = order + SPLINE_ORDER - 1
+ CV_NPIECES(cv) = order
+ CV_SPACING(cv) = (xmax - xmin) / order
+ default:
+ call error (0, "CVINIT: Unknown curve type.")
+ }
+
+ CV_TYPE(cv) = curve_type
+ CV_XMIN(cv) = xmin
+ CV_XMAX(cv) = xmax
+
+ # allocate space for the matrix and vectors
+ call calloc (CV_MATRIX(cv), CV_ORDER(cv)*CV_NCOEFF(cv), TY_REAL)
+ call calloc (CV_CHOFAC(cv), CV_ORDER(cv)*CV_NCOEFF(cv), TY_REAL)
+ call calloc (CV_VECTOR(cv), CV_NCOEFF(cv), TY_REAL)
+ call calloc (CV_COEFF(cv), CV_NCOEFF(cv), TY_REAL)
+
+ # initialize pointer to basis functions to null
+ CV_BASIS(cv) = NULL
+
+ CV_NPTS(cv) = 0
+ CV_YNORM(cv) = 0.
+end
+
+# semi-code for cvaccum
+
+include "curfit.h"
+
+# CVACCUM -- Procedure to add a single point to the data set.
+
+procedure cvaccum (cv, x, y, w, wtflag)
+
+pointer cv # curve descriptor
+real x # x value
+real y # y value
+real w # weight of the data point
+int wtflag # type of weighting desired
+
+begin
+ # calculate the weights
+ switch (wtflag) {
+ case UNIFORM:
+ w = 1.0
+ case NORMAL, SPACING: # problem spacing
+ default:
+ w = 1.0
+ }
+
+ # caculate all non-zero basis functions for a given data point
+ switch (CV_TYPE(cv)) {
+ case CHEBYSHEV:
+ left = 1
+ call chebyshev (cv, x, basis)
+ case LEGENDRE:
+ left = 1
+ call legendre (cv, x, basis)
+ case L2SPLINE4:
+ call l2spline4 (cv, x, left, basis)
+ }
+
+ # accumulate into the matrix
+ leftm1 = left - 1
+ vptr = CV_VECTOR(cv) - 1
+ do i = 1, CV_ORDER(cv) {
+ bw = basis[i] * w
+ jj = leftm1 + i
+ mptr = CV_MATRIX(cv) + jj - 1
+ VECTOR(vptr, jj) = VECTOR(vptr, jj) + bw * y
+ ii = 1
+ do j = i, CV_ORDER(cv) {
+ MATRIX(mptr, ii) = MATRIX(mptr, ii) + basis[j] * bw
+ ii = ii + 1
+ }
+ }
+
+ CV_NPTS(cv) = CV_NPTS(cv) + 1
+ CV_YNORM(cv) = CV_YNORM(cv) + w * y * y
+end
+
+# semi-code for cvreject
+
+include "curfit.h"
+
+# CVREJECT -- Procedure to subtract a single datapoint from the data set
+# to be fitted.
+
+procedure cvreject (cv, x, y, w)
+
+pointer cv # curve fitting image descriptor
+real x # x value
+real y # y value
+real w # weight of the data point
+
+begin
+ # caculate all type non-zero basis functions for a given data point
+ switch (CV_TYPE(cv)) {
+ case CHEBYSHEV:
+ left = 1
+ call chebyshev (cv, x, basis)
+ case LEGENDRE:
+ left = 1
+ call legendre (cv, x, basis)
+ case L2SPLINE4:
+ call l2spline4 (cv, x, left, basis)
+ }
+
+ # subtract the data point from the matrix
+ leftm1 = left - 1
+ vptr = CV_VECTOR(cv) - 1
+ do i = 1, CV_ORDER(cv) {
+ bw = basis[i] * w
+ jj = leftm1 + i
+ mptr = CV_MATRIX(cv) + jj - 1
+ VECTOR(vptr, jj) = VECTOR(vptr, jj) - bw * y
+ ii = 1
+ do j = i, CV_ORDER(cv) {
+ MATRIX(mptr, ii) = MATRIX(mptr, ii) - basis[j] * bw
+ ii = ii + 1
+ }
+ }
+
+ CV_NPTS(cv) = CV_NPTS(cv) - 1
+ CV_NORM(cv) = CV_NORM(cv) - w * y * y
+end
+
+# semi-code for cvsolve
+
+include "curfit.h"
+
+# CVSOLVE -- Procedure to solve a matrix equation of the form Ax = B.
+# The Cholesky factorization of matrix A is calculated in the first
+# step, followed by forward and back substitution to solve for the vector
+# x.
+
+procedure cvsolve (cv, ier)
+
+pointer cv # pointer to the image descriptor structure
+int ier # ier = 0, everything OK
+ # ier = 1, matrix is singular
+
+begin
+ # solve matrix by adapting Deboor's bchfac.f and bchslv.f routines
+ # so that the original matrix and vector are not destroyed
+
+ call chofac (MAT(CV_MATRIX(cv)), CV_ORDER(cv), CV_NCOEFF(cv),
+ CHO(CV_CHOFAC(cv)), ier)
+ call choslv (CHO(CV_CHOFAC(cv)), CV_ORDER(cv), CV_NCOEFF(cv),
+ VECT(CV_VECTOR(cv)), COF(CV_COEFF(cv)))
+end
+
+# semi-code for cvfit
+
+include "curfit.h"
+
+# CVFIT -- Procedure to fit a curve to an array of data points x and y with
+# weights w.
+
+procedure cvfit (x, y, w, npts, wtflag, ier)
+
+real x[npts] # array of abcissa
+real y[npts] # array of ordinates
+real w[npts] # array of weights
+int wtflag # type of weighting
+int ier
+
+begin
+ # calculate weights
+ switch (wtflag) {
+ case UNIFORM:
+ call amovkr (1., w, npts)
+ case SPACING:
+ w[1] = x[2] - x[1] # check for npts > 1
+ do i = 2, npts - 1
+ w[i] = x[i+1] - x[i-1]
+ w[npts] = x[npts] - x[npts-1]
+ case NORMAL:
+ default:
+ call amovkr (1., w, npts)
+ }
+
+ # accumulate data points
+ do i = 1, npts {
+
+ CV_NPTS(cv) = CV_NPTS(cv) + 1
+
+ # calculate the norm of the Y vector
+ CV_YNORM(cv) = CV_YNORM(cv) + w[i] * y[i] * y[i]
+
+ # calculate non zero basis functions
+ switch (CV_TYPE(cv)) {
+ case CHEBYSHEV:
+ left = 1
+ call chebyshev (cv, x, basis)
+ case LEGENDRE:
+ left = 1
+ call legendre (cv, x, basis)
+ case L2SPLINE4:
+ call l2spline4 (cv, x, left, basis)
+ }
+
+ # accumulate the matrix
+ leftm1 = left - 1
+ vptr = CV_VECTOR(cv) - 1
+ do i = 1, CV_ORDER(cv) {
+ bw = basis[i] * w
+ jj = leftm1 + i
+ mptr = CV_MATRIX(cv) + jj - 1
+ VECTOR(vptr, jj) = VECTOR(vptr, jj) + bw * y
+ ii = 1
+ do j = i, CV_ORDER(cv) {
+ MATRIX(mptr, ii) = MATRIX(mptr, ii) + basis[j] * bw
+ ii = ii + 1
+ }
+ }
+ }
+
+ # solve the matrix
+ ier = 0
+ call chofac (MAT(CV_MATRIX(cv)), CV_ORDER(cv), CV_NCOEFF(cv),
+ CHO(CV_CHOFAC(cv)), ier)
+ call choslv (CHO(CV_CHOFAC(cv)), CV_ORDER(cv), CV_NCOEFF(cv),
+ VECT(CV_VECTOR(cv)), COF(CV_COEFF(cv)))
+end
+
+# semi-code for cvrefit
+
+include "curfit.com"
+
+# CV_REFIT -- Procedure to refit the data assuming that the x and w values do
+# not change.
+
+procedure cvrefit (cv, x, y, w, ier)
+
+pointer cv
+real x[ARB]
+real y[ARB]
+real w[ARB]
+int ier
+
+begin
+ # if first call to refit then calculate and store the basis
+ # functions
+
+ vcptr = CV_VECTOR(cv) - 1
+ do i = 1, NCOEFF(cv)
+ VECTOR(vcptr+i) = 0.
+
+ CV_YNORM(cv) = 0.
+ lptr = CV_LEFT(cv) - 1
+ bcptr = CV_BASIS(cv) - CV_NPTS(cv)
+
+ if (CV_BASIS(cv) == NULL) {
+
+ call calloc (CV_BASIS(cv), CV_NPTS(cv)*CV_ORDER(cv), TY_REAL)
+ call calloc (CV_LEFT(cv), CV_NPTS(cv), TY_INT)
+
+ do l = 1, CV_NPTS(cv) {
+ bptr = bcptr + l * CV_NPTS(cv)
+ switch (CV_TYPE(cv)) {
+ case LEGENDRE:
+ LEFT(lptr+l) = 1
+ call legendre (cv, x, BASIS(bptr))
+ case CHEBYSHEV:
+ LEFT(lptr+l) = 1
+ call chebyshev (cv, x, BASIS(bptr))
+ case L2SPLINE4:
+ call l2spline4 (cv, x, LEFT(lptr+l), BASIS(bptr))
+ }
+ }
+ }
+
+ # reset vector to zero
+
+ # accumulate right side of the matrix equation
+ do l = 1, CV_NPTS(cv) {
+
+ CV_YNORM(cv) = CV_YNORM(cv) + w[l] * y[l] * y[l]
+ leftm1 = LEFT(lptr+l) - 1
+ bptr = bcptr + l * CV_NPTS(cv)
+
+ do i = 1, CV_ORDER(cv) {
+ vptr = vcptr + leftm1 + i
+ VECTOR(vptr) = VECTOR(vptr) + BASIS(bptr) * w[l] * y[l]
+ }
+ }
+
+ # solve the matrix
+ call choslv (CHOFAC(CV_CHOFAC(cv)), CV_ORDER(cv), CV_NCOEFF(CV),
+ VECTOR(CV_VECTOR(cv)), COEFF(CV_COEFF(cv)))
+end
+
+# semi-code for cvcoeff
+
+# CVCOEFF -- Procedure to fetch the number and magnitude of the coefficients.
+
+procedure cvcoeff (cv, coeff, ncoeff)
+
+pointer cv # pointer to the curve fitting descriptor
+real coeff[ncoeff] # the coefficients of the fit
+int ncoeff # the number of coefficients
+
+begin
+ ncoeff = CV_NCOEFF(cv)
+ cptr = CV_COEFF(cv) - 1
+ do i = 1, ncoeff
+ coeff[i] = COEFF(cptr, i)
+end
+
+# semi-code for cvvector
+
+include "curfit.h"
+
+# CVVECTOR -- Procedure to evaluate the fitted curve
+
+procedure cvvector (cv, x, npts, yfit)
+
+pointer cv # pointer to the curve descriptor structure
+real x[npts] # data x values
+int npts # number of data points
+real yfit[npts] # the fitted y values
+
+begin
+ do l = 1, npts {
+
+ # calculate the non-zero basis functions
+ switch (CV_TYPE(cv) {
+ case LEGENDRE:
+ left = 1
+ call legendre (cv, x[l], XBASIS(CV_XBASIS(cv)))
+ case CHEBYSHEV:
+ left = 1
+ call chebyshev (cv, x[l], XBASIS(CV_XBASIS(cv)))
+ case L2SPLINE4:
+ call l2spline4 (cv, x[l], left, XBASIS(CV_XBASIS(cv)))
+ }
+
+ sum = 0.0
+ leftm1 = left - 1
+ cptr = CV_COEFF(cv) - 1
+ xptr = CV_XBASIS(cv) - 1
+
+ do i = 1, CV_NCOEFF(cv) {
+ jj = leftm1 + i
+ sum = sum + XBASIS(xptr + i) * COEFF(cptr + jj)
+ }
+ }
+end
+
+# semi-code for cveval
+
+include "curfit.h"
+
+# CVEVAL -- Procedure to evaluate curve at a given x
+
+real procedure cveval (cv, x)
+
+pointer cv # pointer to image descriptor structure
+real x # x value
+
+int left, leftm1, i
+pointer cptr, xptr
+real sum
+
+begin
+ switch (CV_TYPE(cv)) {
+ case CHEBYSHEV:
+ left = 1
+ call chebyshev (cv, x, XBASIS(CV_XBASIS(cv)))
+ case LEGENDRE:
+ left = 1
+ call legendre (cv, x, XBASIS(CV_XBASIS(cv)))
+ case L2SPLINE4:
+ call l2spline4 (cv, x, left, XBASIS(CV_XBASIS(cv)))
+ }
+
+ sum = 0.
+ leftm1 = left - 1
+ cptr = CV_COEFF(cv) - 1
+ xptr = CV_XBASIS(cv) - 1
+ do i = 1, CV_NCOEFF(cv) {
+ jj = leftm1 + i
+ sum = sum + XBASIS(xptr + i) * COEFF(cptr + jj)
+ }
+
+ return (sum)
+end
+
+# semi-code for cverrors
+
+include "curfit.h"
+
+# CVERRORS -- Procedure to calculate the standard deviation of the fit and the
+# standard deviations of the coefficients
+
+procedure cverrors (cv, rms, errors)
+
+pointer cv # curve descriptor
+real rms # standard deviation of data with respect to fit
+real errors[ARB] # errors in coefficients
+
+begin
+ # calculate the variance
+ rms = CV_YNORM(cv)
+ cptr = CV_COEFF(cv) - 1
+ vptr = CV_VECTOR(cv) - 1
+ do i = 1, CV_NCOEFF(cv)
+ rms = rms - COEFF(cptr, i) * VECTOR(vptr, i)
+ rms = rms / (CV_NPTS(cv) - CV_NCOEFF(cv))
+
+ # calculate the standard deviations
+ do i = 1, CV_NCOEFF(cv) {
+ do j = 1, CV_NCOEFF(cv)
+ cov[j] = 0.
+ cov[i] = 1.
+ call choslv (CHO(CV_CHOFAC(cv)), CV_ORDER(cv),
+ CV_NCOEFF(cv), cov, cov)
+ errors[i] = sqrt (cov[i] * rms)
+ }
+
+ rms = sqrt (rms)
+end
+
+# semi-code for CVFREE
+
+# CVFREE -- Procedure to free the curve descriptor
+
+procedure cvfree (cv)
+
+pointer cv
+
+begin
+ call sfree (cv)
+end
+
+include "curfit.h"
+
+# LEGENDRE -- Procedure to calculate the Legendre functions.
+
+procedure legendre (cv, x, basis)
+
+pointer cv
+real x
+real basis[ARB]
+
+begin
+ # normalize to the range x = -1. to 1.
+ xnorm = (2. * x - CV_MAXMIN(cv)) / CV_RANGE(cv)
+
+ b[1] = 1.0
+ if (CV_ORDER(cv) == 1)
+ return
+
+ b[2] = xnorm
+ if (CV_ORDER(cv) == 2)
+ return
+
+ do i = 3, CV_ORDER(cv) {
+ ri = i
+ b[i] = ((2.*ri-3.)*xnorm*b[i-1] - (ri-2.)*b[i-2]) / (ri-1.)
+ }
+end
+
+# CHEBYSHEV -- Procedure to calculate Chebyshev polynomials.
+
+procedure chebyshev (cv, x, basis)
+
+real x
+int order
+real basis[ARB]
+
+begin
+ # normalize to the range -1. to 1.
+ xnorm = (2. * x - CV_MAXMIN(cv)) / CV_RANGE(cv)
+
+ b[1] = 1.
+ if (CV_ORDER(cv) == 1)
+ return
+
+ b[2] = xnorm
+ if (CV_ORDER(cv) == 2)
+ return
+
+ do i = 3, CV_ORDER(cv) {
+ ri = i
+ b[i] = 2.*xnorm*b[i-1] - b[i-2]
+ }
+end
+
+define NPTS_SPLINE 401 # Number of points in the spline lookup table
+define INTERVALS 100 # Number of intervals per spline knot
+
+# L2SPLINE4 -- Procedure to calculate the cubic spline functions
+
+procedure (cv, x, left, basis)
+
+pointer cv
+real x
+int left
+real basis[ARB]
+
+real table[NPTS_SPLINE]
+
+# data table containing the spline
+include "table.dat"
+
+begin
+ xnorm = (x - CV_XMIN(cv)) / CV_SPACING(cv)
+ temp = min (int (xnorm), npieces - 1)
+ left = temp + 1
+ xnorm = xnorm - temp
+
+ near = int ((1. - xnorm + 0.5) * INTERVALS) + 1
+ basis[1] = table[near]
+ near = table[near] + INTERVALS
+ basis[2] = table[near]
+ near = table[near] + INTERVALS
+ basis[3] = table[near]
+ near = table[near] + INTERVALS
+ basis[4] = table[near]
+end
+
+# CHOFAC -- Routine to calculate the Cholesky factorization of a banded
+# matrix.
+
+procedure chofac (matrix, nbands, nrows, matfac, ier)
+
+real matrix[nbands, nrows]
+int nbands
+int nrows
+real matfac[nbands, nrows]
+int ier
+
+begin
+ ier = 0
+
+ if (nrows == 1) {
+ if (matrix[1,1] .gt. 0.)
+ matfac[1,1] = 1./matrix[1,1]
+ return
+ }
+
+
+ # copy matrix into matfac
+ do n = 1, nrows {
+ do j = 1, nbands
+ matfac[j,n] = matrix[j,n]
+ }
+
+ do n = 1, nrows {
+
+ # test to see if matrix is singular
+ if (matfac[1,n] + matrix[1,n] <= matrix[1,n]) {
+ do j = 1, nbands
+ w[j,n] = 0.
+ ier = 1
+ next
+ }
+
+ matfac[1,n] = 1./matfac[1,n]
+ imax = min (nbands - 1, nrows - n)
+ if (imax < 1)
+ next
+
+ jmax = imax
+ do i = 1, imax {
+ ratio = matfac[i+1,n] * matfac[1,n]
+ do j = 1, jmax
+ matfac[j,n+i] = matfac[j,n+i] - matfac[j+i,n] * ratio
+ jmax = jmax - 1
+ matfac[i+1,n] = ratio
+ }
+ }
+end
+
+# CHOSLV -- Solve the matrix whose Cholesky factorization was calculated in
+# CHOFAC.
+
+procedure choslv (matfac, nbands, nrows, vector, coeff)
+
+real matfac[nbands,nrows]
+int nbands
+int nrows
+real vector[nrows]
+real coeff[nrows]
+
+begin
+ if (nrows == 1) {
+ coeff[1] = vector[1] * matfac[1,1]
+ return
+ }
+
+ # copy vector to coefficients
+ do i = 1, nrows
+ coeff[i] = vector[i]
+
+ # forward substitution
+ nbndm1 = nbands - 1
+ do n = 1, nrows {
+ jmax = min (nbndm1, nrows - n)
+ if (jmax < 1)
+ next
+ do j = 1, jmax
+ coeff[j+n] = coeff[j+n] - matfac[j+1,n] * b[n]
+ }
+
+ # back substitution
+ for (n = nrows; n > 0; n = n - 1) {
+ coeff[n] = coeff[n] * matfac[1,n]
+ jmax = min (nbndm1, nrows - 1)
+ if (jmax >= 1) {
+ do j = 1, jmax
+ coeff[n] = coeff[n] - matfac[j+1,n] * coeff[j+n]
+ }
+ }
+
+end
generated by cgit (git 2.34.1) at 2025-09-20 08:41:08 -0400