From 40e5a5811c6ffce9b0974e93cdd927cbcf60c157 Mon Sep 17 00:00:00 2001 From: Joe Hunkeler Date: Tue, 11 Aug 2015 16:51:37 -0400 Subject: Repatch (from linux) of OSX IRAF --- math/gsurfit/gsadd.gx | 181 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 181 insertions(+) create mode 100644 math/gsurfit/gsadd.gx (limited to 'math/gsurfit/gsadd.gx') diff --git a/math/gsurfit/gsadd.gx b/math/gsurfit/gsadd.gx new file mode 100644 index 00000000..516f1b1f --- /dev/null +++ b/math/gsurfit/gsadd.gx @@ -0,0 +1,181 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include +$if (datatype == r) +include "gsurfitdef.h" +$else +include "dgsurfitdef.h" +$endif + +# GSADD -- Procedure to add the fits from two surfaces together. The surfaces +# must be the same type and the fit must cover the same range of data in x +# and y. This is a special function. + +$if (datatype == r) +procedure gsadd (sf1, sf2, sf3) +$else +procedure dgsadd (sf1, sf2, sf3) +$endif + +pointer sf1 # pointer to the first surface +pointer sf2 # pointer to the second surface +pointer sf3 # pointer to the output surface + +int i, order, nmove1, nmove2, nmove3, maxorder1, maxorder2, maxorder3 +pointer ptr1, ptr2, ptr3 +bool fpequal$t() + +begin + # test for NULL surface + if (sf1 == NULL && sf2 == NULL) { + sf3 = NULL + return + } else if (sf1 == NULL) { + $if (datatype == r) + call gscopy (sf2, sf3) + $else + call dgscopy (sf2, sf3) + $endif + return + } else if (sf2 == NULL) { + $if (datatype == r) + call gscopy (sf1, sf3) + $else + call dgscopy (sf1, sf3) + $endif + return + } + + # test that function type is the same + if (GS_TYPE(sf1) != GS_TYPE(sf2)) + call error (0, "GSADD: Incompatable surface types.") + + # test that mins and maxs are the same + if (! fpequal$t (GS_XMIN(sf1), GS_XMIN(sf2))) + call error (0, "GSADD: X ranges not identical.") + if (! fpequal$t (GS_XMAX(sf1), GS_XMAX(sf2))) + call error (0, "GSADD: X ranges not identical.") + if (! fpequal$t (GS_YMIN(sf1), GS_YMIN(sf2))) + call error (0, "GSADD: Y ranges not identical.") + if (! fpequal$t (GS_YMAX(sf1), GS_YMAX(sf2))) + call error (0, "GSADD: Y ranges not identical.") + + # allocate space for the pointer + call calloc (sf3, LEN_GSSTRUCT, TY_STRUCT) + + # copy parameters + GS_TYPE(sf3) = GS_TYPE(sf1) + + switch (GS_TYPE(sf3)) { + case GS_LEGENDRE, GS_CHEBYSHEV, GS_POLYNOMIAL: + GS_NXCOEFF(sf3) = max (GS_NXCOEFF(sf1), GS_NXCOEFF(sf2)) + GS_XORDER(sf3) = max (GS_XORDER(sf1), GS_XORDER(sf2)) + GS_XMIN(sf3) = GS_XMIN(sf1) + GS_XMAX(sf3) = GS_XMAX(sf1) + GS_XRANGE(sf3) = GS_XRANGE(sf1) + GS_XMAXMIN(sf3) = GS_XMAXMIN(sf1) + GS_NYCOEFF(sf3) = max (GS_NYCOEFF(sf1), GS_NYCOEFF(sf2)) + GS_YORDER(sf3) = max (GS_YORDER(sf1), GS_YORDER(sf2)) + GS_YMIN(sf3) = GS_YMIN(sf1) + GS_YMAX(sf3) = GS_YMAX(sf1) + GS_YRANGE(sf3) = GS_YRANGE(sf1) + GS_YMAXMIN(sf3) = GS_YMAXMIN(sf1) + if (GS_XTERMS(sf1) == GS_XTERMS(sf2)) + GS_XTERMS(sf3) = GS_XTERMS(sf1) + else if (GS_XTERMS(sf1) == GS_XFULL || GS_XTERMS(sf2) == GS_XFULL) + GS_XTERMS(sf3) = GS_XFULL + else + GS_XTERMS(sf3) = GS_XHALF + switch (GS_XTERMS(sf3)) { + case GS_XNONE: + GS_NCOEFF(sf3) = GS_NXCOEFF(sf3) + GS_NYCOEFF(sf3) - 1 + case GS_XHALF: + order = min (GS_XORDER(sf3), GS_YORDER(sf3)) + GS_NCOEFF(sf3) = GS_NXCOEFF(sf3) * GS_NYCOEFF(sf3) - order * + (order - 1) / 2 + default: + GS_NCOEFF(sf3) = GS_NXCOEFF(sf3) * GS_NYCOEFF(sf3) + } + default: + call error (0, "GSADD: Unknown curve type.") + } + + # set pointers to NULL + GS_XBASIS(sf3) = NULL + GS_YBASIS(sf3) = NULL + GS_MATRIX(sf3) = NULL + GS_CHOFAC(sf3) = NULL + GS_VECTOR(sf3) = NULL + GS_COEFF(sf3) = NULL + GS_WZ(sf3) = NULL + + # calculate the coefficients + $if (datatype == r) + call calloc (GS_COEFF(sf3), GS_NCOEFF(sf3), TY_REAL) + $else + call calloc (GS_COEFF(sf3), GS_NCOEFF(sf3), TY_DOUBLE) + $endif + + # set up line counters. + maxorder1 = max (GS_XORDER(sf1) + 1, GS_YORDER(sf1) + 1) + maxorder2 = max (GS_XORDER(sf2) + 1, GS_YORDER(sf2) + 1) + maxorder3 = max (GS_XORDER(sf3) + 1, GS_YORDER(sf3) + 1) + + # add in the first surface. + ptr1 = GS_COEFF(sf1) + ptr3 = GS_COEFF(sf3) + nmove1 = GS_NXCOEFF(sf1) + nmove3 = GS_NXCOEFF(sf3) + do i = 1, GS_NYCOEFF(sf1) { + call amov$t (COEFF(ptr1), COEFF(ptr3), nmove1) + ptr1 = ptr1 + nmove1 + ptr3 = ptr3 + nmove3 + switch (GS_XTERMS(sf1)) { + case GS_XNONE: + nmove1 = 1 + case GS_XHALF: + if ((i + GS_XORDER(sf1) + 1) > maxorder1) + nmove1 = nmove1 - 1 + case GS_XFULL: + ; + } + switch (GS_XTERMS(sf3)) { + case GS_XNONE: + nmove3 = 1 + case GS_XHALF: + if ((i + GS_XORDER(sf3) + 1) > maxorder3) + nmove3 = nmove3 - 1 + case GS_XFULL: + ; + } + } + + # add in the second surface. + ptr2 = GS_COEFF(sf2) + ptr3 = GS_COEFF(sf3) + nmove2 = GS_NXCOEFF(sf2) + nmove3 = GS_NXCOEFF(sf3) + do i = 1, GS_NYCOEFF(sf2) { + call aadd$t (COEFF(ptr3), COEFF(ptr2), COEFF(ptr3), nmove2) + ptr2 = ptr2 + nmove2 + ptr3 = ptr3 + nmove3 + switch (GS_XTERMS(sf2)) { + case GS_XNONE: + nmove2 = 1 + case GS_XHALF: + if ((i + GS_XORDER(sf2) + 1) > maxorder2) + nmove2 = nmove2 - 1 + case GS_XFULL: + ; + } + switch (GS_XTERMS(sf3)) { + case GS_XNONE: + nmove3 = 1 + case GS_XHALF: + if ((i + GS_XORDER(sf3) + 1) > maxorder3) + nmove3 = nmove3 - 1 + case GS_XFULL: + ; + } + } +end -- cgit