aboutsummaryrefslogtreecommitdiff
path: root/noao/twodspec/longslit/transform/igsfit
diff options
context:
space:
mode:
Diffstat (limited to 'noao/twodspec/longslit/transform/igsfit')
-rw-r--r--noao/twodspec/longslit/transform/igsfit/Revisions42
-rw-r--r--noao/twodspec/longslit/transform/igsfit/igscolon.x115
-rw-r--r--noao/twodspec/longslit/transform/igsfit/igsdelete.x103
-rw-r--r--noao/twodspec/longslit/transform/igsfit/igsfit.com10
-rw-r--r--noao/twodspec/longslit/transform/igsfit/igsfit.x373
-rw-r--r--noao/twodspec/longslit/transform/igsfit/igsget.x62
-rw-r--r--noao/twodspec/longslit/transform/igsfit/igsgraph.x73
-rw-r--r--noao/twodspec/longslit/transform/igsfit/igsinit.x21
-rw-r--r--noao/twodspec/longslit/transform/igsfit/igsnearest.x51
-rw-r--r--noao/twodspec/longslit/transform/igsfit/igsparams.x23
-rw-r--r--noao/twodspec/longslit/transform/igsfit/igsset.x59
-rw-r--r--noao/twodspec/longslit/transform/igsfit/igssolve.x173
-rw-r--r--noao/twodspec/longslit/transform/igsfit/igsundelete.x107
-rw-r--r--noao/twodspec/longslit/transform/igsfit/mkpkg21
-rw-r--r--noao/twodspec/longslit/transform/igsfit/xgs.x243
15 files changed, 1476 insertions, 0 deletions
diff --git a/noao/twodspec/longslit/transform/igsfit/Revisions b/noao/twodspec/longslit/transform/igsfit/Revisions
new file mode 100644
index 00000000..92b36cca
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/Revisions
@@ -0,0 +1,42 @@
+.help revisions Jun88 noao.twodspec.longslit.transform.igsfit
+.nf
+ igsfit.x
+ igsnearest.x
+ GSCUR was being called with DOUBLE precision values. (12/22/87)
+
+ igsfit.x
+ igscolon.x
+ igsget.x
+ Added colon options to print fit at corners of surface. (8/10/87 Valdes)
+
+ ====
+ V2.5
+ ====
+
+noao$twodspec/longslit/transform/igsfit/*.x
+ Valdes, February 17, 1987
+ 1. GIO changes.
+
+noao$twodspec/longslit/transform/igsfit/igsfit.x
+noao$twodspec/longslit/transform/igsfit/igscolon.x
+ Valdes, January 16, 1987
+ 1. '?' now uses system page facility.
+ 2. Colon command dictionary and switch modified to use macro definitions.
+
+noao$twodspec/longslit/transform/igsfit/igsdelete.x
+noao$twodspec/longslit/transform/igsfit/igsundelete.x
+ Valdes, October 16, 1986
+ 1. Real line type specified in gseti call changed to integer.
+ This caused a crash on AOS/IRAF.
+
+========================================================
+
+From Valdes on Feb 7, 1986:
+
+1. Bug fixed in deleting and undeleting points.
+------
+From Valdes on Jan 3, 1986:
+
+1. Modified IGSFIT to allow zooming on constant x, constant y, constant z,
+and constant feature.
+.endhelp
diff --git a/noao/twodspec/longslit/transform/igsfit/igscolon.x b/noao/twodspec/longslit/transform/igsfit/igscolon.x
new file mode 100644
index 00000000..6847974a
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/igscolon.x
@@ -0,0 +1,115 @@
+include <gset.h>
+
+# List of colon commands
+define CMDS "|show|function|xorder|yorder|corners|"
+
+define SHOW 1 # Show parameters
+define FUNCTION 2 # Set or show function type
+define XORDER 3 # Set or show x order of function
+define YORDER 4 # Set or show y order of function
+define CORNERS 5 # Show corners
+
+# IGS_COLON -- Processes colon commands.
+
+procedure igs_colon (cmdstr, gp, sf)
+
+char cmdstr[ARB] # Command string
+pointer gp # GIO pointer
+pointer sf # Surface pointer
+
+char cmd[SZ_LINE]
+int ncmd, ival
+
+int nscan(), strdic()
+real xgseval()
+
+string funcs "|chebyshev|legendre|"
+
+include "igsfit.com"
+
+begin
+ # Use formated scan to parse the command string.
+ # The first word is the command and it may be minimum match
+ # abbreviated with the list of commands.
+
+ call sscan (cmdstr)
+ call gargwrd (cmd, SZ_LINE)
+ ncmd = strdic (cmd, cmd, SZ_LINE, CMDS)
+
+ switch (ncmd) {
+ case SHOW: # :show - Show the values of the fitting parameters.
+ call gdeactivate (gp, AW_CLEAR)
+ call printf ("function %s\n")
+ call pargstr (function)
+ call printf ("xorder %d\n")
+ call pargi (xorder)
+ call printf ("yorder %d\n")
+ call pargi (yorder)
+ call printf ("Fitted coordinates at the corners of the images:\n")
+ call printf (" (%d, %d) = %g (%d, %d) = %g\n")
+ call pargr (xmin)
+ call pargr (ymin)
+ call pargr (xgseval (sf, xmin, ymin))
+ call pargr (xmax)
+ call pargr (ymin)
+ call pargr (xgseval (sf, xmax, xmin))
+ call printf (" (%d, %d) = %g (%d, %d) = %g\n")
+ call pargr (xmin)
+ call pargr (ymax)
+ call pargr (xgseval (sf, xmin, ymax))
+ call pargr (xmax)
+ call pargr (ymax)
+ call pargr (xgseval (sf, xmax, ymax))
+ call printf ("rms %g\n")
+ call pargr (rms)
+ call greactivate (gp, AW_PAUSE)
+
+ case FUNCTION: # :function - List or set the fitting function.
+ call gargwrd (cmd, SZ_LINE)
+ if (nscan() == 1) {
+ call printf ("function = %s\n")
+ call pargstr (function)
+ } else {
+ if (strdic (cmd, cmd, SZ_LINE, funcs) > 0)
+ call strcpy (cmd, function, SZ_LINE)
+ else
+ call printf ("Unknown or ambiguous function\n")
+ }
+
+ case XORDER: # xorder: List or set the function order.
+ call gargi (ival)
+ if (nscan() == 1) {
+ call printf ("xorder %d\n")
+ call pargi (xorder)
+ } else if (ival < 2)
+ call printf ("xorder must be at least 2\n")
+ else
+ xorder = ival
+
+ case YORDER: # yorder: List or set the function order.
+ call gargi (ival)
+ if (nscan() == 1) {
+ call printf ("yorder %d\n")
+ call pargi (yorder)
+ } else if (ival < 2)
+ call printf ("yorder must be at least 2\n")
+ else
+ yorder = ival
+ case CORNERS: # corners: List coordinates at corners.
+ call printf ("(%d,%d)=%g (%d,%d)=%g (%d,%d)=%g (%d,%d)=%g\n")
+ call pargr (xmin)
+ call pargr (ymin)
+ call pargr (xgseval (sf, xmin, ymin))
+ call pargr (xmax)
+ call pargr (ymin)
+ call pargr (xgseval (sf, xmax, xmin))
+ call pargr (xmin)
+ call pargr (ymax)
+ call pargr (xgseval (sf, xmin, ymax))
+ call pargr (xmax)
+ call pargr (ymax)
+ call pargr (xgseval (sf, xmax, ymax))
+ default:
+ call printf ("Unrecognized or ambiguous command\007")
+ }
+end
diff --git a/noao/twodspec/longslit/transform/igsfit/igsdelete.x b/noao/twodspec/longslit/transform/igsfit/igsdelete.x
new file mode 100644
index 00000000..3de2fb25
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/igsdelete.x
@@ -0,0 +1,103 @@
+include <mach.h>
+include <gset.h>
+include <pkg/gtools.h>
+include <pkg/igsfit.h>
+
+# IGS_NEARESTD -- Nearest point to delete.
+
+int procedure igs_nearestd (gp, ztype, refpt, axis, pts, npts, wx, wy, wcs)
+
+pointer gp # GIO pointer
+int ztype # Zoom type
+int refpt # Reference point
+int axis[2] # Axes
+real pts[npts, ARB] # Data points
+int npts # Number of data points
+real wx, wy # Cursor coordinates
+int wcs # WCS
+
+int i, j, x, y
+real r2, r2min, x0, y0
+
+begin
+ x = axis[1]
+ y = axis[2]
+
+ call gctran (gp, wx, wy, wx, wy, wcs, 0)
+ r2min = MAX_REAL
+ j = 0
+
+ if (IS_INDEFI (ztype)) {
+ do i = 1, npts {
+ if (pts[i,W] == 0.)
+ next
+ call gctran (gp, pts[i, x], pts[i, y], x0, y0, wcs, 0)
+ r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+ if (r2 < r2min) {
+ r2min = r2
+ j = i
+ }
+ }
+ } else {
+ do i = 1, npts {
+ if ((pts[i,ztype] != pts[refpt,ztype]) || (pts[i,W] == 0.))
+ next
+ call gctran (gp, pts[i, x], pts[i, y], x0, y0, wcs, 0)
+ r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+ if (r2 < r2min) {
+ r2min = r2
+ j = i
+ }
+ }
+ }
+
+ return (j)
+end
+
+# IGS_DELETE -- Delete points or subsets.
+
+procedure igs_delete (gp, gt, ztype, refpt, axis, pts, npts, dtype)
+
+pointer gp # GIO pointer
+pointer gt # GTOOLS pointer
+int ztype # Zoom type
+int refpt # Reference point for deletion
+int axis[2] # Axes
+real pts[npts, ARB] # Data points
+int npts # Number of data points
+int dtype # Deletion type
+
+int i, x, y
+real xsize, ysize
+
+real gt_getr()
+
+begin
+ x = axis[1]
+ y = axis[2]
+
+ xsize = gt_getr (gt, GTXSIZE)
+ ysize = gt_getr (gt, GTYSIZE)
+
+ switch (dtype) {
+ case X, Y, Z:
+ do i = 1, npts {
+ if (!IS_INDEFI (ztype))
+ if (pts[i,ztype] != pts[refpt,ztype])
+ next
+ if (pts[i,dtype] != pts[refpt,dtype])
+ next
+ call gseti (gp, G_PMLTYPE, 0)
+ call gmark (gp, pts[i,x], pts[i,y], GM_PLUS, xsize, ysize)
+ call gseti (gp, G_PMLTYPE, 1)
+ call gmark (gp, pts[i,x], pts[i,y], GM_CROSS, xsize, ysize)
+ pts[i,W] = 0.
+ }
+ default:
+ call gseti (gp, G_PMLTYPE, 0)
+ call gmark (gp, pts[refpt,x], pts[refpt,y], GM_PLUS, xsize, ysize)
+ call gseti (gp, G_PMLTYPE, 1)
+ call gmark (gp, pts[refpt,x], pts[refpt,y], GM_CROSS, xsize, ysize)
+ pts[refpt,W] = 0.
+ }
+end
diff --git a/noao/twodspec/longslit/transform/igsfit/igsfit.com b/noao/twodspec/longslit/transform/igsfit/igsfit.com
new file mode 100644
index 00000000..90bf90aa
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/igsfit.com
@@ -0,0 +1,10 @@
+# Common parameters.
+
+char function[SZ_LINE] # Surface function
+int xorder # X order of surface function
+int yorder # Y order of surface function
+real xmin, xmax # X range
+real ymin, ymax # Y range
+real mean, rms # Mean and RMS of fit
+
+common /igscom/ xmin, xmax, ymin, ymax, xorder, yorder, function, mean, rms
diff --git a/noao/twodspec/longslit/transform/igsfit/igsfit.x b/noao/twodspec/longslit/transform/igsfit/igsfit.x
new file mode 100644
index 00000000..14e8e51e
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/igsfit.x
@@ -0,0 +1,373 @@
+include <mach.h>
+include <pkg/gtools.h>
+include <pkg/igsfit.h>
+
+define HELP "noao$lib/scr/igsfit.key"
+define PROMPT "fitcoords surface fitting options"
+
+
+# IGS_FIT1 -- Fit z = f(x, y)
+
+procedure igs_fit1 (sf, gp, gplog, gt, axis, pts, npts, labels, interactive)
+
+pointer sf # GSURFIT pointer
+pointer gp # GIO pointer
+pointer gplog # GIO pointer for plot log
+pointer gt # GTOOLS pointer
+int axis[2] # Axis definitions
+real pts[npts, ARB] # Data
+int npts # Number of pts points
+char labels[SZ_LINE, ARB] # Identification labels
+int interactive # Interactive?
+
+extern igs_solve1()
+
+begin
+ call igs_fit (sf, gp, gplog, gt, axis, pts, npts, labels, interactive,
+ igs_solve1)
+end
+
+
+# IGS_FIT2 -- Fit z = x + f(y)
+
+procedure igs_fit2 (sf, gp, gplog, gt, axis, pts, npts, labels, interactive)
+
+pointer sf # GSURFIT pointer
+pointer gp # GIO pointer
+pointer gplog # GIO pointer for plot log
+pointer gt # GTOOLS pointer
+int axis[2] # Axis definitions
+real pts[npts, ARB] # Data
+int npts # Number of pts points
+char labels[SZ_LINE, ARB] # Identification labels
+int interactive # Interactive?
+
+extern igs_solve2()
+
+begin
+ call igs_fit (sf, gp, gplog, gt, axis, pts, npts, labels, interactive,
+ igs_solve2)
+end
+
+
+# IGS_FIT3 -- Fit z = y + f(x)
+
+procedure igs_fit3 (sf, gp, gplog, gt, axis, pts, npts, labels, interactive)
+
+pointer sf # GSURFIT pointer
+pointer gp # GIO pointer
+pointer gplog # GIO pointer for plot log
+pointer gt # GTOOLS pointer
+int axis[2] # Axis definitions
+real pts[npts, ARB] # Data
+int npts # Number of pts points
+char labels[SZ_LINE, ARB] # Identification labels
+int interactive # Interactive?
+
+extern igs_solve3()
+
+begin
+ call igs_fit (sf, gp, gplog, gt, axis, pts, npts, labels, interactive,
+ igs_solve3)
+end
+
+
+# IGS_FIT -- Interactive surface fitting.
+
+procedure igs_fit (sf, gp, gplog, gt, axis, pts, npts, labels, interactive, igs_solve)
+
+pointer sf # GSURFIT pointer
+pointer gp # GIO pointer
+pointer gplog # GIO pointer for plot log
+pointer gt # GTOOLS pointer
+int axis[2] # Axis definitions
+real pts[npts, ARB] # Data
+int npts # Number of pts points
+char labels[SZ_LINE, ARB] # Identification labels
+int interactive # Interactive?
+extern igs_solve() # Surface solution routine
+
+int i, newgraph, ztype, dtype, refpt, refpt1
+real zval, zval1
+pointer wts
+
+real wx, wy
+int wcs, key
+char cmd[SZ_LINE]
+
+int clgcur(), gt_gcur(), igs_nearest(), igs_nearestd(), igs_nearestu()
+errchk igs_solve
+
+include "igsfit.com"
+
+begin
+ # Compute a solution and set the residuals.
+
+ call igs_solve (sf, pts[1,X], pts[1,Y], pts[1,Z], pts[1,W], npts)
+ call xgsvector (sf, pts[1,X], pts[1,Y], pts[1,S], npts)
+ call asubr (pts[1,Z], pts[1,S], pts[1,R], npts)
+ call aavgr (pts[1,R], npts, mean, rms)
+ call igs_params (gt)
+
+ # Return if not interactive.
+
+ ztype = INDEFI
+ if ((gp == NULL) || (interactive == NO))
+ goto 30
+
+ call malloc (wts, npts, TY_REAL)
+ call amovr (pts[1,W], Memr[wts], npts)
+
+ call igs_graph (gp, gt, ztype, refpt, axis, pts, npts, labels)
+ newgraph = NO
+
+ # Read cursor commands.
+
+10 while (gt_gcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) != EOF) {
+ switch (key) {
+ case '?':
+ # Print help text.
+
+ call gpagefile (gp, HELP, PROMPT)
+
+ case ':':
+ # List or set parameters
+
+ if (cmd[1] == '/')
+ call gt_colon (cmd, gp, gt, newgraph)
+ else
+ call igs_colon (cmd, gp, sf)
+
+ # Set abscissa
+
+ case 'x':
+ call printf ("Select abscissa (x, y, z, s, r): ")
+ if (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF)
+ goto 10
+ call printf ("\n")
+
+ switch (key) {
+ case 'x':
+ i = X
+ case 'y':
+ i = Y
+ case 'z':
+ i = Z
+ case 's':
+ i = S
+ case 'r':
+ i = R
+ default:
+ call printf ("\07\n")
+ goto 10
+ }
+
+ if (axis[1] != i) {
+ axis[1] = i
+ call gt_setr (gt, GTXMIN, INDEF)
+ call gt_setr (gt, GTXMAX, INDEF)
+ }
+
+ # Set ordinate
+
+ case 'y':
+ call printf ("Select ordinate (x, y, z, s, r): ")
+ if(clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF)
+ goto 10
+ call printf ("\n")
+
+ switch (key) {
+ case 'x':
+ i = X
+ case 'y':
+ i = Y
+ case 'z':
+ i = Z
+ case 's':
+ i = S
+ case 'r':
+ i = R
+ default:
+ call printf ("\07\n")
+ goto 10
+ }
+
+ if (axis[2] != i) {
+ axis[2] = i
+ call gt_setr (gt, GTYMIN, INDEF)
+ call gt_setr (gt, GTYMAX, INDEF)
+ }
+
+ case 'r':
+ newgraph = YES
+
+ case 'z':
+ if (IS_INDEFI (ztype)) {
+ refpt = igs_nearest (gp, ztype, refpt, axis, pts, npts, wx,
+ wy, wcs)
+
+ call printf ("Zoom type (x, y, z): ")
+ if (clgcur ("cursor",wx,wy,wcs,key,cmd,SZ_LINE) == EOF)
+ goto 10
+ call printf ("\n")
+
+ switch (key) {
+ case 'x':
+ ztype = X
+ case 'y':
+ ztype = Y
+ case 'z':
+ ztype = Z
+ default:
+ call printf ("\07\n")
+ goto 10
+ }
+
+ newgraph = YES
+ }
+
+ case 'p':
+ if (!IS_INDEFI (ztype)) {
+ ztype = INDEFI
+ newgraph = YES
+ }
+
+ case 'l':
+ if (!IS_INDEFI (ztype)) {
+ refpt1 = 0
+ zval = pts[refpt, ztype]
+ zval1 = -MAX_REAL
+ do i = 1, npts {
+ if ((pts[i,ztype] < zval) && (pts[i,ztype] > zval1)) {
+ refpt1 = i
+ zval1 = pts[refpt1,ztype]
+ }
+ }
+
+ if (refpt1 != 0) {
+ refpt = refpt1
+ newgraph = YES
+ }
+ }
+
+ case 'n':
+ if (!IS_INDEFI (ztype)) {
+ refpt1 = 0
+ zval = pts[refpt, ztype]
+ zval1 = MAX_REAL
+ do i = 1, npts {
+ if ((pts[i,ztype] > zval) && (pts[i,ztype] < zval1)) {
+ refpt1 = i
+ zval1 = pts[refpt1,ztype]
+ }
+ }
+
+ if (refpt1 != 0) {
+ refpt = refpt1
+ newgraph = YES
+ }
+ }
+
+ case 'c':
+ # cursor read
+ i = igs_nearest (gp, ztype, refpt, axis, pts, npts, wx, wy, wcs)
+ call printf ("%g %g %g %g %g %g\n")
+ call pargr (pts[i, X])
+ call pargr (pts[i, Y])
+ call pargr (pts[i, Z])
+ call pargr (pts[i, W])
+ call pargr (pts[i, S])
+ call pargr (pts[i, R])
+
+ case 'd':
+ i = igs_nearestd (gp, ztype, refpt, axis, pts, npts, wx, wy,
+ wcs)
+ if (i == 0)
+ goto 10
+
+ call gscur (gp, real (pts[i,axis[1]]), real (pts[i,axis[2]]))
+
+ call printf ( "Delete 'p'oint or constant 'x', 'y', or 'z': ")
+ if (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF)
+ goto 10
+ call printf ("\n")
+
+ switch (key) {
+ case 'p':
+ dtype = 0
+ case 'x':
+ dtype = X
+ case 'y':
+ dtype = Y
+ case 'z':
+ dtype = Z
+ default:
+ call printf ("\07\n")
+ goto 10
+ }
+
+ call igs_delete (gp, gt, ztype, i, axis, pts, npts, dtype)
+
+ case 'u':
+ i = igs_nearestu (gp, ztype, refpt, axis, pts, npts, wx, wy,
+ wcs)
+ if (i == 0)
+ goto 10
+
+ call gscur (gp, real (pts[i,axis[1]]), real (pts[i,axis[2]]))
+
+ call printf ( "Undelete 'p'oint or constant 'x', 'y', or 'z': ")
+ if (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF)
+ goto 10
+ call printf ("\n")
+
+ switch (key) {
+ case 'p':
+ dtype = 0
+ case 'x':
+ dtype = X
+ case 'y':
+ dtype = Y
+ case 'z':
+ dtype = Z
+ default:
+ call printf ("\07\n")
+ goto 10
+ }
+
+ call igs_undelete (gp, gt, ztype, i, axis, pts, Memr[wts],
+ npts, dtype)
+
+ case 'f':
+ #call printf ("Fitting ...")
+ #call flush (STDOUT)
+ call igs_solve (sf,pts[1,X],pts[1,Y],pts[1,Z],pts[1,W],npts)
+ call xgsvector (sf, pts[1,X], pts[1,Y], pts[1,S], npts)
+ call asubr (pts[1,Z], pts[1,S], pts[1,R], npts)
+ call aavgr (pts[1,R], npts, mean, rms)
+ call igs_params (gt)
+ newgraph = YES
+
+ case 'w':
+ call gt_window (gt, gp, "cursor", newgraph)
+
+ case 'I':
+ call fatal (0, "Interrupt")
+
+ default:
+ # Ring the bell.
+
+ call printf ("\07\n")
+ }
+
+ if (newgraph == YES) {
+ call igs_graph (gp, gt, ztype, refpt, axis, pts, npts, labels)
+ newgraph = NO
+ }
+ }
+
+ call mfree (wts, TY_REAL)
+
+30 call igs_graph (gplog, gt, ztype, refpt, axis, pts, npts, labels)
+
+end
diff --git a/noao/twodspec/longslit/transform/igsfit/igsget.x b/noao/twodspec/longslit/transform/igsfit/igsget.x
new file mode 100644
index 00000000..ccd1fb6c
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/igsget.x
@@ -0,0 +1,62 @@
+include <pkg/igsfit.h>
+
+# IGS_GETI -- Get the value of an integer parameter.
+
+int procedure igs_geti (param)
+
+int param # IGS parameter
+
+include "igsfit.com"
+
+begin
+ switch (param) {
+ case IGS_XORDER:
+ return (xorder)
+ case IGS_YORDER:
+ return (yorder)
+ default:
+ call error (0, "igs_geti: Unknown parameter")
+ }
+end
+
+
+# IGS_GETS -- Get the value of a string parameter.
+
+procedure igs_gets (param, str, maxchar)
+
+int param # IGS parameter
+char str[maxchar] # String
+int maxchar # Maximum number of characters
+
+include "igsfit.com"
+
+begin
+ switch (param) {
+ case IGS_FUNCTION:
+ call strcpy (function, str, maxchar)
+ default:
+ call error (0, "igs_gets: Unknown parameter")
+ }
+end
+
+
+# IGS_GETR -- Get the values of real valued fitting parameters.
+
+real procedure igs_getr (param)
+
+int param # Parameter to be get
+
+include "igsfit.com"
+
+begin
+ switch (param) {
+ case IGS_XMIN:
+ return (xmin)
+ case IGS_XMAX:
+ return (xmax)
+ case IGS_YMIN:
+ return (ymin)
+ case IGS_YMAX:
+ return (ymax)
+ }
+end
diff --git a/noao/twodspec/longslit/transform/igsfit/igsgraph.x b/noao/twodspec/longslit/transform/igsfit/igsgraph.x
new file mode 100644
index 00000000..83eba7e1
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/igsgraph.x
@@ -0,0 +1,73 @@
+include <mach.h>
+include <gset.h>
+include <pkg/gtools.h>
+include <pkg/igsfit.h>
+
+procedure igs_graph (gp, gt, ztype, refpt, axis, pts, npts, labels)
+
+pointer gp # GIO pointer
+pointer gt # GTOOLS pointer
+int ztype # Zoom type
+int refpt # Reference point
+int axis[2] # Axis definitions
+real pts[npts, ARB] # Data
+int npts # Number of pts points
+char labels[SZ_LINE, ARB] # Data labels
+
+int i, x, y
+real xmin, xmax, ymin, ymax, xsize, ysize, gt_getr()
+
+begin
+ if (gp == NULL)
+ return
+
+ x = axis[1]
+ y = axis[2]
+
+ call gt_sets (gt, GTXLABEL, labels[1, x])
+ call gt_sets (gt, GTYLABEL, labels[1, y])
+ xsize = gt_getr (gt, GTXSIZE)
+ ysize = gt_getr (gt, GTYSIZE)
+
+ call gclear (gp)
+
+ if (IS_INDEFI (ztype)) {
+ call gascale (gp, pts[1, x], npts, 1)
+ call gascale (gp, pts[1, y], npts, 2)
+ } else {
+ xmin = MAX_REAL
+ xmax = -MAX_REAL
+ ymin = MAX_REAL
+ ymax = -MAX_REAL
+ do i = 1, npts {
+ if (pts[i,ztype] != pts[refpt,ztype])
+ next
+ xmin = min (xmin, pts[i,x])
+ xmax = max (xmax, pts[i,x])
+ ymin = min (ymin, pts[i,y])
+ ymax = max (ymax, pts[i,y])
+ }
+ call gswind (gp, xmin, xmax, ymin, ymax)
+ }
+
+ call gt_swind (gp, gt)
+ call gt_labax (gp, gt)
+
+ if (IS_INDEFI (ztype)) {
+ do i = 1, npts {
+ if (pts[i,W] == 0.)
+ call gmark (gp, pts[i,x], pts[i,y], GM_CROSS, xsize, ysize)
+ else
+ call gmark (gp, pts[i,x], pts[i,y], GM_PLUS, xsize, ysize)
+ }
+ } else {
+ do i = 1, npts {
+ if (pts[i,ztype] != pts[refpt,ztype])
+ next
+ if (pts[i,W] == 0.)
+ call gmark (gp, pts[i,x], pts[i,y], GM_CROSS, xsize, ysize)
+ else
+ call gmark (gp, pts[i,x], pts[i,y], GM_PLUS, xsize, ysize)
+ }
+ }
+end
diff --git a/noao/twodspec/longslit/transform/igsfit/igsinit.x b/noao/twodspec/longslit/transform/igsfit/igsinit.x
new file mode 100644
index 00000000..f084e7ff
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/igsinit.x
@@ -0,0 +1,21 @@
+include <pkg/igsfit.h>
+
+# IGS_INIT -- Initialize the surface fitting parameters.
+
+procedure igs_init (function, xorder, yorder, xmin, xmax, ymin, ymax)
+
+char function[ARB] # Function
+int xorder # X order
+int yorder # Y order
+real xmin, xmax # X range
+real ymin, ymax # Y range
+
+begin
+ call igs_sets (IGS_FUNCTION, function)
+ call igs_seti (IGS_XORDER, xorder)
+ call igs_seti (IGS_YORDER, yorder)
+ call igs_setr (IGS_XMIN, xmin)
+ call igs_setr (IGS_XMAX, xmax)
+ call igs_setr (IGS_YMIN, ymin)
+ call igs_setr (IGS_YMAX, ymax)
+end
diff --git a/noao/twodspec/longslit/transform/igsfit/igsnearest.x b/noao/twodspec/longslit/transform/igsfit/igsnearest.x
new file mode 100644
index 00000000..69888509
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/igsnearest.x
@@ -0,0 +1,51 @@
+include <mach.h>
+include <gset.h>
+include <pkg/igsfit.h>
+
+int procedure igs_nearest (gp, ztype, refpt, axis, pts, npts, wx, wy, wcs)
+
+pointer gp # GIO pointer
+int ztype # Zoom type
+int refpt # Reference point
+int axis[2] # Axes
+real pts[npts, ARB] # Data points
+int npts # Number of data points
+real wx, wy # Cursor coordinates
+int wcs # WCS
+
+int i, j, x, y
+real r2, r2min, x0, y0
+
+begin
+ x = axis[1]
+ y = axis[2]
+
+ call gctran (gp, wx, wy, wx, wy, wcs, 0)
+ r2min = MAX_REAL
+ j = 0
+
+ if (IS_INDEFI (ztype)) {
+ do i = 1, npts {
+ call gctran (gp, pts[i,x], pts[i,y], x0, y0, wcs, 0)
+ r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+ if (r2 < r2min) {
+ r2min = r2
+ j = i
+ }
+ }
+ } else {
+ do i = 1, npts {
+ if (pts[i,ztype] != pts[refpt,ztype])
+ next
+ call gctran (gp, pts[i,x], pts[i,y], x0, y0, wcs, 0)
+ r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+ if (r2 < r2min) {
+ r2min = r2
+ j = i
+ }
+ }
+ }
+
+ call gscur (gp, real (pts[j,x]), real (pts[j,y]))
+ return (j)
+end
diff --git a/noao/twodspec/longslit/transform/igsfit/igsparams.x b/noao/twodspec/longslit/transform/igsfit/igsparams.x
new file mode 100644
index 00000000..9ecdd422
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/igsparams.x
@@ -0,0 +1,23 @@
+include <pkg/gtools.h>
+
+# IGS_PARAMS -- Set the GTOOLS parameter string.
+
+procedure igs_params (gt)
+
+pointer gt # GTOOLS pointer
+
+pointer params
+
+include "igsfit.com"
+
+begin
+ call malloc (params, SZ_LINE, TY_CHAR)
+ call sprintf (Memc[params], SZ_LINE,
+ "Function = %s, xorder = %d, yorder = %d, rms = %.4g")
+ call pargstr (function)
+ call pargi (xorder)
+ call pargi (yorder)
+ call pargr (rms)
+ call gt_sets (gt, GTPARAMS, Memc[params])
+ call mfree (params, TY_CHAR)
+end
diff --git a/noao/twodspec/longslit/transform/igsfit/igsset.x b/noao/twodspec/longslit/transform/igsfit/igsset.x
new file mode 100644
index 00000000..ea74e8c9
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/igsset.x
@@ -0,0 +1,59 @@
+include <pkg/igsfit.h>
+
+# IGS_SETS -- Set the values of string valued fitting parameters.
+
+procedure igs_sets (param, str)
+
+int param # Parameter to be set
+char str[ARB] # String value
+
+include "igsfit.com"
+
+begin
+ switch (param) {
+ case IGS_FUNCTION:
+ call strcpy (str, function, SZ_LINE)
+ }
+end
+
+
+# IGS_SETI -- Set the values of integer valued fitting parameters.
+
+procedure igs_seti (param, ival)
+
+int param # Parameter to be set
+int ival # Integer value
+
+include "igsfit.com"
+
+begin
+ switch (param) {
+ case IGS_XORDER:
+ xorder = ival
+ case IGS_YORDER:
+ yorder = ival
+ }
+end
+
+
+# IGS_SETR -- Set the values of real valued fitting parameters.
+
+procedure igs_setr (param, rval)
+
+int param # Parameter to be set
+real rval # Real value
+
+include "igsfit.com"
+
+begin
+ switch (param) {
+ case IGS_XMIN:
+ xmin = rval
+ case IGS_XMAX:
+ xmax = rval
+ case IGS_YMIN:
+ ymin = rval
+ case IGS_YMAX:
+ ymax = rval
+ }
+end
diff --git a/noao/twodspec/longslit/transform/igsfit/igssolve.x b/noao/twodspec/longslit/transform/igsfit/igssolve.x
new file mode 100644
index 00000000..a7e39354
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/igssolve.x
@@ -0,0 +1,173 @@
+include <math/gsurfit.h>
+
+
+# IGS_SOLVE1 -- Fit z = f(x, y).
+
+define SFTYPES "|chebyshev|legendre|" # Surface types
+
+procedure igs_solve1 (sf, x, y, z, w, npts)
+
+pointer sf # GSURFIT pointer
+real x[npts] # X points
+real y[npts] # Y points
+real z[npts] # Z points
+real w[npts] # Weights
+int npts # Number of points
+
+int i, nfunc, ix, iy
+pointer sf1, sf2, resids
+
+int strdic()
+
+include "igsfit.com"
+
+begin
+ # Determine the function type.
+
+ nfunc = strdic (function, function, SZ_LINE, SFTYPES)
+
+ # Fit the first surface.
+
+ ix = min (2, xorder)
+ iy = min (2, yorder)
+ call xgsinit (sf1, nfunc, ix, iy, NO, xmin, xmax, ymin, ymax)
+ call xgsfit (sf1, x, y, z, w, npts, WTS_USER, i)
+
+ switch (i) {
+ case SINGULAR:
+ call eprintf ("Singular solution\n")
+ case NO_DEG_FREEDOM:
+ call error (0, "No degrees of freedom")
+ }
+
+ # Evaluate the first surface and fit the residuals.
+
+ call malloc (resids, npts, TY_REAL)
+ call xgsvector (sf1, x, y, Memr[resids], npts)
+ call asubr (z, Memr[resids], Memr[resids], npts)
+
+ call xgsinit (sf2, nfunc, xorder, yorder, YES, xmin,xmax,ymin,ymax)
+ call xgsfit (sf2, x, y, Memr[resids], w, npts, WTS_USER, i)
+
+ switch (i) {
+ case SINGULAR:
+ call eprintf ("Singular solution\n")
+ case NO_DEG_FREEDOM:
+ call error (0, "No degrees of freedom")
+ }
+
+ # Add the two surfaces and free memory.
+
+ call xgsadd (sf1, sf2, sf)
+ call xgsfree (sf1)
+ call xgsfree (sf2)
+ call mfree (resids, TY_REAL)
+end
+
+
+# IGS_SOLVE2 -- Fit z = x + f(y).
+
+
+procedure igs_solve2 (sf, x, y, z, w, npts)
+
+pointer sf # GSURFIT pointer
+real x[npts] # X points
+real y[npts] # Y points
+real z[npts] # Z points
+real w[npts] # Weights
+int npts # Number of points
+
+int i, nfunc
+real a
+pointer sf1
+
+int strdic()
+real xgsgcoeff()
+
+include "igsfit.com"
+
+begin
+ nfunc = strdic (function, function, SZ_LINE, SFTYPES)
+ call xgsinit (sf1, nfunc, 1, yorder, NO, xmin, xmax, ymin, ymax)
+
+ call asubr (z, x, z, npts)
+ call xgsfit (sf1, x, y, z, w, npts, WTS_USER, i)
+ call aaddr (z, x, z, npts)
+
+ switch (i) {
+ case SINGULAR:
+ call eprintf ("Singular solution\n")
+ case NO_DEG_FREEDOM:
+ call error (0, "No degrees of freedom")
+ }
+
+ call xgsfree (sf)
+ call xgsinit (sf, nfunc, 2, yorder, NO, xmin, xmax, ymin, ymax)
+ a = xgsgcoeff (sf1, 1, 1)
+
+ a = a + (xmin + xmax) / 2
+ call xgsscoeff (sf, 1, 1, a)
+
+ a = (xmax - xmin) / 2
+ call xgsscoeff (sf, 2, 1, a)
+
+ do i = 2, yorder {
+ a = xgsgcoeff (sf1, 1, i)
+ call xgsscoeff (sf, 1, i, a)
+ }
+
+ call xgsfree (sf1)
+end
+
+# IGS_SOLVE3 -- Fit z = y + f(x).
+
+procedure igs_solve3 (sf, x, y, z, w, npts)
+
+pointer sf # GSURFIT pointer
+real x[npts] # X points
+real y[npts] # Y points
+real z[npts] # Z points
+real w[npts] # Weights
+int npts # Number of points
+
+int i, nfunc
+real a
+pointer sf1
+
+int strdic()
+real xgsgcoeff()
+
+include "igsfit.com"
+
+begin
+ nfunc = strdic (function, function, SZ_LINE, SFTYPES)
+ call xgsinit (sf1, nfunc, xorder, 1, NO, xmin, xmax, ymin, ymax)
+
+ call asubr (z, y, z, npts)
+ call xgsfit (sf1, x, y, z, w, npts, WTS_USER, i)
+ call aaddr (z, y, z, npts)
+
+ switch (i) {
+ case SINGULAR:
+ call eprintf ("Singular solution\n")
+ case NO_DEG_FREEDOM:
+ call error (0, "No degrees of freedom")
+ }
+
+ call xgsfree (sf)
+ call xgsinit (sf, nfunc, xorder, 2, NO, xmin, xmax, ymin, ymax)
+ a = xgsgcoeff (sf1, 1, 1)
+
+ a = a + (ymin + ymax) / 2
+ call xgsscoeff (sf, 1, 1, a)
+
+ a = (ymax - ymin) / 2
+ call xgsscoeff (sf, 1, 2, a)
+
+ do i = 2, xorder {
+ a = xgsgcoeff (sf1, i, 1)
+ call xgsscoeff (sf, i, 1, a)
+ }
+
+ call xgsfree (sf1)
+end
diff --git a/noao/twodspec/longslit/transform/igsfit/igsundelete.x b/noao/twodspec/longslit/transform/igsfit/igsundelete.x
new file mode 100644
index 00000000..dc7b802e
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/igsundelete.x
@@ -0,0 +1,107 @@
+include <mach.h>
+include <gset.h>
+include <pkg/gtools.h>
+include <pkg/igsfit.h>
+
+int procedure igs_nearestu (gp, ztype, refpt, axis, pts, npts, wx, wy, wcs)
+
+pointer gp # GIO pointer
+int ztype # Zoom type
+int refpt # Reference point
+int axis[2] # Axes
+real pts[npts, ARB] # Data points
+int npts # Number of data points
+real wx, wy # Cursor coordinates
+int wcs # WCS
+
+int i, j, x, y
+real r2, r2min, x0, y0
+
+begin
+ x = axis[1]
+ y = axis[2]
+
+ call gctran (gp, wx, wy, wx, wy, wcs, 0)
+ r2min = MAX_REAL
+ j = 0
+
+ if (IS_INDEFI (ztype)) {
+ do i = 1, npts {
+ if (pts[i,W] != 0.)
+ next
+ call gctran (gp, pts[i, x], pts[i, y], x0, y0, wcs, 0)
+ r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+ if (r2 < r2min) {
+ r2min = r2
+ j = i
+ }
+ }
+ } else {
+ do i = 1, npts {
+ if ((pts[i,ztype] != pts[refpt,ztype]) || (pts[i,W] != 0.))
+ next
+ call gctran (gp, pts[i, x], pts[i, y], x0, y0, wcs, 0)
+ r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+ if (r2 < r2min) {
+ r2min = r2
+ j = i
+ }
+ }
+ }
+
+ return (j)
+end
+
+
+# IGS_UNDELETE - Undelete point or subset.
+
+procedure igs_undelete (gp, gt, ztype, refpt, axis, pts, wts, npts, dtype)
+
+pointer gp # GIO pointer
+pointer gt # GTOOLS pointer
+int ztype # Zoom type
+int refpt # Reference point for undeletion
+int axis[2] # Axes
+real pts[npts, ARB] # Data points
+real wts[npts] # Original weights
+int npts # Number of data points
+int dtype # Undeletion type
+
+int i, x, y
+real xsize, ysize
+
+real gt_getr()
+
+begin
+ x = axis[1]
+ y = axis[2]
+
+ xsize = gt_getr (gt, GTXSIZE)
+ ysize = gt_getr (gt, GTYSIZE)
+
+ switch (dtype) {
+ case X, Y, Z:
+ do i = 1, npts {
+ if (!IS_INDEFI (ztype))
+ if (pts[refpt,ztype] != pts[i,ztype])
+ next
+ if (pts[refpt,dtype] != pts[i,dtype])
+ next
+ call gseti (gp, G_PMLTYPE, 0)
+ call gmark (gp, pts[i,x], pts[i,y], GM_CROSS, xsize, ysize)
+ call gseti (gp, G_PMLTYPE, 1)
+ call gmark (gp, pts[i,x], pts[i,y], GM_PLUS, xsize, ysize)
+ if (wts[i] == 0)
+ wts[i] = 1
+ pts[i,W] = wts[i]
+ }
+ default:
+ call gseti (gp, G_PMLTYPE, 0)
+ call gmark (gp, pts[refpt,x], pts[refpt,y], GM_CROSS, xsize, ysize)
+ call gseti (gp, G_PMLTYPE, 1)
+ call gmark (gp, pts[refpt,x], pts[refpt,y], GM_PLUS, xsize, ysize)
+ if (wts[refpt] == 0)
+ wts[refpt] = 1
+ pts[refpt,W] = wts[refpt]
+ }
+end
diff --git a/noao/twodspec/longslit/transform/igsfit/mkpkg b/noao/twodspec/longslit/transform/igsfit/mkpkg
new file mode 100644
index 00000000..ac5a6ca9
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/mkpkg
@@ -0,0 +1,21 @@
+# Interactive General Surface Fitting Package
+
+$checkout libpkg.a ../../
+$update libpkg.a
+$checkin libpkg.a ../../
+$exit
+
+libpkg.a:
+ igscolon.x igsfit.com <gset.h>
+ igsdelete.x <gset.h> <mach.h> <pkg/gtools.h> <pkg/igsfit.h>
+ igsfit.x igsfit.com <mach.h> <pkg/gtools.h> <pkg/igsfit.h>
+ igsget.x igsfit.com <pkg/igsfit.h>
+ igsgraph.x <gset.h> <mach.h> <pkg/gtools.h> <pkg/igsfit.h>
+ igsinit.x <pkg/igsfit.h>
+ igsnearest.x <gset.h> <mach.h> <pkg/igsfit.h>
+ igsparams.x igsfit.com <pkg/gtools.h>
+ igsset.x igsfit.com <pkg/igsfit.h>
+ igssolve.x igsfit.com <math/gsurfit.h>
+ igsundelete.x <gset.h> <mach.h> <pkg/gtools.h> <pkg/igsfit.h>
+ xgs.x <math/gsurfit.h>
+ ;
diff --git a/noao/twodspec/longslit/transform/igsfit/xgs.x b/noao/twodspec/longslit/transform/igsfit/xgs.x
new file mode 100644
index 00000000..7d2ea331
--- /dev/null
+++ b/noao/twodspec/longslit/transform/igsfit/xgs.x
@@ -0,0 +1,243 @@
+include <math/gsurfit.h>
+
+# XGS -- These routines provide an interface between real input data and
+# the double precision surface fitting. Rather than make the input data
+# be double precision we only want the internal surface fitting arithmetic
+# to be double. But the surface fitting package only provides real
+# arithmetic for real input and double precision arithmetic for double
+# precision input. Hence these interfaces. Note that the save and restore
+# functions use double precision.
+
+# XGSINIT -- Procedure to initialize the surface descriptor.
+
+procedure xgsinit (sf, surface_type, xorder, yorder, xterms, xmin, xmax,
+ ymin, ymax)
+
+pointer sf # surface descriptor
+int surface_type # type of surface to be fitted
+int xorder # x order of surface to be fit
+int yorder # y order of surface to be fit
+int xterms # presence of cross terms
+real xmin # minimum value of x
+real xmax # maximum value of x
+real ymin # minimum value of y
+real ymax # maximum value of y
+
+begin
+ call dgsinit (sf, surface_type, xorder, yorder, xterms, double (xmin),
+ double (xmax), double (ymin), double (ymax))
+end
+
+
+# XGSFIT -- Procedure to solve the normal equations for a surface.
+
+procedure xgsfit (sf, x, y, z, w, npts, wtflag, ier)
+
+pointer sf # surface descriptor
+real x[npts] # array of x values
+real y[npts] # array of y values
+real z[npts] # data array
+real w[npts] # array of weights
+int npts # number of data points
+int wtflag # type of weighting
+int ier # ier = OK, everything OK
+ # ier = SINGULAR, matrix is singular, 1 or more
+ # coefficients are 0.
+ # ier = NO_DEG_FREEDOM, too few points to solve matrix
+
+pointer sp, xd, yd, zd, wd
+errchk salloc
+
+begin
+ call smark (sp)
+ call salloc (xd, npts, TY_DOUBLE)
+ call salloc (yd, npts, TY_DOUBLE)
+ call salloc (zd, npts, TY_DOUBLE)
+ call salloc (wd, npts, TY_DOUBLE)
+ call achtrd (x, Memd[xd], npts)
+ call achtrd (y, Memd[yd], npts)
+ call achtrd (z, Memd[zd], npts)
+ call achtrd (w, Memd[wd], npts)
+ call dgsfit (sf, Memd[xd], Memd[yd], Memd[zd], Memd[wd], npts,
+ wtflag, ier)
+ call sfree (sp)
+end
+
+
+# XGSVECTOR -- Procedure to evaluate the fitted surface at an array of points.
+
+procedure xgsvector (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
+
+pointer sp, xd, yd, zd
+errchk salloc
+
+begin
+ call smark (sp)
+ call salloc (xd, npts, TY_DOUBLE)
+ call salloc (yd, npts, TY_DOUBLE)
+ call salloc (zd, npts, TY_DOUBLE)
+ call achtrd (x, Memd[xd], npts)
+ call achtrd (y, Memd[yd], npts)
+ call dgsvector (sf, Memd[xd], Memd[yd], Memd[zd], npts)
+ call achtdr (Memd[zd], zfit, npts)
+ call sfree (sp)
+end
+
+
+# XGSEVAL -- Procedure to evaluate the fitted surface at a single point.
+
+real procedure xgseval (sf, x, y)
+
+pointer sf # pointer to surface descriptor structure
+real x # x value
+real y # y value
+
+double dgseval()
+
+begin
+ return (real (dgseval (sf, double (x), double (y))))
+end
+
+
+# XGSADD -- Procedure to add the fits from two surfaces together.
+
+procedure xgsadd (sf1, sf2, sf3)
+
+pointer sf1 # pointer to the first surface
+pointer sf2 # pointer to the second surface
+pointer sf3 # pointer to the output surface
+
+begin
+ call dgsadd (sf1, sf2, sf3)
+end
+
+
+# XGSFREE -- Procedure to free the surface descriptor
+
+procedure xgsfree (sf)
+
+pointer sf # the surface descriptor
+
+begin
+ call dgsfree (sf)
+end
+
+
+# XGSGCOEFF -- Procedure to fetch a particular coefficient.
+
+real procedure xgsgcoeff (sf, xorder, yorder)
+
+pointer sf # pointer to the surface fitting descriptor
+int xorder # X order of desired coefficent
+int yorder # Y order of desired coefficent
+
+double dgsgcoeff()
+
+begin
+ return (real (dgsgcoeff (sf, xorder, yorder)))
+end
+
+
+# XGSSCOEFF -- Procedure to set a particular coefficient.
+
+procedure xgsscoeff (sf, xorder, yorder, coeff)
+
+pointer sf # pointer to the surface fitting descriptor
+int xorder # X order of desired coefficent
+int yorder # Y order of desired coefficent
+real coeff # Coefficient value
+
+begin
+ call dgsscoeff (sf, xorder, yorder, double (coeff))
+end
+
+
+# XGSGETR -- Procedure to fetch a real gsurfit parameter
+
+real procedure xgsgetr (sf, parameter)
+
+pointer sf # pointer to the surface fit
+int parameter # parameter to be fetched
+
+double dgsgetd()
+
+begin
+ return (real (dgsgetd (sf, parameter)))
+end
+
+
+# XGSGETI -- Procedure to fetch an integer parameter
+
+int procedure xgsgeti (sf, parameter)
+
+pointer sf # pointer to the surface fit
+int parameter # integer parameter
+
+int dgsgeti()
+
+begin
+ return (dgsgeti (sf, parameter))
+end
+
+
+# XGSSAVE -- Procedure to save the surface fit for later use by the
+# evaluate routines.
+#
+# NOTE THAT THIS USES DOUBLE PRECISION FOR THE COEFFICIENTS.
+
+procedure xgssave (sf, fit)
+
+pointer sf # pointer to the surface descriptor
+double fit[ARB] # array for storing fit
+
+begin
+ call dgssave (sf, fit)
+end
+
+
+# XGSRESTORE -- Procedure to restore the surface fit stored by GSSAVE
+# to the surface descriptor for use by the evaluating routines.
+#
+# NOTE THAT THIS USES DOUBLE PRECISION FOR THE COEFFICIENTS.
+
+procedure xgsrestore (sf, fit)
+
+pointer sf # surface descriptor
+double fit[ARB] # array containing the surface parameters and
+
+begin
+ call dgsrestore (sf, fit)
+end
+
+
+# XGSDER -- Procedure to calculate a new surface which is a derivative of
+# the previous surface
+
+procedure xgsder (sf1, x, y, zfit, npts, nxd, nyd)
+
+pointer sf1 # pointer to the previous surface
+real x[npts] # x values
+real y[npts] # y values
+real zfit[npts] # fitted values
+int npts # number of points
+int nxd, nyd # order of the derivatives in x and y
+
+pointer sp, xd, yd, zd
+
+begin
+ call smark (sp)
+ call salloc (xd, npts, TY_DOUBLE)
+ call salloc (yd, npts, TY_DOUBLE)
+ call salloc (zd, npts, TY_DOUBLE)
+ call achtrd (x, Memd[xd], npts)
+ call achtrd (y, Memd[yd], npts)
+ call dgsder (sf1, Memd[xd], Memd[yd], Memd[zd], npts, nxd, nyd)
+ call achtdr (Memd[zd], zfit, npts)
+ call sfree (sp)
+end