aboutsummaryrefslogtreecommitdiff
path: root/pkg/images/lib/geograph.x
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/images/lib/geograph.x')
-rw-r--r--pkg/images/lib/geograph.x1740
1 files changed, 1740 insertions, 0 deletions
diff --git a/pkg/images/lib/geograph.x b/pkg/images/lib/geograph.x
new file mode 100644
index 00000000..6597311a
--- /dev/null
+++ b/pkg/images/lib/geograph.x
@@ -0,0 +1,1740 @@
+# 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
+
+
+
+# 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
+
+
+
+
+
+
+# GEO_1DELETE -- Delete a point from the fit.
+
+procedure geo_1deleter (gd, xin, yin, wts, userwts, npts, wx, wy, delete)
+
+pointer gd #I pointer to graphics descriptor
+real xin[ARB] #I x array
+real yin[ARB] #I y array
+real wts[ARB] #I array of weights
+real 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] <= real(0.0))
+ next
+ call gctran (gd, xin[i], yin[i], x0, y0, 1, 0)
+ r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+ if (r2 < r2min) {
+ r2min = r2
+ j = i
+ }
+ }
+
+ # Mark point and set weights to 0.
+ if (j != 0) {
+ call gscur (gd, xin[j], yin[j])
+ call gmark (gd, xin[j], yin[j], GM_CROSS, 2., 2.)
+ wts[j] = real(0.0)
+ }
+
+ } else {
+
+ # Search for the nearest deleted point.
+ do i = 1, npts {
+ if (wts[i] > real(0.0))
+ next
+ call gctran (gd, xin[i], yin[i], x0, y0, 1, 0)
+ r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+ if (r2 < r2min) {
+ r2min = r2
+ j = i
+ }
+ }
+
+ # Erase cross and remark with a plus.
+ if (j != 0) {
+ 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.)
+ wts[j] = userwts[j]
+ }
+ }
+end
+
+
+# GEO_2DELETE -- Delete the residuals.
+
+procedure geo_2deleter (gd, x, resid, wts, userwts, npts, wx, wy, delete)
+
+pointer gd #I pointer to graphics descriptor
+real x[ARB] #I reference x values
+real resid[ARB] #I residuals
+real wts[ARB] #I weight array
+real 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] <= real(0.0))
+ next
+ call gctran (gd, x[i], resid[i], x0, y0, 1, 0)
+ 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) {
+ call gscur (gd, x[j], resid[j])
+ call gmark (gd, x[j], resid[j], GM_CROSS, 2., 2.)
+ wts[j] = real(0.0)
+ }
+
+ } else {
+
+ # Find the nearest deleted point.
+ do i = 1, npts {
+ if (wts[i] > real(0.0))
+ next
+ call gctran (gd, x[i], resid[i], x0, y0, 1, 0)
+ 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) {
+ 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.)
+ 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_1graphr (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
+real xref[ARB] #I x reference values
+real yref[ARB] #I y reference values
+real xin[ARB] #I x values
+real yin[ARB] #I y values
+real wts[ARB] #I array of weights
+int npts #I number of points
+
+int i, j
+
+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)
+ call gascale (gd, xin, npts, 1)
+ call gascale (gd, yin, npts, 2)
+ call gt_swind (gd, gt)
+ call gtlabax (gd, gt)
+
+ # Mark the data and deleted points.
+ do i = 1, npts {
+ if (wts[i] == real(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.)
+ }
+
+ call gflush (gd)
+ }
+
+ # Mark the rejected points.
+ do i = 1, GM_NREJECT(fit) {
+ j = Memi[GM_REJ(fit)+i-1]
+ call gmark (gd, xin[j], yin[j], GM_CIRCLE, 2., 2.)
+ }
+
+ 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_2graphr (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
+real x[ARB] #I x reference values
+real resid[ARB] #I residual
+real wts[ARB] #I array of weights
+int npts #I number of points
+
+int i, j
+pointer sp, zero
+
+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.
+ call gascale (gd, x, npts, 1)
+ call gascale (gd, resid, npts, 2)
+ call gt_swind (gd, gt)
+ call gtlabax (gd, gt)
+
+ call gpline (gd, x, Memr[zero], npts)
+ }
+
+ # Graph residuals and mark deleted points.
+ if (GG_OVERPLOT(gfit) == NO || GG_NEWFUNCTION(gfit) == YES) {
+ do i = 1, npts {
+ if (wts[i] == real(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.)
+ }
+ }
+
+ # plot rejected points
+ if (GM_NREJECT(fit) > 0) {
+ do i = 1, GM_NREJECT(fit) {
+ j = Memi[GM_REJ(fit)+i-1]
+ call gmark (gd, x[j], resid[j], GM_CIRCLE, 2., 2.)
+ }
+ }
+
+ # 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_conxyr (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
+real xint, yint, dx, dy
+
+begin
+ # allocate temporary space
+ call smark (sp)
+ call salloc (xtemp, NGRAPH, TY_REAL)
+ call salloc (ytemp, NGRAPH, TY_REAL)
+ call salloc (xfit1, NGRAPH, TY_REAL)
+ call salloc (yfit1, NGRAPH, TY_REAL)
+ call salloc (xfit2, NGRAPH, TY_REAL)
+ call salloc (yfit2, NGRAPH, TY_REAL)
+
+ # 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.
+ Memr[ytemp] = GM_YMIN(fit)
+ do i = 2, NGRAPH
+ Memr[ytemp+i-1] = Memr[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 amovkr (xint, Memr[xtemp], NGRAPH)
+
+ # X fit.
+ call gsvector (sx1, Memr[xtemp], Memr[ytemp], Memr[xfit1],
+ NGRAPH)
+ if (sx2 != NULL) {
+ call gsvector (sx2, Memr[xtemp], Memr[ytemp], Memr[xfit2],
+ NGRAPH)
+ call aaddr (Memr[xfit1], Memr[xfit2], Memr[xfit1], NGRAPH)
+ }
+
+ # Y fit.
+ call gsvector (sy1, Memr[xtemp], Memr[ytemp], Memr[yfit1],
+ NGRAPH)
+ if (sy2 != NULL) {
+ call gsvector (sy2, Memr[xtemp], Memr[ytemp], Memr[yfit2],
+ NGRAPH)
+ call aaddr (Memr[yfit1], Memr[yfit2], Memr[yfit1], NGRAPH)
+ }
+
+ # Plot line of constant x.
+ call gpline (gd, Memr[xfit1], Memr[yfit1], NGRAPH)
+
+ # 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.
+ Memr[xtemp] = GM_XMIN(fit)
+ do i = 2, NGRAPH
+ Memr[xtemp+i-1] = Memr[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 amovkr (yint, Memr[ytemp], NGRAPH)
+
+ # X fit.
+ call gsvector (sx1, Memr[xtemp], Memr[ytemp], Memr[xfit1],
+ NGRAPH)
+ if (sx2 != NULL) {
+ call gsvector (sx2, Memr[xtemp], Memr[ytemp], Memr[xfit2],
+ NGRAPH)
+ call aaddr (Memr[xfit1], Memr[xfit2], Memr[xfit1], NGRAPH)
+ }
+
+
+ # Y fit.
+ call gsvector (sy1, Memr[xtemp], Memr[ytemp], Memr[yfit1],
+ NGRAPH)
+ if (sy2 != NULL) {
+ call gsvector (sy2, Memr[xtemp], Memr[ytemp], Memr[yfit2],
+ NGRAPH)
+ call aaddr (Memr[yfit1], Memr[yfit2], Memr[yfit1], NGRAPH)
+ }
+
+ # Plot line of constant y.
+ call gpline (gd, Memr[xfit1], Memr[yfit1], NGRAPH)
+
+ # 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_lxyr (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
+real xref[ARB] #I x reference values
+real yref[ARB] #I y reference values
+real xin[ARB] #I x input values
+real 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
+real x0, y0, r2, r2min
+real delta, deltax, deltay
+real gseval()
+
+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 {
+ call gctran (gd, xin[i], yin[i], x0, y0, 1, 0)
+ 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_REAL)
+ call salloc (ytemp, NGRAPH, TY_REAL)
+ call salloc (xfit1, NGRAPH, TY_REAL)
+ call salloc (yfit1, NGRAPH, TY_REAL)
+ call salloc (xfit2, NGRAPH, TY_REAL)
+ call salloc (yfit2, NGRAPH, TY_REAL)
+
+ # Compute the deltas.
+ 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])
+
+ # Set up line of constant x.
+ call amovkr (xref[j], Memr[xtemp], NGRAPH)
+ delta = (GM_YMAX(fit) - GM_YMIN(fit)) / (NGRAPH - 1)
+ Memr[ytemp] = GM_YMIN(fit)
+ do i = 2, NGRAPH
+ Memr[ytemp+i-1] = Memr[ytemp+i-2] + delta
+
+ # X solution.
+ call gsvector (sx1, Memr[xtemp], Memr[ytemp], Memr[xfit1],
+ NGRAPH)
+ if (sx2 != NULL) {
+ call gsvector (sx2, Memr[xtemp], Memr[ytemp], Memr[xfit2],
+ NGRAPH)
+ call aaddr (Memr[xfit1], Memr[xfit2], Memr[xfit1], NGRAPH)
+ }
+ call aaddkr (Memr[xfit1], deltax, Memr[xfit1], NGRAPH)
+
+ # Y solution.
+ call gsvector (sy1, Memr[xtemp], Memr[ytemp], Memr[yfit1],
+ NGRAPH)
+ if (sy2 != NULL) {
+ call gsvector (sy2, Memr[xtemp], Memr[ytemp], Memr[yfit2],
+ NGRAPH)
+ call aaddr (Memr[yfit1], Memr[yfit2], Memr[yfit1], NGRAPH)
+ }
+ call aaddkr (Memr[yfit1], deltay, Memr[yfit1], NGRAPH)
+
+ # Plot line of constant x.
+ call gpline (gd, Memr[xfit1], Memr[yfit1], NGRAPH)
+ call gflush (gd)
+
+ # Set up line of constant y.
+ call amovkr (yref[j], Memr[ytemp], NGRAPH)
+ delta = (GM_XMAX(fit) - GM_XMIN(fit)) / (NGRAPH - 1)
+ Memr[xtemp] = GM_XMIN(fit)
+ do i = 2, NGRAPH
+ Memr[xtemp+i-1] = Memr[xtemp+i-2] + delta
+
+ # X fit.
+ call gsvector (sx1, Memr[xtemp], Memr[ytemp], Memr[xfit1],
+ NGRAPH)
+ if (sx2 != NULL) {
+ call gsvector (sx2, Memr[xtemp], Memr[ytemp], Memr[xfit2],
+ NGRAPH)
+ call aaddr (Memr[xfit1], Memr[xfit2], Memr[xfit1], NGRAPH)
+ }
+ call aaddkr (Memr[xfit1], deltax, Memr[xfit1], NGRAPH)
+
+ # Y fit.
+ call gsvector (sy1, Memr[xtemp], Memr[ytemp], Memr[yfit1],
+ NGRAPH)
+ if (sy2 != NULL) {
+ call gsvector (sy2, Memr[xtemp], Memr[ytemp], Memr[yfit2],
+ NGRAPH)
+ call aaddr (Memr[yfit1], Memr[yfit2], Memr[yfit1], NGRAPH)
+ }
+ call aaddkr (Memr[yfit1], deltay, Memr[yfit1], NGRAPH)
+
+ # Plot line of constant y.
+ call gpline (gd, Memr[xfit1], Memr[yfit1], NGRAPH)
+ 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_gcoeffr (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
+real xshift #O output x shift
+real yshift #O output y shift
+real a #O output x coefficient of x fit
+real b #O output y coefficient of x fit
+real c #O output x coefficient of y fit
+real d #O output y coefficient of y fit
+
+int nxxcoeff, nxycoeff, nyxcoeff, nyycoeff
+pointer sp, xcoeff, ycoeff
+real xxrange, xyrange, xxmaxmin, xymaxmin
+real yxrange, yyrange, yxmaxmin, yymaxmin
+
+int gsgeti()
+real gsgetr()
+
+begin
+ # Allocate working space.
+ call smark (sp)
+ call salloc (xcoeff, gsgeti (sx, GSNCOEFF), TY_REAL)
+ call salloc (ycoeff, gsgeti (sy, GSNCOEFF), TY_REAL)
+
+ # Get coefficients and numbers of coefficients.
+ call gscoeff (sx, Memr[xcoeff], nxxcoeff)
+ call gscoeff (sy, Memr[ycoeff], nyycoeff)
+ nxxcoeff = gsgeti (sx, GSNXCOEFF)
+ nxycoeff = gsgeti (sx, GSNYCOEFF)
+ nyxcoeff = gsgeti (sy, GSNXCOEFF)
+ nyycoeff = gsgeti (sy, GSNYCOEFF)
+
+ # Get the data range.
+ 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 {
+ xxrange = real(1.0)
+ xxmaxmin = real(0.0)
+ xyrange = real(1.0)
+ xymaxmin = real(0.0)
+ }
+
+ 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 {
+ yxrange = real(1.0)
+ yxmaxmin = real(0.0)
+ yyrange = real(1.0)
+ yymaxmin = real(0.0)
+ }
+
+ # Get the shifts.
+ xshift = Memr[xcoeff] + Memr[xcoeff+1] * xxmaxmin / xxrange +
+ Memr[xcoeff+2] * xymaxmin / xyrange
+ yshift = Memr[ycoeff] + Memr[ycoeff+1] * yxmaxmin / yxrange +
+ Memr[ycoeff+2] * yymaxmin / yyrange
+
+ # Get the rotation and scaling parameters and correct for normalization.
+ if (nxxcoeff > 1)
+ a = Memr[xcoeff+1] / xxrange
+ else
+ a = real(0.0)
+ if (nxycoeff > 1)
+ b = Memr[xcoeff+nxxcoeff] / xyrange
+ else
+ b = real(0.0)
+ if (nyxcoeff > 1)
+ c = Memr[ycoeff+1] / yxrange
+ else
+ c = real(0.0)
+ if (nyycoeff > 1)
+ d = Memr[ycoeff+nyxcoeff] / yyrange
+ else
+ d = real(0.0)
+
+ call sfree (sp)
+end
+
+
+
+# GEO_1DELETE -- Delete a point from the fit.
+
+procedure geo_1deleted (gd, xin, yin, wts, userwts, npts, wx, wy, delete)
+
+pointer gd #I pointer to graphics descriptor
+double xin[ARB] #I x array
+double yin[ARB] #I y array
+double wts[ARB] #I array of weights
+double 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] <= double(0.0))
+ next
+ call gctran (gd, real (xin[i]), real (yin[i]), x0, y0, 1, 0)
+ r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+ if (r2 < r2min) {
+ r2min = r2
+ j = i
+ }
+ }
+
+ # Mark point and set weights to 0.
+ if (j != 0) {
+ call gscur (gd, real(xin[j]), real(yin[j]))
+ call gmark (gd, real(xin[j]), real(yin[j]), GM_CROSS, 2., 2.)
+ wts[j] = double(0.0)
+ }
+
+ } else {
+
+ # Search for the nearest deleted point.
+ do i = 1, npts {
+ if (wts[i] > double(0.0))
+ next
+ call gctran (gd, real(xin[i]), real(yin[i]), x0, y0, 1, 0)
+ r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+ if (r2 < r2min) {
+ r2min = r2
+ j = i
+ }
+ }
+
+ # Erase cross and remark with a plus.
+ if (j != 0) {
+ 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.)
+ wts[j] = userwts[j]
+ }
+ }
+end
+
+
+# GEO_2DELETE -- Delete the residuals.
+
+procedure geo_2deleted (gd, x, resid, wts, userwts, npts, wx, wy, delete)
+
+pointer gd #I pointer to graphics descriptor
+double x[ARB] #I reference x values
+double resid[ARB] #I residuals
+double wts[ARB] #I weight array
+double 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] <= double(0.0))
+ next
+ call gctran (gd, real(x[i]), real(resid[i]), x0, y0, 1, 0)
+ 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) {
+ call gscur (gd, real(x[j]), real(resid[j]))
+ call gmark (gd, real(x[j]), real(resid[j]), GM_CROSS, 2., 2.)
+ wts[j] = double(0.0)
+ }
+
+ } else {
+
+ # Find the nearest deleted point.
+ do i = 1, npts {
+ if (wts[i] > double(0.0))
+ next
+ call gctran (gd, real(x[i]), real(resid[i]), x0, y0, 1, 0)
+ 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) {
+ 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.)
+ 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_1graphd (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
+double xref[ARB] #I x reference values
+double yref[ARB] #I y reference values
+double xin[ARB] #I x values
+double yin[ARB] #I y values
+double wts[ARB] #I array of weights
+int npts #I number of points
+
+int i, j
+pointer sp, rxin, ryin
+
+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)
+ 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)
+ call gt_swind (gd, gt)
+ call gtlabax (gd, gt)
+
+ # Mark the data and deleted points.
+ do i = 1, npts {
+ if (wts[i] == double(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.)
+ }
+
+ call gflush (gd)
+ }
+
+ # Mark the rejected points.
+ do i = 1, GM_NREJECT(fit) {
+ j = Memi[GM_REJ(fit)+i-1]
+ call gmark (gd, real(xin[j]), real(yin[j]), GM_CIRCLE, 2., 2.)
+ }
+
+ 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_2graphd (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
+double x[ARB] #I x reference values
+double resid[ARB] #I residual
+double wts[ARB] #I array of weights
+int npts #I number of points
+
+int i, j
+pointer sp, zero
+pointer rxin, ryin
+
+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.
+ 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)
+ call gt_swind (gd, gt)
+ call gtlabax (gd, gt)
+
+ call gpline (gd, Memr[rxin], Memr[zero], npts)
+ }
+
+ # Graph residuals and mark deleted points.
+ if (GG_OVERPLOT(gfit) == NO || GG_NEWFUNCTION(gfit) == YES) {
+ do i = 1, npts {
+ if (wts[i] == double(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.)
+ }
+ }
+
+ # plot rejected points
+ if (GM_NREJECT(fit) > 0) {
+ do i = 1, GM_NREJECT(fit) {
+ j = Memi[GM_REJ(fit)+i-1]
+ call gmark (gd, Memr[rxin+j-1], Memr[ryin+j-1], GM_CIRCLE,
+ 2., 2.)
+ }
+ }
+
+ # 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_conxyd (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
+pointer xbuf, ybuf
+double xint, yint, dx, dy
+
+begin
+ # allocate temporary space
+ call smark (sp)
+ call salloc (xtemp, NGRAPH, TY_DOUBLE)
+ call salloc (ytemp, NGRAPH, TY_DOUBLE)
+ call salloc (xfit1, NGRAPH, TY_DOUBLE)
+ call salloc (yfit1, NGRAPH, TY_DOUBLE)
+ call salloc (xfit2, NGRAPH, TY_DOUBLE)
+ call salloc (yfit2, NGRAPH, TY_DOUBLE)
+ call salloc (xbuf, NGRAPH, TY_REAL)
+ call salloc (ybuf, NGRAPH, TY_REAL)
+
+ # 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.
+ Memd[ytemp] = GM_YMIN(fit)
+ do i = 2, NGRAPH
+ Memd[ytemp+i-1] = Memd[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 amovkd (xint, Memd[xtemp], NGRAPH)
+
+ # X fit.
+ call dgsvector (sx1, Memd[xtemp], Memd[ytemp], Memd[xfit1],
+ NGRAPH)
+ if (sx2 != NULL) {
+ call dgsvector (sx2, Memd[xtemp], Memd[ytemp], Memd[xfit2],
+ NGRAPH)
+ call aaddd (Memd[xfit1], Memd[xfit2], Memd[xfit1], NGRAPH)
+ }
+
+ # Y fit.
+ call dgsvector (sy1, Memd[xtemp], Memd[ytemp], Memd[yfit1],
+ NGRAPH)
+ if (sy2 != NULL) {
+ call dgsvector (sy2, Memd[xtemp], Memd[ytemp], Memd[yfit2],
+ NGRAPH)
+ call aaddd (Memd[yfit1], Memd[yfit2], Memd[yfit1], NGRAPH)
+ }
+
+ # Plot line of constant x.
+ call achtdr (Memd[xfit1], Memr[xbuf], NGRAPH)
+ call achtdr (Memd[yfit1], Memr[ybuf], NGRAPH)
+ call gpline (gd, Memr[xbuf], Memr[ybuf], NGRAPH)
+
+ # 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.
+ Memd[xtemp] = GM_XMIN(fit)
+ do i = 2, NGRAPH
+ Memd[xtemp+i-1] = Memd[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 amovkd (yint, Memd[ytemp], NGRAPH)
+
+ # X fit.
+ call dgsvector (sx1, Memd[xtemp], Memd[ytemp], Memd[xfit1],
+ NGRAPH)
+ if (sx2 != NULL) {
+ call dgsvector (sx2, Memd[xtemp], Memd[ytemp], Memd[xfit2],
+ NGRAPH)
+ call aaddd (Memd[xfit1], Memd[xfit2], Memd[xfit1], NGRAPH)
+ }
+
+
+ # Y fit.
+ call dgsvector (sy1, Memd[xtemp], Memd[ytemp], Memd[yfit1],
+ NGRAPH)
+ if (sy2 != NULL) {
+ call dgsvector (sy2, Memd[xtemp], Memd[ytemp], Memd[yfit2],
+ NGRAPH)
+ call aaddd (Memd[yfit1], Memd[yfit2], Memd[yfit1], NGRAPH)
+ }
+
+ # Plot line of constant y.
+ call achtdr (Memd[xfit1], Memr[xbuf], NGRAPH)
+ call achtdr (Memd[yfit1], Memr[ybuf], NGRAPH)
+ call gpline (gd, Memr[xbuf], Memr[ybuf], NGRAPH)
+
+ # 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_lxyd (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
+double xref[ARB] #I x reference values
+double yref[ARB] #I y reference values
+double xin[ARB] #I x input values
+double 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
+pointer xbuf, ybuf
+real x0, y0, r2, r2min
+double delta, deltax, deltay
+double dgseval()
+
+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 {
+ call gctran (gd, real(xin[i]), real(yin[i]), x0, y0, 1, 0)
+ 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_DOUBLE)
+ call salloc (ytemp, NGRAPH, TY_DOUBLE)
+ call salloc (xfit1, NGRAPH, TY_DOUBLE)
+ call salloc (yfit1, NGRAPH, TY_DOUBLE)
+ call salloc (xfit2, NGRAPH, TY_DOUBLE)
+ call salloc (yfit2, NGRAPH, TY_DOUBLE)
+ call salloc (xbuf, NGRAPH, TY_REAL)
+ call salloc (ybuf, NGRAPH, TY_REAL)
+
+ # Compute the deltas.
+ 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])
+
+ # Set up line of constant x.
+ call amovkd (xref[j], Memd[xtemp], NGRAPH)
+ delta = (GM_YMAX(fit) - GM_YMIN(fit)) / (NGRAPH - 1)
+ Memd[ytemp] = GM_YMIN(fit)
+ do i = 2, NGRAPH
+ Memd[ytemp+i-1] = Memd[ytemp+i-2] + delta
+
+ # X solution.
+ call dgsvector (sx1, Memd[xtemp], Memd[ytemp], Memd[xfit1],
+ NGRAPH)
+ if (sx2 != NULL) {
+ call dgsvector (sx2, Memd[xtemp], Memd[ytemp], Memd[xfit2],
+ NGRAPH)
+ call aaddd (Memd[xfit1], Memd[xfit2], Memd[xfit1], NGRAPH)
+ }
+ call aaddkd (Memd[xfit1], deltax, Memd[xfit1], NGRAPH)
+
+ # Y solution.
+ call dgsvector (sy1, Memd[xtemp], Memd[ytemp], Memd[yfit1],
+ NGRAPH)
+ if (sy2 != NULL) {
+ call dgsvector (sy2, Memd[xtemp], Memd[ytemp], Memd[yfit2],
+ NGRAPH)
+ call aaddd (Memd[yfit1], Memd[yfit2], Memd[yfit1], NGRAPH)
+ }
+ call aaddkd (Memd[yfit1], deltay, Memd[yfit1], NGRAPH)
+
+ # Plot line of constant x.
+ call achtdr (Memd[xfit1], Memr[xbuf], NGRAPH)
+ call achtdr (Memd[yfit1], Memr[ybuf], NGRAPH)
+ call gpline (gd, Memr[xbuf], Memr[ybuf], NGRAPH)
+ call gflush (gd)
+
+ # Set up line of constant y.
+ call amovkd (yref[j], Memd[ytemp], NGRAPH)
+ delta = (GM_XMAX(fit) - GM_XMIN(fit)) / (NGRAPH - 1)
+ Memd[xtemp] = GM_XMIN(fit)
+ do i = 2, NGRAPH
+ Memd[xtemp+i-1] = Memd[xtemp+i-2] + delta
+
+ # X fit.
+ call dgsvector (sx1, Memd[xtemp], Memd[ytemp], Memd[xfit1],
+ NGRAPH)
+ if (sx2 != NULL) {
+ call dgsvector (sx2, Memd[xtemp], Memd[ytemp], Memd[xfit2],
+ NGRAPH)
+ call aaddd (Memd[xfit1], Memd[xfit2], Memd[xfit1], NGRAPH)
+ }
+ call aaddkd (Memd[xfit1], deltax, Memd[xfit1], NGRAPH)
+
+ # Y fit.
+ call dgsvector (sy1, Memd[xtemp], Memd[ytemp], Memd[yfit1],
+ NGRAPH)
+ if (sy2 != NULL) {
+ call dgsvector (sy2, Memd[xtemp], Memd[ytemp], Memd[yfit2],
+ NGRAPH)
+ call aaddd (Memd[yfit1], Memd[yfit2], Memd[yfit1], NGRAPH)
+ }
+ call aaddkd (Memd[yfit1], deltay, Memd[yfit1], NGRAPH)
+
+ # Plot line of constant y.
+ call achtdr (Memd[xfit1], Memr[xbuf], NGRAPH)
+ call achtdr (Memd[yfit1], Memr[ybuf], NGRAPH)
+ call gpline (gd, Memr[xbuf], Memr[ybuf], NGRAPH)
+ 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_gcoeffd (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
+double xshift #O output x shift
+double yshift #O output y shift
+double a #O output x coefficient of x fit
+double b #O output y coefficient of x fit
+double c #O output x coefficient of y fit
+double d #O output y coefficient of y fit
+
+int nxxcoeff, nxycoeff, nyxcoeff, nyycoeff
+pointer sp, xcoeff, ycoeff
+double xxrange, xyrange, xxmaxmin, xymaxmin
+double yxrange, yyrange, yxmaxmin, yymaxmin
+
+int dgsgeti()
+double dgsgetd()
+
+begin
+ # Allocate working space.
+ call smark (sp)
+ call salloc (xcoeff, dgsgeti (sx, GSNCOEFF), TY_DOUBLE)
+ call salloc (ycoeff, dgsgeti (sy, GSNCOEFF), TY_DOUBLE)
+
+ # Get coefficients and numbers of coefficients.
+ call dgscoeff (sx, Memd[xcoeff], nxxcoeff)
+ call dgscoeff (sy, Memd[ycoeff], nyycoeff)
+ nxxcoeff = dgsgeti (sx, GSNXCOEFF)
+ nxycoeff = dgsgeti (sx, GSNYCOEFF)
+ nyxcoeff = dgsgeti (sy, GSNXCOEFF)
+ nyycoeff = dgsgeti (sy, GSNYCOEFF)
+
+ # Get the data range.
+ 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
+ } else {
+ xxrange = double(1.0)
+ xxmaxmin = double(0.0)
+ xyrange = double(1.0)
+ xymaxmin = double(0.0)
+ }
+
+ 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
+ } else {
+ yxrange = double(1.0)
+ yxmaxmin = double(0.0)
+ yyrange = double(1.0)
+ yymaxmin = double(0.0)
+ }
+
+ # Get the shifts.
+ xshift = Memd[xcoeff] + Memd[xcoeff+1] * xxmaxmin / xxrange +
+ Memd[xcoeff+2] * xymaxmin / xyrange
+ yshift = Memd[ycoeff] + Memd[ycoeff+1] * yxmaxmin / yxrange +
+ Memd[ycoeff+2] * yymaxmin / yyrange
+
+ # Get the rotation and scaling parameters and correct for normalization.
+ if (nxxcoeff > 1)
+ a = Memd[xcoeff+1] / xxrange
+ else
+ a = double(0.0)
+ if (nxycoeff > 1)
+ b = Memd[xcoeff+nxxcoeff] / xyrange
+ else
+ b = double(0.0)
+ if (nyxcoeff > 1)
+ c = Memd[ycoeff+1] / yxrange
+ else
+ c = double(0.0)
+ if (nyycoeff > 1)
+ d = Memd[ycoeff+nyxcoeff] / yyrange
+ else
+ d = double(0.0)
+
+ call sfree (sp)
+end
+
+
generated by cgit (git 2.34.1) at 2025-09-20 05:25:32 -0400