From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- math/gsurfit/gsder.gx | 264 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 264 insertions(+) create mode 100644 math/gsurfit/gsder.gx (limited to 'math/gsurfit/gsder.gx') diff --git a/math/gsurfit/gsder.gx b/math/gsurfit/gsder.gx new file mode 100644 index 00000000..e0ee95bd --- /dev/null +++ b/math/gsurfit/gsder.gx @@ -0,0 +1,264 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include +$if (datatype == r) +include "gsurfitdef.h" +$else +include "dgsurfitdef.h" +$endif + +# GSDER -- Procedure to calculate a new surface which is a derivative of +# the previous surface + +$if (datatype == r) +procedure gsder (sf1, x, y, zfit, npts, nxd, nyd) +$else +procedure dgsder (sf1, x, y, zfit, npts, nxd, nyd) +$endif + +pointer sf1 # pointer to the previous surface +PIXEL x[npts] # x values +PIXEL y[npts] # y values +PIXEL zfit[npts] # fitted values +int npts # number of points +int nxd, nyd # order of the derivatives in x and y + +PIXEL norm +int ncoeff, nxder, nyder, i, j +int order, maxorder1, maxorder2, nmove1, nmove2 +pointer sf2, sp, coeff, ptr1, ptr2 + +begin + if (sf1 == NULL) + return + + if (nxd < 0 || nyd < 0) + call error (0, "GSDER: Order of derivatives cannot be < 0") + + if (nxd == 0 && nyd == 0) { + $if (datatype == r) + call gsvector (sf1, x, y, zfit, npts) + $else + call dgsvector (sf1, x, y, zfit, npts) + $endif + return + } + + # allocate space for new surface + call calloc (sf2, LEN_GSSTRUCT, TY_STRUCT) + + # check the order of the derivatives and return 0 if the order is + # high + nxder = min (nxd, GS_NXCOEFF(sf1)) + nyder = min (nyd, GS_NYCOEFF(sf1)) + if (nxder >= GS_NXCOEFF(sf1) && nyder >= GS_NYCOEFF(sf1)) + call amovk$t (PIXEL(0.0), zfit, npts) + + # set up new surface + GS_TYPE(sf2) = GS_TYPE(sf1) + + # set the derivative surface parameters + switch (GS_TYPE(sf2)) { + case GS_LEGENDRE, GS_CHEBYSHEV, GS_POLYNOMIAL: + + GS_XTERMS(sf2) = GS_XTERMS(sf1) + + # find the order of the new surface + switch (GS_XTERMS(sf2)) { + case GS_XNONE: + if (nxder > 0 && nyder > 0) { + GS_NXCOEFF(sf2) = 1 + GS_XORDER(sf2) = 1 + GS_NYCOEFF(sf2) = 1 + GS_YORDER(sf2) = 1 + GS_NCOEFF(sf2) = 1 + } else if (nxder > 0) { + GS_NXCOEFF(sf2) = max (1, GS_NXCOEFF(sf1) - nxder) + GS_XORDER(sf2) = max (1, GS_NXCOEFF(sf1) - nxder) + GS_NYCOEFF(sf2) = 1 + GS_YORDER(sf2) = 1 + GS_NCOEFF(sf2) = GS_NXCOEFF(sf2) + } else if (nyder > 0) { + GS_NXCOEFF(sf2) = 1 + GS_XORDER(sf2) = 1 + GS_NYCOEFF(sf2) = max (1, GS_NYCOEFF(sf1) - nyder) + GS_YORDER(sf2) = max (1, GS_NYCOEFF(sf1) - nyder) + GS_NCOEFF(sf2) = GS_NYCOEFF(sf2) + } + + case GS_XHALF: + if ((nxder >= GS_NXCOEFF(sf1)) || (nyder >= GS_NYCOEFF(sf1)) || + (nxder + nyder) >= max (GS_NXCOEFF(sf1), + GS_NYCOEFF(sf1))) { + GS_NXCOEFF(sf2) = 1 + GS_XORDER(sf2) = 1 + GS_NYCOEFF(sf2) = 1 + GS_YORDER(sf2) = 1 + GS_NCOEFF(sf2) = 1 + } else { + maxorder1 = max (GS_XORDER(sf1) + 1, GS_YORDER(sf1) + 1) + order = max (1, min (maxorder1 - 1 - nyder - nxder, + GS_NXCOEFF(sf1) - nxder)) + GS_NXCOEFF(sf2) = order + GS_XORDER(sf2) = order + order = max (1, min (maxorder1 - 1 - nyder - nxder, + GS_NYCOEFF(sf1) - nyder)) + GS_NYCOEFF(sf2) = order + GS_YORDER(sf2) = order + order = min (GS_XORDER(sf2), GS_YORDER(sf2)) + GS_NCOEFF(sf2) = GS_NXCOEFF(sf2) * GS_NYCOEFF(sf2) - + order * (order - 1) / 2 + } + + default: + if (nxder >= GS_NXCOEFF(sf1) || nyder >= GS_NYCOEFF(sf1)) { + GS_NXCOEFF(sf2) = 1 + GS_XORDER(sf2) = 1 + GS_NYCOEFF(sf2) = 1 + GS_YORDER(sf2) = 1 + GS_NCOEFF(sf2) = 1 + } else { + GS_NXCOEFF(sf2) = max (1, GS_NXCOEFF(sf1) - nxder) + GS_XORDER(sf2) = max (1, GS_XORDER(sf1) - nxder) + GS_NYCOEFF(sf2) = max (1, GS_NYCOEFF(sf1) - nyder) + GS_YORDER(sf2) = max (1, GS_YORDER(sf1) - nyder) + GS_NCOEFF(sf2) = GS_NXCOEFF(sf2) * GS_NYCOEFF(sf2) + } + } + + # define the data limits + GS_XMIN(sf2) = GS_XMIN(sf1) + GS_XMAX(sf2) = GS_XMAX(sf1) + GS_XRANGE(sf2) = GS_XRANGE(sf1) + GS_XMAXMIN(sf2) = GS_XMAXMIN(sf1) + GS_YMIN(sf2) = GS_YMIN(sf1) + GS_YMAX(sf2) = GS_YMAX(sf1) + GS_YRANGE(sf2) = GS_YRANGE(sf1) + GS_YMAXMIN(sf2) = GS_YMAXMIN(sf1) + + default: + call error (0, "GSDER: Unknown surface type.") + } + + # set remaining surface pointers to NULL + GS_XBASIS(sf2) = NULL + GS_YBASIS(sf2) = NULL + GS_MATRIX(sf2) = NULL + GS_CHOFAC(sf2) = NULL + GS_VECTOR(sf2) = NULL + GS_COEFF(sf2) = NULL + GS_WZ(sf2) = NULL + + # allocate space for coefficients + call calloc (GS_COEFF(sf2), GS_NCOEFF(sf2), TY_PIXEL) + + # get coefficients + call smark (sp) + call salloc (coeff, GS_NCOEFF(sf1), TY_PIXEL) + $if (datatype == r) + call gscoeff (sf1, Mem$t[coeff], ncoeff) + $else + call dgscoeff (sf1, Mem$t[coeff], ncoeff) + $endif + + # compute the new coefficients + switch (GS_XTERMS(sf2)) { + case GS_XFULL: + if (nxder >= GS_NXCOEFF(sf1) || nyder >= GS_NYCOEFF(sf1)) + COEFF(GS_COEFF(sf2)) = 0. + else { + ptr2 = GS_COEFF(sf2) + (GS_NYCOEFF(sf2) - 1) * GS_NXCOEFF(sf2) + ptr1 = coeff + (GS_NYCOEFF(sf1) - 1) * GS_NXCOEFF(sf1) + do i = GS_NYCOEFF(sf1), nyder + 1, -1 { + call amov$t (Mem$t[ptr1+nxder], COEFF(ptr2), + GS_NXCOEFF(sf2)) + ptr2 = ptr2 - GS_NXCOEFF(sf2) + ptr1 = ptr1 - GS_NXCOEFF(sf1) + } + } + + case GS_XHALF: + if ((nxder >= GS_NXCOEFF(sf1)) || (nyder >= GS_NYCOEFF(sf1)) || + (nxder + nyder) >= max (GS_NXCOEFF(sf1), GS_NYCOEFF(sf1))) + COEFF(GS_COEFF(sf2)) = 0. + else { + maxorder1 = max (GS_XORDER(sf1) + 1, GS_YORDER(sf1) + 1) + maxorder2 = max (GS_XORDER(sf2) + 1, GS_YORDER(sf2) + 1) + ptr2 = GS_COEFF(sf2) + GS_NCOEFF(sf2) + ptr1 = coeff + GS_NCOEFF(sf1) + do i = GS_NYCOEFF(sf1), nyder + 1, -1 { + nmove1 = max (0, min (maxorder1 - i, GS_NXCOEFF(sf1))) + nmove2 = max (0, min (maxorder2 - i + nyder, + GS_NXCOEFF(sf2))) + ptr1 = ptr1 - nmove1 + ptr2 = ptr2 - nmove2 + call amov$t (Mem$t[ptr1+nxder], COEFF(ptr2), nmove2) + } + } + + default: + if (nxder > 0 && nyder > 0) + COEFF(GS_COEFF(sf2)) = 0. + else if (nxder > 0) { + if (nxder >= GS_NXCOEFF(sf1)) + COEFF(GS_COEFF(sf2)) = 0. + else { + ptr1 = coeff + ptr2 = GS_COEFF(sf2) + GS_NCOEFF(sf2) - 1 + do j = GS_NXCOEFF(sf1), nxder + 1, -1 { + COEFF(ptr2) = Mem$t[ptr1+j-1] + ptr2 = ptr2 - 1 + } + } + } else if (nyder > 0) { + if (nyder >= GS_NYCOEFF(sf1)) + COEFF(GS_COEFF(sf2)) = 0. + else { + ptr1 = coeff + GS_NCOEFF(sf1) - 1 + ptr2 = GS_COEFF(sf2) + do i = GS_NYCOEFF(sf1), nyder + 1, -1 + ptr1 = ptr1 - 1 + call amov$t (Mem$t[ptr1+1], COEFF(ptr2), GS_NCOEFF(sf2)) + } + } + } + + # evaluate the derivatives + switch (GS_TYPE(sf2)) { + case GS_POLYNOMIAL: + call $tgs_derpoly (COEFF(GS_COEFF(sf2)), x, y, zfit, npts, + GS_XTERMS(sf2), GS_XORDER(sf2), GS_YORDER(sf2), nxder, + nyder, GS_XMAXMIN(sf2), GS_XRANGE(sf2), GS_YMAXMIN(sf2), + GS_YRANGE(sf2)) + + case GS_CHEBYSHEV: + call $tgs_dercheb (COEFF(GS_COEFF(sf2)), x, y, zfit, npts, + GS_XTERMS(sf2), GS_XORDER(sf2), GS_YORDER(sf2), nxder, + nyder, GS_XMAXMIN(sf2), GS_XRANGE(sf2), GS_YMAXMIN(sf2), + GS_YRANGE(sf2)) + + case GS_LEGENDRE: + call $tgs_derleg (COEFF(GS_COEFF(sf2)), x, y, zfit, npts, + GS_XTERMS(sf2), GS_XORDER(sf2), GS_YORDER(sf2), nxder, + nyder, GS_XMAXMIN(sf2), GS_XRANGE(sf2), GS_YMAXMIN(sf2), + GS_YRANGE(sf2)) + + default: + call error (0, "GSVECTOR: Unknown surface type.") + } + + # Normalize. + if (GS_TYPE(sf2) != GS_POLYNOMIAL) { + norm = (2. / (GS_XMAX(sf2) - GS_XMIN(sf2))) ** nxder * (2. / + (GS_YMAX(sf2) - GS_YMIN(sf2))) ** nyder + call amulk$t (zfit, norm, zfit, npts) + } + + # free the space + $if (datatype == r) + call gsfree (sf2) + $else + call dgsfree (sf2) + $endif + call sfree (sp) +end -- cgit