aboutsummaryrefslogtreecommitdiff
path: root/pkg/images/lib/geograph.gx
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/images/lib/geograph.gx')
-rw-r--r--pkg/images/lib/geograph.gx1379
1 files changed, 1379 insertions, 0 deletions
diff --git a/pkg/images/lib/geograph.gx b/pkg/images/lib/geograph.gx
new file mode 100644
index 00000000..5c42de24
--- /dev/null
+++ b/pkg/images/lib/geograph.gx
@@ -0,0 +1,1379 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <math/gsurfit.h>
+include <pkg/gtools.h>
+include <mach.h>
+include <math.h>
+include <gset.h>
+include "geomap.h"
+include "geogmap.h"
+
+define MAX_PARAMS (10 * SZ_LINE)
+define NINTERVALS 5
+define NGRAPH 100
+
+$for (r)
+
+# GEO_LABEL -- Annotate the plot.
+
+procedure geo_label (plot_type, gt, fit)
+
+int plot_type #I type of plot
+pointer gt #I gtools descriptor
+pointer fit #I geomap fit parameters
+
+int npts
+pointer sp, params, xtermlab, ytermlab
+real xrms, yrms, rej
+int strlen(), rg_wrdstr()
+
+begin
+ call smark (sp)
+ call salloc (params, MAX_PARAMS, TY_CHAR)
+ call salloc (xtermlab, SZ_FNAME, TY_CHAR)
+ call salloc (ytermlab, SZ_FNAME, TY_CHAR)
+
+ npts = max (0, GM_NPTS(fit) - GM_NWTS0(fit))
+ xrms = max (0.0d0, GM_XRMS(fit))
+ yrms = max (0.0d0, GM_YRMS(fit))
+ if (npts > 1) {
+ xrms = sqrt (xrms / (npts - 1))
+ yrms = sqrt (yrms / (npts - 1))
+ } else {
+ xrms = 0.0
+ yrms = 0.0
+ }
+ if (IS_INDEFD(GM_REJECT(fit)))
+ rej = INDEFR
+ else if (GM_REJECT(fit) > MAX_REAL)
+ rej = INDEFR
+ else
+ rej = GM_REJECT(fit)
+
+ # Print data parameters.
+ if (GM_PROJECTION(fit) == GM_NONE)
+ call sprintf (Memc[params], MAX_PARAMS,
+ "GEOMAP: function = %s npts = %d reject = %g nrej = %d\n")
+ else
+ call sprintf (Memc[params], MAX_PARAMS,
+ "CCMAP: function = %s npts = %d reject = %g nrej = %d\n")
+
+ switch (GM_FUNCTION(fit)) {
+ case GS_LEGENDRE:
+ call pargstr ("legendre")
+ case GS_CHEBYSHEV:
+ call pargstr ("chebyshev")
+ case GS_POLYNOMIAL:
+ call pargstr ("polynomial")
+ }
+ call pargi (GM_NPTS(fit))
+ call pargr (rej)
+ call pargi (GM_NWTS0(fit))
+
+ # Print fit parameters.
+ switch (plot_type) {
+ case FIT:
+
+ if (rg_wrdstr ((GM_XXTERMS(fit) + 1), Memc[xtermlab], SZ_FNAME,
+ GM_XFUNCS) <= 0)
+ call strcpy ("none", Memc[xtermlab], SZ_FNAME)
+ if (rg_wrdstr ((GM_YXTERMS(fit) + 1), Memc[ytermlab], SZ_FNAME,
+ GM_XFUNCS) <= 0)
+ call strcpy ("none", Memc[ytermlab], SZ_FNAME)
+ if (GM_PROJECTION(fit) == GM_NONE)
+ call sprintf (Memc[params+strlen(Memc[params])], MAX_PARAMS,
+ "X fit: xorder = %d yorder = %d xterms = %s stdev = %8.3g\n")
+ else
+ call sprintf (Memc[params+strlen(Memc[params])], MAX_PARAMS,
+ "XI fit: xorder = %d yorder = %d xterms = %s stdev = %8.3g arcsec\n")
+ call pargi (GM_XXORDER(fit))
+ call pargi (GM_XYORDER(fit))
+ call pargstr (Memc[xtermlab])
+ call pargr (xrms)
+
+ if (GM_PROJECTION(fit) == GM_NONE)
+ call sprintf (Memc[params+strlen(Memc[params])], MAX_PARAMS,
+ "Y fit: xorder = %d yorder = %d xterms = %s stdev = %8.3g\n")
+ else
+ call sprintf (Memc[params+strlen(Memc[params])], MAX_PARAMS,
+ "ETA fit: xorder = %d yorder = %d xterms = %s stdev = %8.3g arcsec\n")
+ call pargi (GM_YXORDER(fit))
+ call pargi (GM_YYORDER(fit))
+ call pargstr (Memc[ytermlab])
+ call pargr (yrms)
+
+ case XXRESID, XYRESID:
+
+ if (rg_wrdstr ((GM_XXTERMS(fit) + 1), Memc[xtermlab], SZ_FNAME,
+ GM_XFUNCS) <= 0)
+ call strcpy ("none", Memc[xtermlab], SZ_FNAME)
+ if (GM_PROJECTION(fit) == GM_NONE)
+ call sprintf (Memc[params+strlen(Memc[params])], MAX_PARAMS,
+ "X fit: xorder = %d yorder = %d xterms = %s rms = %8.3g\n")
+ else
+ call sprintf (Memc[params+strlen(Memc[params])], MAX_PARAMS,
+ "XI fit: xorder = %d yorder = %d xterms = %s rms = %8.3g arcsec\n")
+ call pargi (GM_XXORDER(fit))
+ call pargi (GM_XYORDER(fit))
+ call pargstr (Memc[xtermlab])
+ call pargr (xrms)
+
+ case YXRESID, YYRESID:
+
+ if (rg_wrdstr ((GM_YXTERMS(fit) + 1), Memc[ytermlab], SZ_FNAME,
+ GM_XFUNCS) <= 0)
+ call strcpy ("none", Memc[ytermlab], SZ_FNAME)
+ if (GM_PROJECTION(fit) == GM_NONE)
+ call sprintf (Memc[params+strlen(Memc[params])], MAX_PARAMS,
+ "Y fit: xorder = %d yorder = %d xterms = %s rms = %8.3g\n")
+ else
+ call sprintf (Memc[params+strlen(Memc[params])], MAX_PARAMS,
+ "ETA fit: xorder = %d yorder = %d xterms = %s rms = %8.3g arcsec\n")
+ call pargi (GM_YXORDER(fit))
+ call pargi (GM_YYORDER(fit))
+ call pargstr (Memc[ytermlab])
+ call pargr (yrms)
+
+ default:
+
+ # do nothing gracefully
+ }
+
+ call gt_sets (gt, GTPARAMS, Memc[params])
+
+ call sfree (sp)
+end
+
+
+# GEO_GTSET -- Write title and labels.
+
+procedure geo_gtset (plot_type, gt, fit)
+
+int plot_type #I plot type
+pointer gt #I plot descriptor
+pointer fit #I fit descriptor
+
+char str[SZ_LINE]
+int nchars
+int gstrcpy()
+
+begin
+ nchars = gstrcpy (GM_RECORD(fit), str, SZ_LINE)
+
+ switch (plot_type) {
+ case FIT:
+
+ if (GM_PROJECTION(fit) == GM_NONE)
+ call strcpy (": Coordinate Transformation", str[nchars+1],
+ SZ_LINE)
+ else
+ call strcpy (": Celestial Coordinate Transformation",
+ str[nchars+1], SZ_LINE)
+ call gt_sets (gt, GTTITLE, str)
+ if (GM_PROJECTION(fit) == GM_NONE) {
+ call gt_sets (gt, GTXLABEL, "X (in units)")
+ call gt_sets (gt, GTYLABEL, "Y (in units)")
+ } else {
+ call gt_sets (gt, GTXLABEL, "XI (arcsec)")
+ call gt_sets (gt, GTYLABEL, "ETA (arcsec)")
+ }
+
+ case XXRESID:
+
+ if (GM_PROJECTION(fit) == GM_NONE)
+ call strcpy (": X fit Residuals", str[nchars+1], SZ_LINE)
+ else
+ call strcpy (": XI fit Residuals", str[nchars+1], SZ_LINE)
+ call gt_sets (gt, GTTITLE, str)
+ if (GM_PROJECTION(fit) == GM_NONE) {
+ call gt_sets (gt, GTXLABEL, "X (ref units)")
+ call gt_sets (gt, GTYLABEL, "X Residuals (in units)")
+ } else {
+ call gt_sets (gt, GTXLABEL, "X (pixels)")
+ call gt_sets (gt, GTYLABEL, "XI Residuals (arcsec)")
+ }
+
+ case XYRESID:
+
+ if (GM_PROJECTION(fit) == GM_NONE)
+ call strcpy (": X fit Residuals", str[nchars+1], SZ_LINE)
+ else
+ call strcpy (": XI fit Residuals", str[nchars+1], SZ_LINE)
+ call gt_sets (gt, GTTITLE, str)
+ if (GM_PROJECTION(fit) == GM_NONE) {
+ call gt_sets (gt, GTXLABEL, "Y (ref units)")
+ call gt_sets (gt, GTYLABEL, "X Residuals (in units)")
+ } else {
+ call gt_sets (gt, GTXLABEL, "Y (pixels)")
+ call gt_sets (gt, GTYLABEL, "XI Residuals (arcsec)")
+ }
+
+ case YXRESID:
+
+ if (GM_PROJECTION(fit) == GM_NONE)
+ call strcpy (": Y fit Residuals", str[nchars+1], SZ_LINE)
+ else
+ call strcpy (": ETA fit Residuals", str[nchars+1], SZ_LINE)
+ call gt_sets (gt, GTTITLE, str)
+ if (GM_PROJECTION(fit) == GM_NONE) {
+ call gt_sets (gt, GTXLABEL, "X (ref units)")
+ call gt_sets (gt, GTYLABEL, "Y (Residuals (in units)")
+ } else {
+ call gt_sets (gt, GTXLABEL, "X (pixels)")
+ call gt_sets (gt, GTYLABEL, "ETA Residuals (arcsec)")
+ }
+
+ case YYRESID:
+
+ if (GM_PROJECTION(fit) == GM_NONE)
+ call strcpy (": Y fit Residuals", str[nchars+1], SZ_LINE)
+ else
+ call strcpy (": ETA fit Residuals", str[nchars+1], SZ_LINE)
+ call gt_sets (gt, GTTITLE, str)
+ if (GM_PROJECTION(fit) == GM_NONE) {
+ call gt_sets (gt, GTXLABEL, "Y (ref units)")
+ call gt_sets (gt, GTYLABEL, "Y Residuals (in units)")
+ } else {
+ call gt_sets (gt, GTXLABEL, "Y (pixels)")
+ call gt_sets (gt, GTYLABEL, "ETA Residuals (arcsec)")
+ }
+
+ default:
+
+ # do nothing gracefully
+ }
+end
+
+
+# GEO_COLON -- Process the colon commands.
+
+procedure geo_colon (gd, fit, gfit, cmdstr, newgraph)
+
+pointer gd #I graphics stream
+pointer fit #I pointer to fit structure
+pointer gfit #I pointer to the gfit structure
+char cmdstr[ARB] #I command string
+int newgraph #I plot new graph
+
+int ncmd, ival
+pointer sp, str, cmd
+real rval
+int nscan(), strdic(), rg_wrdstr()
+
+begin
+ call smark (sp)
+ call salloc (cmd, SZ_LINE, TY_CHAR)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ call sscan (cmdstr)
+ call gargwrd (Memc[cmd], SZ_LINE)
+ if (nscan() == 0) {
+ call sfree (sp)
+ return
+ }
+
+ ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, GM_CMDS)
+ switch (ncmd) {
+ case GMCMD_SHOW:
+ call gdeactivate (gd, AW_CLEAR)
+ call printf ("Current Fitting Parameters\n\n")
+ if (GM_PROJECTION(fit) != GM_NONE) {
+ if (rg_wrdstr (GM_PROJECTION(fit), Memc[str], SZ_FNAME,
+ GM_PROJLIST) <= 0)
+ ;
+ call printf ("\tprojection = %s\n")
+ call pargstr (Memc[str])
+ call printf ("\tlngref = %h\n")
+ call pargd (GM_XREFPT(fit))
+ call printf ("\tlatref = %h\n")
+ call pargd (GM_YREFPT(fit))
+ }
+ if (rg_wrdstr (GM_FIT(fit), Memc[str], SZ_FNAME,
+ GM_GEOMETRIES) <= 0)
+ call strcpy ("general", Memc[str], SZ_FNAME)
+ call printf ("\tfitgeometry = %s\n")
+ call pargstr (Memc[str])
+ if (rg_wrdstr (GM_FUNCTION(fit), Memc[str], SZ_FNAME,
+ GM_FUNCS) <= 0)
+ call strcpy ("polynomial", Memc[str], SZ_FNAME)
+ call printf ("\tfunction = %s\n")
+ Call pargstr (Memc[str])
+ call printf ("\txxorder = %d\n")
+ call pargi (GM_XXORDER(fit))
+ call printf ("\txyorder = %d\n")
+ call pargi (GM_XYORDER(fit))
+ if (rg_wrdstr ((GM_XXTERMS(fit) + 1), Memc[str], SZ_FNAME,
+ GM_XFUNCS) <= 0)
+ call strcpy ("none", Memc[str], SZ_FNAME)
+ call printf ("\txxterms = %s\n")
+ call pargstr (Memc[str])
+ call printf ("\tyxorder = %d\n")
+ call pargi (GM_YXORDER(fit))
+ call printf ("\tyyorder = %d\n")
+ call pargi (GM_YYORDER(fit))
+ if (rg_wrdstr ((GM_YXTERMS(fit) + 1), Memc[str], SZ_FNAME,
+ GM_XFUNCS) <= 0)
+ call strcpy ("none", Memc[str], SZ_FNAME)
+ call printf ("\tyxterms = %s\n")
+ call pargstr (Memc[str])
+ if (IS_INDEFD(GM_REJECT(fit)))
+ rval = INDEFR
+ else if (GM_REJECT(fit) > MAX_REAL)
+ rval = INDEFR
+ else
+ rval = GM_REJECT(fit)
+ call printf ("\treject = %f\n")
+ call pargr (rval)
+ call greactivate (gd, AW_PAUSE)
+
+ case GMCMD_PROJECTION:
+ if (rg_wrdstr (GM_PROJECTION(fit), Memc[str], SZ_FNAME,
+ GM_PROJLIST) <= 0)
+ call strcpy ("INDEF", Memc[str], SZ_FNAME)
+ call printf ("projection = %s\n")
+ call pargstr (Memc[str])
+
+ case GMCMD_REFPOINT:
+ call printf ("lngref = %h latref = %h\n")
+ call pargd (GM_XREFPT(fit))
+ call pargd (GM_YREFPT(fit))
+
+ case GMCMD_GEOMETRY:
+ call gargwrd (Memc[cmd], SZ_LINE)
+ if (nscan () == 1) {
+ if (rg_wrdstr (GM_FIT(fit), Memc[str], SZ_FNAME,
+ GM_GEOMETRIES) <= 0)
+ call strcpy ("general", Memc[str], SZ_FNAME)
+ call printf ("fitgeometry = %s\n")
+ call pargstr (Memc[str])
+ } else {
+ ival = strdic (Memc[cmd], Memc[cmd], SZ_LINE, GM_GEOMETRIES)
+ if (ival > 0) {
+ GM_FIT(fit) = ival
+ GG_NEWFUNCTION(gfit) = YES
+ GG_FITERROR(gfit) = NO
+ }
+ }
+
+ case GMCMD_FUNCTION:
+ call gargwrd (Memc[cmd], SZ_LINE)
+ if (nscan () == 1) {
+ if (rg_wrdstr (GM_FUNCTION(fit), Memc[str], SZ_FNAME,
+ GM_FUNCS) <= 0)
+ call strcpy ("polynomial", Memc[str], SZ_FNAME)
+ call printf ("function = %s\n")
+ call pargstr (Memc[str])
+ } else {
+ ival = strdic (Memc[cmd], Memc[cmd], SZ_LINE, GM_FUNCS)
+ if (ival > 0) {
+ GM_FUNCTION(fit) = ival
+ GG_NEWFUNCTION(gfit) = YES
+ GG_FITERROR(gfit) = NO
+ }
+ }
+
+ case GMCMD_ORDER:
+ call gargi (ival)
+ if (nscan () == 1) {
+ call printf (
+ "xxorder = %d xyorder = %d yxorder = %d yyorder = %d\n")
+ call pargi (GM_XXORDER(fit))
+ call pargi (GM_XYORDER(fit))
+ call pargi (GM_YXORDER(fit))
+ call pargi (GM_YYORDER(fit))
+ } else {
+ GM_XXORDER(fit) = max (ival, 2)
+ GM_XYORDER(fit) = max (ival, 2)
+ GM_YXORDER(fit) = max (ival, 2)
+ GM_YYORDER(fit) = max (ival, 2)
+ GG_NEWFUNCTION(gfit) = YES
+ GG_FITERROR(gfit) = NO
+ }
+
+ case GMCMD_XXORDER:
+ call gargi (ival)
+ if (nscan () == 1) {
+ call printf ("xxorder = %d\n")
+ call pargi (GM_XXORDER(fit))
+ } else {
+ GM_XXORDER(fit) = max (ival, 2)
+ GG_NEWFUNCTION(gfit) = YES
+ GG_FITERROR(gfit) = NO
+ }
+
+ case GMCMD_XYORDER:
+ call gargi (ival)
+ if (nscan () == 1) {
+ call printf ("xyorder = %d\n")
+ call pargi (GM_XYORDER(fit))
+ } else {
+ GM_XYORDER(fit) = max (ival,2)
+ GG_NEWFUNCTION(gfit) = YES
+ GG_FITERROR(gfit) = NO
+ }
+
+ case GMCMD_YXORDER:
+ call gargi (ival)
+ if (nscan () == 1) {
+ call printf ("yxorder = %d\n")
+ call pargi (GM_YXORDER(fit))
+ } else {
+ GM_YXORDER(fit) = max (ival, 2)
+ GG_NEWFUNCTION(gfit) = YES
+ GG_FITERROR(gfit) = NO
+ }
+
+ case GMCMD_YYORDER:
+ call gargi (ival)
+ if (nscan () == 1) {
+ call printf ("yyorder = %d\n")
+ call pargi (GM_YYORDER(fit))
+ } else {
+ GM_YYORDER(fit) = max (ival, 2)
+ GG_NEWFUNCTION(gfit) = YES
+ GG_FITERROR(gfit) = NO
+ }
+
+ case GMCMD_XXTERMS:
+ call gargwrd (Memc[cmd], SZ_LINE)
+ if (nscan () == 1) {
+ if (rg_wrdstr ((GM_XXTERMS(fit) + 1), Memc[str], SZ_FNAME,
+ GM_XFUNCS) <= 0)
+ call strcpy ("none", Memc[str], SZ_FNAME)
+ call printf ("xxterms = %s\n")
+ call pargstr (Memc[str])
+ } else {
+ ival = strdic (Memc[cmd], Memc[cmd], SZ_LINE, GM_XFUNCS)
+ if (ival > 0) {
+ GM_XXTERMS(fit) = ival - 1
+ GG_NEWFUNCTION(gfit) = YES
+ GG_FITERROR(gfit) = NO
+ }
+ }
+
+ case GMCMD_YXTERMS:
+ call gargwrd (Memc[cmd], SZ_LINE)
+ if (nscan () == 1) {
+ if (rg_wrdstr ((GM_YXTERMS(fit) + 1), Memc[str], SZ_FNAME,
+ GM_XFUNCS) <= 0)
+ call strcpy ("none", Memc[str], SZ_FNAME)
+ call printf ("yxterms = %s\n")
+ call pargstr (Memc[str])
+ } else {
+ ival = strdic (Memc[cmd], Memc[cmd], SZ_LINE, GM_XFUNCS)
+ if (ival > 0) {
+ GM_YXTERMS(fit) = ival - 1
+ GG_NEWFUNCTION(gfit) = YES
+ GG_FITERROR(gfit) = NO
+ }
+ }
+
+ case GMCMD_REJECT:
+ call gargr (rval)
+ if (nscan() == 1) {
+ if (IS_INDEFD(GM_REJECT(fit)))
+ rval = INDEFR
+ else if (GM_REJECT(fit) > MAX_REAL)
+ rval = INDEFR
+ else
+ rval = GM_REJECT(fit)
+ call printf ("reject = %f\n")
+ call pargr (rval)
+ } else {
+ GM_REJECT(fit) = rval
+ GG_NEWFUNCTION(gfit) = YES
+ GG_FITERROR(gfit) = NO
+ }
+
+ case GMCMD_MAXITER:
+ call gargi (ival)
+ if (nscan() == 1) {
+ call printf ("maxiter = %d\n")
+ call pargi (GM_MAXITER(fit))
+ } else {
+ GM_MAXITER(fit) = ival
+ GG_NEWFUNCTION(gfit) = YES
+ GG_FITERROR(gfit) = NO
+ }
+
+ }
+
+ call sfree (sp)
+end
+
+$endfor
+
+
+$for (rd)
+
+# GEO_1DELETE -- Delete a point from the fit.
+
+procedure geo_1delete$t (gd, xin, yin, wts, userwts, npts, wx, wy, delete)
+
+pointer gd #I pointer to graphics descriptor
+PIXEL xin[ARB] #I x array
+PIXEL yin[ARB] #I y array
+PIXEL wts[ARB] #I array of weights
+PIXEL userwts[ARB] #I array of user weights
+int npts #I number of points
+real wx, wy #I world coordinates
+int delete #I delete points ?
+
+int i, j, pmltype
+real r2min, r2, x0, y0
+int gstati()
+
+begin
+ call gctran (gd, wx, wy, wx, wy, 1, 0)
+ r2min = MAX_REAL
+ j = 0
+
+ if (delete == YES) {
+
+ # Search for nearest point that has not been deleted.
+ do i = 1, npts {
+ if (wts[i] <= PIXEL(0.0))
+ next
+$if (datatype == r)
+ call gctran (gd, xin[i], yin[i], x0, y0, 1, 0)
+$else
+ call gctran (gd, real (xin[i]), real (yin[i]), x0, y0, 1, 0)
+$endif
+ r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+ if (r2 < r2min) {
+ r2min = r2
+ j = i
+ }
+ }
+
+ # Mark point and set weights to 0.
+ if (j != 0) {
+$if (datatype == r)
+ call gscur (gd, xin[j], yin[j])
+ call gmark (gd, xin[j], yin[j], GM_CROSS, 2., 2.)
+$else
+ call gscur (gd, real(xin[j]), real(yin[j]))
+ call gmark (gd, real(xin[j]), real(yin[j]), GM_CROSS, 2., 2.)
+$endif
+ wts[j] = PIXEL(0.0)
+ }
+
+ } else {
+
+ # Search for the nearest deleted point.
+ do i = 1, npts {
+ if (wts[i] > PIXEL(0.0))
+ next
+$if (datatype == r)
+ call gctran (gd, xin[i], yin[i], x0, y0, 1, 0)
+$else
+ call gctran (gd, real(xin[i]), real(yin[i]), x0, y0, 1, 0)
+$endif
+ r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+ if (r2 < r2min) {
+ r2min = r2
+ j = i
+ }
+ }
+
+ # Erase cross and remark with a plus.
+ if (j != 0) {
+$if (datatype == r)
+ call gscur (gd, xin[j], yin[j])
+ pmltype = gstati (gd, G_PMLTYPE)
+ call gseti (gd, G_PMLTYPE, 0)
+ call gmark (gd, xin[j], yin[j], GM_CROSS, 2., 2.)
+ call gseti (gd, G_PMLTYPE, pmltype)
+ call gmark (gd, xin[j], yin[j], GM_PLUS, 2., 2.)
+$else
+ call gscur (gd, real(xin[j]), real(yin[j]))
+ pmltype = gstati (gd, G_PMLTYPE)
+ call gseti (gd, G_PMLTYPE, 0)
+ call gmark (gd, real(xin[j]), real(yin[j]), GM_CROSS, 2., 2.)
+ call gseti (gd, G_PMLTYPE, pmltype)
+ call gmark (gd, real(xin[j]), real(yin[j]), GM_PLUS, 2., 2.)
+$endif
+ wts[j] = userwts[j]
+ }
+ }
+end
+
+
+# GEO_2DELETE -- Delete the residuals.
+
+procedure geo_2delete$t (gd, x, resid, wts, userwts, npts, wx, wy, delete)
+
+pointer gd #I pointer to graphics descriptor
+PIXEL x[ARB] #I reference x values
+PIXEL resid[ARB] #I residuals
+PIXEL wts[ARB] #I weight array
+PIXEL userwts[ARB] #I user weight array
+int npts #I number of points
+real wx #I world x coordinate
+real wy #I world y coordinate
+int delete #I delete point
+
+int i, j, pmltype
+real r2, r2min, x0, y0
+int gstati()
+
+begin
+ # Delete the point.
+ call gctran (gd, wx, wy, wx, wy, 1, 0)
+ r2min = MAX_REAL
+ j = 0
+
+ # Delete or add a point.
+ if (delete == YES) {
+
+ # Find the nearest undeleted point.
+ do i = 1, npts {
+ if (wts[i] <= PIXEL(0.0))
+ next
+$if (datatype == r)
+ call gctran (gd, x[i], resid[i], x0, y0, 1, 0)
+$else
+ call gctran (gd, real(x[i]), real(resid[i]), x0, y0, 1, 0)
+$endif
+ r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+ if (r2 < r2min) {
+ r2min = r2
+ j = i
+ }
+ }
+
+ # Mark the point with a cross and set weight to zero.
+ if (j != 0) {
+$if (datatype == r)
+ call gscur (gd, x[j], resid[j])
+ call gmark (gd, x[j], resid[j], GM_CROSS, 2., 2.)
+$else
+ call gscur (gd, real(x[j]), real(resid[j]))
+ call gmark (gd, real(x[j]), real(resid[j]), GM_CROSS, 2., 2.)
+$endif
+ wts[j] = PIXEL(0.0)
+ }
+
+ } else {
+
+ # Find the nearest deleted point.
+ do i = 1, npts {
+ if (wts[i] > PIXEL(0.0))
+ next
+$if (datatype == r)
+ call gctran (gd, x[i], resid[i], x0, y0, 1, 0)
+$else
+ call gctran (gd, real(x[i]), real(resid[i]), x0, y0, 1, 0)
+$endif
+ r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+ if (r2 < r2min) {
+ r2min = r2
+ j = i
+ }
+ }
+
+ # Erase the cross and remark with a plus.
+ if (j != 0) {
+$if (datatype == r)
+ call gscur (gd, x[j], resid[j])
+ pmltype = gstati (gd, G_PMLTYPE)
+ call gseti (gd, G_PMLTYPE, 0)
+ call gmark (gd, x[j], resid[j], GM_CROSS, 2., 2.)
+ call gseti (gd, G_PMLTYPE, pmltype)
+ call gmark (gd, x[j], resid[j], GM_PLUS, 2., 2.)
+$else
+ call gscur (gd, real(x[j]), real(resid[j]))
+ pmltype = gstati (gd, G_PMLTYPE)
+ call gseti (gd, G_PMLTYPE, 0)
+ call gmark (gd, real(x[j]), real(resid[j]), GM_CROSS, 2., 2.)
+ call gseti (gd, G_PMLTYPE, pmltype)
+ call gmark (gd, real(x[j]), real(resid[j]), GM_PLUS, 2., 2.)
+$endif
+ wts[j] = userwts[j]
+ }
+ }
+end
+
+
+# GEO_1GRAPH - Procedure to graph the distribution of the data in the x-y
+# plane. Rejected points are marked by a ' ' and deleted points are marked
+# by a ' '. The shift in position of the data points are indicated by
+# vectors. Sample fits of constant x and y are marked on the plots.
+
+procedure geo_1graph$t (gd, gt, fit, gfit, xref, yref, xin, yin, wts, npts)
+
+pointer gd #I pointer to the graphics device
+pointer gt #I pointer to the plot descriptor
+pointer fit #I pointer to the geofit structure
+pointer gfit #I pointer to the plot structure
+PIXEL xref[ARB] #I x reference values
+PIXEL yref[ARB] #I y reference values
+PIXEL xin[ARB] #I x values
+PIXEL yin[ARB] #I y values
+PIXEL wts[ARB] #I array of weights
+int npts #I number of points
+
+int i, j
+$if (datatype == d)
+pointer sp, rxin, ryin
+$endif
+
+begin
+ # If previous plot different type don't overplot.
+ if (GG_PLOTTYPE(gfit) != FIT)
+ GG_OVERPLOT(gfit) = NO
+
+ # If not overplottting start new plot.
+ if (GG_OVERPLOT(gfit) == NO) {
+
+ # Set scale and axes.
+ call gclear (gd)
+$if (datatype == r)
+ call gascale (gd, xin, npts, 1)
+ call gascale (gd, yin, npts, 2)
+$else
+ call smark (sp)
+ call salloc (rxin, npts, TY_REAL)
+ call salloc (ryin, npts, TY_REAL)
+ call achtdr (xin, Memr[rxin], npts)
+ call achtdr (yin, Memr[ryin], npts)
+ call gascale (gd, Memr[rxin], npts, 1)
+ call gascale (gd, Memr[ryin], npts, 2)
+ call sfree (sp)
+$endif
+ call gt_swind (gd, gt)
+ call gtlabax (gd, gt)
+
+ # Mark the data and deleted points.
+ do i = 1, npts {
+$if (datatype == r)
+ if (wts[i] == PIXEL(0.0))
+ call gmark (gd, xin[i], yin[i], GM_CROSS, 2., 2.)
+ else
+ call gmark (gd, xin[i], yin[i], GM_PLUS, 2., 2.)
+$else
+ if (wts[i] == PIXEL(0.0))
+ call gmark (gd, real(xin[i]), real(yin[i]), GM_CROSS,
+ 2., 2.)
+ else
+ call gmark (gd, real(xin[i]), real(yin[i]), GM_PLUS,
+ 2., 2.)
+$endif
+ }
+
+ call gflush (gd)
+ }
+
+ # Mark the rejected points.
+ do i = 1, GM_NREJECT(fit) {
+ j = Memi[GM_REJ(fit)+i-1]
+$if (datatype == r)
+ call gmark (gd, xin[j], yin[j], GM_CIRCLE, 2., 2.)
+$else
+ call gmark (gd, real(xin[j]), real(yin[j]), GM_CIRCLE, 2., 2.)
+$endif
+ }
+
+ call gflush (gd)
+
+ # Reset the status flags
+ GG_OVERPLOT(gfit) = NO
+end
+
+
+# GEO_2GRAPH -- Graph the x and y fit residuals versus x or y .
+
+procedure geo_2graph$t (gd, gt, fit, gfit, x, resid, wts, npts)
+
+pointer gd #I pointer to the graphics device
+pointer gt #I pointer to the plot descriptor
+pointer fit #I pointer to geomap structure
+pointer gfit #I pointer to the plot structure
+PIXEL x[ARB] #I x reference values
+PIXEL resid[ARB] #I residual
+PIXEL wts[ARB] #I array of weights
+int npts #I number of points
+
+int i, j
+pointer sp, zero
+$if (datatype == d)
+pointer rxin, ryin
+$endif
+
+begin
+ # Allocate space.
+ call smark (sp)
+ call salloc (zero, npts, TY_REAL)
+ call amovkr (0.0, Memr[zero], npts)
+
+ # Calculate the residuals.
+ if (GG_PLOTTYPE(gfit) == FIT)
+ GG_OVERPLOT(gfit) = NO
+
+ if (GG_OVERPLOT(gfit) == NO) {
+
+ call gclear (gd)
+
+ # Set scale and axes.
+$if (datatype == r)
+ call gascale (gd, x, npts, 1)
+ call gascale (gd, resid, npts, 2)
+$else
+ call salloc (rxin, npts, TY_REAL)
+ call salloc (ryin, npts, TY_REAL)
+ call achtdr (x, Memr[rxin], npts)
+ call achtdr (resid, Memr[ryin], npts)
+ call gascale (gd, Memr[rxin], npts, 1)
+ call gascale (gd, Memr[ryin], npts, 2)
+$endif
+ call gt_swind (gd, gt)
+ call gtlabax (gd, gt)
+
+$if (datatype == r)
+ call gpline (gd, x, Memr[zero], npts)
+$else
+ call gpline (gd, Memr[rxin], Memr[zero], npts)
+$endif
+ }
+
+ # Graph residuals and mark deleted points.
+ if (GG_OVERPLOT(gfit) == NO || GG_NEWFUNCTION(gfit) == YES) {
+ do i = 1, npts {
+$if (datatype == r)
+ if (wts[i] == PIXEL(0.0))
+ call gmark (gd, x[i], resid[i], GM_CROSS, 2., 2.)
+ else
+ call gmark (gd, x[i], resid[i], GM_PLUS, 2., 2.)
+$else
+ if (wts[i] == PIXEL(0.0))
+ call gmark (gd, Memr[rxin+i-1], Memr[ryin+i-1],
+ GM_CROSS, 2., 2.)
+ else
+ call gmark (gd, Memr[rxin+i-1], Memr[ryin+i-1],
+ GM_PLUS, 2., 2.)
+$endif
+ }
+ }
+
+ # plot rejected points
+ if (GM_NREJECT(fit) > 0) {
+ do i = 1, GM_NREJECT(fit) {
+ j = Memi[GM_REJ(fit)+i-1]
+$if (datatype == r)
+ call gmark (gd, x[j], resid[j], GM_CIRCLE, 2., 2.)
+$else
+ call gmark (gd, Memr[rxin+j-1], Memr[ryin+j-1], GM_CIRCLE,
+ 2., 2.)
+$endif
+ }
+ }
+
+ # Reset the status flag.
+ GG_OVERPLOT(gfit) = NO
+
+ call gflush (gd)
+ call sfree (sp)
+end
+
+
+# GEO_CONXY -- Plot a set of default lines of xref = const and yref = const.
+
+procedure geo_conxy$t (gd, fit, sx1, sy1, sx2, sy2)
+
+pointer gd #I graphics file descriptor
+pointer fit #I fit descriptor
+pointer sx1, sy1 #I pointer to the linear x and y surface fits
+pointer sx2, sy2 #I pointer to the linear x and y surface fits
+
+int i
+pointer sp, xtemp, ytemp, xfit1, yfit1, xfit2, yfit2
+$if (datatype == d)
+pointer xbuf, ybuf
+$endif
+PIXEL xint, yint, dx, dy
+
+begin
+ # allocate temporary space
+ call smark (sp)
+ call salloc (xtemp, NGRAPH, TY_PIXEL)
+ call salloc (ytemp, NGRAPH, TY_PIXEL)
+ call salloc (xfit1, NGRAPH, TY_PIXEL)
+ call salloc (yfit1, NGRAPH, TY_PIXEL)
+ call salloc (xfit2, NGRAPH, TY_PIXEL)
+ call salloc (yfit2, NGRAPH, TY_PIXEL)
+$if (datatype == d)
+ call salloc (xbuf, NGRAPH, TY_REAL)
+ call salloc (ybuf, NGRAPH, TY_REAL)
+$endif
+
+ # Calculate intervals in x and y.
+ dx = (GM_XMAX(fit) - GM_XMIN(fit)) / NINTERVALS
+ dy = (GM_YMAX(fit) - GM_YMIN(fit)) / (NGRAPH - 1)
+
+ # Set up an array of y values.
+ Mem$t[ytemp] = GM_YMIN(fit)
+ do i = 2, NGRAPH
+ Mem$t[ytemp+i-1] = Mem$t[ytemp+i-2] + dy
+
+ # Mark lines of constant x.
+ xint = GM_XMIN(fit)
+ for (i = 1; i <= NINTERVALS + 1; i = i + 1) {
+
+ # Set the x value.
+ call amovk$t (xint, Mem$t[xtemp], NGRAPH)
+
+ # X fit.
+$if (datatype == r)
+ call gsvector (sx1, Mem$t[xtemp], Mem$t[ytemp], Mem$t[xfit1],
+ NGRAPH)
+$else
+ call dgsvector (sx1, Mem$t[xtemp], Mem$t[ytemp], Mem$t[xfit1],
+ NGRAPH)
+$endif
+ if (sx2 != NULL) {
+$if (datatype == r)
+ call gsvector (sx2, Mem$t[xtemp], Mem$t[ytemp], Mem$t[xfit2],
+ NGRAPH)
+$else
+ call dgsvector (sx2, Mem$t[xtemp], Mem$t[ytemp], Mem$t[xfit2],
+ NGRAPH)
+$endif
+ call aadd$t (Mem$t[xfit1], Mem$t[xfit2], Mem$t[xfit1], NGRAPH)
+ }
+
+ # Y fit.
+$if (datatype == r)
+ call gsvector (sy1, Mem$t[xtemp], Mem$t[ytemp], Mem$t[yfit1],
+ NGRAPH)
+$else
+ call dgsvector (sy1, Mem$t[xtemp], Mem$t[ytemp], Mem$t[yfit1],
+ NGRAPH)
+$endif
+ if (sy2 != NULL) {
+$if (datatype == r)
+ call gsvector (sy2, Mem$t[xtemp], Mem$t[ytemp], Mem$t[yfit2],
+ NGRAPH)
+$else
+ call dgsvector (sy2, Mem$t[xtemp], Mem$t[ytemp], Mem$t[yfit2],
+ NGRAPH)
+$endif
+ call aadd$t (Mem$t[yfit1], Mem$t[yfit2], Mem$t[yfit1], NGRAPH)
+ }
+
+ # Plot line of constant x.
+$if (datatype == r)
+ call gpline (gd, Memr[xfit1], Memr[yfit1], NGRAPH)
+$else
+ call achtdr (Memd[xfit1], Memr[xbuf], NGRAPH)
+ call achtdr (Memd[yfit1], Memr[ybuf], NGRAPH)
+ call gpline (gd, Memr[xbuf], Memr[ybuf], NGRAPH)
+$endif
+
+ # Update the x value.
+ xint = xint + dx
+ }
+
+ call gflush (gd)
+
+ # Calculate x and y intervals.
+ dx = (GM_XMAX(fit) - GM_XMIN(fit)) / (NGRAPH - 1)
+ dy = (GM_YMAX(fit) - GM_YMIN(fit)) / NINTERVALS
+
+ # Set up array of x values.
+ Mem$t[xtemp] = GM_XMIN(fit)
+ do i = 2, NGRAPH
+ Mem$t[xtemp+i-1] = Mem$t[xtemp+i-2] + dx
+
+ # Mark lines of constant y.
+ yint = GM_YMIN(fit)
+ for (i = 1; i <= NINTERVALS + 1; i = i + 1) {
+
+ # set the y value
+ call amovk$t (yint, Mem$t[ytemp], NGRAPH)
+
+ # X fit.
+$if (datatype == r)
+ call gsvector (sx1, Mem$t[xtemp], Mem$t[ytemp], Mem$t[xfit1],
+ NGRAPH)
+$else
+ call dgsvector (sx1, Mem$t[xtemp], Mem$t[ytemp], Mem$t[xfit1],
+ NGRAPH)
+$endif
+ if (sx2 != NULL) {
+$if (datatype == r)
+ call gsvector (sx2, Mem$t[xtemp], Mem$t[ytemp], Mem$t[xfit2],
+ NGRAPH)
+$else
+ call dgsvector (sx2, Mem$t[xtemp], Mem$t[ytemp], Mem$t[xfit2],
+ NGRAPH)
+$endif
+ call aadd$t (Mem$t[xfit1], Mem$t[xfit2], Mem$t[xfit1], NGRAPH)
+ }
+
+
+ # Y fit.
+$if (datatype == r)
+ call gsvector (sy1, Mem$t[xtemp], Mem$t[ytemp], Mem$t[yfit1],
+ NGRAPH)
+$else
+ call dgsvector (sy1, Mem$t[xtemp], Mem$t[ytemp], Mem$t[yfit1],
+ NGRAPH)
+$endif
+ if (sy2 != NULL) {
+$if (datatype == r)
+ call gsvector (sy2, Mem$t[xtemp], Mem$t[ytemp], Mem$t[yfit2],
+ NGRAPH)
+$else
+ call dgsvector (sy2, Mem$t[xtemp], Mem$t[ytemp], Mem$t[yfit2],
+ NGRAPH)
+$endif
+ call aadd$t (Mem$t[yfit1], Mem$t[yfit2], Mem$t[yfit1], NGRAPH)
+ }
+
+ # Plot line of constant y.
+$if (datatype == r)
+ call gpline (gd, Memr[xfit1], Memr[yfit1], NGRAPH)
+$else
+ call achtdr (Memd[xfit1], Memr[xbuf], NGRAPH)
+ call achtdr (Memd[yfit1], Memr[ybuf], NGRAPH)
+ call gpline (gd, Memr[xbuf], Memr[ybuf], NGRAPH)
+$endif
+
+ # Update the y value.
+ yint = yint + dy
+ }
+
+ call gflush (gd)
+
+ call sfree (sp)
+end
+
+
+# GEO_LXY -- Draw a line of constant x-y.
+
+procedure geo_lxy$t (gd, fit, sx1, sy1, sx2, sy2, xref, yref, xin, yin, npts,
+ wx, wy)
+
+pointer gd #I pointer to graphics descriptor
+pointer fit #I pointer to the fit parameters
+pointer sx1 #I pointer to the linear x fit
+pointer sy1 #I pointer to the linear y fit
+pointer sx2 #I pointer to the higher order x fit
+pointer sy2 #I pointer to the higher order y fit
+PIXEL xref[ARB] #I x reference values
+PIXEL yref[ARB] #I y reference values
+PIXEL xin[ARB] #I x input values
+PIXEL yin[ARB] #I y input values
+int npts #I number of data points
+real wx, wy #I x and y world coordinates
+
+int i, j
+pointer sp, xtemp, ytemp, xfit1, yfit1, xfit2, yfit2
+$if (datatype == d)
+pointer xbuf, ybuf
+$endif
+real x0, y0, r2, r2min
+PIXEL delta, deltax, deltay
+$if (datatype == r)
+real gseval()
+$else
+double dgseval()
+$endif
+
+begin
+ # Transform world coordinates.
+ call gctran (gd, wx, wy, wx, wy, 1, 0)
+ r2min = MAX_REAL
+ j = 0
+
+ # Find the nearest data point.
+ do i = 1, npts {
+$if (datatype == r)
+ call gctran (gd, xin[i], yin[i], x0, y0, 1, 0)
+$else
+ call gctran (gd, real(xin[i]), real(yin[i]), x0, y0, 1, 0)
+$endif
+ r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+ if (r2 < r2min) {
+ r2min = r2
+ j = i
+ }
+ }
+
+ # Fit the line
+ if (j != 0) {
+
+ # Allocate temporary space.
+ call smark (sp)
+ call salloc (xtemp, NGRAPH, TY_PIXEL)
+ call salloc (ytemp, NGRAPH, TY_PIXEL)
+ call salloc (xfit1, NGRAPH, TY_PIXEL)
+ call salloc (yfit1, NGRAPH, TY_PIXEL)
+ call salloc (xfit2, NGRAPH, TY_PIXEL)
+ call salloc (yfit2, NGRAPH, TY_PIXEL)
+$if (datatype == d)
+ call salloc (xbuf, NGRAPH, TY_REAL)
+ call salloc (ybuf, NGRAPH, TY_REAL)
+$endif
+
+ # Compute the deltas.
+$if (datatype == r)
+ deltax = xin[j] - gseval (sx1, xref[j], yref[j])
+ if (sx2 != NULL)
+ deltax = deltax - gseval (sx2, xref[j], yref[j])
+ deltay = yin[j] - gseval (sy1, xref[j], yref[j])
+ if (sy2 != NULL)
+ deltay = deltay - gseval (sy2, xref[j], yref[j])
+$else
+ deltax = xin[j] - dgseval (sx1, xref[j], yref[j])
+ if (sx2 != NULL)
+ deltax = deltax - dgseval (sx2, xref[j], yref[j])
+ deltay = yin[j] - dgseval (sy1, xref[j], yref[j])
+ if (sy2 != NULL)
+ deltay = deltay - dgseval (sy2, xref[j], yref[j])
+$endif
+
+ # Set up line of constant x.
+ call amovk$t (xref[j], Mem$t[xtemp], NGRAPH)
+ delta = (GM_YMAX(fit) - GM_YMIN(fit)) / (NGRAPH - 1)
+ Mem$t[ytemp] = GM_YMIN(fit)
+ do i = 2, NGRAPH
+ Mem$t[ytemp+i-1] = Mem$t[ytemp+i-2] + delta
+
+ # X solution.
+$if (datatype == r)
+ call gsvector (sx1, Mem$t[xtemp], Mem$t[ytemp], Mem$t[xfit1],
+ NGRAPH)
+$else
+ call dgsvector (sx1, Mem$t[xtemp], Mem$t[ytemp], Mem$t[xfit1],
+ NGRAPH)
+$endif
+ if (sx2 != NULL) {
+$if (datatype == r)
+ call gsvector (sx2, Mem$t[xtemp], Mem$t[ytemp], Mem$t[xfit2],
+ NGRAPH)
+$else
+ call dgsvector (sx2, Mem$t[xtemp], Mem$t[ytemp], Mem$t[xfit2],
+ NGRAPH)
+$endif
+ call aadd$t (Mem$t[xfit1], Mem$t[xfit2], Mem$t[xfit1], NGRAPH)
+ }
+ call aaddk$t (Mem$t[xfit1], deltax, Mem$t[xfit1], NGRAPH)
+
+ # Y solution.
+$if (datatype == r)
+ call gsvector (sy1, Mem$t[xtemp], Mem$t[ytemp], Mem$t[yfit1],
+ NGRAPH)
+$else
+ call dgsvector (sy1, Mem$t[xtemp], Mem$t[ytemp], Mem$t[yfit1],
+ NGRAPH)
+$endif
+ if (sy2 != NULL) {
+$if (datatype == r)
+ call gsvector (sy2, Mem$t[xtemp], Mem$t[ytemp], Mem$t[yfit2],
+ NGRAPH)
+$else
+ call dgsvector (sy2, Mem$t[xtemp], Mem$t[ytemp], Mem$t[yfit2],
+ NGRAPH)
+$endif
+ call aadd$t (Mem$t[yfit1], Mem$t[yfit2], Mem$t[yfit1], NGRAPH)
+ }
+ call aaddk$t (Mem$t[yfit1], deltay, Mem$t[yfit1], NGRAPH)
+
+ # Plot line of constant x.
+$if (datatype == r)
+ call gpline (gd, Memr[xfit1], Memr[yfit1], NGRAPH)
+$else
+ call achtdr (Memd[xfit1], Memr[xbuf], NGRAPH)
+ call achtdr (Memd[yfit1], Memr[ybuf], NGRAPH)
+ call gpline (gd, Memr[xbuf], Memr[ybuf], NGRAPH)
+$endif
+ call gflush (gd)
+
+ # Set up line of constant y.
+ call amovk$t (yref[j], Mem$t[ytemp], NGRAPH)
+ delta = (GM_XMAX(fit) - GM_XMIN(fit)) / (NGRAPH - 1)
+ Mem$t[xtemp] = GM_XMIN(fit)
+ do i = 2, NGRAPH
+ Mem$t[xtemp+i-1] = Mem$t[xtemp+i-2] + delta
+
+ # X fit.
+$if (datatype == r)
+ call gsvector (sx1, Mem$t[xtemp], Mem$t[ytemp], Mem$t[xfit1],
+ NGRAPH)
+$else
+ call dgsvector (sx1, Mem$t[xtemp], Mem$t[ytemp], Mem$t[xfit1],
+ NGRAPH)
+$endif
+ if (sx2 != NULL) {
+$if (datatype == r)
+ call gsvector (sx2, Mem$t[xtemp], Mem$t[ytemp], Mem$t[xfit2],
+ NGRAPH)
+$else
+ call dgsvector (sx2, Mem$t[xtemp], Mem$t[ytemp], Mem$t[xfit2],
+ NGRAPH)
+$endif
+ call aadd$t (Mem$t[xfit1], Mem$t[xfit2], Mem$t[xfit1], NGRAPH)
+ }
+ call aaddk$t (Mem$t[xfit1], deltax, Mem$t[xfit1], NGRAPH)
+
+ # Y fit.
+$if (datatype == r)
+ call gsvector (sy1, Mem$t[xtemp], Mem$t[ytemp], Mem$t[yfit1],
+ NGRAPH)
+$else
+ call dgsvector (sy1, Mem$t[xtemp], Mem$t[ytemp], Mem$t[yfit1],
+ NGRAPH)
+$endif
+ if (sy2 != NULL) {
+$if (datatype == r)
+ call gsvector (sy2, Mem$t[xtemp], Mem$t[ytemp], Mem$t[yfit2],
+ NGRAPH)
+$else
+ call dgsvector (sy2, Mem$t[xtemp], Mem$t[ytemp], Mem$t[yfit2],
+ NGRAPH)
+$endif
+ call aadd$t (Mem$t[yfit1], Mem$t[yfit2], Mem$t[yfit1], NGRAPH)
+ }
+ call aaddk$t (Mem$t[yfit1], deltay, Mem$t[yfit1], NGRAPH)
+
+ # Plot line of constant y.
+$if (datatype == r)
+ call gpline (gd, Memr[xfit1], Memr[yfit1], NGRAPH)
+$else
+ call achtdr (Memd[xfit1], Memr[xbuf], NGRAPH)
+ call achtdr (Memd[yfit1], Memr[ybuf], NGRAPH)
+ call gpline (gd, Memr[xbuf], Memr[ybuf], NGRAPH)
+$endif
+ call gflush (gd)
+
+ # Free space.
+ call sfree (sp)
+ }
+end
+
+
+# GEO_GCOEFF -- Print the coefficents of the linear portion of the
+# fit, xshift, yshift,
+
+procedure geo_gcoeff$t (sx, sy, xshift, yshift, a, b, c, d)
+
+pointer sx #I pointer to the x surface fit
+pointer sy #I pointer to the y surface fit
+PIXEL xshift #O output x shift
+PIXEL yshift #O output y shift
+PIXEL a #O output x coefficient of x fit
+PIXEL b #O output y coefficient of x fit
+PIXEL c #O output x coefficient of y fit
+PIXEL d #O output y coefficient of y fit
+
+int nxxcoeff, nxycoeff, nyxcoeff, nyycoeff
+pointer sp, xcoeff, ycoeff
+PIXEL xxrange, xyrange, xxmaxmin, xymaxmin
+PIXEL yxrange, yyrange, yxmaxmin, yymaxmin
+
+$if (datatype == r)
+int gsgeti()
+real gsgetr()
+$else
+int dgsgeti()
+double dgsgetd()
+$endif
+
+begin
+ # Allocate working space.
+ call smark (sp)
+$if (datatype == r)
+ call salloc (xcoeff, gsgeti (sx, GSNCOEFF), TY_PIXEL)
+ call salloc (ycoeff, gsgeti (sy, GSNCOEFF), TY_PIXEL)
+$else
+ call salloc (xcoeff, dgsgeti (sx, GSNCOEFF), TY_PIXEL)
+ call salloc (ycoeff, dgsgeti (sy, GSNCOEFF), TY_PIXEL)
+$endif
+
+ # Get coefficients and numbers of coefficients.
+$if (datatype == r)
+ call gscoeff (sx, Mem$t[xcoeff], nxxcoeff)
+ call gscoeff (sy, Mem$t[ycoeff], nyycoeff)
+ nxxcoeff = gsgeti (sx, GSNXCOEFF)
+ nxycoeff = gsgeti (sx, GSNYCOEFF)
+ nyxcoeff = gsgeti (sy, GSNXCOEFF)
+ nyycoeff = gsgeti (sy, GSNYCOEFF)
+$else
+ call dgscoeff (sx, Mem$t[xcoeff], nxxcoeff)
+ call dgscoeff (sy, Mem$t[ycoeff], nyycoeff)
+ nxxcoeff = dgsgeti (sx, GSNXCOEFF)
+ nxycoeff = dgsgeti (sx, GSNYCOEFF)
+ nyxcoeff = dgsgeti (sy, GSNXCOEFF)
+ nyycoeff = dgsgeti (sy, GSNYCOEFF)
+$endif
+
+ # Get the data range.
+$if (datatype == r)
+ if (gsgeti (sx, GSTYPE) != GS_POLYNOMIAL) {
+ xxrange = (gsgetr (sx, GSXMAX) - gsgetr (sx, GSXMIN)) / 2.0
+ xxmaxmin = - (gsgetr (sx, GSXMAX) + gsgetr (sx, GSXMIN)) / 2.0
+ xyrange = (gsgetr (sx, GSYMAX) - gsgetr (sx, GSYMIN)) / 2.0
+ xymaxmin = - (gsgetr (sx, GSYMAX) + gsgetr (sx, GSYMIN)) / 2.0
+$else
+ if (dgsgeti (sx, GSTYPE) != GS_POLYNOMIAL) {
+ xxrange = (dgsgetd (sx, GSXMAX) - dgsgetd (sx, GSXMIN)) / 2.0d0
+ xxmaxmin = - (dgsgetd (sx, GSXMAX) + dgsgetd (sx, GSXMIN)) / 2.0d0
+ xyrange = (dgsgetd (sx, GSYMAX) - dgsgetd (sx, GSYMIN)) / 2.0d0
+ xymaxmin = - (dgsgetd (sx, GSYMAX) + dgsgetd (sx, GSYMIN)) / 2.0d0
+$endif
+ } else {
+ xxrange = PIXEL(1.0)
+ xxmaxmin = PIXEL(0.0)
+ xyrange = PIXEL(1.0)
+ xymaxmin = PIXEL(0.0)
+ }
+
+$if (datatype == r)
+ if (gsgeti (sy, GSTYPE) != GS_POLYNOMIAL) {
+ yxrange = (gsgetr (sy, GSXMAX) - gsgetr (sy, GSXMIN)) / 2.0
+ yxmaxmin = - (gsgetr (sy, GSXMAX) + gsgetr (sy, GSXMIN)) / 2.0
+ yyrange = (gsgetr (sy, GSYMAX) - gsgetr (sy, GSYMIN)) / 2.0
+ yymaxmin = - (gsgetr (sy, GSYMAX) + gsgetr (sy, GSYMIN)) / 2.0
+$else
+ if (dgsgeti (sy, GSTYPE) != GS_POLYNOMIAL) {
+ yxrange = (dgsgetd (sy, GSXMAX) - dgsgetd (sy, GSXMIN)) / 2.0d0
+ yxmaxmin = - (dgsgetd (sy, GSXMAX) + dgsgetd (sy, GSXMIN)) / 2.0d0
+ yyrange = (dgsgetd (sy, GSYMAX) - dgsgetd (sy, GSYMIN)) / 2.0d0
+ yymaxmin = - (dgsgetd (sy, GSYMAX) + dgsgetd (sy, GSYMIN)) / 2.0d0
+$endif
+ } else {
+ yxrange = PIXEL(1.0)
+ yxmaxmin = PIXEL(0.0)
+ yyrange = PIXEL(1.0)
+ yymaxmin = PIXEL(0.0)
+ }
+
+ # Get the shifts.
+ xshift = Mem$t[xcoeff] + Mem$t[xcoeff+1] * xxmaxmin / xxrange +
+ Mem$t[xcoeff+2] * xymaxmin / xyrange
+ yshift = Mem$t[ycoeff] + Mem$t[ycoeff+1] * yxmaxmin / yxrange +
+ Mem$t[ycoeff+2] * yymaxmin / yyrange
+
+ # Get the rotation and scaling parameters and correct for normalization.
+ if (nxxcoeff > 1)
+ a = Mem$t[xcoeff+1] / xxrange
+ else
+ a = PIXEL(0.0)
+ if (nxycoeff > 1)
+ b = Mem$t[xcoeff+nxxcoeff] / xyrange
+ else
+ b = PIXEL(0.0)
+ if (nyxcoeff > 1)
+ c = Mem$t[ycoeff+1] / yxrange
+ else
+ c = PIXEL(0.0)
+ if (nyycoeff > 1)
+ d = Mem$t[ycoeff+nyxcoeff] / yyrange
+ else
+ d = PIXEL(0.0)
+
+ call sfree (sp)
+end
+
+$endfor
generated by cgit (git 2.34.1) at 2025-09-20 05:12:48 -0400