diff options
Diffstat (limited to 'math/surfit')
41 files changed, 3663 insertions, 0 deletions
diff --git a/math/surfit/doc/iscoeff.hlp b/math/surfit/doc/iscoeff.hlp new file mode 100644 index 00000000..ea0b292c --- /dev/null +++ b/math/surfit/doc/iscoeff.hlp @@ -0,0 +1,63 @@ +.help iscoeff Apr85 "Surfit Package" +.ih +NAME +iscoeff - get the number and values of the surface coefficients +.ih +SYNOPSIS +iscoeff (sf, coeff, ncoeff) + +.nf +pointer sf # surface descriptor +real coeff[ncoeff] # coefficient array +int ncoeff # number of coefficients +.fi +.ih +ARGUMENTS +.ls sf +Pointer to the surface descriptor. +.le +.ls coeff +Array of coefficients. +.le +.ls ncoeff +The number of coefficients. +.le + +.ih +DESCRIPTION +ISCOEFF fetches the coefficient array and the number of coefficients from +the surface descriptor structure. A 1-D array ncoeff elements long +is required to hold the coefficients where ncoeff is defined as follows: + +.nf + SF_LEGENDRE: nxcoeff = xorder + nycoeff = yorder + ncoeff = xorder * yorder xterms = yes + ncoeff = nycoeff + nxcoeff - 1 xterms = no + + SF_CHEBYSHEV: nxcoeff = xorder + nycoeff = yorder + ncoeff = xorder * yorder xterms = yes + ncoeff = nycoeff + nxcoeff - 1 xterms = no + + SF_SPLINE3: nxcoeff = xorder + 3 + nycoeff = yorder + 3 + ncoeff = (xorder + 3) * (yorder + 3) + + SF_SPLINE1: nxcoeff = xorder + 1 + nycoeff = yorder + 1 + ncoeff = (xorder + 1) * (yorder + 1) + +.fi + +The coefficient of basis function B(i,x) * B(j,y) will be stored in +element (i - 1) * +nycoeff + j of the array coeff if the xterms parameter was set +to yes by ISINIT. Otherwise the nycoeff y-term coefficients will be output +first followed by the (nxcoeff - 1) x-term coefficients. + +.ih +NOTES +.ih +SEE ALSO +.endhelp diff --git a/math/surfit/doc/iseval.hlp b/math/surfit/doc/iseval.hlp new file mode 100644 index 00000000..53853cfe --- /dev/null +++ b/math/surfit/doc/iseval.hlp @@ -0,0 +1,34 @@ +.help iseval Apr85 "Surfit Package" +.ih +NAME +iseval -- evaluate the fitted surface at x and y +.ih +SYNOPSIS +y = iseval (sf, x, y) + +.nf +pointer sf # surface descriptor +real x # x value, 1 <= x <= ncols +real y # y value, 1 <= y <= nlines +.fi +.ih +ARGUMENTS +.ls sf +Pointer to the surface descriptor structure. +.le +.ls x, y +X and y values at which the surface is to be evaluated. +.le +.ih +DESCRIPTION +Evaluate the surface at the specified value of x and y. ISEVAL is a real +valued function which returns the fitted value. +.ih +NOTES +ISEVAL uses the coefficient array stored in the surface descriptor structure. +Checking for out of bounds x and y values is the responsibility of the calling +program. +.ih +SEE ALSO +isvector +.endhelp diff --git a/math/surfit/doc/isfree.hlp b/math/surfit/doc/isfree.hlp new file mode 100644 index 00000000..04f52b5f --- /dev/null +++ b/math/surfit/doc/isfree.hlp @@ -0,0 +1,27 @@ +.help isfree Apr85 "Surfit Package" +.ih +NAME +isfree -- free the surface grid descriptor structure +.ih +SYNOPSIS +isfree (sf) + +.nf +pointer sf # surface descriptor +.fi +.ih +ARGUMENTS +.ls sf +Pointer to the surface descriptor structure. +.le +.ih +DESCRIPTION +Frees the surface descriptor structure and all the vectors and arrays used +by the numerical routines. +.ih +NOTES +ISFREE should be called after the surface fit is complete. +.ih +SEE ALSO +isinit +.endhelp diff --git a/math/surfit/doc/isinit.hlp b/math/surfit/doc/isinit.hlp new file mode 100644 index 00000000..e0f2123c --- /dev/null +++ b/math/surfit/doc/isinit.hlp @@ -0,0 +1,61 @@ +.help isinit Apr85 "Surfit Package" +.ih +NAME +isinit -- initialize surface grid descriptor +.ih +SYNOPSIS +include <math/surfit.h> + +.nf +isinit (sf, surface_type, xorder, yorder, xterms, ncols, nlines) +.fi + +.nf +pointer sf # surface descriptor +int surface_type # surface function +int xorder # order of function in x +int yorder # order of function in y +int xterms # include cross-terms? (YES/NO) +int ncols # number of columns in the surface grid +int nlines # number of lines in the surface grid +.fi +.ih +ARGUMENTS +.ls sf +Pointer to the surface descriptor structure. +.le +.ls surface_type +Fitting function. Permitted values are SF_LEGENDRE and SF_CHEBYSHEV for +the Legendre and Chebyshev polynomials and SF_SPLINE1 and SF_SPLINE3 +for the linear and bicubic splines. +.le +.ls xorder, yorder +Order of the polynomial to be fit in x and y or the number of spline pieces to +be fit in x and y. The orders must be greater than or equal to 1. +.le +.ls xterms +Include cross-terms? If xterms = YES coefficients are fit to terms containing +the cross products of x and y polynomials. Xterms defaults to YES for the +spline functions. +.le +.ls ncols +The number of columns in the surface grid. The surface is assumed to lie +on a rectangular grid such that 1 <= x <= ncols. +.le +.ls nlines +The number of lines in the surface to be grid. The surface is assumed to lie +on a rectangular grid such that 1 <= y <= nlines. +.le +.ih +DESCRIPTION +ISINIT allocates space for the surface descriptor and the arrays and vectors +used by the numerical routines. It initializes all arrays and vectors to zero, +calculates and stores the basis functions in x and y +and returns the surface descriptor to the calling routine. +.ih +NOTES +ISINIT must be the first SURFIT routine called. +.ih +SEE ALSO +isfree +.endhelp diff --git a/math/surfit/doc/islaccum.hlp b/math/surfit/doc/islaccum.hlp new file mode 100644 index 00000000..48481cf0 --- /dev/null +++ b/math/surfit/doc/islaccum.hlp @@ -0,0 +1,60 @@ +.help islaccum Apr85 "Surfit Package" +.ih +NAME +islaccum -- accumulate a surface grid line into the fit +.ih +SYNOPSIS +include <math/surfit.h> + +islaccum (sf, cols, lineno, z, w, npts, wtflag) + +.nf +pointer sf # surface grid descriptor +int cols[npts] # column values, 1 <= col[i] <= ncols +int lineno # number of line, 1 <= lineno <= nlines +real z[npts] # surface values on lineno at cols +real w[npts] # array of weights +int npts # number of surface values, npts <= ncols +int wtflag # type of weighting desired +.fi +.ih +ARGUMENTS +.ls sf +Pointer to the surface grid descriptor structure. +.le +.ls cols +The column numbers of surface grid points on line lineno to be added to the +dataset. +.le +.ls lineno +The line number of the surface grid line to be added to the data set. +.le +.ls z +The data values. +.le +.ls w +The array of weights for the data points. +.le +.ls npts +The number of surface grid values on line number lineno to be added to the +data set. +.le +.ls wtflag +Type of weighting. The options are SF_USER and SF_UNIFORM. If wtflag +equals SF_USER the weight for each data point is supplied by the user +and points with zero-valued weights are not included in the fit. +If wtflag equals SF_UNIFORM the routine sets the weights to 1. +.le +.ih +DESCRIPTION +ISLACCUM computes the contribution of each data point to the normal +equations in x, sums that contribution into the appropriate arrays and +vectors and stores these intermediate results for use by ISLSOLVE. +.ih +NOTES +Checking for out of bounds col values and INDEF valued data points is +the responsibility of the calling program. +.ih +SEE ALSO +isinit, islfit, islrefit, islsolve, islzero, issolve, isfree +.endhelp diff --git a/math/surfit/doc/islfit.hlp b/math/surfit/doc/islfit.hlp new file mode 100644 index 00000000..5adf0a59 --- /dev/null +++ b/math/surfit/doc/islfit.hlp @@ -0,0 +1,67 @@ +.help islfit Apr85 "Surfit Package" +.ih +NAME +islfit -- fit a surface grid line +.ih +SYNOPSIS +include <math/surfit.h> + +islfit (sf, cols, lineno, z, w, npts, wtflag, ier) + +.nf +pointer sf # surface grid descriptor +real cols[npts] # array of column numbers, 1 <= cols[i] <= ncols +int lineno # number of surface grid line to be added +real z[npts] # data values +real w[npts] # weight array +int npts # number of data points, npts <= ncols +int wtflag # type of weighting +int ier # error code +.fi +.ih +ARGUMENTS +.ls sf +Pointer to the surface grid descriptor structure. +.le +.ls cols +The column numbers of surface grid points to be added to the dataset. +.le +.ls lineno +The line number of the surface grid line to be added to the dataset. +.le +.ls z +Array of data values. +.le +.ls w +Array of weights. +.le +.ls npts +Number of data points. +.le +.ls wtflag +Type of weighting. The options are SF_USER and SF_UNIFORM. If wtflag = +SF_USER individual weights for each data point are supplied by the calling +program and points with zero-valued weights are not included in the fit. +If wtflag = SF_UNIFORM, all weights are assigned values of 1. +.le +.ls ier +Error code for the fit. The options are OK, SINGULAR and NO_DEG_FREEDOM. +If ier = SINGULAR, the numerical routines will compute a solution but one +or more of the coefficients will be zero. If ier = NO_DEG_FREEDOM there +were too few data points to solve the matrix equations and the routine +returns without fitting the data. +.le +.ih +DESCRIPTION +ISLFIT zeroes the appropriate arrays and vectors, +computes the contribution of each data point to the normal equations +in x and accumulates it into the appropriate array and vector elements. The +x coefficients are stored for later use by ISSOLVE. +.ih +NOTES +Checking for out of bounds col values and INDEF values pixels is the +responsibility of the user. +.ih +SEE ALSO +isinit, islrefit, islaccum, islsolve, islzero, issolve, isfree +.endhelp diff --git a/math/surfit/doc/islrefit.hlp b/math/surfit/doc/islrefit.hlp new file mode 100644 index 00000000..c34476d6 --- /dev/null +++ b/math/surfit/doc/islrefit.hlp @@ -0,0 +1,51 @@ +.help islrefit Aug85 "Gsurfit Package" +.ih +NAME +islrefit -- refit surface grid line assuming old cols, w vector +.ih +SYNOPSIS +include <math/surfit.h> + +islrefit (sf, cols, lineno, z, w) + +.nf +pointer sf # surface grid descriptor +int cols[ARB] # columns to be fit, 1 <= cols[i] <= ncols +int lineno # line number of surface grid line to be fit +real z[ARB] # array of data values +real w[ARB] # array of weights +.fi +.ih +ARGUMENTS +.ls sf +Pointer to the surface grid descriptor structure. +.le +.ls cols +The column numbers of the surface grid line to be added to the dataset. +.le +.ls lineno +The number of the surface grid line to be added to the dataset. +.le +.ls z +Array of data values. +.le +.ls w +Array of weights. +.le +.ih +DESCRIPTION +In some applications the cols, and w values remain unchanged from fit +to fit and only the z values vary. In this case it is redundant to +reaccumulate the matrix and perform the Cholesky factorization for each +succeeding surface grid line. ISLREFIT +zeros and reaccumulates the vector on the right hand side of the matrix +equation and performs the forward and back substitution phase to fit for +a new coefficient vector. +.ih +NOTES +It is the responsibility of the calling program to check for out of bounds +column values and INDEF valued pixels. +.ih +SEE ALSO +isinit, islfit, islaccum, islsolve, islzero, issolve, isfree +.endhelp diff --git a/math/surfit/doc/islsolve.hlp b/math/surfit/doc/islsolve.hlp new file mode 100644 index 00000000..752cbc19 --- /dev/null +++ b/math/surfit/doc/islsolve.hlp @@ -0,0 +1,45 @@ +.help islsolve Apr85 "Surfit Package" +.ih +NAME +islsolve -- solve the equations for a surface grid line +.ih +SYNOPSIS +include <math/surfit.h> + +islsolve (sf, lineno, ier) + +.nf +pointer sf # surface grid descriptor +int lineno # surface grid line number +int ier # error code +.fi +.ih +ARGUMENTS +.ls sf +Pointer to the surface grid descriptor structure. +.le +.ls lineno +The number of the surface grid line to be added to the dataset. +.le +.ls ier +Error code returned by the fitting routines. The options are OK, SINGULAR, +and NO_DEG_FREEDOM. If ier = SINGULAR the matrix is singular, ISLSOLVE +will compute a solution to the normal equations but one or more of the +coefficients will be zero. If ier = NO_DEG_FREEDOM, too few data points +exist for a reasonable solution to be computed and ISLSOLVE returns without +fitting the data. +.le +.ih +DESCRIPTION +ISLSOLVE computes the Cholesky factorization of the data matrix and +solves for the coefficients +of the x fitting function by forward and back substitution. +An error code is +returned by ISLSOLVE if it is unable to solve the normal equations as +formulated. +.ih +NOTES +.ih +SEE ALSO +isinit, islfit, islrefit, islaccum, islzero, issolve, isfree +.endhelp diff --git a/math/surfit/doc/islzero.hlp b/math/surfit/doc/islzero.hlp new file mode 100644 index 00000000..cbdb9515 --- /dev/null +++ b/math/surfit/doc/islzero.hlp @@ -0,0 +1,32 @@ +.help islzero Apr85 "Surfit Package" +.ih +NAME +islzero -- set up for a new surface grid line fit +.ih +SYNOPSIS +islzero (sf, lineno) + +.nf +pointer sf # surface grid descriptor +int lineno # surface grid line number +.fi +.ih +ARGUMENTS +.ls sf +Pointer to the surface grid descriptor structure. +.le +.ls lineno +Line number of the surface grid line fit to be zeroed. +.le +.ih +DESCRIPTION +ISLZERO zeros the appropriate arrays and vectors for the surface grid line +number line number. The line can then be refit by calls to ISLACCUM and +ISLSOLVE. Alternatively a single call to ISLFIT can perform the combined +functions of ISLZERO, ISLACCUM and ISLSOLVE. +.ih +NOTES +.ih +SEE ALSO +isinit, islfit, islrefit, islaccum, islsolve, issolve +.endhelp diff --git a/math/surfit/doc/isreplace.hlp b/math/surfit/doc/isreplace.hlp new file mode 100644 index 00000000..d248e5f8 --- /dev/null +++ b/math/surfit/doc/isreplace.hlp @@ -0,0 +1,33 @@ +.help isreplace Apr1985 "Surfit Package" +.ih +NAME +isreplace -- restore a surface fit saved by issave +.ih +SYNOPSIS +isreplace (sf, fit) + +.nf +pointer sf # surface grid descriptor +real fit[ARB] # fit array +.fi + +.ih +ARGUMENTS + +.ls sf +The pointer to the surface grid descriptor. +.le +.ls fit +The array containing the fit parameters and coefficients. +.le + +.ih +DESCRIPTION + +ISREPLACE restores the fit saved by ISSAVE for use by ISEVAL or ISVECTOR. +.ih +NOTES +.ih +SEE ALSO +issave +.endhelp diff --git a/math/surfit/doc/isresolve.hlp b/math/surfit/doc/isresolve.hlp new file mode 100644 index 00000000..afed8ed8 --- /dev/null +++ b/math/surfit/doc/isresolve.hlp @@ -0,0 +1,45 @@ +.help isresolve Aug85 "Surfit Package" +.ih +NAME +isresolve -- resolve the surface, same lines array +.ih +SYNOPSIS +include <math/surfit.h> + +isresolve (sf, lines, ier) + +.nf +pointer sf # surface grid descriptor +int lines[nlines] # surface grid line numbers, 1 <= lines[i] <= nlines +int ier # error code +.fi +.ih +ARGUMENTS +.ls sf +Pointer to the surface grid descriptor structure. +.le +.ls nlines +The number of surface grid lines to be fit. +.le +.ls ier +Error code returned by the fitting routines. The options are OK, SINGULAR, +and NO_DEG_FREEDOM. If ier = SINGULAR the matrix is singular, ISSOLVE +will compute a solution to the normal equations but one or more of the +coefficients will be zero. If ier = NO_DEG_FREEDOM, too few data points +exist for a reasonable solution to be computed. ISSOLVE returns without +fitting the data. +.le +.ih +DESCRIPTION +In some applications it is necessary to refit the same surface several +times with the same lines array. In this case it is inefficient to reaccumulate +the matrices and calculate the Cholesky factorization for each fit. +ISERESOLVE reaccumulates only the right side of the equation. +.ih +NOTES +It is the responsibility of the calling program to check for out of bounds +lines values. +.ih +SEE ALSO +issolve, isinit, slzero, islfit, islrefit, islaccum, islresolve, isfree +.endhelp diff --git a/math/surfit/doc/issave.hlp b/math/surfit/doc/issave.hlp new file mode 100644 index 00000000..6ec97cd6 --- /dev/null +++ b/math/surfit/doc/issave.hlp @@ -0,0 +1,55 @@ +.help issave Apr85 "Surfit Package" +.ih +NAME +issave - save the surface fit +.ih +SYNOPSIS +issave (sf, fit) + +.nf +pointer sf # surface descriptor +real fit[ARB] # fit array +.fi +.ih +ARGUMENTS +.ls sf +Pointer to the surface grid descriptor. +.le +.ls fit +Array for storing the fit. The fit array must be at least SAVE_COEFF + ncoeff +elements long. +.le + +.ih +DESCRIPTION +ISSAVE stores the surface parameters and coefficient array +in the 1-D reall array fit. Fit must be at least ncoeff + SAVE_COEFF +elements long where ncoeff is defined as follows. + +.nf + SF_LEGENDRE: nxcoeff = xorder + nycoeff = yorder + ncoeff = xorder * yorder xterms = yes + ncoeff = nycoeff + nxcoeff - 1 xterms = no + + SF_CHEBYSHEV: nxcoeff = xorder + nycoeff = yorder + ncoeff = xorder * yorder xterms = yes + ncoeff = nycoeff + nxcoeff - 1 xterms = no + + SF_SPLINE3: nxcoeff = xorder + 3 + nycoeff = yorder + 3 + ncoeff = (xorder + 3) * (yorder + 3) + + SF_SPLINE1: nxcoeff = xorder + 1 + nycoeff = yorder + 1 + ncoeff = (xorder + 1) * (yorder + 1) + +.fi + +.ih +NOTES +.ih +SEE ALSO +isreplace +.endhelp diff --git a/math/surfit/doc/issolve.hlp b/math/surfit/doc/issolve.hlp new file mode 100644 index 00000000..312d69a3 --- /dev/null +++ b/math/surfit/doc/issolve.hlp @@ -0,0 +1,49 @@ +.help issolve Aug85 "Surfit Package" +.ih +NAME +issolve -- solve the surface +.ih +SYNOPSIS +include <math/surfit.h> + +issolve (sf, lines, nlines, ier) + +.nf +pointer sf # surface grid descriptor +int lines[nlines] # surface grid line numbers, 1 <= lines[i] <= nlines +int nlines # number of lines +int ier # error code +.fi +.ih +ARGUMENTS +.ls sf +Pointer to the surface grid descriptor structure. +.le +.ls lines +The line number of surface grid lines to be included in the fit. +.le +.ls nlines +The number of surface grid lines to be fit. +.le +.ls ier +Error code returned by the fitting routines. The options are OK, SINGULAR, +and NO_DEG_FREEDOM. If ier = SINGULAR the matrix is singular, ISSOLVE +will compute a solution to the normal equations but one or more of the +coefficients will be zero. If ier = NO_DEG_FREEDOM, too few data points +exist for a reasonable solution to be computed. ISSOLVE returns without +fitting the data. +.le +.ih +DESCRIPTION +ISSOLVE solves for the surface coefficients. +An error code is +returned by ISSOLVE if it is unable to solve the normal equations as +formulated. +.ih +NOTES +It is the responsibility of the calling program to check for out of bounds +lines values. +.ih +SEE ALSO +isinit, slzero, islfit, islrefit, islaccum, islresolve, isfree +.endhelp diff --git a/math/surfit/doc/isvector.hlp b/math/surfit/doc/isvector.hlp new file mode 100644 index 00000000..fabbdcc7 --- /dev/null +++ b/math/surfit/doc/isvector.hlp @@ -0,0 +1,41 @@ +.help isvector Aug85 "Surfit Package" +.ih +NAME +isvector -- evaluate the fitted surface at a set of points +.ih +SYNOPSIS +isvector (sf, x, y, zfit, npts) + +.nf +pointer sf # surface grid descriptor +real x[npts] # x array, 1 <= x[i] <= ncols +real y[npts] # y array, 1 <= y[i] <= nlines +real zfit[npts] # data values +int npts # number of data points +.fi +.ih +ARGUMENTS +.ls sf +Pointer to the surface grid descriptor structure. +.le +.ls x, y +Arrays of x and y values. +.le +.ls zfit +Array of fitted values. +.le +.ls npts +The number of points to be fit. +.le +.ih +DESCRIPTION +Fit the surface to an array of x and y values. ISVECTOR uses the coefficients +stored in the surface descriptor structure. +.ih +NOTES +Checking for out of bounds x and y values is the responsibility of the +calling program. +.ih +SEE ALSO +iseval +.endhelp diff --git a/math/surfit/doc/iszero.hlp b/math/surfit/doc/iszero.hlp new file mode 100644 index 00000000..4c338807 --- /dev/null +++ b/math/surfit/doc/iszero.hlp @@ -0,0 +1,26 @@ +.help iszero Aug85 "Surfit Package" +.ih +NAME +iszero -- set up for a new surface fit +.ih +SYNOPSIS +iszero (sf) + +.nf +pointer sf # surface grid descriptor +.fi +.ih +ARGUMENTS +.ls sf +Pointer to the surface grid descriptor structure. +.le +.ih +DESCRIPTION +ISZERO zeros all matrices and vectors used by the fitting program. It +should be used only if one wishes to refit the surface from scratch. +.ih +NOTES +.ih +SEE ALSO +isinit, islzero, islaccum, islfit, islrefit, islsolve, issolve, isfree +.endhelp diff --git a/math/surfit/doc/surfit.hd b/math/surfit/doc/surfit.hd new file mode 100644 index 00000000..18e7ec1a --- /dev/null +++ b/math/surfit/doc/surfit.hd @@ -0,0 +1,19 @@ +# HELP Directory for SURFIT package + +$surfit = "math$surfit/" + +iscoeff hlp = iscoeff.hlp, src = surfit$iscoeff.x +iseval hlp = iseval.hlp, src = surfit$iseval.x +isfree hlp = isfree.hlp, src = surfit$isfree.x +isinit hlp = isinit.hlp, src = surfit$isinit.x +islaccum hlp = islaccum.hlp, src = surfit$islaccum.x +islfit hlp = islfit.hlp, src = surfit$islfit.x +islrefit hlp = islrefit.hlp, src = surfit$islrefit.x +islsolve hlp = islsolve.hlp, src = surfit$islsolve.x +islzero hlp = islzero.hlp, src = surfit$islzero.x +isreplace hlp = isreplace.hlp, src = surfit$isreplace.x +isresolve hlp = isresolve.hlp, src = surfit$isresolve.x +issave hlp = issave.hlp, src = surfit$issave.x +issolve hlp = issolve.hlp, src = surfit$issolve.x +isvector hlp = isvector.hlp, src = surfit$isvector.x +iszero hlp = iszero.hlp, src = surfit$iszero.x diff --git a/math/surfit/doc/surfit.hlp b/math/surfit/doc/surfit.hlp new file mode 100644 index 00000000..a7c347cc --- /dev/null +++ b/math/surfit/doc/surfit.hlp @@ -0,0 +1,157 @@ +.help surfit Aug84 "Math Package" + +.ih +NAME +surfit -- set of routines for fitting gridded surface data + + +.ih +SYNOPSIS + +The surface is defined in the following way. The x and y values are +assumed to be integers and lie in the range 0 <= x <= ncols + 1 and +0 <= y <= nlines + 1. The package prefix is for image surface fit. +Package routines with a prefix followed by l operate on individual +image lines only. + +.nf + isinit (sf, surf_type, xorder, yorder, xterms, ncols, nlines) + iszero (sf) + islzero (sf, lineno) + islaccum (sf, cols, lineno, z, w, ncols, wtflag) + islsolve (sf, lineno, ier) + islfit (sf, cols, lineno, z, w, ncols, wtflag, ier) + islrefit (sf, cols, lineno, z, w, ier) + issolve (sf, lines, nlines, ier) + isresolve (sf, lines, ier) + z = iseval (sf, x, y) + isvector (sf, x, y, zfit, npts) + issave (sf, surface) + isrestore (sf, surface) + isfree (sf) +.fi + +.ih +DESCRIPTION + +The SURFIT package provides a set of routines for fitting surfaces to functions +linear in their coefficients using the tensor product method and least squares +techniques. The basic numerical +technique employed is the solution of the normal equations by the Cholesky +method. + +.ih +NOTES + +The following series of steps illustrates the use of the package. + +.ls 4 +.ls (1) +Include an include <math/surfit.h> statement in the calling program to make the +SURFIT package definitions available to the user program. +.le +.ls (2) +Call ISINIT to initialize the surface fitting parameters. +.le +.ls (3) +Call ISLACCUM to select a weighting function and accumulate data points +for each line into the appropriate arrays and vectors. Call ISLSOLVE +to compute the coefficients in x for each line. The coefficients are +stored inside SURFIT. Repeat this procedure for each line. ISLACCUM +and SFLSOLVE can be combined by a call to SFLFIT. If the x values +and weights remain the same from line to line ISLREFIT can be called. +.le +.ls (4) +Call ISSOLVE to solve for the surface coefficients. +.le +.ls (5) +Call ISEVAL or ISVECTOR to evaluate the fitted surface at the x and y +value(s) of interest. +.le +.ls (6) +Call ISFREE to release the space allocated for the fit. +.le +.le + +.ih +EXAMPLES + +.nf +Example 1: Fit a 2nd order Lengendre polynomial in x and y to an image +and output the fitted image + +include <imhdr.h> +include <math/surfit.h> + + + old = immap (old_image, READ_ONLY, 0) + new = immap (new_image, NEW_COPY, 0) + + ncols = IM_LEN(old, 1) + nlines = IM_LEN(old, 2) + + # initialize surface fit + call isinit (sf, LEGENDRE, 3, 3, YES, ncols, nlines) + + # allocate space for lines, columns and weight arrays + call malloc (cols, ncols, TY_INT) + call malloc (lines, nlines, TY_INT) + call malloc (weight, ncols, TY_REAL) + + # initialize lines and columns arrays + do i = 1, ncols + Memi[cols - 1 + i] = i + do i = 1, nlines + Memi[lines - 1 + i] = i + + # fit the surface in x line by line + call amovkl (long(1), v, IM_MAXDIM) + do i = 1, nlines { + if (imgnlr (old, inbuf, v) == EOF) + call error (0, "Error reading image.") + if (i == 1) { + call islfit (sf, Memi[cols], i, Memr[inbuf], Memr[weight], + ncols, SF_WTSUNIFORM, ier) + if (ier != OK) + ... + } else + call islrefit (sf, Memi[cols], i, Memr[inbuf], Memr[weight], + ier) + } + + # solve for surface coefficients + call issolve (sf, Memi[lines], nlines, ier) + + # free space used in fitting arrays + call mfree (cols, TY_INT) + call mfree (lines, TY_INT) + call mfree (weight, TY_REAL) + + # allocate space for x and y arrays + call malloc (x, ncols, TY_REAL) + call malloc (y, ncols, TY_REAL) + + # intialize z array + do i = 1, ncols + Memr[x - 1 + i] = real (i) + + # create fitted image + call amovkl (long(10, v, IM_MAXDIM) + do i = 1, nlines { + if (impnlr (new, outbuf, v) == EOF) + call error (0, "Error writing image.") + call amovkr (real (i), Memr[y], ncols) + call isvector (sf, Memr[x], Memr[y], Memr[outbuf], ncols) + } + + # close files and cleanup + call mfree (x, TY_REAL) + call mfree (y, TY_REAL + + call isfree (sf) + + call imunmap (old) + call imunmap (new) + +.fi +.endhelp diff --git a/math/surfit/doc/surfit.men b/math/surfit/doc/surfit.men new file mode 100644 index 00000000..4a3b3258 --- /dev/null +++ b/math/surfit/doc/surfit.men @@ -0,0 +1,15 @@ + isinit - initialize surface descriptor + iscoeff - get coefficients of surface + iseval - evaluate the surface at a point + isfree - free surface descriptor + islaccum - accumulate surface at a given line number + islfit - fit the surface at a given line number + islrefit - refit the surface at a given line number, same columns + islsolve - solve the surface at a given line number + islzero - zero the fit for a given line number + issolve - solve the surface + isreplace - restore the surface fit + isresolve - resolve surface, same lines + issave - save the surface fit + isvector - fit surface at a set of (x,y) values + iszero - zero the surface fit diff --git a/math/surfit/doc/surfit.spc b/math/surfit/doc/surfit.spc new file mode 100644 index 00000000..07f2b934 --- /dev/null +++ b/math/surfit/doc/surfit.spc @@ -0,0 +1,500 @@ +.help surfit Aug84 "Math Package" +.CE +Specifications for the Surface Fitting Package +.CE +Lindsey Davis +.CE +August 1984 +.CE +First Draft + +.sh +1. Introduction + +The SURFIT package provides a set of routines for fitting surfaces to functions +linear in their coefficients using least squares techniques. The basic numerical +technique employed is the solution of the normal equations by the Cholesky +method. + +.sh +2. Requirements + +.ls (1) +The package shall take as input a surface or section of +a surface. +The data are assumed to be on a rectangular grid and to have the following +ranges, 0 <= x <= ncols + 1 and 0 <= y <= nlines + 1. +The package routines assume that data values equal to INDEF have +been removed from the data set or replaced with appropriate interpolated +values prior to entering the package routines. It is not necessary that +all the data be present to perform the fit. +.le +.ls (2) +The package shall perform the following operations: +.ls o +Determine the coefficients of the fitting function by solving the normal +equations. The surface fitting function is selected at run time from the +following list: (1) SF_LEGENDRE, Legendre polynomials in x and y, +(2) SF_CHEBYSHEV, Chebyshev polynomials in x and y, (3) SF_SPLINE3, bicubic +spline with break points evenly spaced in x and y. The calling sequence +must be invariant to the form of the fitting function. +.le +.ls o +Set an error code if the numerical routines are unable to fit the specified +function. +.le +.ls o +Optionally output the values of the coefficients. The coefficients will be +stored internal to the SURFIT package. However in some applications it is +the coefficients which are of primary interest. A package routine shall +exist to extract the coefficients from the curve descriptor structure. +.le +.ls o +Evaluate the surface at arbitrary value(s) of x and y. The evaluating +routines shall use the calculated coefficients and the user supplied +x values(s). +.le +.ls o +Be capable of storing the fitted surface parameters and coefficients in a +user supplied array and restoring the saved fit for later use by the +evaluating routines. +.le +.ls o +Input surface parameters and surface coefficients derived external to +the surfit package into the surfit format for use by the evaluate +routines. +.le +.le +.ls (3) +The program shall perform a weighted fit to the surface using a user +supplied weight array and weight flag. The weighting options will be +SF_WTSUSER and SF_WTSUNIFORM. In SF_WTSUNIFORM mode the package routines +set all the weights to 1. In SF_WTSUSER mode the package routines apply +user supplied weights to the individual data points. +.le +.ls (4) +The input data set, output coefficients, error and fitted z arrays are +single precision real quantities. All package arithmetic shall be performed +in single precision. +.le + +.sh +3. Specifications + +.sh +3.1. List of Routines + + +The surface is defined in the following way. The x and y values are +assumed to be integers and lie in the range 0 <= x <= ncols + 1 and +0 <= y <= nlines + 1. The package prefix is is for image surface fit. +Package routines with a prefix followed by l operate on individual +image lines only. + +.nf + z = f (col, line) +.fi + +.nf + isinit (sf, surf_type, xorder, yorder, xterms, ncols, nlines) + iszero (sf) + islzero (sf, lineno) + islaccum (sf, cols, lineno, z, w, ncols, wtflag) + islsolve (sf, lineno, ier) + islfit (sf, cols, lineno, z, w, ncols, wtflag, ier) + islrefit (sf, cols, lineno, z, w, ier) + issolve (sf, lines, nlines, ier) + isresolve (sf, lines, ier) + z = iseval (sf, x, y) + isvector (sf, x, y, zfit, npts) + issave (sf, surface) + isrestore (sf, surface) + isfree (sf) +.fi + +.sh +3.2. Algorithms + +.sh +3.2.1. Polynomial Basis Functions + +The approximating function is assumed to be of the following form, + +.nf + f(x,y) = sum (a[i,j] * F(i,x) * F(j,y)) i = 1...nxcoeff j = 1...nycoeff +.fi + +where the F(i,x) are the polynomial basis functions containing terms of order +x**(i-1) and the a(i,j) are the coefficients. In order to avoid a very +ill-conditioned linear system for moderate or large i, the Legendre and +Chebyshev polynomials were chosen for the basis functions. The Chebyshev +and Legendre polynomials are orthogonal over -1. <= x,y <= 1. The data x and +y values are normalized to this regime using the number of lines and columns +supplied by the user. For each data point the basis +functions are calculated using the following recursion relations. The +cross terms are optional. + +Legendre series: +.nf + + F(1,x) = 1. + F(2,x) = x + F(i,x) = [(2*i-1)*x*F(i-1,x)-(i-2)*F(i-2,x)]/(i-1) + F(1,y) = 1. + F(2,y) = y + F(j,y) = [(2*j-1)*y*F(j-1,y)-(j-2)*F(j-2,y)]/(j-1) + +.fi +Chebyshev series: +.nf + + F(1,x) = 1. + F(2,x) = x + F(i,x) = 2.*x*F(i-1,x)-F(i-2,x) + F(1,y) = 1. + F(2,y) = y + F(j,y) = 2.*y*F(j-1,y)-F(j-2,y) +.fi + + +.sh +3.2.2. Bicubic Cardinal B-spline + +The approximating function is of the form + +.nf + f(x,y) = sum (x(i,j) * F(i,x) * F(j,y)) i=1...nxcoeff j=1...nycoeff +.fi + +where the basis functions, F(i,x), are the cubic cardinal B-splines +(Prenter 1975). The user supplies the number of columns and lines and the +number of polynomial pieces in x and y, nxpieces and nypieces, to be fit +to the data set. The number of bicubic spline coefficients, ncoeff, will be + +.nf + nxcoeff = (nxpieces + 3) + nycoeff = (nypieces + 3) + ncoeff = nxcoeff * nycoeff +.fi + +The cardinal B-spline is stored in a lookup table. For each x and y the +appropriate break point is selected and the four non-zero B-splines are +calculated by nearest neighbour interpolation in the lookup table. + +.sh +3.2.3. Method of Solution + +.sh +3.2.3.1. Full Surface Fit + +The normal equations are accumulated in the following form. + +.nf + c * x = b + c[i,j] = (B(i,x,y), B(j,x,y)) + b[j] = (B(j,x,y), S(x,y)) +.fi + +B(i,x,y) is the ith basis function at x and y, S(x,y) is the surface to +be approximated and the inner product of two functions G and H is +given by + +.nf + (G,H) = sum (w[i] * G(x[i],y[i]) * H(x[i],y[i])) i = 1...npts +.fi + +Since the matrix c is symmetric and positive semi-definite it may +be solved by the Cholesky method. +The coefficient matrix c can be written as + +.nf + c = l * d * l-transpose +.fi + +where l is a unit lower triangular matrix and d is the diagonal of c +(de Boor 1978). Near zero pivots are handled in the following way. +At the nth elimination step the current value of the nth diagonal +element is compared with the original nth diagonal element. If the +diagonal element has been reduced by one computer word length, the +entire nth row is declared linearly dependent on the previous n-1 +rows and x(n) = 0. + +The triangular system + +.nf + l * w = b +.fi + +is solved for w (forward substitution), the vector d ** (-1) * w is +computed and the triangular system + +.nf + l-transpose * x = d ** (-1) * w +.fi + +solved for the coefficients, x (backward substitution). + + +.sh +3.2.3.2. Line by Line Fit + +Each line of the surface can be represented by the following equation. + +.nf + S(x,yo) = a[i](yo) * B[i](x) i = 1...nxcoeff for each yo +.fi + +The normal equations for each image line are formed + +.nf + c * a = b + c[i,j] = (B(i,x), B(j,x)) + b[j] = (B(j,x), S(x,yo)) +.fi + +and solved for a as described in the previous section. +After fitting the entire image matrix a has nxcoeff columns and +nlines rows. + +For each column i in a the normal equations +are formed and solved for the c[i,j] + +.nf + c * x = b + c[j,k] = (B(j,y), B(k,x)) + b[i,j] = (a[i](y), B(j,y)) +.fi + +.sh +3.2.4. Number of Operations + +It is worth while to calculate the number of operations reqired to compute +the surface fit by the full surface method versus the line by line method. +The number of operations required to accumulate the normal equations +and solve the set is the following + +.nf + nop = npts * order ** 2 / 2 + order ** 3 /6 +.fi + +where the first term is the number of operations required to accumulate the +matrix and the second is the number required to solve the system. + +.nf +Example 1: full surface, 3 by 3 polynomial no cross terms, 512 by 512 image + + norder = 5 + npts = 512 * 512 + nops(accum) = 3,276,800 + nops(solve) = 21 + +Example 2: full surface, 30 by 30 spline, 512 by 512 image + + norder = 900 + npts = 512 * 512 + nops(accum) = 1.062 * 10 ** 11 + nops(solve) = 1.216 * 10 ** 8 + +Example 3: line by line method, 3 by 3 polynomial, 512 by 512 image + + step 1 solve for a coefficients + norder = 3 + npts = 512 + nlines = 512 + nops(accum) = 1,179,648 + nops(solve) = 2304 + + step 2 solve for c coefficients + norder = 3 + npts = 512 + nlines = 3 + nops(accum) = 6912 + nops(solve) = 14 + +Example 4: line by line method, 30 by 30 spline, 512 by 512 image + + step 1 solve coefficients + norder = 30 + npts = 512 + nlines = 512 + nops(accum) = 117,964,800 + nops(solve) = 2,304,000 + + step 2 solve for c coefficients + norder = 30 + npts = 512 + nlines = 30 + nops(accum) = 6,912,000 + nops(solve) = 135,000 +.fi + +The figures for the line by line method are worst case numbers. If the +x values remain the same from line to line then the coefficient matrix +only has to be accumulated and inverted twice. +For the bicubic spline function the number of operations is significantly +reduced by taking advantage of the banded nature of the matrices. + +.sh +4. Usage + +.sh +4.1. User Notes + +The following series of steps illustrates the use of the package. + +.ls 4 +.ls (1) +Include an include <surfit.h> statement in the calling program to make the +SURFIT package definitions available to the user program. +.le +.ls (2) +Call SFINIT to initialize the surface fitting parameters. +.le +.ls (3) +Call SFLACCUM to select a weighting function and accumulate data points +for each line into the appropriate arrays and vectors. Call SFLSOLVE +to compute the coefficients in x for each line. The coefficents are +stored inside SURFIT. Repeat this procedure for each line. SFLACCUM +and SFLSOLVE can be combined by a call to SFLFIT. If the x values +and weights remain the same from line to line SFLREFIT can be called. +.le +.ls (4) +Call SFSOLVE to solve for the surface coefficients. +.le +.ls (5) +Call SFEVAL or SFVECTOR to evaluate the fitted surface at the x and y +value(s) of interest. +.le +.ls (6) +Call SFFREE to release the space allocated for the fit. +.le +.le + +.sh +4.2. Examples + +.nf +Example 1: Fit a 2nd order Lengendre polynomial in x and y to an image +and output the fitted image + +include <imhdr.h> +include <math/surfit.h> + + + old = immap (old_image, READ_ONLY, 0) + new = immap (new_image, NEW_COPY, 0) + + ncols = IM_LEN(old, 1) + nlines = IM_LEN(old, 2) + + # initialize surface fit + call isinit (sf, LEGENDRE, 3, 3, YES, ncols, nlines) + + # allocate space for lines, columns and weight arrays + call malloc (cols, ncols, TY_INT) + call malloc (lines, nlines, TY_INT) + call malloc (weight, ncols, TY_REAL) + + # initialize lines and columns arrays + do i = 1, ncols + Memi[cols - 1 + i] = i + do i = 1, nlines + Memi[lines - 1 + i] = i + + # fit the surface in x line by line + call amovkl (long(1), v, IM_MAXDIM) + do i = 1, nlines { + if (imgnlr (old, inbuf, v) == EOF) + call error (0, "Error reading image.") + if (i == 1) { + call islfit (sf, Memi[cols], i, Memr[inbuf], Memr[weight], + ncols, SF_WTSUNIFORM, ier) + if (ier != OK) + ... + } else + call sflrefit (sf, Memi[cols], i, Memr[inbuf], Memr[weight], + ier) + } + + # solve for surface coefficients + call issolve (sf, Memi[lines], nlines, ier) + + # free space used in fitting arrays + call mfree (cols, TY_INT) + call mfree (lines, TY_INT) + call mfree (weight, TY_REAL) + + # allocate space for x and y arrays + call malloc (x, ncols, TY_REAL) + call malloc (y, ncols, TY_REAL) + + # intialize z array + do i = 1, ncols + Memr[x - 1 + i] = real (i) + + # create fitted image + call amovkl (long(10, v, IM_MAXDIM) + do i = 1, nlines { + if (impnlr (new, outbuf, v) == EOF) + call error (0, "Error writing image.") + call amovkr (real (i), Memr[y], ncols) + call isvector (sf, Memr[x], Memr[y], Memr[outbuf], ncols) + } + + # close files and cleanup + call mfree (x, TY_REAL) + call mfree (y, TY_REAL + + call isfree (sf) + + call imunmap (old) + call imunmap (new) + +.fi +.sh +5. Detailed Design + +.sh +5.1. Surface Descriptor Structure + +To be written when specifications are finalised. + +.sh +5.2. Storage Requirements + +The minimum storage requirements for the full surface fit method +assuming that points are to be rejected from the fit without revaluating +the matrix are the following. + +.nf + (nxcoeff*nycoeff) ** 2 -- coefficient array plus triangularized matrix + nxcoeff*nycoeff -- right side + nxcoeff*nycoeff -- solution vector +.fi + +For a 30 by 30 spline roughly 811,800 storage units are required. +For a 3 by 3 polynomial the requirements are 675 storage units. The +requirements are roughly half when rejection is disabled. + +The minimum storage requirements for the line by line method are + +.nf + nx*nx*2 -- The x matrix and its factorization + nx*nlines -- The x coefficient matrix + ny*ny*2 -- The y matrix and its factorization + nx*ny -- The solution +.fi + + +.sh +6. References + +.ls (1) +Carl de Boor, "A Practical Guide to Splines", 1978, Springer-Verlag New York +Inc. +.le +.ls (2) +P.M. Prenter, "Splines and Variational Methods", 1975, John Wiley and Sons +Inc. +.le +.endhelp diff --git a/math/surfit/iscoeff.x b/math/surfit/iscoeff.x new file mode 100644 index 00000000..f656e7e6 --- /dev/null +++ b/math/surfit/iscoeff.x @@ -0,0 +1,37 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "surfitdef.h" + +# ISCOEFF -- Procedure to fetch the number and magnitude of the coefficients. +# If the surface type is a bicubic spline or polynomials with cross terms +# then the number of coefficients is SF_NXCOEFF(sf) * SF_NYCOEFF(sf) and +# the coefficient of B(i,x) * B(j,y) will be stored in element +# (i - 1) * SF_NYCOEFF(sf) + j of the array coeff. Otherwise the number +# of coefficients will be (SF_NXCOEFF(sf) + SF_NYCOEFF(sf) - 1), and +# the SF_NYCOEFF(sf) y coefficients will be output first followed by +# the (SF_NXCOEFF(sf) - 1) x coefficients. + +procedure iscoeff (sf, coeff, ncoeff) + +pointer sf # pointer to the surface fitting descriptor +real coeff[ARB] # the coefficients of the fit +int ncoeff # the number of coefficients + +int i +pointer cptr + +begin + # calculate the number of coefficients + if (SF_XTERMS(sf) == NO) { + ncoeff = SF_NXCOEFF(sf) + SF_NYCOEFF(sf) - 1 + call amovr (COEFF(SF_COEFF(sf)), coeff, SF_NYCOEFF(sf)) + cptr = SF_COEFF(sf) + SF_NYCOEFF(sf) + do i = SF_NYCOEFF(sf) + 1, ncoeff { + coeff[i] = COEFF(cptr) + cptr = cptr + SF_NYCOEFF(sf) + } + } else { + ncoeff = SF_NXCOEFF(sf) * SF_NYCOEFF(sf) + call amovr (COEFF(SF_COEFF(sf)), coeff, ncoeff) + } +end diff --git a/math/surfit/iseval.x b/math/surfit/iseval.x new file mode 100644 index 00000000..08ef5f89 --- /dev/null +++ b/math/surfit/iseval.x @@ -0,0 +1,92 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <math/surfit.h> +include "surfitdef.h" + +# ISEVAL -- Procedure to evaluate the fitted surface at a single point. +# The SF_NXCOEFF(sf) by SF_NYCOEFF(sf) coefficients are stored in the +# SF_NYCOEFF(sf) by SF_NXCOEFF(sf) matrix COEFF. The j-th element of the ith +# row of COEFF contains the coefficient of the i-th basis function in x and +# the j-th basis function in y. + +real procedure iseval (sf, x, y) + +pointer sf # pointer to surface descriptor structure +real x # x value +real y # y value + +real sum, accum +int i, k, leftx, lefty, yorder +pointer sp, xb, xzb, yb, yzb, czptr + +begin + # allocate space for the basis functions + call smark (sp) + call salloc (xb, SF_XORDER(sf), MEM_TYPE) + xzb = xb - 1 + call salloc (yb, SF_YORDER(sf), MEM_TYPE) + yzb = yb - 1 + + # calculate the basis functions + switch (SF_TYPE(sf)) { + case SF_CHEBYSHEV: + leftx = 0 + lefty = 0 + czptr = SF_COEFF(sf) - 1 + call sf_b1cheb (x, SF_XORDER(sf), SF_XMAXMIN(sf), SF_XRANGE(sf), + XBS(xb)) + call sf_b1cheb (y, SF_YORDER(sf), SF_YMAXMIN(sf), SF_YRANGE(sf), + YBS(yb)) + + case SF_LEGENDRE: + leftx = 0 + lefty = 0 + czptr = SF_COEFF(sf) - 1 + call sf_b1leg (x, SF_XORDER(sf), SF_XMAXMIN(sf), SF_XRANGE(sf), + XBS(xb)) + call sf_b1leg (y, SF_YORDER(sf), SF_YMAXMIN(sf), SF_YRANGE(sf), + YBS(yb)) + + case SF_SPLINE3: + call sf_b1spline3 (x, SF_NXPIECES(sf), -SF_XMIN(sf), + SF_XSPACING(sf), XBS(xb), leftx) + call sf_b1spline3 (y, SF_NYPIECES(sf), -SF_YMIN(sf), + SF_YSPACING(sf), YBS(yb), lefty) + czptr = SF_COEFF(sf) - 1 + lefty + leftx * SF_NYCOEFF(sf) + + case SF_SPLINE1: + call sf_b1spline1 (x, SF_NXPIECES(sf), -SF_XMIN(sf), + SF_XSPACING(sf), XBS(xb), leftx) + call sf_b1spline1 (y, SF_NYPIECES(sf), -SF_YMIN(sf), + SF_YSPACING(sf), YBS(yb), lefty) + czptr = SF_COEFF(sf) - 1 + lefty + leftx * SF_NYCOEFF(sf) + } + + # initialize accumulator + # basis functions + sum = 0. + + # loop over y basis functions + yorder = SF_YORDER(sf) + do i = 1, SF_XORDER(sf) { + + # loop over the x basis functions + accum = 0. + do k = 1, yorder { + accum = accum + COEFF(czptr+k) * YBS(yzb+k) + } + accum = accum * XBS(xzb+i) + sum = sum + accum + + # elements of COEFF where neither k = 1 or i = 1 + # are not calculated if SF_XTERMS(sf) = NO + if (SF_XTERMS(sf) == NO) + yorder = 1 + + czptr = czptr + SF_NYCOEFF(sf) + } + + call sfree (sp) + + return (sum) +end diff --git a/math/surfit/isfree.x b/math/surfit/isfree.x new file mode 100644 index 00000000..2eb25891 --- /dev/null +++ b/math/surfit/isfree.x @@ -0,0 +1,45 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "surfitdef.h" + +# ISFREE -- Procedure to free the surface descriptor + +procedure isfree (sf) + +pointer sf # pointer to the surface descriptor +errchk mfree + +begin + # free arrays in memory + + # first the basis functions + if (SF_XBASIS(sf) != NULL) + call mfree (SF_XBASIS(sf), MEM_TYPE) + if (SF_XLEFT(sf) != NULL) + call mfree (SF_XLEFT(sf), TY_INT) + if (SF_YBASIS(sf) != NULL) + call mfree (SF_YBASIS(sf), MEM_TYPE) + if (SF_YLEFT(sf) != NULL) + call mfree (SF_YLEFT(sf), TY_INT) + + # next the x and y matrices + if (SF_XMATRIX(sf) != NULL) + call mfree (SF_XMATRIX(sf), MEM_TYPE) + if (SF_YMATRIX(sf) != NULL) + call mfree (SF_YMATRIX(sf), MEM_TYPE) + + # last the coefficient matrices + if (SF_XCOEFF(sf) != NULL) + call mfree (SF_XCOEFF(sf), MEM_TYPE) + if (SF_COEFF(sf) != NULL) + call mfree (SF_COEFF(sf), MEM_TYPE) + + if (SF_WZ(sf) != NULL) + call mfree (SF_WZ(sf), MEM_TYPE) + if (SF_TLEFT(sf) != NULL) + call mfree (SF_TLEFT(sf), TY_INT) + + # free surface descriptor + if (sf != NULL) + call mfree (sf, TY_STRUCT) +end diff --git a/math/surfit/isinit.x b/math/surfit/isinit.x new file mode 100644 index 00000000..ade35b47 --- /dev/null +++ b/math/surfit/isinit.x @@ -0,0 +1,167 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <math/surfit.h> +include "surfitdef.h" + +# ISINIT -- Procedure to set up the surface descriptor. + +procedure isinit (sf, surf_type, xorder, yorder, xterms, ncols, nlines) + +pointer sf # pointer to surface descriptor structure +int surf_type # type of surface to be fitted +int xorder # x order of surface to be fit, or in the case of the + # spline the number of polynomial pieces in x to be fit +int yorder # y order of surface to be fit, or in the case of the + # spline the number of polynomial pieces in y to be fit +int xterms # cross terms for polynomials? +int ncols # number of columns in the surface +int nlines # number of lines in the surface + +int i +pointer x, y +pointer sp +errchk malloc, calloc + +begin + # allocate space for the surface descriptor + call malloc (sf, LEN_SFSTRUCT, TY_STRUCT) + + if (xorder < 1 || yorder < 1) + call error (0, "SFLINIT: Illegal order.") + + if (ncols < 1) + call error (0, "SFLINIT: x data range is 0.") + if (nlines < 1) + call error (0, "SFLINIT: y data range is 0.") + + switch (surf_type) { + + case SF_CHEBYSHEV, SF_LEGENDRE: + SF_NXCOEFF(sf) = xorder + SF_XORDER(sf) = xorder + SF_XRANGE(sf) = 2. / real (ncols + 1) + SF_XMAXMIN(sf) = - real (ncols + 1) / 2. + SF_XMIN(sf) = 0. + SF_XMAX(sf) = real (ncols + 1) + SF_NYCOEFF(sf) = yorder + SF_YORDER(sf) = yorder + SF_YRANGE(sf) = 2. / real (nlines + 1) + SF_YMAXMIN(sf) = - real (nlines + 1) / 2. + SF_YMIN(sf) = 0. + SF_YMAX(sf) = real (nlines + 1) + SF_XTERMS(sf) = xterms + + case SF_SPLINE3: + SF_NXCOEFF(sf) = (xorder + SPLINE3_ORDER - 1) + SF_XORDER(sf) = SPLINE3_ORDER + SF_NXPIECES(sf) = xorder - 1 + SF_XSPACING(sf) = xorder / real (ncols + 1) + SF_NYCOEFF(sf) = (yorder + SPLINE3_ORDER - 1) + SF_YORDER(sf) = SPLINE3_ORDER + SF_NYPIECES(sf) = yorder - 1 + SF_YSPACING(sf) = yorder / real (nlines + 1) + SF_XMIN(sf) = 0. + SF_XMAX(sf) = real (ncols + 1) + SF_YMIN(sf) = 0. + SF_YMAX(sf) = real (nlines + 1) + SF_XTERMS(sf) = YES + + case SF_SPLINE1: + SF_NXCOEFF(sf) = (xorder + SPLINE1_ORDER - 1) + SF_XORDER(sf) = SPLINE1_ORDER + SF_NXPIECES(sf) = xorder - 1 + SF_XSPACING(sf) = xorder / real (ncols + 1) + SF_NYCOEFF(sf) = (yorder + SPLINE1_ORDER - 1) + SF_YORDER(sf) = SPLINE1_ORDER + SF_NYPIECES(sf) = yorder - 1 + SF_YSPACING(sf) = yorder / real (nlines + 1) + SF_XMIN(sf) = 0. + SF_XMAX(sf) = real (ncols + 1) + SF_YMIN(sf) = 0. + SF_YMAX(sf) = real (nlines + 1) + SF_XTERMS(sf) = YES + + default: + call error (0, "SFINIT: Unknown surface type.") + } + + SF_TYPE(sf) = surf_type + SF_NLINES(sf) = nlines + SF_NCOLS(sf) = ncols + + # allocate space for the matrix and vectors + call calloc (SF_XBASIS(sf), SF_XORDER(sf) * SF_NCOLS(sf), + MEM_TYPE) + call calloc (SF_YBASIS(sf), SF_YORDER(sf) * SF_NLINES(sf), + MEM_TYPE) + call calloc (SF_XMATRIX(sf), SF_XORDER(sf) * SF_NXCOEFF(sf), MEM_TYPE) + call calloc (SF_XCOEFF(sf), SF_NLINES(sf) * SF_NXCOEFF(sf), MEM_TYPE) + call calloc (SF_YMATRIX(sf), SF_YORDER(sf) * SF_NYCOEFF(sf), MEM_TYPE) + call calloc (SF_COEFF(sf), SF_NXCOEFF(sf) * SF_NYCOEFF(sf), MEM_TYPE) + + # allocate temporary space + call smark (sp) + call salloc (x, SF_NCOLS(sf), MEM_TYPE) + call salloc (y, SF_NLINES(sf), MEM_TYPE) + + # calculate all possible x basis functions and store + do i = 1, SF_NCOLS(sf) + Memr[x+i-1] = i + + switch (SF_TYPE(sf)) { + case SF_LEGENDRE: + SF_XLEFT(sf) = NULL + call sf_bleg (Memr[x], SF_NCOLS(sf), SF_XORDER(sf), SF_XMAXMIN(sf), + SF_XRANGE(sf), XBASIS(SF_XBASIS(sf))) + + case SF_CHEBYSHEV: + SF_XLEFT(sf) = NULL + call sf_bcheb (Memr[x], SF_NCOLS(sf), SF_XORDER(sf), SF_XMAXMIN(sf), + SF_XRANGE(sf), XBASIS(SF_XBASIS(sf))) + + case SF_SPLINE3: + call calloc (SF_XLEFT(sf), SF_NCOLS(sf), TY_INT) + call sf_bspline3 (Memr[x], SF_NCOLS(sf), SF_NXPIECES(sf), + -SF_XMIN(sf), SF_XSPACING(sf), XBASIS(SF_XBASIS(sf)), + XLEFT(SF_XLEFT(sf))) + + case SF_SPLINE1: + call calloc (SF_XLEFT(sf), SF_NCOLS(sf), TY_INT) + call sf_bspline1 (Memr[x], SF_NCOLS(sf), SF_NXPIECES(sf), + -SF_XMIN(sf), SF_XSPACING(sf), XBASIS(SF_XBASIS(sf)), + XLEFT(SF_XLEFT(sf))) + } + + # calculate all possible y basis functions and store + do i = 1, SF_NLINES(sf) + Memr[y+i-1] = i + + switch (SF_TYPE(sf)) { + case SF_LEGENDRE: + SF_YLEFT(sf) = NULL + call sf_bleg (Memr[y], SF_NLINES(sf), SF_YORDER(sf), + SF_YMAXMIN(sf), SF_YRANGE(sf), YBASIS(SF_YBASIS(sf))) + + case SF_CHEBYSHEV: + SF_YLEFT(sf) = NULL + call sf_bcheb (Memr[y], SF_NLINES(sf), SF_YORDER(sf), + SF_YMAXMIN(sf), SF_YRANGE(sf), YBASIS(SF_YBASIS(sf))) + + case SF_SPLINE3: + call calloc (SF_YLEFT(sf), SF_NLINES(sf), TY_INT) + call sf_bspline3 (Memr[y], SF_NLINES(sf), SF_NYPIECES(sf), + -SF_YMIN(sf), SF_YSPACING(sf), YBASIS(SF_YBASIS(sf)), + YLEFT(SF_YLEFT(sf))) + + case SF_SPLINE1: + call calloc (SF_YLEFT(sf), SF_NLINES(sf), TY_INT) + call sf_bspline1 (Memr[y], SF_NLINES(sf), SF_NYPIECES(sf), + -SF_YMIN(sf), SF_YSPACING(sf), YBASIS(SF_YBASIS(sf)), + YLEFT(SF_YLEFT(sf))) + } + + SF_WZ(sf) = NULL + SF_TLEFT(sf) = NULL + + call sfree (sp) +end diff --git a/math/surfit/islaccum.x b/math/surfit/islaccum.x new file mode 100644 index 00000000..cc69323a --- /dev/null +++ b/math/surfit/islaccum.x @@ -0,0 +1,117 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <math/surfit.h> +include "surfitdef.h" + +# ISLACCUM -- Procedure to add points on a line to the data set. +# The inner products of the non-zero x basis functions are stored +# in the SF_XORDER(sf) by SF_NXCOEFF(sf) matrix XMATRIX. The +# main diagonal is stored in the first row of XMATRIX. Successive +# non-zero diagonals are stored in the succeeding rows. This method +# of storage is particularly efficient for the large symmetric +# banded matrices produced during spline fits. The inner products +# of the data ordinates and the non-zero x basis functions are +# caculated and stored in the lineno-th row of the SF_NXCOEFF(sf) +# by SF_NLINES(sf) matrix XCOEFF. + +procedure islaccum (sf, cols, lineno, z, w, ncols, wtflag) + +pointer sf # pointer to surface descriptor +int cols[ncols] # column values +int lineno # lineno of data being accumulated +real z[ncols] # surface values on lineno at cols +real w[ncols] # weight of the data points +int ncols # number of data points +int wtflag # type of weighting desired + +int i, ii, j, k +pointer xbzptr, xbptr +pointer xlzptr +pointer xmzptr, xmindex +pointer xczptr, xcindex +pointer bw, rows, left +pointer sp + +begin + # count the number of points + SF_NXPTS(sf) = SF_NXPTS(sf) + ncols + + # calculate the weights, default is uniform weighting + switch (wtflag) { + case SF_UNIFORM: + call amovkr (1.0, w, ncols) + case SF_USER: + # do not alter weights + default: + call amovkr (1.0, w, ncols) + } + + # set up temporary storage + call smark (sp) + call salloc (bw, ncols, TY_REAL) + call salloc (left, ncols, TY_INT) + call salloc (rows, ncols, TY_INT) + + # set up the pointers + xbzptr = SF_XBASIS(sf) - 1 + xmzptr = SF_XMATRIX(sf) + xczptr = SF_XCOEFF(sf) + (lineno - 1) * SF_NXCOEFF(sf) - 1 + + # accumulate the line + switch (SF_TYPE(sf)) { + case SF_LEGENDRE, SF_CHEBYSHEV: + + do i = 1, SF_XORDER(sf) { + do j = 1, ncols + Memr[bw+j-1] = w[j] * XBASIS(xbzptr+cols[j]) + xcindex = xczptr + i + do j = 1, ncols + XCOEFF(xcindex) = XCOEFF(xcindex) + Memr[bw+j-1] * z[j] + xbptr = xbzptr + ii = 0 + do k = i, SF_XORDER(sf) { + xmindex = xmzptr + ii + do j = 1, ncols + XMATRIX(xmindex) = XMATRIX(xmindex) + Memr[bw+j-1] * + XBASIS(xbptr+cols[j]) + ii = ii + 1 + xbptr = xbptr + SF_NCOLS(sf) + } + xbzptr = xbzptr + SF_NCOLS(sf) + xmzptr = xmzptr + SF_XORDER(sf) + } + + case SF_SPLINE3, SF_SPLINE1: + + xlzptr = SF_XLEFT(sf) - 1 + do j = 1, ncols + Memi[left+j-1] = XLEFT(xlzptr+cols[j]) + call amulki (Memi[left], SF_XORDER(sf), Memi[rows], ncols) + call aaddki (Memi[rows], SF_XMATRIX(sf), Memi[rows], ncols) + call aaddki (Memi[left], xczptr, Memi[left], ncols) + + do i = 1, SF_XORDER(sf) { + do j = 1, ncols { + Memr[bw+j-1] = w[j] * XBASIS(xbzptr+cols[j]) + xcindex = Memi[left+j-1] + i + XCOEFF(xcindex) = XCOEFF(xcindex) + Memr[bw+j-1] * z[j] + } + xbptr = xbzptr + ii = 0 + do k = i, SF_XORDER(sf) { + do j = 1, ncols { + xmindex = Memi[rows+j-1] + ii + XMATRIX(xmindex) = XMATRIX(xmindex) + Memr[bw+j-1] * + XBASIS(xbptr+cols[j]) + } + ii = ii + 1 + xbptr = xbptr + SF_NCOLS(sf) + } + xbzptr = xbzptr + SF_NCOLS(sf) + call aaddki (Memi[rows], SF_XORDER(sf), Memi[rows], ncols) + } + } + + # release space + call sfree (sp) +end diff --git a/math/surfit/islfit.x b/math/surfit/islfit.x new file mode 100644 index 00000000..7b60c361 --- /dev/null +++ b/math/surfit/islfit.x @@ -0,0 +1,150 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <math/surfit.h> +include "surfitdef.h" + +# ISLFIT -- Procedure to fit a single line of a surface. The inner products +# of the x basis functions are calculated and accumulated into the +# SF_XORDER(sf) by SF_NXCOEFF(sf) matrix XMATRIX. The main diagonal is stored +# in the first row of XMATRIX followed by the non-zero lower diagonals. This +# method of storage is very efficient for the large symmetric banded matrices +# generated by spline fits. The Cholesky factorization of XMATRIX is calculated +# and stored in XMATRIX destroying the original data. The inner products +# of the x basis functions and the data ordinates are stored in the lineno-th +# row of the SF_NXCOEFF(sf) by SF_NLINES(sf) matrix XCOEFF. The x coefficients +# for each line are calculated by forward and back substitution and +# stored in the lineno-th line of XCOEFF destroying the original data. + +procedure islfit (sf, cols, lineno, z, w, ncols, wtflag, ier) + +pointer sf # pointer to the surface descriptor +int cols[ncols] # array of columns +int lineno # lineno +real z[ncols] # the surface values +real w[ncols] # array of weights +int ncols # the number of columns +int wtflag # type of weighting +int ier # error codes + +int i, ii, j, k +pointer xbzptr, xbptr +pointer xmzptr, xmindex +pointer xczptr, xcindex +pointer xlzptr +pointer left, rows, bw +pointer sp +real adotr() + +begin + # index the pointers + xbzptr = SF_XBASIS(sf) - 1 + xmzptr = SF_XMATRIX(sf) + xczptr = SF_XCOEFF(sf) + (lineno - 1) * SF_NXCOEFF(sf) - 1 + + # zero the accumulators + call aclrr (XMATRIX(SF_XMATRIX(sf)), SF_NXCOEFF(sf) * SF_XORDER(sf)) + call aclrr (XCOEFF(xczptr + 1), SF_NXCOEFF(sf)) + + # free space used by previous islrefit calls + if (SF_WZ(sf) != NULL) + call mfree (SF_WZ(sf), MEM_TYPE) + if (SF_TLEFT(sf) != NULL) + call mfree (SF_TLEFT(sf), TY_INT) + + # reset number of points + SF_NXPTS(sf) = ncols + + # calculate the weights, default is uniform weighting + switch (wtflag) { + case SF_UNIFORM: + call amovkr (1.0, w, ncols) + case SF_USER: + # do not alter weights + default: + call amovkr (1.0, w, ncols) + } + + # allocate temporary space + call smark (sp) + call salloc (bw, ncols, TY_REAL) + + switch (SF_TYPE(sf)) { + case SF_LEGENDRE, SF_CHEBYSHEV: + + do i = 1, SF_XORDER(sf) { + do j = 1, ncols + Memr[bw+j-1] = w[j] * XBASIS(xbzptr+cols[j]) + + xcindex = xczptr + i + XCOEFF(xcindex) = XCOEFF(xcindex) + adotr (Memr[bw], z, ncols) + + xbptr = xbzptr + ii = 0 + do k = i, SF_XORDER(sf) { + xmindex = xmzptr + ii + do j = 1, ncols + XMATRIX(xmindex) = XMATRIX(xmindex) + Memr[bw+j-1] * + XBASIS(xbptr+cols[j]) + ii = ii + 1 + xbptr = xbptr + SF_NCOLS(sf) + } + + xbzptr = xbzptr + SF_NCOLS(sf) + xmzptr = xmzptr + SF_XORDER(sf) + } + + case SF_SPLINE3, SF_SPLINE1: + xlzptr = SF_XLEFT(sf) - 1 + + call salloc (left, ncols, TY_INT) + call salloc (rows, ncols, TY_INT) + + do j = 1, ncols + Memi[left+j-1] = XLEFT(xlzptr+cols[j]) + call amulki (Memi[left], SF_XORDER(sf), Memi[rows], ncols) + call aaddki (Memi[rows], SF_XMATRIX(sf), Memi[rows], ncols) + call aaddki (Memi[left], xczptr, Memi[left], ncols) + + do i = 1, SF_XORDER(sf) { + do j = 1, ncols { + Memr[bw+j-1] = w[j] * XBASIS(xbzptr+cols[j]) + xcindex = Memi[left+j-1] + i + XCOEFF(xcindex) = XCOEFF(xcindex) + Memr[bw+j-1] * z[j] + } + + xbptr = xbzptr + ii = 0 + do k = i, SF_XORDER(sf) { + do j = 1, ncols { + xmindex = Memi[rows+j-1] + ii + XMATRIX(xmindex) = XMATRIX(xmindex) + Memr[bw+j-1] * + XBASIS(xbptr+cols[j]) + } + ii = ii + 1 + xbptr = xbptr + SF_NCOLS(sf) + } + + xbzptr = xbzptr + SF_NCOLS(sf) + call aaddki (Memi[rows], SF_XORDER(sf), Memi[rows], ncols) + } + } + + # release space + call sfree (sp) + + # return if not enough data points + ier = OK + if ((SF_NXPTS(sf) - SF_NXCOEFF(sf)) < 0) { + ier = NO_DEG_FREEDOM + return + } + + # calculate the Cholesky factorization of XMATRIX + call sfchofac (XMATRIX(SF_XMATRIX(sf)), SF_XORDER(sf), SF_NXCOEFF(sf), + XMATRIX(SF_XMATRIX(sf)), ier) + + # calculate the x coefficients for lineno-th image line and place in the + # lineno-th row of XCOEFF + call sfchoslv (XMATRIX(SF_XMATRIX(sf)), SF_XORDER(sf), SF_NXCOEFF(sf), + XCOEFF(xczptr + 1), XCOEFF(xczptr + 1)) +end diff --git a/math/surfit/islrefit.x b/math/surfit/islrefit.x new file mode 100644 index 00000000..b0430e5b --- /dev/null +++ b/math/surfit/islrefit.x @@ -0,0 +1,74 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <math/surfit.h> +include "surfitdef.h" + +# ISLREFIT -- Procedure to refit the data assuming that the cols and w +# arrays do not change. SIFLREFIT assumes that the Cholesky factorization +# of XMATRIX is stored in XMATRIX. The inner products of the x basis +# functions and the data ordinates are accumulated into the lineno-th +# row of the SF_NXCOEFF(sf) by SF_NLINES(sf) matrix XCOEFF. The coefficients +# for line number lineno are calculated by forward and back +# substitution and placed in the lineno-th row of XCOEFF replacing the +# original data. + +procedure islrefit (sf, cols, lineno, z, w) + +pointer sf # pointer to surface descriptor +int cols[ARB] # columns to be fit +int lineno # line number +real z[ARB] # surface values +real w[ARB] # weight values + +int i, j +pointer xbzptr, xczptr, xcindex, xlzptr + +begin + # set pointers + xbzptr = SF_XBASIS(sf) - 1 + xczptr = SF_XCOEFF(sf) + (lineno - 1) * SF_NXCOEFF(sf) - 1 + xlzptr = SF_XLEFT(sf) - 1 + + # reset lineno-th row of the x coefficient matrix + call aclrr (XCOEFF(xczptr+1), SF_NXCOEFF(sf)) + + if (SF_WZ(sf) == NULL) + call malloc (SF_WZ(sf), SF_NXPTS(sf), MEM_TYPE) + + # calculate new right sides + call amulr (w, z, Memr[SF_WZ(sf)], SF_NXPTS(sf)) + + switch (SF_TYPE(sf)) { + case SF_LEGENDRE, SF_CHEBYSHEV: + + do i = 1, SF_XORDER(sf) { + xcindex = xczptr + i + do j = 1, SF_NXPTS(sf) + XCOEFF(xcindex) = XCOEFF(xcindex) + Memr[SF_WZ(sf)+j-1] * + XBASIS(xbzptr+cols[j]) + xbzptr = xbzptr + SF_NCOLS(sf) + } + + case SF_SPLINE3, SF_SPLINE1: + + if (SF_TLEFT(sf) == NULL) + call malloc (SF_TLEFT(sf), SF_NXPTS(sf), TY_INT) + + do i = 1, SF_NXPTS(sf) + Memi[SF_TLEFT(sf)+i-1] = XLEFT(xlzptr+cols[i]) + xczptr + + do i = 1, SF_XORDER(sf) { + do j = 1, SF_NXPTS(sf) { + xcindex = Memi[SF_TLEFT(sf)+j-1] + i + XCOEFF(xcindex) = XCOEFF(xcindex) + Memr[SF_WZ(sf)+j-1] * + XBASIS(xbzptr+cols[j]) + } + xbzptr = xbzptr + SF_NCOLS(sf) + } + + } + + # solve for the new x coefficients for line number lineno + call sfchoslv (XMATRIX(SF_XMATRIX(sf)), SF_XORDER(sf), SF_NXCOEFF(sf), + XCOEFF(xczptr+1), XCOEFF(xczptr+1)) +end diff --git a/math/surfit/islsolve.x b/math/surfit/islsolve.x new file mode 100644 index 00000000..11d0d313 --- /dev/null +++ b/math/surfit/islsolve.x @@ -0,0 +1,48 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <math/surfit.h> +include "surfitdef.h" + +# ISLSOLVE -- Procedure to solve for the x coefficients of image line +# number lineno. The inner products of the x basis functions are assumed +# to be in the SF_XORDER(sf) by SF_NXCOEFF(sf) array XMATRIX, +# while the inner products of the basis functions and +# the data ordinated for line number lineno are assumed to be in the +# lineno-th row of the SF_NXCOEFF(sf) by SF_NLINES(sf) matrix XCOEFF. +# The Cholesky factorization of XMATRIX is calculated and placed +# in XMATRIX overwriting the original data. The x coefficients for +# line number lineno are calculated and placed in the lineno-th row +# of XCOEFF replacing the original data. + +procedure islsolve (sf, lineno, ier) + +pointer sf # pointer to the surface descriptor structure +int lineno # line being fitted in x +int ier # ier = 0, everything OK + # ier = 1, matrix is singular + # ier = 2, no degree of freedom + +pointer xcptr + +begin + # return if there are insuffucient points to solve the matrix + ier = OK + if ((SF_NXPTS(sf) - SF_NXCOEFF(sf)) < 0 ) { + ier = NO_DEG_FREEDOM + return + } + + # calculate the Cholesky factorization of the x matrix and store + # separately for possible use by SFLREFIT + + call sfchofac (XMATRIX(SF_XMATRIX(sf)), SF_XORDER(sf), SF_NXCOEFF(sf), + XMATRIX(SF_XMATRIX(sf)), ier) + + # solve for the x coefficients for line lineno assuming the + # data are in row lineno of xcoeff, the solution is placed + # on top of the data + + xcptr = SF_XCOEFF(sf) + (lineno - 1) * SF_NXCOEFF(sf) + call sfchoslv (XMATRIX(SF_XMATRIX(sf)), SF_XORDER(sf), SF_NXCOEFF(sf), + XCOEFF(xcptr), XCOEFF(xcptr)) +end diff --git a/math/surfit/islzero.x b/math/surfit/islzero.x new file mode 100644 index 00000000..03034f01 --- /dev/null +++ b/math/surfit/islzero.x @@ -0,0 +1,25 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "surfitdef.h" + +# ISLZERO -- Procedure to zero the accumulators for line number lineno. + +procedure islzero (sf, lineno) + +pointer sf # pointer to the surface descriptor +int lineno # line number +pointer xcptr + +begin + SF_NXPTS(sf) = 0 + + call aclrr (XMATRIX(SF_XMATRIX(sf)), + SF_XORDER(sf) * SF_NXCOEFF(sf)) + xcptr = SF_XCOEFF(sf) + (lineno - 1) * SF_NXCOEFF(sf) + call aclrr (XCOEFF(xcptr), SF_NXCOEFF(sf)) + + if (SF_WZ(sf) != NULL) + call mfree (SF_WZ(sf), MEM_TYPE) + if (SF_TLEFT(sf) != NULL) + call mfree (SF_TLEFT(sf), TY_INT) +end diff --git a/math/surfit/isreplace.x b/math/surfit/isreplace.x new file mode 100644 index 00000000..1ba69479 --- /dev/null +++ b/math/surfit/isreplace.x @@ -0,0 +1,114 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <math/surfit.h> +include "surfitdef.h" + +# ISREPLACE -- Procedure to restore the surface fit stored by SIFSAVE +# to the surface descriptor for use by the evaluating routines. The +# surface parameters, surface type, xorder (or number of polynomial +# pieces in x), yorder (or number of polynomial pieces in y), xterms, +# number of columns and number of lines, are stored in the first +# six elements of the real array fit, followed by the SF_NYCOEFF(sf) * +# SF_NYCOEFF(sf) surface coefficients. The coefficient of B(i,x) * B(j,y) +# is stored in element number 6 + (i - 1) * SF_NYCOEFF(sf) + j of the +# array fit. + +procedure isreplace (sf, fit) + +pointer sf # surface descriptor +real fit[ARB] # array containing the surface parameters and + # coefficients + +int surface_type, xorder, yorder, ncols, nlines + +begin + # allocate space for the surface descriptor + call calloc (sf, LEN_SFSTRUCT, TY_STRUCT) + + xorder = nint (SF_SAVEXORDER(fit)) + if (xorder < 1) + call error (0, "SFRESTORE: Illegal x order.") + yorder = nint (SF_SAVEYORDER(fit)) + if (yorder < 1) + call error (0, "SFRESTORE: Illegal y order.") + + ncols = nint (SF_SAVENCOLS(fit)) + if (ncols < 1) + call error (0, "SFRESTORE: Illegal x range.") + nlines = nint (SF_SAVENLINES(fit)) + if (nlines < 1) + call error (0, "SFRESTORE: Illegal y range.") + + # set surface type dependent surface descriptor parameters + surface_type = nint (SF_SAVETYPE(fit)) + + switch (surface_type) { + case SF_LEGENDRE, SF_CHEBYSHEV: + SF_NXCOEFF(sf) = xorder + SF_XORDER(sf) = xorder + SF_XRANGE(sf) = 2. / real (ncols + 1) + SF_XMAXMIN(sf) = - real (ncols + 1) / 2. + SF_XMIN(sf) = 0. + SF_XMAX(sf) = real (ncols + 1) + SF_NYCOEFF(sf) = yorder + SF_YORDER(sf) = yorder + SF_YRANGE(sf) = 2. / real (nlines + 1) + SF_YMAXMIN(sf) = - real (nlines + 1) / 2. + SF_YMIN(sf) = 0. + SF_YMAX(sf) = real (nlines + 1) + SF_XTERMS(sf) = SF_SAVEXTERMS(fit) + + case SF_SPLINE3: + SF_NXCOEFF(sf) = (xorder + SPLINE3_ORDER - 1) + SF_XORDER(sf) = SPLINE3_ORDER + SF_NXPIECES(sf) = xorder - 1 + SF_XSPACING(sf) = xorder / real (ncols + 1) + SF_XMIN(sf) = 0. + SF_XMAX(sf) = real (ncols + 1) + SF_NYCOEFF(sf) = (yorder + SPLINE3_ORDER - 1) + SF_YORDER(sf) = SPLINE3_ORDER + SF_NYPIECES(sf) = yorder - 1 + SF_YSPACING(sf) = yorder / real (nlines + 1) + SF_YMIN(sf) = 0. + SF_YMAX(sf) = real (nlines + 1) + SF_XTERMS(sf) = YES + + case SF_SPLINE1: + SF_NXCOEFF(sf) = (xorder + SPLINE1_ORDER - 1) + SF_XORDER(sf) = SPLINE1_ORDER + SF_NXPIECES(sf) = xorder - 1 + SF_XSPACING(sf) = xorder / real (ncols + 1) + SF_XMIN(sf) = 0. + SF_XMAX(sf) = real (ncols + 1) + SF_NYCOEFF(sf) = (yorder + SPLINE1_ORDER - 1) + SF_YORDER(sf) = SPLINE1_ORDER + SF_NYPIECES(sf) = yorder - 1 + SF_YSPACING(sf) = yorder / real (nlines + 1) + SF_YMIN(sf) = 0. + SF_YMAX(sf) = real (nlines + 1) + SF_XTERMS(sf) = YES + + default: + call error (0, "SFRESTORE: Unknown surface type.") + } + + # set remaining curve parameters + SF_TYPE(sf) = surface_type + SF_NLINES(sf) = nlines + SF_NCOLS(sf) = ncols + + # allocate space for the coefficient array + SF_XBASIS(sf) = NULL + SF_YBASIS(sf) = NULL + SF_XMATRIX(sf) = NULL + SF_YMATRIX(sf) = NULL + SF_XCOEFF(sf) = NULL + SF_WZ(sf) = NULL + SF_TLEFT(sf) = NULL + + call calloc (SF_COEFF(sf), SF_NXCOEFF(sf) * SF_NYCOEFF(sf), MEM_TYPE) + + # restore coefficient array + call amovr (fit[SF_SAVECOEFF+1], COEFF(SF_COEFF(sf)), SF_NYCOEFF(sf) * + SF_NXCOEFF(sf)) +end diff --git a/math/surfit/isresolve.x b/math/surfit/isresolve.x new file mode 100644 index 00000000..7af6a8da --- /dev/null +++ b/math/surfit/isresolve.x @@ -0,0 +1,127 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <math/surfit.h> +include "surfitdef.h" + +# ISRESOLVE -- Procedure to solve the x coefficient matrix for the surface +# coefficients assuming that the lines array is unchanged since the last +# call to SIFSOLVE. The Cholesky factorization of YMATRIX is assumed to be +# stored in YMATRIX. The inner product of the y basis functions and i-th +# column of XCOEFF containing the i-th x coefficients for each line are +# calculated and stored in the i-th row of the SF_NYCOEFF(sf) by +# SF_NXCOEFF(sf) array COEFF. Each of the SF_NXCOEFF(sf) rows of COEFF +# is solved to determine the SF_NYCOEFF(sf) by SF_NXCOEFF(sf) surface +# coefficients. After a call to SIFSOLVE the coefficient for the i-th +# y basis function and the j-th x coefficient is found in the j-th row +# and i-th column of COEFF. + +procedure isresolve (sf, lines, ier) + +pointer sf # pointer to the surface descriptor structure +int lines[ARB] # line numbers included in the fit +int ier # error code + +int i, j, k, nxcoeff +pointer ybzptr +pointer ylzptr +pointer xczptr, xcptr, xcindex +pointer czptr, cptr +pointer left, tleft +pointer sp + +begin + # define pointers + ybzptr = SF_YBASIS(sf) - 1 + xczptr = SF_XCOEFF(sf) - SF_NXCOEFF(sf) - 1 + czptr = SF_COEFF(sf) - 1 + + # zero out coefficient matrix + call aclrr (COEFF(SF_COEFF(sf)), SF_NXCOEFF(sf) * SF_NYCOEFF(sf)) + + switch (SF_TYPE(sf)) { + case SF_LEGENDRE, SF_CHEBYSHEV: + + # loop over the y basis functions + nxcoeff = SF_NXCOEFF(sf) + do i = 1, SF_YORDER(sf) { + cptr = czptr + i + do k = 1, nxcoeff { + xcptr = xczptr + k + do j = 1, SF_NYPTS(sf) { + xcindex = xcptr + lines[j] * SF_NXCOEFF(sf) + COEFF(cptr) = COEFF(cptr) + + YBASIS(ybzptr+lines[j]) * XCOEFF(xcindex) + } + cptr = cptr + SF_NYCOEFF(sf) + } + + ybzptr = ybzptr + SF_NYPTS(sf) + + # if SF_XTERMS(sf) = NO do not accumulate elements + # of COEFF where i != 1 and k != 1 + if (SF_XTERMS(sf) == NO) + nxcoeff = 1 + } + + case SF_SPLINE3, SF_SPLINE1: + call smark (sp) + call salloc (left, SF_NYPTS(sf), TY_INT) + call salloc (tleft, SF_NYPTS(sf), TY_INT) + + ylzptr = SF_YLEFT(sf) - 1 + do j = 1, SF_NYPTS(sf) + Memi[left+j-1] = YLEFT(ylzptr+lines[j]) + call aaddki (Memi[left], czptr, Memi[left], SF_NYPTS(sf)) + + nxcoeff = SF_NXCOEFF(sf) + do i = 1, SF_YORDER(sf) { + call aaddki (Memi[left], i, Memi[tleft], SF_NYPTS(sf)) + do k = 1, nxcoeff { + xcptr = xczptr + k + do j = 1, SF_NYPTS(sf) { + cptr = Memi[tleft+j-1] + xcindex = xcptr + lines[j] * SF_NXCOEFF(sf) + COEFF(cptr) = COEFF(cptr) + YBASIS(ybzptr+lines[j]) * + XCOEFF(xcindex) + } + call aaddki (Memi[tleft], SF_NYCOEFF(sf), Memi[tleft], + SF_NYPTS(sf)) + } + + ybzptr = ybzptr + SF_NYPTS(sf) + } + + call sfree (sp) + } + + # return if not enough points + ier = OK + if ((SF_NYPTS(sf) - SF_NYCOEFF(sf)) < 0) { + ier = NO_DEG_FREEDOM + return + } + + if (SF_XTERMS(sf) == YES) { + + # solve for the nxcoeff right sides + cptr = SF_COEFF(sf) + do i = 1, SF_NXCOEFF(sf) { + call sfchoslv (YMATRIX(SF_YMATRIX(sf)), SF_YORDER(sf), + SF_NYCOEFF(sf), COEFF(cptr), COEFF(cptr)) + cptr = cptr + SF_NYCOEFF(sf) + } + + } else { + + # solve for the y coefficients + cptr = SF_NYCOEFF(sf) + call sfchoslv (YMATRIX(SF_YMATRIX(sf)), SF_YORDER(sf), + SF_NYCOEFF(sf), COEFF(cptr), COEFF(cptr)) + + # solve for the x coefficients + do i = 2, SF_NXCOEFF(sf) { + cptr = czptr + SF_NYCOEFF(sf) + COEFF(cptr) = COEFF(cptr) / SF_NYPTS(sf) + } + } +end diff --git a/math/surfit/issave.x b/math/surfit/issave.x new file mode 100644 index 00000000..5a3ecab8 --- /dev/null +++ b/math/surfit/issave.x @@ -0,0 +1,44 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <math/surfit.h> +include "surfitdef.h" + +# ISSAVE -- Procedure to save the surface fit for later use by the +# evaluate routines. After a call to SIFSAVE the first six elements +# of fit contain the surface type, xorder (or number of polynomial pieces +# in x), yorder (or the number of polynomial pieces in y), xterms, ncols +# and nlines. The remaining spaces are filled by the SF_NYCOEFF(sf) * +# SF_NXCOEFF(sf) surface coefficients. The coefficient of B(i,x) * B(j,y) +# is located in element number 6 + (i - 1) * SF_NYCOEFF(sf) + j of the +# array fit where i <= SF_NXCOEFF(sf) and j <= SF_NYCOEFF(sf). + +procedure issave (sf, fit) + +pointer sf # pointer to the surface descriptor +real fit[ARB] # array for storing fit + +begin + # get the surface parameters + + # order is surface type dependent + switch (SF_TYPE(sf)) { + case SF_LEGENDRE, SF_CHEBYSHEV: + SF_SAVEXORDER(fit) = SF_XORDER(sf) + SF_SAVEYORDER(fit) = SF_YORDER(sf) + case SF_SPLINE3, SF_SPLINE1: + SF_SAVEXORDER(fit) = SF_NXPIECES(sf) + 1 + SF_SAVEYORDER(fit) = SF_NYPIECES(sf) + 1 + default: + call error (0, "SIFSAVE: Unknown surface type.") + } + + # save remaining parameters + SF_SAVETYPE(fit) = SF_TYPE(sf) + SF_SAVENLINES(fit) = SF_NLINES(sf) + SF_SAVENCOLS(fit) = SF_NCOLS(sf) + SF_SAVEXTERMS(fit) = SF_XTERMS(sf) + + # save the coefficients + call amovr (COEFF(SF_COEFF(sf)), fit[SF_SAVECOEFF+1], SF_NXCOEFF(sf) * + SF_NYCOEFF(sf)) +end diff --git a/math/surfit/issolve.x b/math/surfit/issolve.x new file mode 100644 index 00000000..7790ef60 --- /dev/null +++ b/math/surfit/issolve.x @@ -0,0 +1,169 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <math/surfit.h> +include "surfitdef.h" + +# ISSOLVE -- Procedure to solve the x coefficient matrix for the surface +# coefficients. The inner products of the y basis functions are accumulated +# and stored in the SF_YORDER(sf) by SF_NYCOEFF(sf) array YMATRIX. The +# main diagonal of YMATRIX is stored in the first row of YMATRIX followed +# by the remaining non-zero lower diagonals. The Cholesky factorization +# of YMATRIX is calculated and stored on top of YMATRIX destroying the +# original data. The inner products of the y basis functions and +# and the i-th column of the SF_NXCOEFF(sf) by SF_NLINES(sf) matrix XCOEFF +# containing the i-th x coefficients for each line are calculated and +# placed in the i-th row of the SF_NYCOEFF(sf) by SF_NXCOEFF(sf) matrix +# COEFF. Each of the SF_NXCOEFF(sf) rows of COEFF is solved to determine +# the SF_NXCOEFF(sf) by SF_NYCOEFF(sf) surface coefficients. After a +# call to SIFSOLVE the coefficient of the i-th x basis function and the +# j-th y basis function will be found in the j-th column and i-th row +# of COEFF. + +procedure issolve (sf, lines, nlines, ier) + +pointer sf # pointer to the curve descriptor structure +int lines[ARB] # line numbers included in the fit +int nlines # number of lines fit +int ier # error code + +int i, ii, j, k, nxcoeff +pointer ybzptr, ybptr +pointer ylzptr +pointer ymzptr, ymindex +pointer xczptr, xcptr, xcindex +pointer czptr, cptr +pointer left, tleft, rows +pointer sp + +begin + # define pointers + ybzptr = SF_YBASIS(sf) - 1 + ymzptr = SF_YMATRIX(sf) + xczptr = SF_XCOEFF(sf) - SF_NXCOEFF(sf) - 1 + czptr = SF_COEFF(sf) - 1 + + # zero out coefficient matrix and the y coefficient matrix + call aclrr (YMATRIX(SF_YMATRIX(sf)), SF_YORDER(sf) * SF_NYCOEFF(sf)) + call aclrr (COEFF(SF_COEFF(sf)), SF_NXCOEFF(sf) * SF_NYCOEFF(sf)) + + # increment the number of points + SF_NYPTS(sf) = nlines + + switch (SF_TYPE(sf)) { + case SF_LEGENDRE, SF_CHEBYSHEV: + + # accumulate the y value in the y matrix + nxcoeff = SF_NXCOEFF(sf) + do i = 1, SF_YORDER(sf) { + cptr = czptr + i + do k = 1, nxcoeff { + xcptr = xczptr + k + do j = 1, nlines { + xcindex = xcptr + lines[j] * SF_NXCOEFF(sf) + COEFF(cptr) = COEFF(cptr) + + YBASIS(ybzptr+lines[j]) * XCOEFF(xcindex) + } + cptr = cptr + SF_NYCOEFF(sf) + } + ii = 0 + ybptr = ybzptr + do k = i, SF_YORDER(sf) { + ymindex = ymzptr + ii + do j = 1, nlines + YMATRIX(ymindex) = YMATRIX(ymindex) + + YBASIS(ybzptr+lines[j]) * + YBASIS(ybptr+lines[j]) + ii = ii + 1 + ybptr = ybptr + SF_NLINES(sf) + } + + if (SF_XTERMS(sf) == NO) + nxcoeff = 1 + + ybzptr = ybzptr + SF_NLINES(sf) + ymzptr = ymzptr + SF_YORDER(sf) + } + + case SF_SPLINE3, SF_SPLINE1: + + call smark (sp) + call salloc (left, nlines, TY_INT) + call salloc (tleft, nlines, TY_INT) + call salloc (rows, nlines, TY_INT) + + ylzptr = SF_YLEFT(sf) - 1 + do j = 1, nlines + Memi[left+j-1] = YLEFT(ylzptr+lines[j]) + call amulki (Memi[left], SF_YORDER(sf), Memi[rows], nlines) + call aaddki (Memi[rows], SF_YMATRIX(sf), Memi[rows], nlines) + call aaddki (Memi[left], czptr, Memi[left], nlines) + + # accumulate the y value in the y matrix + nxcoeff = SF_NXCOEFF(sf) + do i = 1, SF_YORDER(sf) { + call aaddki (Memi[left], i, Memi[tleft], nlines) + do k = 1, nxcoeff { + xcptr = xczptr + k + do j = 1, nlines { + cptr = Memi[tleft+j-1] + xcindex = xcptr + lines[j] * SF_NXCOEFF(sf) + COEFF(cptr) = COEFF(cptr) + YBASIS(ybzptr+lines[j]) * + XCOEFF(xcindex) + } + call aaddki (Memi[tleft], SF_NYCOEFF(sf), Memi[tleft], + nlines) + } + ii = 0 + ybptr = ybzptr + do k = i, SF_YORDER(sf) { + do j = 1, nlines { + ymindex = Memi[rows+j-1] + ii + YMATRIX(ymindex) = YMATRIX(ymindex) + + YBASIS(ybzptr+lines[j]) * + YBASIS(ybptr+lines[j]) + } + ii = ii + 1 + ybptr = ybptr + SF_NLINES(sf) + } + + ybzptr = ybzptr + SF_NLINES(sf) + call aaddki (Memi[rows], SF_YORDER(sf), Memi[rows], nlines) + } + + call sfree (sp) + + } + + # return if not enough points + ier = OK + if ((SF_NYPTS(sf) - SF_NYCOEFF(sf)) < 0) { + ier = NO_DEG_FREEDOM + return + } + + # calculate the Cholesky factorization of the y matrix + call sfchofac (YMATRIX(SF_YMATRIX(sf)), SF_YORDER(sf), SF_NYCOEFF(sf), + YMATRIX(SF_YMATRIX(sf)), ier) + + if (SF_XTERMS(sf) == YES) { + + # solve for the nxcoeff right sides + cptr = SF_COEFF(sf) + do i = 1, SF_NXCOEFF(sf) { + call sfchoslv (YMATRIX(SF_YMATRIX(sf)), SF_YORDER(sf), + SF_NYCOEFF(sf), COEFF(cptr), COEFF(cptr)) + cptr = cptr + SF_NYCOEFF(sf) + } + + } else { + + cptr = SF_COEFF(sf) + call sfchoslv (YMATRIX(SF_YMATRIX(sf)), SF_YORDER(sf), + SF_NYCOEFF(sf), COEFF(cptr), COEFF(cptr)) + + do i = 2, SF_NXCOEFF(sf) { + cptr = cptr + SF_NYCOEFF(sf) + COEFF(cptr) = COEFF(cptr) / SF_NYPTS(sf) + } + } +end diff --git a/math/surfit/isvector.x b/math/surfit/isvector.x new file mode 100644 index 00000000..023d3f4d --- /dev/null +++ b/math/surfit/isvector.x @@ -0,0 +1,76 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <math/surfit.h> +include "surfitdef.h" + +# ISVECTOR -- Procedure to evaluate the fitted surface at an array of points. +# The SF_NXCOEFF(sf) by SF_NYCOEFF(sf) coefficients are stored in the +# SF_NYCOEFF(sf) by SF_NXCOEFF(sf) matrix COEFF. The j-th element of the ith +# row of COEFF contains the coefficient of the i-th basis function in x and +# the j-th basis function in y. + +procedure isvector (sf, x, y, zfit, npts) + +pointer sf # pointer to surface descriptor structure +real x[ARB] # x value +real y[ARB] # y value +real zfit[ARB] # fits surface values +int npts # number of data points + +int i +pointer xcoeff, cptr, sp + +begin + # evaluate the surface along the vector + switch (SF_TYPE(sf)) { + case SF_CHEBYSHEV: + if (SF_XORDER(sf) == 1) { + call cv_evcheb (COEFF(SF_COEFF(sf)), y, zfit, npts, + SF_YORDER(sf), SF_YMAXMIN(sf), SF_YRANGE(sf)) + } else if (SF_YORDER(sf) == 1) { + call smark (sp) + call salloc (xcoeff, SF_NXCOEFF(sf), MEM_TYPE) + cptr = SF_COEFF(sf) + do i = 1, SF_NXCOEFF(sf) { + Memr[xcoeff+i-1] = COEFF(cptr) + cptr = cptr + SF_NYCOEFF(sf) + } + call cv_evcheb (Memr[xcoeff], x, zfit, npts, + SF_XORDER(sf), SF_XMAXMIN(sf), SF_XRANGE(sf)) + call sfree (sp) + } else + call sf_evcheb (COEFF(SF_COEFF(sf)), x, y, zfit, npts, + SF_XTERMS(sf), SF_XORDER(sf), SF_YORDER(sf), SF_XMAXMIN(sf), + SF_XRANGE(sf), SF_YMAXMIN(sf), SF_YRANGE(sf)) + + case SF_LEGENDRE: + if (SF_XORDER(sf) == 1) { + call cv_evleg (COEFF(SF_COEFF(sf)), y, zfit, npts, + SF_YORDER(sf), SF_YMAXMIN(sf), SF_YRANGE(sf)) + } else if (SF_YORDER(sf) == 1) { + call smark (sp) + call salloc (xcoeff, SF_NXCOEFF(sf), MEM_TYPE) + cptr = SF_COEFF(sf) + do i = 1, SF_NXCOEFF(sf) { + Memr[xcoeff+i-1] = COEFF(cptr) + cptr = cptr + SF_NYCOEFF(sf) + } + call cv_evcheb (Memr[xcoeff], x, zfit, npts, + SF_XORDER(sf), SF_XMAXMIN(sf), SF_XRANGE(sf)) + call sfree (sp) + } else + call sf_evleg (COEFF(SF_COEFF(sf)), x, y, zfit, npts, + SF_XTERMS(sf), SF_XORDER(sf), SF_YORDER(sf), SF_XMAXMIN(sf), + SF_XRANGE(sf), SF_YMAXMIN(sf), SF_YRANGE(sf)) + + case SF_SPLINE3: + call sf_evspline3 (COEFF(SF_COEFF(sf)), x, y, zfit, npts, + SF_NXPIECES(sf), SF_NYPIECES(sf), -SF_XMIN(sf), SF_XSPACING(sf), + -SF_YMIN(sf), SF_YSPACING(sf)) + + case SF_SPLINE1: + call sf_evspline1 (COEFF(SF_COEFF(sf)), x, y, zfit, npts, + SF_NXPIECES(sf), SF_NYPIECES(sf), -SF_XMIN(sf), SF_XSPACING(sf), + -SF_YMIN(sf), SF_YSPACING(sf)) + } +end diff --git a/math/surfit/iszero.x b/math/surfit/iszero.x new file mode 100644 index 00000000..d453d1fd --- /dev/null +++ b/math/surfit/iszero.x @@ -0,0 +1,26 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "surfitdef.h" + +# ISZERO -- Procedure to zero the accumulators for line number lineno. + +procedure iszero (sf) + +pointer sf # pointer to the surface descriptor + +begin + SF_NXPTS(sf) = 0 + SF_NYPTS(sf) = 0 + + call aclrr (XMATRIX(SF_XMATRIX(sf)), + SF_XORDER(sf) * SF_NXCOEFF(sf)) + call aclrr (YMATRIX(SF_YMATRIX(sf)), + SF_YORDER(sf) * SF_NYCOEFF(sf)) + call aclrr (XCOEFF(SF_XCOEFF(sf)), SF_NXCOEFF(sf) * SF_NLINES(sf)) + call aclrr (COEFF(SF_COEFF(sf)), SF_NXCOEFF(sf) * SF_NYCOEFF(sf)) + + if (SF_WZ(sf) != NULL) + call mfree (SF_WZ(sf), MEM_TYPE) + if (SF_TLEFT(sf) != NULL) + call mfree (SF_TLEFT(sf), TY_INT) +end diff --git a/math/surfit/mkpkg b/math/surfit/mkpkg new file mode 100644 index 00000000..58948e29 --- /dev/null +++ b/math/surfit/mkpkg @@ -0,0 +1,29 @@ +# Surface fitting tools library. + +$checkout libsurfit.a lib$ +$update libsurfit.a +$checkin libsurfit.a lib$ +$exit + +libsurfit.a: + iscoeff.x surfitdef.h + iseval.x surfitdef.h <math/surfit.h> + isfree.x surfitdef.h + isinit.x surfitdef.h <math/surfit.h> + islaccum.x surfitdef.h <math/surfit.h> + islfit.x surfitdef.h <math/surfit.h> + islrefit.x surfitdef.h <math/surfit.h> + islsolve.x surfitdef.h <math/surfit.h> + islzero.x surfitdef.h + isreplace.x surfitdef.h <math/surfit.h> + isresolve.x surfitdef.h <math/surfit.h> + issave.x surfitdef.h <math/surfit.h> + issolve.x surfitdef.h <math/surfit.h> + isvector.x surfitdef.h <math/surfit.h> + iszero.x surfitdef.h + sf_b1eval.x + sf_beval.x + sf_f1deval.x + sf_feval.x + sfchomat.x surfitdef.h <mach.h> <math/surfit.h> + ; diff --git a/math/surfit/sf_b1eval.x b/math/surfit/sf_b1eval.x new file mode 100644 index 00000000..d07006fc --- /dev/null +++ b/math/surfit/sf_b1eval.x @@ -0,0 +1,108 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +# SF_B1LEG -- Procedure to evaluate all the non-zero Legendrefunctions for +# a single point and given order. + +procedure sf_b1leg (x, order, k1, k2, basis) + +real x # array of data points +int order # order of polynomial, order = 1, constant +real k1, k2 # normalizing constants +real basis[ARB] # basis functions + +int i +real ri, xnorm + +begin + basis[1] = 1. + if (order == 1) + return + + xnorm = (x + k1) * k2 + basis[2] = xnorm + if (order == 2) + return + + do i = 3, order { + ri = i + basis[i] = ((2. * ri - 3.) * xnorm * basis[i-1] - + (ri - 2.) * basis[i-2]) / (ri - 1.) + } +end + + +# SF_B1CHEB -- Procedure to evaluate all the non zero Chebyshev function +# for a given x and order. + +procedure sf_b1cheb (x, order, k1, k2, basis) + +real x # number of data points +int order # order of polynomial, 1 is a constant +real k1, k2 # normalizing constants +real basis[ARB] # array of basis functions + +int i +real xnorm + +begin + basis[1] = 1. + if (order == 1) + return + + xnorm = (x + k1) * k2 + basis[2] = xnorm + if (order == 2) + return + + do i = 3, order + basis[i] = 2. * xnorm * basis[i-1] - basis[i-2] +end + + +# SF_B1SPLINE1 -- Evaluate all the non-zero spline1 functions for a +# single point. + +procedure sf_b1spline1 (x, npieces, k1, k2, basis, left) + +real x # set of data points +int npieces # number of polynomial pieces minus 1 +real k1, k2 # normalizing constants +real basis[ARB] # basis functions +int left # index of the appropriate spline functions + +real xnorm + +begin + xnorm = (x + k1) * k2 + left = min (int (xnorm), npieces) + + basis[2] = xnorm - left + basis[1] = 1. - basis[2] +end + + +# SF_B1SPLINE3 -- Procedure to evaluate all the non-zero basis functions +# for a cubic spline. + +procedure sf_b1spline3 (x, npieces, k1, k2, basis, left) + +real x # array of data points +int npieces # number of polynomial pieces +real k1, k2 # normalizing constants +real basis[ARB] # array of basis functions +int left # array of indices for first non-zero spline + +real sx, tx + +begin + sx = (x + k1) * k2 + left = min (int (sx), npieces) + + sx = sx - left + tx = 1. - sx + + basis[1] = tx * tx * tx + basis[2] = 1. + tx * (3. + tx * (3. - 3. * tx)) + basis[3] = 1. + sx * (3. + sx * (3. - 3. * sx)) + basis[4] = sx * sx * sx +end diff --git a/math/surfit/sf_beval.x b/math/surfit/sf_beval.x new file mode 100644 index 00000000..301967f9 --- /dev/null +++ b/math/surfit/sf_beval.x @@ -0,0 +1,143 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +# SF_BCHEB -- Procedure to evaluate all the non-zero Chebyshev functions for +# a set of points and given order. + +procedure sf_bcheb (x, npts, order, k1, k2, basis) + +real x[npts] # array of data points +int npts # number of points +int order # order of polynomial, order = 1, constant +real k1, k2 # normalizing constants +real basis[ARB] # basis functions + +int k, bptr + +begin + bptr = 1 + do k = 1, order { + + if (k == 1) + call amovkr (1., basis, npts) + else if (k == 2) + call altar (x, basis[bptr], npts, k1, k2) + else { + call amulr (basis[1+npts], basis[bptr-npts], basis[bptr], + npts) + call amulkr (basis[bptr], 2., basis[bptr], npts) + call asubr (basis[bptr], basis[bptr-2*npts], basis[bptr], npts) + } + + bptr = bptr + npts + } +end + + +# SF_BLEG -- Procedure to evaluate all the non zero Legendre function +# for a given order and set of points. + +procedure sf_bleg (x, npts, order, k1, k2, basis) + +real x[npts] # number of data points +int npts # number of points +int order # order of polynomial, 1 is a constant +real k1, k2 # normalizing constants +real basis[ARB] # array of basis functions + +int k, bptr +real ri, ri1, ri2 + +begin + bptr = 1 + + do k = 1, order { + if (k == 1) + call amovkr (1., basis, npts) + else if (k == 2) + call altar (x, basis[bptr], npts, k1, k2) + else { + ri = k + ri1 = (2. * ri - 3.) / (ri - 1.) + ri2 = - (ri - 2.) / (ri - 1.) + call amulr (basis[1+npts], basis[bptr-npts], basis[bptr], + npts) + call awsur (basis[bptr], basis[bptr-2*npts], + basis[bptr], npts, ri1, ri2) + } + + bptr = bptr + npts + } +end + + +# SF_BSPLINE1 -- Evaluate all the non-zero spline1 functions for a set +# of points. + +procedure sf_bspline1 (x, npts, npieces, k1, k2, basis, left) + +real x[npts] # set of data points +int npts # number of points +int npieces # number of polynomial pieces minus 1 +real k1, k2 # normalizing constants +real basis[ARB] # basis functions +int left[ARB] # indices of the appropriate spline functions + +int k + +begin + call altar (x, basis[1+npts], npts, k1, k2) + call achtri (basis[1+npts], left, npts) + call aminki (left, npieces, left, npts) + + do k = 1, npts { + basis[npts+k] = basis[npts+k] - left[k] + basis[k] = 1. - basis[npts+k] + } +end + + +# SF_BSPLINE3 -- Procedure to evaluate all the non-zero basis functions +# for a cubic spline. + +procedure sf_bspline3 (x, npts, npieces, k1, k2, basis, left) + +real x[npts] # array of data points +int npts # number of data points +int npieces # number of polynomial pieces minus 1 +real k1, k2 # normalizing constants +real basis[ARB] # array of basis functions +int left[ARB] # array of indices for first non-zero spline + +int i +pointer sp, sx, tx + +begin + # allocate space + call smark (sp) + call salloc (sx, npts, TY_REAL) + call salloc (tx, npts, TY_REAL) + + # calculate the index of the first non-zero coeff + call altar (x, Memr[sx], npts, k1, k2) + call achtri (Memr[sx], left, npts) + call aminki (left, npieces, left, npts) + + # normalize x to 0 to 1 + do i = 1, npts { + Memr[sx+i-1] = Memr[sx+i-1] - left[i] + Memr[tx+i-1] = 1. - Memr[sx+i-1] + } + + # calculate the basis function + call apowkr (Memr[tx], 3, basis, npts) + do i = 1, npts { + basis[npts+i] = 1. + Memr[tx+i-1] * (3. + Memr[tx+i-1] * (3. - + 3. * Memr[tx+i-1])) + basis[2*npts+i] = 1. + Memr[sx+i-1] * (3. + Memr[sx+i-1] * (3. - + 3. * Memr[sx+i-1])) + } + call apowkr (Memr[sx], 3, basis[1+3*npts], npts) + + # release space + call sfree (sp) +end diff --git a/math/surfit/sf_f1deval.x b/math/surfit/sf_f1deval.x new file mode 100644 index 00000000..01cb98d0 --- /dev/null +++ b/math/surfit/sf_f1deval.x @@ -0,0 +1,233 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +# CV_EVCHEB -- Procedure to evaluate a Chebyshev polynomial assuming that +# the coefficients have been calculated. + +procedure cv_evcheb (coeff, x, yfit, npts, order, k1, k2) + +real coeff[ARB] # EV array of coefficients +real x[npts] # x values of points to be evaluated +real yfit[npts] # the fitted points +int npts # number of points to be evaluated +int order # order of the polynomial, 1 = constant +real k1, k2 # normalizing constants + +int i +pointer sx, pn, pnm1, pnm2 +pointer sp +real c1, c2 + +begin + # fit a constant + call amovkr (coeff[1], yfit, npts) + if (order == 1) + return + + # fit a linear function + c1 = k2 * coeff[2] + c2 = c1 * k1 + coeff[1] + call altmr (x, yfit, npts, c1, c2) + if (order == 2) + return + + # allocate temporary space + call smark (sp) + call salloc (sx, npts, TY_REAL) + call salloc (pn, npts, TY_REAL) + call salloc (pnm1, npts, TY_REAL) + call salloc (pnm2, npts, TY_REAL) + + # a higher order polynomial + call amovkr (1., Memr[pnm2], npts) + call altar (x, Memr[sx], npts, k1, k2) + call amovr (Memr[sx], Memr[pnm1], npts) + call amulkr (Memr[sx], 2., Memr[sx], npts) + do i = 3, order { + call amulr (Memr[sx], Memr[pnm1], Memr[pn], npts) + call asubr (Memr[pn], Memr[pnm2], Memr[pn], npts) + if (i < order) { + call amovr (Memr[pnm1], Memr[pnm2], npts) + call amovr (Memr[pn], Memr[pnm1], npts) + } + call amulkr (Memr[pn], coeff[i], Memr[pn], npts) + call aaddr (yfit, Memr[pn], yfit, npts) + } + + # free temporary space + call sfree (sp) + +end + + +# CV_EVLEG -- Procedure to evaluate a Legendre polynomial assuming that +# the coefficients have been calculated. + +procedure cv_evleg (coeff, x, yfit, npts, order, k1, k2) + +real coeff[ARB] # EV array of coefficients +real x[npts] # x values of points to be evaluated +real yfit[npts] # the fitted points +int npts # number of data points +int order # order of the polynomial, 1 = constant +real k1, k2 # normalizing constants + +int i +pointer sx, pn, pnm1, pnm2 +pointer sp +real ri, ri1, ri2 + +begin + # fit a constant + call amovkr (coeff[1], yfit, npts) + if (order == 1) + return + + # fit a linear function + ri1 = k2 * coeff[2] + ri2 = ri1 * k1 + coeff[1] + call altmr (x, yfit, npts, ri1, ri2) + if (order == 2) + return + + # allocate temporary space + call smark (sp) + call salloc (sx, npts, TY_REAL) + call salloc (pn, npts, TY_REAL) + call salloc (pnm1, npts, TY_REAL) + call salloc (pnm2, npts, TY_REAL) + + # a higher order polynomial + call amovkr (1., Memr[pnm2], npts) + call altar (x, Memr[sx], npts, k1, k2) + call amovr (Memr[sx], Memr[pnm1], npts) + do i = 3, order { + ri = i + ri1 = (2. * ri - 3.) / (ri - 1.) + ri2 = - (ri - 2.) / (ri - 1.) + call amulr (Memr[sx], Memr[pnm1], Memr[pn], npts) + call awsur (Memr[pn], Memr[pnm2], Memr[pn], npts, ri1, ri2) + if (i < order) { + call amovr (Memr[pnm1], Memr[pnm2], npts) + call amovr (Memr[pn], Memr[pnm1], npts) + } + call amulkr (Memr[pn], coeff[i], Memr[pn], npts) + call aaddr (yfit, Memr[pn], yfit, npts) + } + + # free temporary space + call sfree (sp) + +end + + +# CV_EVSPLINE1 -- Procedure to evaluate a piecewise linear spline function +# assuming that the coefficients have been calculated. + +procedure cv_evspline1 (coeff, x, yfit, npts, npieces, k1, k2) + +real coeff[ARB] # array of coefficients +real x[npts] # array of x values +real yfit[npts] # array of fitted values +int npts # number of data points +int npieces # number of fitted points minus 1 +real k1, k2 # normalizing constants + +int j +pointer sx, tx, index +pointer sp + +begin + # allocate the required space + call smark (sp) + call salloc (sx, npts, TY_REAL) + call salloc (tx, npts, TY_REAL) + call salloc (index, npts, TY_INT) + + # calculate the index of the first non-zero coefficient + # for each point + call altar (x, Memr[sx], npts, k1, k2) + call achtri (Memr[sx], Memi[index], npts) + call aminki (Memi[index], npieces, Memi[index], npts) + + # transform sx to range 0 to 1 + do j = 1, npts { + Memr[sx+j-1] = Memr[sx+j-1] - Memi[index+j-1] + Memr[tx+j-1] = 1. - Memr[sx+j-1] + } + + # calculate yfit using the two non-zero basis function + call aclrr (yfit, npts) + do j = 1, npts + yfit[j] = Memr[tx+j-1] * coeff[1+Memi[index+j-1]] + + Memr[sx+j-1] * coeff[2+Memi[index+j-1]] + + # free space + call sfree (sp) + +end + + +# CV_EVSPLINE3 -- Procedure to evaluate the cubic spline assuming that +# the coefficients of the fit are known. + +procedure cv_evspline3 (coeff, x, yfit, npts, npieces, k1, k2) + +real coeff[ARB] # array of coeffcients +real x[npts] # array of x values +real yfit[npts] # array of fitted values +int npts # number of data points +int npieces # number of polynomial pieces +real k1, k2 # normalizing constants + +int i, j +pointer sx, tx, temp, index, sp + +begin + + # allocate the required space + call smark (sp) + call salloc (sx, npts, TY_REAL) + call salloc (tx, npts, TY_REAL) + call salloc (temp, npts, TY_REAL) + call salloc (index, npts, TY_INT) + + # calculate to which coefficients the x values contribute to + call altar (x, Memr[sx], npts, k1, k2) + call achtri (Memr[sx], Memi[index], npts) + call aminki (Memi[index], npieces, Memi[index], npts) + + # transform sx to range 0 to 1 + do j = 1, npts { + Memr[sx+j-1] = Memr[sx+j-1] - Memi[index+j-1] + Memr[tx+j-1] = 1. - Memr[sx+j-1] + } + + # calculate yfit using the four non-zero basis function + call aclrr (yfit, npts) + do i = 1, 4 { + + switch (i) { + case 1: + call apowkr (Memr[tx], 3, Memr[temp], npts) + case 2: + do j = 1, npts { + Memr[temp+j-1] = 1. + Memr[tx+j-1] * (3. + Memr[tx+j-1] * + (3. - 3. * Memr[tx+j-1])) + } + case 3: + do j = 1, npts { + Memr[temp+j-1] = 1. + Memr[sx+j-1] * (3. + Memr[sx+j-1] * + (3. - 3. * Memr[sx+j-1])) + } + case 4: + call apowkr (Memr[sx], 3, Memr[temp], npts) + } + + do j = 1, npts + Memr[temp+j-1] = Memr[temp+j-1] * coeff[i+Memi[index+j-1]] + call aaddr (yfit, Memr[temp], yfit, npts) + } + + # free space + call sfree (sp) +end 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 diff --git a/math/surfit/sfchomat.x b/math/surfit/sfchomat.x new file mode 100644 index 00000000..419fe777 --- /dev/null +++ b/math/surfit/sfchomat.x @@ -0,0 +1,105 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <math/surfit.h> +include "surfitdef.h" +include <mach.h> + +# SFCHOFAC -- Routine to calculate the Cholesky factorization of a +# symmetric, positive semi-definite banded matrix. This routines was +# adapted from the bchfac.f routine described in "A Practical Guide +# to Splines", Carl de Boor (1978). + +procedure sfchofac (matrix, nbands, nrows, matfac, ier) + +VAR_TYPE matrix[nbands, nrows] # data matrix +int nbands # number of bands +int nrows # number of rows +VAR_TYPE matfac[nbands, nrows] # Cholesky factorization +int ier # error code + +int i, n, j, imax, jmax +VAR_TYPE ratio + +begin + if (nrows == 1) { + if (matrix[1,1] > 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]) <= DELTA) { + do j = 1, nbands + matfac[j,n] = 0. + ier = SINGULAR + 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 + + +# SFCHOSLV -- Solve the matrix whose Cholesky factorization was calculated in +# SFCHOFAC for the coefficients. This routine was adapted from bchslv.f +# described in "A Practical Guide to Splines", by Carl de Boor (1978). + +procedure sfchoslv (matfac, nbands, nrows, vector, coeff) + +VAR_TYPE matfac[nbands,nrows] # Cholesky factorization +int nbands # number of bands +int nrows # number of rows +VAR_TYPE vector[nrows] # right side of matrix equation +VAR_TYPE coeff[nrows] # coefficients + +int i, n, j, jmax, nbndm1 + +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) { + do j = 1, jmax + coeff[j+n] = coeff[j+n] - matfac[j+1,n] * coeff[n] + } + } + + # back substitution + for (n = nrows; n >= 1; n = n - 1) { + coeff[n] = coeff[n] * matfac[1,n] + jmax = min (nbndm1, nrows - n) + if (jmax >= 1) { + do j = 1, jmax + coeff[n] = coeff[n] - matfac[j+1,n] * coeff[j+n] + } + } +end diff --git a/math/surfit/surfitdef.h b/math/surfit/surfitdef.h new file mode 100644 index 00000000..846109ea --- /dev/null +++ b/math/surfit/surfitdef.h @@ -0,0 +1,74 @@ +# set up the curve fitting structure + +define LEN_SFSTRUCT 35 + +define SF_TYPE Memi[$1] # Type of curve to be fitted +define SF_NXCOEFF Memi[$1+1] # Number of coefficients +define SF_XORDER Memi[$1+2] # Order of the fit in x +define SF_NXPIECES Memi[$1+3] # Number of x polynomial pieces - 1 +define SF_NCOLS Memi[$1+4] # Maximum x value +define SF_NXPTS Memi[$1+5] # Number of points in x +define SF_NYCOEFF Memi[$1+6] # Number of y coefficients +define SF_YORDER Memi[$1+7] # Order of the fit in y +define SF_NYPIECES Memi[$1+8] # Number of y polynomial pieces - 1 +define SF_NLINES Memi[$1+9] # Minimum x value +define SF_NYPTS Memi[$1+10] # Number of y points +define SF_XTERMS Memi[$1+11] # cross terms? + +define SF_XBASIS Memi[$1+12] # Pointer to the x basis functions +define SF_XLEFT Memi[$1+13] # Indices to x basis functions, spline +define SF_YBASIS Memi[$1+14] # Pointer to the y basis functions +define SF_YLEFT Memi[$1+15] # Indices to y basis functions, spline +define SF_XMATRIX Memi[$1+16] # Pointer to x data matrix +define SF_YMATRIX Memi[$1+17] # Pointer to y data matrix +define SF_XCOEFF Memi[$1+18] # X coefficient matrix +define SF_COEFF Memi[$1+19] # Pointer to coefficient vector + +define SF_XMIN Memr[P2R($1+20)] # Min x value +define SF_XMAX Memr[P2R($1+21)] # Max x value +define SF_XRANGE Memr[P2R($1+22)] # 2. / (xmax - xmin), polynomials +define SF_XMAXMIN Memr[P2R($1+23)] # - (xmax + xmin) / 2., polynomials +define SF_XSPACING Memr[P2R($1+24)] # order / (xmax - xmin), splines +define SF_YMIN Memr[P2R($1+25)] # Min y value +define SF_YMAX Memr[P2R($1+26)] # Max y value +define SF_YRANGE Memr[P2R($1+27)] # 2. / (ymax - ymin), polynomials +define SF_YMAXMIN Memr[P2R($1+28)] # - (ymax + ymin) / 2., polynomials +define SF_YSPACING Memr[P2R($1+29)] # order / (ymax - ymin), splines + +define SF_WZ Memi[$1+30] +define SF_TLEFT Memi[$1+31] + +# matrix and vector element definitions + +define XBASIS Memr[P2P($1)] # +define XBS Memr[P2P($1)] # +define YBASIS Memr[P2P($1)] # +define YBS Memr[P2P($1)] # +define XMATRIX Memr[P2P($1)] # +define XCHOFAC Memr[P2P($1)] # +define YMATRIX Memr[P2P($1)] # +define XCOEFF Memr[P2P($1)] # +define COEFF Memr[P2P($1)] # +define XLEFT Memi[$1] # +define YLEFT Memi[$1] # + +# structure definitions for the save restore functions + +define SF_SAVETYPE $1[1] +define SF_SAVEXORDER $1[2] +define SF_SAVEYORDER $1[3] +define SF_SAVEXTERMS $1[4] +define SF_SAVENCOLS $1[5] +define SF_SAVENLINES $1[6] +define SF_SAVECOEFF 6 + +# miscellaneous + +define SPLINE3_ORDER 4 +define SPLINE1_ORDER 2 + +# data type + +define MEM_TYPE TY_REAL +define VAR_TYPE real +define DELTA EPSILON |