aboutsummaryrefslogtreecommitdiff
path: root/math/gsurfit/zzdebug.x
diff options
context:
space:
mode:
Diffstat (limited to 'math/gsurfit/zzdebug.x')
-rw-r--r--math/gsurfit/zzdebug.x348
1 files changed, 348 insertions, 0 deletions
diff --git a/math/gsurfit/zzdebug.x b/math/gsurfit/zzdebug.x
new file mode 100644
index 00000000..522da747
--- /dev/null
+++ b/math/gsurfit/zzdebug.x
@@ -0,0 +1,348 @@
+task test = t_test
+
+include <math/gsurfit.h>
+
+procedure t_test()
+
+int i, j, k
+int xorder, yorder, xterms, stype
+int ncoeff, maxorder, xincr, npts, stype1, ier
+pointer gs, ags, sgs
+double dx, dy, const, accum, sum, rms1, rms2
+double x[121], y[121], z[121], w[121], zfit[121], coeff[121], save[121]
+int clgeti()
+double dgseval(), dgsgcoeff()
+
+begin
+ # Generate x and y grid.
+ dy = -1.0d0
+ npts = 0
+ do i = 1, 9 {
+ dx = -1.0d0
+ do j = 1, 9 {
+ x[npts+1] = dx
+ y[npts+1] = dy
+ npts = npts + 1
+ dx = dx + 0.25d0
+ }
+ dy = dy + 0.25d0
+ }
+
+ stype = clgeti ("stype")
+ xorder = clgeti ("xorder")
+ yorder = clgeti ("yorder")
+ xterms = clgeti ("xterms")
+ call printf ("\n\nSURFACE: %d XORDER: %d YORDER: %d XTERMS: %d\n")
+ call pargi (stype)
+ call pargi (xorder)
+ call pargi (yorder)
+ call pargi (xterms)
+
+ # Generate data
+ if (stype > 3) {
+ switch (xterms) {
+ case GS_XNONE:
+
+ do i = 1, npts {
+ sum = 0.0d0
+ do j = 2, yorder
+ sum = sum + j * y[i] ** (j - 1)
+ do j = 2, xorder
+ sum = sum + j * x[i] ** (j - 1)
+ z[i] = sum
+ }
+
+ case GS_XHALF:
+
+ maxorder = max (xorder + 1, yorder + 1)
+ do i = 1, npts {
+ sum = 0.0d0
+ xincr = xorder
+ do j = 1, yorder {
+ const = j * y[i] ** (j - 1)
+ accum= 0.0d0
+ do k = 1, xincr {
+ if (j > 1 || k > 1)
+ accum = accum + k * x[i] ** (k - 1)
+ }
+ sum = sum + const * accum
+ if ((j + xorder + 1) > maxorder)
+ xincr = xincr - 1
+ }
+ z[i] = sum
+ }
+
+ case GS_XFULL:
+
+ do i = 1, npts {
+ sum = 0.0d0
+ do j = 1, yorder {
+ const = j * y[i] ** (j - 1)
+ accum = 0.0d0
+ do k = 1, xorder {
+ if (j > 1 || k > 1)
+ accum = accum + k * x[i] ** (k - 1)
+ }
+ sum = sum + const * accum
+ }
+ z[i] = sum
+ }
+ }
+
+ stype1 = stype - 3
+ } else {
+ switch (xterms) {
+ case GS_XNONE:
+
+ do i = 1, npts {
+ sum = 0.0d0
+ do j = 2, yorder
+ sum = sum + j * y[i] ** (j - 1)
+ do j = 1, xorder
+ sum = sum + j * x[i] ** (j - 1)
+ z[i] = sum
+ }
+
+ case GS_XHALF:
+
+ maxorder = max (xorder + 1, yorder + 1)
+ do i = 1, npts {
+ sum = 0.0d0
+ xincr = xorder
+ do j = 1, yorder {
+ const = j * y[i] ** (j - 1)
+ accum= 0.0d0
+ do k = 1, xincr {
+ accum = accum + k * x[i] ** (k - 1)
+ }
+ sum = sum + const * accum
+ if ((j + xorder + 1) > maxorder)
+ xincr = xincr - 1
+ }
+ z[i] = sum
+ }
+
+ case GS_XFULL:
+
+ do i = 1, npts {
+ sum = 0.0d0
+ do j = 1, yorder {
+ const = j * y[i] ** (j - 1)
+ accum = 0.0d0
+ do k = 1, xorder {
+ accum = accum + k * x[i] ** (k - 1)
+ }
+ sum = sum + const * accum
+ }
+ z[i] = sum
+ }
+ }
+
+ stype1 = stype
+ }
+
+ # Print out the data.
+ call printf ("\nXIN:\n")
+ do i = 1, npts {
+ call printf ("%6.3f ")
+ call pargd (x[i])
+ if (mod (i, 9) == 0)
+ call printf ("\n")
+ }
+ call printf ("\n")
+
+ call printf ("\nYIN:\n")
+ do i = 1, npts {
+ call printf ("%6.3f ")
+ call pargd (y[i])
+ if (mod (i, 9) == 0)
+ call printf ("\n")
+ }
+ call printf ("\n")
+
+ call printf ("\nZIN:\n")
+ do i = 1, npts {
+ call printf ("%6.3f ")
+ call pargd (z[i])
+ if (mod (i, 9) == 0)
+ call printf ("\n")
+ }
+ call printf ("\n")
+
+ # Fit surface.
+ call dgsinit (gs, stype1, xorder, yorder, xterms, -1.0d0, 1.0d0,
+ -1.0d0, 1.0d0)
+ if (stype > 3) {
+ call dgsset (gs, GSXREF, 0d0)
+ call dgsset (gs, GSYREF, 0d0)
+ call dgsset (gs, GSZREF, 0d0)
+ }
+ call dgsfit (gs, x, y, z, w, npts, WTS_UNIFORM, ier)
+ call printf ("\nFIT ERROR CODE: %d\n")
+ call pargi (ier)
+
+ # Evaluate the fit and its rms.
+ call dgsvector (gs, x, y, zfit, npts)
+ call printf ("\nZFIT:\n")
+ do i = 1, npts {
+ call printf ("%6.3f ")
+ call pargd (zfit[i])
+ if (mod (i, 9) == 0)
+ call printf ("\n")
+ }
+ call printf ("\n")
+ rms1 = 0.0d0
+ do i = 1, npts
+ rms1 = rms1 + (z[i] - zfit[i]) ** 2
+ rms1 = sqrt (rms1 / (npts - 1))
+ rms2 = 0.0d0
+ do i = 1, npts
+ rms2 = rms2 + (z[i] - dgseval (gs, x[i], y[i])) ** 2
+ rms2 = sqrt (rms2 / (npts - 1))
+ #call printf ("\nRMS: vector = %0.14g point = %0.14g\n\n")
+ call printf ("\nRMS: vector = %0.4f point = %0.4f\n\n")
+ call pargd (rms1)
+ call pargd (rms2)
+
+ # Print the coefficients.
+ call dgscoeff (gs, coeff, ncoeff)
+ call printf ("GSFIT coeff:\n")
+ call printf ("first %0.14g %0.14g\n")
+ call pargd (dgsgcoeff (gs, 1, 1))
+ call pargd (dgsgcoeff (gs, xorder, 1))
+ do i = 1, ncoeff {
+ call printf ("%d %0.14g\n")
+ call pargi (i)
+ call pargd (coeff[i])
+ }
+ call printf ("last %0.14g %0.14g\n")
+ call pargd (dgsgcoeff (gs, 1, yorder))
+ call pargd (dgsgcoeff (gs, xorder, yorder))
+ call printf ("\n")
+
+ call dgsfree (gs)
+ return
+
+ # Evaluate the first derivatives.
+ call dgsder (gs, x, y, zfit, npts, 1, 0)
+ call printf ("\nZDER: 1 0\n")
+ do i = 1, npts {
+ call printf ("%0.7g ")
+ call pargd (zfit[i])
+ if (mod (i, 9) == 0)
+ call printf ("\n")
+ }
+ call printf ("\n")
+
+ call dgsder (gs, x, y, zfit, npts, 0, 1)
+ call printf ("\nZDER: 0 1\n")
+ do i = 1, npts {
+ call printf ("%0.7g ")
+ call pargd (zfit[i])
+ if (mod (i, 9) == 0)
+ call printf ("\n")
+ }
+ call printf ("\n")
+
+ call dgsder (gs, x, y, zfit, npts, 1, 1)
+ call printf ("\nZDER: 1 1\n")
+ do i = 1, npts {
+ call printf ("%0.7g ")
+ call pargd (zfit[i])
+ if (mod (i, 9) == 0)
+ call printf ("\n")
+ }
+ call printf ("\n")
+
+ # Refit the surface point by point.
+ call dgszero (gs)
+ do i = 1, npts {
+ call dgsaccum (gs, x[i], y[i], z[i], w[i], WTS_UNIFORM)
+ }
+ if (stype > 3)
+ call dgssolve1 (gs, ier)
+ else
+ call dgssolve (gs, ier)
+ call printf ("\nACCUM FIT ERROR CODE: %d\n")
+ call pargi (ier)
+ call dgsrej (gs, x[1], y[1], z[1], w[1], WTS_UNIFORM)
+ call dgsrej (gs, x[npts], y[npts], z[npts], w[npts], WTS_UNIFORM)
+ call dgsaccum (gs, x[1], y[1], z[1], w[1], WTS_UNIFORM)
+ call dgsaccum (gs, x[npts], y[npts], z[npts], w[npts], WTS_UNIFORM)
+ call dgssolve (gs, ier)
+ call printf ("\nREJ FIT ERROR CODE: %d\n")
+ call pargi (ier)
+
+ call dgscoeff (gs, coeff, ncoeff)
+ call printf ("GSACCUM coeff:\n")
+ call printf ("first %0.14g %0.14g\n")
+ call pargd (dgsgcoeff (gs, 1, 1))
+ call pargd (dgsgcoeff (gs, xorder, 1))
+ do i = 1, ncoeff {
+ call printf ("%d %0.14g\n")
+ call pargi (i)
+ call pargd (coeff[i])
+ }
+ call printf ("last %0.14g %0.14g\n")
+ call pargd (dgsgcoeff (gs, 1, yorder))
+ call pargd (dgsgcoeff (gs, xorder, yorder))
+ call printf ("\n")
+
+ # Save and restore.
+ call dgssave (gs, save)
+ call dgsfree (gs)
+ call dgsrestore (gs, save)
+
+ call dgscoeff (gs, coeff, ncoeff)
+ call printf ("RESTORE coeff:\n")
+ call printf ("first %0.14g %0.14g\n")
+ call pargd (dgsgcoeff (gs, 1, 1))
+ call pargd (dgsgcoeff (gs, xorder, 1))
+ do i = 1, ncoeff {
+ call printf ("%d %0.14g\n")
+ call pargi (i)
+ call pargd (coeff[i])
+ }
+ call printf ("last %0.14g %0.14g\n")
+ call pargd (dgsgcoeff (gs, 1, yorder))
+ call pargd (dgsgcoeff (gs, xorder, yorder))
+ call printf ("\n")
+
+ # Add two surfaces.
+ call dgsadd (gs, gs, ags)
+ call dgscoeff (ags, coeff, ncoeff)
+ call printf ("GSADD coeff:\n")
+ call printf ("first %0.14g %0.14g\n")
+ call pargd (dgsgcoeff (ags, 1, 1))
+ call pargd (dgsgcoeff (ags, xorder, 1))
+ do i = 1, ncoeff {
+ call printf ("%d %0.14g\n")
+ call pargi (i)
+ call pargd (coeff[i])
+ }
+ call printf ("last %0.14g %0.14g\n")
+ call pargd (dgsgcoeff (ags, 1, yorder))
+ call pargd (dgsgcoeff (ags, xorder, yorder))
+ call printf ("\n")
+
+ # Subtract two surfaces.
+ call dgssub (gs, gs, sgs)
+ call dgscoeff (sgs, coeff, ncoeff)
+ call printf ("GSSUB coeff:\n")
+ call printf ("first %0.14g %0.14g\n")
+ call pargd (dgsgcoeff (sgs, 1, 1))
+ call pargd (dgsgcoeff (sgs, xorder, 1))
+ do i = 1, ncoeff {
+ call printf ("%d %0.14g\n")
+ call pargi (i)
+ call pargd (coeff[i])
+ }
+ call printf ("last %0.14g %0.14g\n")
+ call pargd (dgsgcoeff (sgs, 1, yorder))
+ call pargd (dgsgcoeff (sgs, xorder, yorder))
+ call printf ("\n")
+
+ call dgsfree (gs)
+ call dgsfree (ags)
+ call dgsfree (sgs)
+end
generated by cgit (git 2.34.1) at 2025-09-20 08:49:55 -0400