diff options
Diffstat (limited to 'math/curfit/cv_bevald.x')
| -rw-r--r-- | math/curfit/cv_bevald.x | 147 |
1 files changed, 147 insertions, 0 deletions
diff --git a/math/curfit/cv_bevald.x b/math/curfit/cv_bevald.x new file mode 100644 index 00000000..7d7f6e44 --- /dev/null +++ b/math/curfit/cv_bevald.x @@ -0,0 +1,147 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +# CV_BCHEB -- Procedure to evaluate all the non-zero Chebyshev functions for +# a set of points and given order. + +procedure dcv_bcheb (x, npts, order, k1, k2, basis) + +double x[npts] # array of data points +int npts # number of points +int order # order of polynomial, order = 1, constant +double k1, k2 # normalizing constants +double basis[ARB] # basis functions + +int k, bptr + +begin + bptr = 1 + do k = 1, order { + if (k == 1) + call amovkd (double(1.0), basis, npts) + else if (k == 2) + call altad (x, basis[bptr], npts, k1, k2) + else { + call amuld (basis[1+npts], basis[bptr-npts], basis[bptr], + npts) + call amulkd (basis[bptr], double(2.0), basis[bptr], npts) + call asubd (basis[bptr], basis[bptr-2*npts], basis[bptr], npts) + } + bptr = bptr + npts + } +end + + +# CV_BLEG -- Procedure to evaluate all the non zero Legendre function +# for a given order and set of points. + +procedure dcv_bleg (x, npts, order, k1, k2, basis) + +double x[npts] # number of data points +int npts # number of points +int order # order of polynomial, 1 is a constant +double k1, k2 # normalizing constants +double basis[ARB] # array of basis functions + +int k, bptr +double ri, ri1, ri2 + +begin + bptr = 1 + do k = 1, order { + if (k == 1) + call amovkd (double(1.0), basis, npts) + else if (k == 2) + call altad (x, basis[bptr], npts, k1, k2) + else { + ri = k + ri1 = (double(2.0) * ri - double(3.0)) / (ri - double(1.0)) + ri2 = - (ri - double(2.0)) / (ri - double(1.0)) + call amuld (basis[1+npts], basis[bptr-npts], basis[bptr], + npts) + call awsud (basis[bptr], basis[bptr-2*npts], + basis[bptr], npts, ri1, ri2) + } + bptr = bptr + npts + } +end + + +# CV_BSPLINE1 -- Evaluate all the non-zero spline1 functions for a set +# of points. + +procedure dcv_bspline1 (x, npts, npieces, k1, k2, basis, left) + +double x[npts] # set of data points +int npts # number of points +int npieces # number of polynomial pieces minus 1 +double k1, k2 # normalizing constants +double basis[ARB] # basis functions +int left[ARB] # indices of the appropriate spline functions + +int k + +begin + call altad (x, basis[1+npts], npts, k1, k2) + call achtdi (basis[1+npts], left, npts) + call aminki (left, npieces, left, npts) + + do k = 1, npts { + basis[npts+k] = max (double(0.0), min (double(1.0), + basis[npts+k] - left[k])) + basis[k] = max (double(0.0), min (double(1.0), double(1.0) - + basis[npts+k])) + } +end + + +# CV_BSPLINE3 -- Procedure to evaluate all the non-zero basis functions +# for a cubic spline. + +procedure dcv_bspline3 (x, npts, npieces, k1, k2, basis, left) + +double x[npts] # array of data points +int npts # number of data points +int npieces # number of polynomial pieces minus 1 +double k1, k2 # normalizing constants +double basis[ARB] # array of basis functions +int left[ARB] # array of indices for first non-zero spline + +int i +pointer sp, sx, tx +double dsx, dtx + +begin + # allocate space + call smark (sp) + call salloc (sx, npts, TY_DOUBLE) + call salloc (tx, npts, TY_DOUBLE) + + # calculate the index of the first non-zero coeff + call altad (x, Memd[sx], npts, k1, k2) + call achtdi (Memd[sx], left, npts) + call aminki (left, npieces, left, npts) + + do i = 1, npts { + Memd[sx+i-1] = max (double(0.0), min (double(1.0), + Memd[sx+i-1] - left[i])) + Memd[tx+i-1] = max (double(0.0), min (double(1.0), double(1.0) - + Memd[sx+i-1])) + } + + # calculate the basis function + #call apowk$t (Mem$t[tx], 3, basis, npts) + do i = 1, npts { + dsx = Memd[sx+i-1] + dtx = Memd[tx+i-1] + basis[i] = dtx * dtx * dtx + basis[npts+i] = double(1.0) + dtx * (double(3.0) + dtx * + (double(3.0) - double(3.0) * dtx)) + basis[2*npts+i] = double(1.0) + dsx * (double(3.0) + dsx * + (double(3.0) - double(3.0) * dsx)) + basis[3*npts+i] = dsx * dsx * dsx + } + #call apowk$t (Mem$t[sx], 3, basis[1+3*npts], npts) + + # release space + call sfree (sp) +end |