From 40e5a5811c6ffce9b0974e93cdd927cbcf60c157 Mon Sep 17 00:00:00 2001 From: Joe Hunkeler Date: Tue, 11 Aug 2015 16:51:37 -0400 Subject: Repatch (from linux) of OSX IRAF --- noao/twodspec/longslit/transform/igsfit/Revisions | 42 +++ noao/twodspec/longslit/transform/igsfit/igscolon.x | 115 +++++++ .../twodspec/longslit/transform/igsfit/igsdelete.x | 103 ++++++ noao/twodspec/longslit/transform/igsfit/igsfit.com | 10 + noao/twodspec/longslit/transform/igsfit/igsfit.x | 373 +++++++++++++++++++++ noao/twodspec/longslit/transform/igsfit/igsget.x | 62 ++++ noao/twodspec/longslit/transform/igsfit/igsgraph.x | 73 ++++ noao/twodspec/longslit/transform/igsfit/igsinit.x | 21 ++ .../longslit/transform/igsfit/igsnearest.x | 51 +++ .../twodspec/longslit/transform/igsfit/igsparams.x | 23 ++ noao/twodspec/longslit/transform/igsfit/igsset.x | 59 ++++ noao/twodspec/longslit/transform/igsfit/igssolve.x | 173 ++++++++++ .../longslit/transform/igsfit/igsundelete.x | 107 ++++++ noao/twodspec/longslit/transform/igsfit/mkpkg | 21 ++ noao/twodspec/longslit/transform/igsfit/xgs.x | 243 ++++++++++++++ 15 files changed, 1476 insertions(+) create mode 100644 noao/twodspec/longslit/transform/igsfit/Revisions create mode 100644 noao/twodspec/longslit/transform/igsfit/igscolon.x create mode 100644 noao/twodspec/longslit/transform/igsfit/igsdelete.x create mode 100644 noao/twodspec/longslit/transform/igsfit/igsfit.com create mode 100644 noao/twodspec/longslit/transform/igsfit/igsfit.x create mode 100644 noao/twodspec/longslit/transform/igsfit/igsget.x create mode 100644 noao/twodspec/longslit/transform/igsfit/igsgraph.x create mode 100644 noao/twodspec/longslit/transform/igsfit/igsinit.x create mode 100644 noao/twodspec/longslit/transform/igsfit/igsnearest.x create mode 100644 noao/twodspec/longslit/transform/igsfit/igsparams.x create mode 100644 noao/twodspec/longslit/transform/igsfit/igsset.x create mode 100644 noao/twodspec/longslit/transform/igsfit/igssolve.x create mode 100644 noao/twodspec/longslit/transform/igsfit/igsundelete.x create mode 100644 noao/twodspec/longslit/transform/igsfit/mkpkg create mode 100644 noao/twodspec/longslit/transform/igsfit/xgs.x (limited to 'noao/twodspec/longslit/transform/igsfit') 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 + +# 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 +include +include +include + +# 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 +include +include + +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 + +# 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 +include +include +include + +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 + +# 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 +include +include + +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 + +# 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 + +# 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 + + +# 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 +include +include +include + +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 + igsdelete.x + igsfit.x igsfit.com + igsget.x igsfit.com + igsgraph.x + igsinit.x + igsnearest.x + igsparams.x igsfit.com + igsset.x igsfit.com + igssolve.x igsfit.com + igsundelete.x + xgs.x + ; 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 + +# 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 -- cgit