From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- noao/onedspec/t_specplot.x | 2030 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 2030 insertions(+) create mode 100644 noao/onedspec/t_specplot.x (limited to 'noao/onedspec/t_specplot.x') diff --git a/noao/onedspec/t_specplot.x b/noao/onedspec/t_specplot.x new file mode 100644 index 00000000..414b275c --- /dev/null +++ b/noao/onedspec/t_specplot.x @@ -0,0 +1,2030 @@ +include +include +include +include +include +include +include +include +include "specplot.h" + +# Define the help information. +define HELP "noao$onedspec/specplot.key" +define PROMPT "specplot options" + + +# T_SPECPLOT -- Plot multiple spectra in a variety of formats and layouts. +# The spectra may be individually scaled and offset in intensity, shifted +# and scaled in wavelength, and plotted in uniform steps. The plotting +# type may be symbols or lines. The spectra may be labeled. Each spectrum +# is read into memory in a structre defined in "specplot.h". An array +# of structures is then manipulated. Each line of two dimensional images +# are treated as separate spectra. + +procedure t_specplot () + +pointer list # List of input spectra +real step # Initial separation step +int labels # Labeling mode +real fraction # Fraction of minimum step +bool yscale # Draw y scale? + +bool wscale +int i, j, n, fd, nspec, wcs, key, redraw +real wx, wy, wx1, wy1, wx2, wy2 +pointer stack, units, cmd, sp, sh, spsave, sps, gp, gt + +bool clgetb() +int clgwrd(), clgcur() +int open(), imtgetim(), getline(), scan(), nscan() +int stridxs(), nowhite(), btoi(), gt_geti() +real clgetr() +pointer sp_nearest(), imtopenp(), gopen(), gt_init() +errchk sp_gdata, un_changer + +define nospec_ 99 + +begin + call smark (stack) + call salloc (units, SZ_LINE, TY_CHAR) + call salloc (cmd, SZ_LINE, TY_CHAR) + call calloc (sps, 100, TY_POINTER) + spsave = NULL + + # Read the input spectrum list into an array of structures. + i = 0 + nspec = 0 + list = imtopenp ("spectra") + call clgstr ("units", Memc[units], SZ_LINE) + if (nowhite (Memc[units], Memc[cmd], SZ_LINE) == 0) + call strcpy ("display", Memc[units], SZ_LINE) + while (imtgetim (list, Memc[cmd], SZ_FNAME) != EOF) { + iferr (call sp_gdata (Memc[cmd], Memc[units], i, sps, nspec)) + call erract (EA_WARN) + } + call imtclose (list) + + + # Set the layout of the spectra. + step = clgetr ("step") + fraction = clgetr ("fraction") + if (clgetb ("autolayout")) { + if (clgetb ("autoscale")) + call sp_autolayout (Memi[sps], nspec, true, fraction, step) + else + call sp_autolayout (Memi[sps], nspec, false, fraction, step) + } + call sp_scale (Memi[sps], nspec, step) + + # Get optional user labels from a file and set the label type. + call clgstr ("ulabels", Memc[cmd], SZ_FNAME) + ifnoerr (fd = open (Memc[cmd], READ_ONLY, TEXT_FILE)) { + do i = 1, nspec { + sp = Memi[sps+i-1] + if (getline (fd, Memc[cmd]) != EOF) + call strcpy (Memc[cmd], SP_ULABEL(sp), SP_SZULABEL) + else + SP_ULABEL(sp) = EOS + j = stridxs ("\n", SP_ULABEL(sp)) + if (j > 0) + call strcpy (SP_ULABEL(sp), SP_ULABEL(sp), j-1) + } + call close (fd) + } + labels = clgwrd ("labels", Memc[cmd], SZ_FNAME, LABELS) + call sp_labels (Memi[sps], nspec, labels) + + # Initialize the graphics + call clgstr ("graphics", Memc[cmd], SZ_FNAME) + gp = gopen (Memc[cmd], NEW_FILE, STDGRAPH) + + gt = gt_init () + call gt_seti (gt, GTSYSID, btoi (clgetb ("sysid"))) + call clgstr ("title", Memc[cmd], SZ_LINE) + call gt_sets (gt, GTTITLE, Memc[cmd]) + call clgstr ("xlabel", Memc[cmd], SZ_LINE) + if (Memc[cmd] != EOS) { + call gt_sets (gt, GTXLABEL, Memc[cmd]) + call gt_sets (gt, GTXUNITS, "") + } else if (nspec > 0) { + if (UN_LABEL(UN(SP_SH(Memi[sps]))) != EOS) { + call gt_sets (gt, GTXLABEL, UN_LABEL(UN(SP_SH(Memi[sps])))) + call gt_sets (gt, GTXUNITS, UN_UNITS(UN(SP_SH(Memi[sps])))) + } else { + call gt_sets (gt, GTXLABEL, LABEL(SP_SH(Memi[sps]))) + call gt_sets (gt, GTXUNITS, UNITS(SP_SH(Memi[sps]))) + } + } + call clgstr ("ylabel", Memc[cmd], SZ_LINE) + call gt_sets (gt, GTYLABEL, Memc[cmd]) + wx = clgetr ("xmin") + call gt_setr (gt, GTXMIN, wx) + wx = clgetr ("xmax") + call gt_setr (gt, GTXMAX, wx) + wx = clgetr ("ymin") + call gt_setr (gt, GTYMIN, wx) + wx = clgetr ("ymax") + call gt_setr (gt, GTYMAX, wx) + wscale = true + yscale = clgetb ("yscale") + #if (!scale) + # call gseti (gp, G_YDRAWTICKS, NO) + + # Draw the graph on the first pass and then read the cursor. + key = 'r' + repeat { + switch (key) { + case '?': # Page help summary + call gpagefile (gp, HELP, PROMPT) + case ':': # Process colon commands + if (Memc[cmd] == '/') + call gt_colon (Memc[cmd], gp, gt, redraw) + else { + i = sp_nearest (gp, wx, wy, key, Memc[cmd], Memi[sps], + nspec) + call sp_colon (Memc[cmd], gp, gt, Memi[sps], nspec, + Memc[units], labels, i, step, fraction, redraw) + if (nspec == 0) { + redraw = NO + goto nospec_ + } + } + case 'a', 'i': # Append or insert a new spectrum + i = sp_nearest (gp, wx, wy, key, Memc[cmd], Memi[sps], nspec) + if (key == 'i') + i = max (0, i - 1) + call printf ("Spectrum: ") + call flush (STDOUT) + if (scan() != EOF) { + call gargwrd (Memc[cmd], SZ_LINE) + if (nscan() == 1) { + iferr { + call sp_gdata (Memc[cmd], Memc[units], + i, sps, nspec) + call sp_labels (Memi[sps], nspec, labels) + call sp_scale (Memi[sps], nspec, step) + redraw = YES + } then + call erract (EA_WARN) + } + } + case 'd': # Delete a spectrum + if (nspec == 0) + goto nospec_ + + i = sp_nearest (gp, wx, wy, key, Memc[cmd], Memi[sps], nspec) + sp = Memi[sps+i-1] + call sp_ptype (SP_PTYPE(sp), SP_COLOR(sp), YES, gp, gt) + call gt_plot (gp, gt, SP_X(sp), SP_Y(sp), SP_NPTS(sp)) + call gline (gp, SP_X(sp), SP_Y(sp), SP_X(sp), SP_Y(SP)) + call sp_delete (i, sps, nspec) + call sp_labels (Memi[sps], nspec, labels) + call sp_scale (Memi[sps], nspec, step) + if (spsave != NULL) + call sp_free (spsave) + spsave = sp +# redraw = YES + case 'e': # Undelete a spectrum + if (spsave != NULL) { + i = sp_nearest (gp, wx, wy, key, Memc[cmd], Memi[sps], + nspec) + i = max (0, i - 1) + call sp_add (spsave, i, sps, nspec) + call sp_labels (Memi[sps], nspec, labels) + call sp_scale (Memi[sps], nspec, step) + spsave = NULL + redraw = YES + } + case 'f': # Toggle wavelength scale + if (wscale) { + call gt_sets (gt, GTXLABEL, "Pixels") + call gt_sets (gt, GTXUNITS, "") + wscale = false + } else { + if (nspec > 0) { + sp = Memi[sps] + sh = SP_SH(sp) + call gt_sets (gt, GTXLABEL, UN_LABEL(UN(sh))) + call gt_sets (gt, GTXUNITS, UN_UNITS(UN(sh))) + } + wscale = true + } + redraw = YES + case 'l', 'p': # Mark label position and enter label. + if (nspec == 0) + goto nospec_ + + i = sp_nearest (gp, wx, wy, key, Memc[cmd], Memi[sps], nspec) + sp = Memi[sps+i-1] + call printf ( + "Spectrum %d: Mark position for label ('q' to cancel)") + call pargi (SP_INDEX(sp)) + i = clgcur ("cursor", wx, wy, wcs, j, Memc[cmd], SZ_LINE) + if (j != 'q') { + call ggwind (gp, wx1, wx2, wy1, wy2) + wx2 = wx2 - wx1 + wy2 = wy2 - wy1 + SP_XLPOS(sp) = (wx - wx1) / wx2 + SP_YLPOS(sp) = (wy - SP_MEAN(sp)) / wy2 + + if (key == 'l') { + call printf ("Spectrum %d: Label = ") + call pargi (SP_INDEX(sp)) + call flush (STDOUT) + if (scan() != EOF) { + call gargstr (SP_ULABEL(sp), SP_SZULABEL) + j = stridxs ("\n", SP_ULABEL(sp)) + if (j > 0) + call strcpy (SP_ULABEL(sp), SP_ULABEL(sp), j-1) + call strcpy (SP_ULABEL(sp), SP_LABEL(sp), + SP_SZLABEL) + } + } + call gtext (gp, wx, wy, SP_LABEL(sp), "") + } + call printf ("\n") + case 'o': # Reorder the spectra to eliminate gaps. + if (nspec == 0) + goto nospec_ + + do i = 1, nspec { + sp = Memi[sps+i-1] + if (SP_INDEX(sp) != i) { + SP_INDEX (sp) = i + redraw = YES + } + } + if (redraw == YES) { + call sp_labels (Memi[sps], nspec, labels) + call sp_scale (Memi[sps], nspec, step) + } + case 'q', 'I': # Quit or interrupt + break + case 'r': # Redraw the current graph + redraw = YES + case 's': # Shift the spectrum nearest the cursor + if (nspec == 0) + goto nospec_ + + i = sp_nearest (gp, wx, wy, key, Memc[cmd], Memi[sps], nspec) + sp = Memi[sps+i-1] + call printf ( "Shift spectrum %d: (q, r, s, t, x, y, z)") + call pargi (SP_INDEX(sp)) + while (clgcur ("cursor", wx1, wy1, wcs, key, Memc[cmd], + SZ_LINE) != EOF) { + switch (key) { + case 's': + if (wy != SP_OFFSET(sp)) { + call sp_ptype (SP_PTYPE(sp), SP_COLOR(sp), + YES, gp, gt) + call gt_plot (gp, gt, SP_X(sp), SP_Y(sp), + SP_NPTS(sp)) + call gline (gp, SP_X(sp), SP_Y(sp), SP_X(sp), + SP_Y(sp)) + SP_SCALE(sp) = SP_SCALE(sp) * + (wy1 - SP_OFFSET(sp)) / (wy - SP_OFFSET(sp)) + call sp_scale (sp, 1, step) + call sp_ptype (SP_PTYPE(sp), SP_COLOR(sp), + NO, gp, gt) + call gt_plot (gp, gt, SP_X(sp), SP_Y(sp), + SP_NPTS(sp)) + wy = wy1 + } + case 't': + if (wy != SP_OFFSET(sp)) { + call sp_ptype (SP_PTYPE(sp), SP_COLOR(sp), + YES, gp, gt) + call gt_plot (gp, gt, SP_X(sp), SP_Y(sp), + SP_NPTS(sp)) + call gline (gp, SP_X(sp), SP_Y(sp), SP_X(sp), + SP_Y(sp)) + if (UN_CLASS(UN(SP_SH(sp))) == UN_VEL) + SP_XOFFSET(sp) = SP_XOFFSET(sp) + wx1 - wx + else + SP_XSCALE(sp) = SP_XSCALE(sp) * wx1 / wx + SP_SCALE(sp) = SP_SCALE(sp) * + (wy1 - SP_OFFSET(sp)) / (wy - SP_OFFSET(sp)) + call sp_scale (sp, 1, step) + call sp_ptype (SP_PTYPE(sp), SP_COLOR(sp), + NO, gp, gt) + call gt_plot (gp, gt, SP_X(sp), SP_Y(sp), + SP_NPTS(sp)) + wx = wx1 + wy = wy1 + } + case 'x': + call sp_ptype (SP_PTYPE(sp), SP_COLOR(sp), + YES, gp, gt) + call gt_plot (gp, gt, SP_X(sp), SP_Y(sp), SP_NPTS(sp)) + call gline (gp, SP_X(sp), SP_Y(sp), SP_X(sp), SP_Y(SP)) + if (UN_CLASS(UN(SP_SH(sp))) == UN_VEL) + SP_XOFFSET(sp) = SP_XOFFSET(sp) + wx1 - wx + else + SP_XSCALE(sp) = SP_XSCALE(sp) * wx1 / wx + call sp_scale (sp, 1, step) + call sp_ptype (SP_PTYPE(sp), SP_COLOR(sp), + NO, gp, gt) + call gt_plot (gp, gt, SP_X(sp), SP_Y(sp), SP_NPTS(sp)) + wx = wx1 + case 'y': + call sp_ptype (SP_PTYPE(sp), SP_COLOR(sp), + YES, gp, gt) + call gt_plot (gp, gt, SP_X(sp), SP_Y(sp), SP_NPTS(sp)) + call gline (gp, SP_X(sp), SP_Y(sp), SP_X(sp), SP_Y(SP)) + SP_OFFSET(sp) = SP_OFFSET(sp) + wy1 - wy + call sp_scale (sp, 1, step) + call sp_ptype (SP_PTYPE(sp), SP_COLOR(sp), + NO, gp, gt) + call gt_plot (gp, gt, SP_X(sp), SP_Y(sp), SP_NPTS(sp)) + wy = wy1 + case 'z': + call sp_ptype (SP_PTYPE(sp), SP_COLOR(sp), + YES, gp, gt) + call gt_plot (gp, gt, SP_X(sp), SP_Y(sp), SP_NPTS(sp)) + call gline (gp, SP_X(sp), SP_Y(sp), SP_X(sp), SP_Y(SP)) + if (UN_CLASS(UN(SP_SH(sp))) == UN_VEL) + SP_XOFFSET(sp) = SP_XOFFSET(sp) + wx1 - wx + else + SP_XSCALE(sp) = SP_XSCALE(sp) * wx1 / wx + SP_OFFSET(sp) = SP_OFFSET(sp) + wy1 - wy + call sp_scale (sp, 1, step) + call sp_ptype (SP_PTYPE(sp), SP_COLOR(sp), + NO, gp, gt) + call gt_plot (gp, gt, SP_X(sp), SP_Y(sp), SP_NPTS(sp)) + wx = wx1 + wy = wy1 + case 'r': + if (gt_geti (gt, GTSYSID) == YES) { + call sprintf (Memc[cmd], SZ_LINE, + "Separation step = %g") + call pargr (step) + call gt_sets (gt, GTPARAMS, Memc[cmd]) + } else + call gt_sets (gt, GTPARAMS, "") + call sp_plot (gp, gt, Memi[sps], nspec, wscale, yscale) + case 'q': + break + } + call printf ( "Shift spectrum %d: (q, r, s, t, x, y, z)") + call pargi (SP_INDEX(sp)) + } + call printf ("\n") + case 't': # Set a wavelength scale using the cursor. + if (nspec == 0) + goto nospec_ + + i = sp_nearest (gp, wx, wy, key, Memc[cmd], Memi[sps], nspec) + sp = Memi[sps+i-1] + call printf ("X coordinate (%g): ") + call pargr (wx) + call flush (STDOUT) + if (scan() != EOF) { + call gargr (wy) + if (nscan() == 0) + wy = wx + } else + wy = wx + call printf ("Mark another position") + i = clgcur ("cursor", wx1, wy1, wcs, key, Memc[cmd], SZ_LINE) + call printf ("X coordinate (%g): ") + call pargr (wx1) + call flush (STDOUT) + if (scan() != EOF) { + call gargr (wy1) + if (nscan() == 0) + wy1 = wx1 + } else + wy1 = wx1 + if (wx != wx1) { + n = SP_NPTS(sp) - 1 + sh = SP_PX(sp) - 1 + if (SP_WPC(sp) > 0.) { + for (i=1; in && wx>Memr[sh+i]; i=i+1) + ; + for (j=1; j>n && wx1>Memr[sh+j]; j=j+1) + ; + } + wx = i + (wx - Memr[sh+i]) / (Memr[sh+i+1] - Memr[sh+i]) + wx1 = j + (wx1 - Memr[sh+j]) / (Memr[sh+j+1] - Memr[sh+j]) + SP_WPC(sp) = (wy - wy1) / (wx - wx1) + SP_W0(sp) = wy - SP_WPC(sp) * (wx - 1) + call sp_linear (sp) + call sp_scale (sp, 1, step) + redraw = YES + } + case 'u': # Set a wavelength point using the cursor. + if (nspec == 0) + goto nospec_ + + i = sp_nearest (gp, wx, wy, key, Memc[cmd], Memi[sps], nspec) + sp = Memi[sps+i-1] + call printf ("X coordinate (%g): ") + call pargr (wx) + call flush (STDOUT) + if (scan() != EOF) { + call gargr (wx1) + if (nscan() == 1) { + SP_XOFFSET(sp) = SP_XOFFSET(sp) + wx1 - wx + call sp_scale (sp, 1, step) + redraw = YES + } + } + case 'v': # Change to velocity scale + if (nspec == 0) + goto nospec_ + + iferr { + do i = 1, nspec { + sp = Memi[sps+i-1] + sh = SP_SH(sp) + if (i == 1) { + call un_changer (UN(sh), "angstroms", wx, 1, NO) + call sprintf (Memc[units], SZ_LINE, + "km/s %g angstroms") + call pargr (wx) + call un_changer (UN(sh), Memc[units], Memr[SX(sh)], + SN(sh), YES) + call gt_sets (gt, GTXLABEL, UN_LABEL(UN(sh))) + call gt_sets (gt, GTXUNITS, UN_UNITS(UN(sh))) + redraw = YES + } else + call un_changer (UN(sh), Memc[units], Memr[SX(sh)], + SN(sh), YES) + SP_W0(sp) = Memr[SX(sh)] + SP_WPC(sp) = (Memr[SX(sh)+SN(sh)-1] - Memr[SX(sh)]) / + (SN(sh) - 1) + SP_XSCALE(sp) = 1. + SP_XOFFSET(sp) = 0. + call sp_scale (sp, 1, step) + } + } then + call erract (EA_WARN) + case 'w': # Window the graph + call gt_window (gt, gp, "cursor", redraw) + case 'x': # No layout + if (nspec == 0) + goto nospec_ + + do i = 1, nspec { + sp = Memi[sps+i-1] + SP_SCALE(sp) = 1. + SP_OFFSET(sp) = 0. + } + call sp_scale (Memi[sps], nspec, step) + redraw = YES + case 'y': # Layout the spectra offsets to common mean + if (nspec == 0) + goto nospec_ + + call sp_autolayout (Memi[sps], nspec, false, fraction, step) + call sp_scale (Memi[sps], nspec, step) + redraw = YES + case 'z': # Layout the spectra scaled to common mean + if (nspec == 0) + goto nospec_ + + call sp_autolayout (Memi[sps], nspec, true, fraction, step) + call sp_scale (Memi[sps], nspec, step) + redraw = YES + default: + call printf ("\007") + } + + # Redraw the graph as needed. + if (redraw == YES) { + if (gt_geti (gt, GTSYSID) == YES) { + call sprintf (Memc[cmd], SZ_LINE, "Separation step = %g") + call pargr (step) + call gt_sets (gt, GTPARAMS, Memc[cmd]) + } else + call gt_sets (gt, GTPARAMS, "") + call sp_plot (gp, gt, Memi[sps], nspec, wscale, yscale) + redraw = NO + } +nospec_ + if (nspec == 0) + call printf ("No spectra defined\007") + + } until (clgcur ("cursor", wx, wy, wcs, key, Memc[cmd], SZ_LINE) == EOF) + + + call clgstr ("logfile", Memc[cmd], SZ_LINE) + if (nowhite (Memc[cmd], Memc[cmd], SZ_LINE) > 0) + iferr (call sp_vshow (Memc[cmd], NULL, Memi[sps], nspec, step)) + call erract (EA_WARN) + + # Close the graphics device and free memory. + call gclose (gp) + call gt_free (gt) + + if (nspec > 0) { + do i = 1, nspec + call sp_free (Memi[sps+i-1]) + } + if (spsave != NULL) + call sp_free (spsave) + call mfree (sps, TY_POINTER) + call sfree (stack) +end + + +# SP_SCALE -- Scale the spectra. This uses the wavelength scale and intensity +# scale parameters defined for each spectrum and adds the intensity offset. + +procedure sp_scale (sps, nspec, step) + +pointer sps[ARB] # Spectrum structures +int nspec # Number of spectra +real step # Final step + +int i, npts +real scale, offset +pointer sp, sh + +begin + do i = 1, nspec { + sp = sps[i] + sh = SP_SH(sp) + npts = SP_NPTS(sp) + + scale = SP_XSCALE(sp) + offset = SP_XOFFSET(sp) + call altmr (Memr[SX(sh)], SP_X(sp), npts, scale, offset) + + scale = SP_SCALE(sp) + offset = SP_OFFSET(sp) + (SP_INDEX(sp) - 1) * step + call altmr (Memr[SY(sh)], SP_Y(sp), npts, scale, offset) + + SP_MEAN(sp) = SP_OMEAN(sp) * scale + offset + SP_MIN(sp) = SP_OMIN(sp) * scale + offset + SP_MAX(sp) = SP_OMAX(sp) * scale + offset + } +end + + +# SP_AUTOLAYOUT -- Apply an automatic layout algorithm in which the spectra +# are scaled or offset to a common mean and a separation step is computed +# to provide a specified degree of overlap between the nearest spectra. + +procedure sp_autolayout (sps, nspec, autoscale, fraction, step) + +pointer sps[ARB] # Spectrum structures +int nspec # Number of spectra +bool autoscale # Scale spectra to common mean? +real fraction # Fraction to adjust step +real step # Final step + +int i +real a, b, scale, offset +pointer sp + +begin + if (nspec < 2) + return + + # Scale to the lowest indexed spectrum (usually 1). + sp = sps[1] + scale = SP_SCALE(sp) + offset = SP_OFFSET(sp) + a = SP_OMEAN(sp) + + # If desired use a multiplicative scaling to a common mean. + # If the mean is <= 0 then use offset to common mean. + + if (autoscale) { + do i = 2, nspec { + sp = sps[i] + if (a * SP_OMEAN(sp) > 0.) { + SP_SCALE(sp) = a / SP_OMEAN(sp) * scale + SP_OFFSET(sp) = offset + } else { + SP_SCALE(sp) = scale + SP_OFFSET(sp) = (a - SP_OMEAN(sp)) * scale + offset + } + } + + # Otherwise use an offset scaling to a common mean. + } else { + do i = 2, nspec { + sp = sps[i] + SP_SCALE(sp) = scale + SP_OFFSET(sp) = (a - SP_OMEAN(sp)) * scale + offset + } + } + + # Compute the minimum step which just separates the maximum of + # one spectrum from the minimum of the next spectrum. A degree + # of overlap can be set using the fraction parameter. + + step = -MAX_REAL + do i = 2, nspec { + sp = sps[i-1] + a = SP_OMAX(sp) * SP_SCALE(sp) + SP_OFFSET(sp) + sp = sps[i] + b = SP_OMIN(sp) * SP_SCALE(sp) + SP_OFFSET(sp) + step = max (step, a - b) + } + step = fraction * step +end + + +# SP_PLOT -- Determine the range of all the data and then make a plot with +# specified labels. The GTOOLS procedures are used to allow user adjustment. + +procedure sp_plot (gp, gt, sps, nspec, wscale, yscale) + +pointer gp # GIO pointer +pointer gt # GTOOLS pointer +pointer sps[ARB] # Spectrum structures +int nspec # Number of spectra +bool wscale # Draw in world coordinates? +bool yscale # Draw Y scale? + +int i, n +real x, y, xmin, xmax, ymin, ymax +pointer sp, pix + +begin + # Set the default limits from the data. + xmin = MAX_REAL + xmax = -MAX_REAL + ymin = MAX_REAL + ymax = -MAX_REAL + n = 0 + do i = 1, nspec { + sp = sps[i] + if (wscale) { + xmin = min (xmin, SP_X(sp), Memr[SP_PX(sp)+SP_NPTS(sp)-1]) + xmax = max (xmax, SP_X(sp), Memr[SP_PX(sp)+SP_NPTS(sp)-1]) + } else { + n = max (n, SP_NPTS(sp)) + xmin = 1 + xmax = n + } + ymin = min (ymin, SP_MIN(sp)) + ymax = max (ymax, SP_MAX(sp)) + } + + if (xmin > xmax) { + xmin = 0. + xmax = 1. + } + if (ymin > ymax) { + ymin = 0. + ymax = 1. + } + + # Draw the axes with GTOOLS limits override. + #call gframe (gp) + call gclear (gp) + if (!yscale) + call gseti (gp, G_YDRAWTICKS, NO) + call gswind (gp, xmin, xmax, ymin, ymax) + call gt_swind (gp, gt) + call gt_labax (gp, gt) + + # The label positions are based on the limits of the graph. + call ggwind (gp, xmin, xmax, ymin, ymax) + xmax = xmax - xmin + ymax = ymax - ymin + + if (!wscale) { + call malloc (pix, n, TY_REAL) + do i = 1, n + Memr[pix+i-1] = i + } + + do i = 1, nspec { + sp = sps[i] + call sp_ptype (SP_PTYPE(sp), SP_COLOR(sp), NO, gp, gt) + if (wscale) + call gt_plot (gp, gt, SP_X(sp), SP_Y(sp), SP_NPTS(sp)) + else + call gt_plot (gp, gt, Memr[pix], SP_Y(sp), SP_NPTS(sp)) + x = SP_XLPOS(sp) * xmax + xmin + y = SP_YLPOS(sp) * ymax + SP_MEAN(sp) + call gtext (gp, x, y, SP_LABEL(sp), "") + } + + if (!wscale) + call mfree (pix, TY_REAL) +end + + +# SP_PTYPE -- Decode the plotting type and set the GTOOLS structure. + +procedure sp_ptype (ptype, color, erase, gp, gt) + +char ptype[ARB] # Plotting type string +int color # Color +int erase # Erase plot? +pointer gp # GIO pointer +pointer gt # GTOOLS pointer + +int i, j, ctoi() +pointer sp, gttype + +begin + call smark (sp) + call salloc (gttype, SZ_LINE, TY_CHAR) + call gt_gets (gt, GTTYPE, Memc[gttype], SZ_LINE) + + i = 1 + if (ctoi (ptype, i, j) > 0) { + if (j < 0) + call gt_sets (gt, GTTYPE, "histogram") + else + call gt_sets (gt, GTTYPE, "line") + if (erase == YES) + call gt_seti (gt, GTLINE, 0) + else + call gt_seti (gt, GTLINE, abs(j)) + } else { + call gt_sets (gt, GTTYPE, "mark") + call gt_sets (gt, GTMARK, ptype) + if (erase == YES) + call gseti (gp, G_PMLTYPE, 0) + else + call gseti (gp, G_PMLTYPE, 1) + } + call gt_seti (gt, GTCOLOR, color) + + call sfree (sp) +end + + +# List of colon commands. +define CMDS "|show|vshow|step|fraction|move|shift|w0|wpc|velocity|redshift\ + |offset|scale|xlpos|ylpos|label|ulabel|ptype|units|color|" + +define SHOW 1 # Show +define VSHOW 2 # Verbose show +define STEP 3 # Separation step +define FRACTION 4 # Fraction for autolayout +define MOVE 5 # Move spectrum index +define SHIFT 6 # Shift spectrum indices +define WZP 7 # Wavelength zero point +define WPC 8 # Wavelength per channel +define VELOCITY 9 # Radial velocity +define REDSHIFT 10 # Redshift +define OFFSET 11 # Intensity offset +define SCALE 12 # Intensity scale +define XLPOS 13 # X label position +define YLPOS 14 # Y label position +define LABEL 15 # Type of labels +define ULABEL 16 # User label +define PTYPE 17 # Plot type +define UNITS 18 # Plot units +define COLOR 19 # Color + +# SP_COLON -- Interpret colon commands. + +procedure sp_colon (cmdstr, gp, gt, sps, nspec, units, labels, current, step, + fraction, redraw) + +char cmdstr[ARB] # Colon command +pointer gp # GIO pointer (used for paging screen) +pointer gt # GTOOLS pointer +pointer sps[ARB] # Array of spectra structures +int nspec # Number of spectra +char units[SZ_LINE] # Units string +int labels # Label type +int current # Current spectrum element (0 if not defined) +real step # Separation step +real fraction # Fraction for autolayout +int redraw # Redraw graph + +int i, j, index, ncmd +real rval +pointer stack, cmd, sp, sh, un1, un2 + +int nscan(), strdic(), ctoi(), stridxs() +pointer un_open() + +define done_ 10 + +begin + call smark (stack) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Scan the command string and get the first word. + call sscan (cmdstr) + call gargwrd (Memc[cmd], SZ_LINE) + + # Parse the optional spectrum index. Moving the the end of string. + # Set the spectrum element to 0 if a non-numeric index is specified. + # If an index number is given find the appropriate element and print + # an error if the spectrum index is not defined. + i = stridxs ("[", Memc[cmd]) + j = 0 + if (i > 0) { + Memc[cmd+i-1] = EOS + current = 0 + + i = i + 1 + if (ctoi (Memc[cmd], i, index) > 0) { + for (i=1; (i<=nspec)&&(SP_INDEX(sps[i])!=index); i=i+1) + ; + + current = i + if (current > nspec) { + call printf ("Spectrum %d not defined") + call pargi (index) + call sfree (stack) + return + } + } + j = current + } + + # Parse the command. Print the command if unknown. + ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, CMDS) + + switch (ncmd) { + case SHOW: # show spectrum parameters + call gargwrd (Memc[cmd], SZ_LINE) + if (nscan() == 1) + call sp_show ("STDOUT", gp, sps, nspec, step) + else + iferr (call sp_show (Memc[cmd], NULL, sps, nspec, step)) + call erract (EA_WARN) + case VSHOW: # show spectrum parameters + call gargwrd (Memc[cmd], SZ_LINE) + if (nscan() == 1) + call sp_vshow ("STDOUT", gp, sps, nspec, step) + else + iferr (call sp_vshow (Memc[cmd], NULL, sps, nspec, step)) + call erract (EA_WARN) + case STEP: # set or show step + call gargr (rval) + if (nscan() == 1) { + call printf ("step %g") + call pargr (step) + } else { + step = rval + call sp_scale (sps, nspec, step) + redraw = YES + } + case FRACTION: # set or show autolayout fraction + call gargr (rval) + if (nscan() == 1) { + call printf ("fraction %g") + call pargr (fraction) + } else + fraction = rval + case MOVE: # Move spectrum by index + call gargi (index) + if (nscan() > 1) { + if (current > 0) { + sp = sps[current] + if (index != SP_INDEX(sp)) { + SP_INDEX(sp) = index + + for (i=current; iSP_INDEX(sps[i])); i=i+1) + ; + for (j=nspec; j>i; j=j-1) + sps[j] = sps[j-1] + sps[i] = sp + current = i + + for (j=i; j 1) { + if (current > 0) { + if (j > 0) { + for (i=current; i<=nspec; i=i+1) { + sp = sps[i] + SP_INDEX(sp) = SP_INDEX(sp) + j + call sp_scale (sp, 1, step) + call sp_labels (sp, 1, labels) + redraw = YES + } + } else if (j < 0) { + for (i=current; i>0; i=i-1) { + sp = sps[i] + SP_INDEX(sp) = SP_INDEX(sp) + j + call sp_scale (sp, 1, step) + call sp_labels (sp, 1, labels) + redraw = YES + } + } + } else { + for (i=1; i<=nspec; i=i+1) { + sp = sps[i] + SP_INDEX(sp) = SP_INDEX(sp) + j + call sp_scale (sp, 1, step) + call sp_labels (sp, 1, labels) + redraw = YES + } + } + } + case WZP: # set or show zero point wavelength + call gargr (rval) + if (current > 0) { + sp = sps[current] + if (nscan() == 1) { + call printf ("w0[%d] %g") + call pargi (SP_INDEX(sp)) + call pargr (SP_W0(sp)*SP_XSCALE(sp)+SP_XOFFSET(sp)) + } else { + SP_XOFFSET(sp) = rval - SP_W0(sp) * SP_XSCALE(sp) + call sp_scale (sp, 1, step) + redraw = YES + } + } else { + if (nscan() == 1) { + call printf ("w0:") + do i = 1, nspec { + sp = sps[i] + call printf (" %d=%g") + call pargi (SP_INDEX(sp)) + call pargr (SP_W0(sp)*SP_XSCALE(sp)+SP_XOFFSET(sp)) + } + } else { + do i = 1, nspec { + sp = sps[i] + SP_XOFFSET(sp) = rval - SP_W0(sp) * SP_XSCALE(sp) + call sp_scale (sp, 1, step) + redraw = YES + } + } + } + case WPC: # set or show wavelength per channel + call gargr (rval) + if (current > 0) { + sp = sps[current] + if (nscan() == 1) { + call printf ("wpc[%d] %g") + call pargi (SP_INDEX(sp)) + call pargr (SP_WPC(sp)*SP_XSCALE(sp)) + } else { + SP_WPC(sp) = rval + call sp_linear (sp) + call sp_scale (sp, 1, step) + redraw = YES + } + } else { + if (nscan() == 1) { + call printf ("wpc:") + do i = 1, nspec { + sp = sps[i] + call printf (" %d=%g") + call pargi (SP_INDEX(sp)) + call pargr (SP_WPC(sp)*SP_XSCALE(sp)) + } + } else { + do i = 1, nspec { + sp = sps[i] + SP_WPC(sp) = rval + call sp_linear (sp) + call sp_scale (sp, 1, step) + redraw = YES + } + } + } + case VELOCITY: # set or show radial velocity + if (nspec < 0) + goto done_ + call gargr (rval) + un1 = UN(SP_SH(sps[1])) + if (UN_CLASS(un1) == UN_VEL) { + if (current > 0) { + sp = sps[current] + if (nscan() == 1) { + call printf ("velocity[%d] %g") + call pargi (SP_INDEX(sp)) + call pargr (SP_XOFFSET(sp)) + } else { + SP_XOFFSET(sp) = rval + call sp_scale (sp, 1, step) + redraw = YES + } + } else { + if (nscan() == 1) { + call printf ("velocity:") + do i = 1, nspec { + sp = sps[i] + call printf (" %d=%g") + call pargi (SP_INDEX(sp)) + call pargr (SP_XOFFSET(sp)) + } + } else { + do i = 1, nspec { + sp = sps[i] + SP_XOFFSET(sp) = rval + call sp_scale (sp, 1, step) + redraw = YES + } + } + } + } else if (UN_CLASS(un1) != UN_UNKNOWN) { + if (current > 0) { + sp = sps[current] + call sprintf (Memc[cmd], SZ_LINE, "km/s %g %s") + call pargr (SP_W0(sp)) + call pargstr (UN_UNITS(un1)) + if (nscan() == 1) { + if (SP_XSCALE(sp) != 1.) { + rval = SP_W0(sp) * SP_XSCALE(sp) + call un_changer (un1, Memc[cmd], rval, 1, NO) + } else + rval = 0. + call printf ("velocity[%d] %g") + call pargi (SP_INDEX(sp)) + call pargr (rval) + } else { + un2 = un_open (Memc[cmd]) + call un_ctranr (un2, un1, rval, rval, 1) + call un_close (un2) + SP_XSCALE(sp) = rval / SP_W0(sp) + call sp_scale (sp, 1, step) + redraw = YES + } + } else { + if (nscan() == 1) { + call printf ("velocity:") + do i = 1, nspec { + sp = sps[i] + if (SP_XSCALE(sp) != 1.) { + call sprintf (Memc[cmd], SZ_LINE, "km/s %g %s") + call pargr (SP_W0(sp)) + call pargstr (UN_UNITS(un1)) + rval = SP_W0(sp) * SP_XSCALE(sp) + call un_changer (un1, Memc[cmd], rval, 1, NO) + } else + rval = 0. + call printf (" %d=%g") + call pargi (SP_INDEX(sp)) + call pargr (rval) + } + } else { + do i = 1, nspec { + sp = sps[i] + call sprintf (Memc[cmd], SZ_LINE, "km/s %g %s") + call pargr (SP_W0(sp)) + call pargstr (UN_UNITS(un1)) + un2 = un_open (Memc[cmd]) + call un_ctranr (un2, un1, rval, rval, 1) + call un_close (un1) + SP_XSCALE(sp) = rval / SP_W0(sp) + call sp_scale (sps[i], 1, step) + redraw = YES + } + } + } + } + case REDSHIFT: # set or show redshift + if (nspec < 0) + goto done_ + call gargr (rval) + un1 = UN(SP_SH(sps[1])) + if (UN_CLASS(un1) == UN_VEL) { + if (current > 0) { + sp = sps[current] + if (nscan() == 1) { + call printf ("redshift[%d] %g") + call pargi (SP_INDEX(sp)) + call pargr (SP_XOFFSET(sp)/UN_SCALE(un1)) + } else { + SP_XOFFSET(sp) = rval * UN_SCALE(un1) + call sp_scale (sp, 1, step) + redraw = YES + } + } else { + if (nscan() == 1) { + call printf ("redshift:") + do i = 1, nspec { + sp = sps[i] + call printf (" %d=%g") + call pargi (SP_INDEX(sp)) + call pargr (SP_XOFFSET(sp)/UN_SCALE(un1)) + } + } else { + do i = 1, nspec { + SP_XOFFSET(sp) = rval * UN_SCALE(un1) + call sp_scale (sps[i], 1, step) + redraw = YES + } + } + } + } else if (UN_CLASS(un1) == UN_WAVE) { + if (current > 0) { + sp = sps[current] + if (nscan() == 1) { + call printf ("redshift[%d] %g") + call pargi (SP_INDEX(sp)) + call pargr (SP_XSCALE(sp)-1) + } else { + rval = 1. + rval + SP_XSCALE(sp) = rval + call sp_scale (sp, 1, step) + redraw = YES + } + } else { + if (nscan() == 1) { + call printf ("redshift:") + do i = 1, nspec { + sp = sps[i] + call printf (" %d=%g") + call pargi (SP_INDEX(sp)) + call pargr (SP_XSCALE(sp)-1) + } + } else { + rval = 1. + rval + do i = 1, nspec { + SP_XSCALE(sps[i]) = rval + call sp_scale (sps[i], 1, step) + redraw = YES + } + } + } + } else if (UN_CLASS(un1) == UN_FREQ || UN_CLASS(un1) == UN_ENERGY) { + if (current > 0) { + sp = sps[current] + if (nscan() == 1) { + call printf ("redshift[%d] %g") + call pargi (SP_INDEX(sp)) + call pargr (1./SP_XSCALE(sp)-1) + } else { + rval = 1. / (1. + rval) + SP_XSCALE(sp) = rval + call sp_scale (sp, 1, step) + redraw = YES + } + } else { + if (nscan() == 1) { + call printf ("redshift:") + do i = 1, nspec { + sp = sps[i] + call printf (" %d=%g") + call pargi (SP_INDEX(sp)) + call pargr (1./SP_XSCALE(sp)-1) + } + } else { + rval = 1./ (1. + rval) + do i = 1, nspec { + SP_XSCALE(sps[i]) = rval + call sp_scale (sps[i], 1, step) + redraw = YES + } + } + } + } + case OFFSET: # set or show intensity offset + call gargr (rval) + if (current > 0) { + sp = sps[current] + if (nscan() == 1) { + call printf ("offset[%d] %g") + call pargi (SP_INDEX(sp)) + call pargr (SP_OFFSET(sp)) + } else { + SP_OFFSET(sp) = rval + call sp_scale (sp, 1, step) + redraw = YES + } + } else { + if (nscan() == 1) { + call printf ("offset:") + do i = 1, nspec { + sp = sps[i] + call printf (" %d=%g") + call pargi (SP_INDEX(sp)) + call pargr (SP_OFFSET(sp)) + } + } else { + do i = 1, nspec { + SP_OFFSET(sps[i]) = rval + call sp_scale (sps[i], 1, step) + redraw = YES + } + } + } + case SCALE: # set or show intensity scale + call gargr (rval) + if (current > 0) { + sp = sps[current] + if (nscan() == 1) { + call printf ("scale[%d] %g") + call pargi (SP_INDEX(sp)) + call pargr (SP_SCALE(sp)) + } else { + SP_SCALE(sp) = rval + call sp_scale (sp, 1, step) + redraw = YES + } + } else { + if (nscan() == 1) { + call printf ("scale:") + do i = 1, nspec { + sp = sps[i] + call printf (" %d=%g") + call pargi (SP_INDEX(sp)) + call pargr (SP_SCALE(sp)) + } + } else { + do i = 1, nspec { + SP_SCALE(sps[i]) = rval + call sp_scale (sps[i], 1, step) + redraw = YES + } + } + } + case XLPOS: # set or show X label position + call gargr (rval) + if (current > 0) { + sp = sps[current] + if (nscan() == 1) { + call printf ("xlpos[%d] %g") + call pargi (SP_INDEX(sp)) + call pargr (SP_XLPOS(sp)) + } else { + SP_XLPOS(sp) = rval + redraw = YES + } + } else { + if (nscan() == 1) { + call printf ("xlpos:") + do i = 1, nspec { + sp = sps[i] + call printf (" %d=%g") + call pargi (SP_INDEX(sp)) + call pargr (SP_XLPOS(sp)) + } + } else { + do i = 1, nspec { + SP_XLPOS(sps[i]) = rval + redraw = YES + } + } + } + case YLPOS: # set or show Y label position + call gargr (rval) + if (current > 0) { + sp = sps[current] + if (nscan() == 1) { + call printf ("ylpos[%d] %g") + call pargi (SP_INDEX(sp)) + call pargr (SP_YLPOS(sp)) + } else { + SP_YLPOS(sp) = rval + redraw = YES + } + } else { + if (nscan() == 1) { + call printf ("ylpos:") + do i = 1, nspec { + sp = sps[i] + call printf (" %d=%g") + call pargi (SP_INDEX(sp)) + call pargr (SP_YLPOS(sp)) + } + } else { + do i = 1, nspec { + SP_YLPOS(sps[i]) = rval + redraw = YES + } + } + } + case LABEL: # Set or show label type + call gargwrd (Memc[cmd], SZ_LINE) + if (nscan() == 1) { + switch (labels) { + case LABEL_NONE: + call printf ("labels none") + case LABEL_IMNAME: + call printf ("labels imname") + case LABEL_IMTITLE: + call printf ("labels imtitle") + case LABEL_INDEX: + call printf ("labels index") + case LABEL_USER: + call printf ("labels user") + } + } else { + ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, LABELS) + if (ncmd == 0) { + call printf ("Unknown label type: %s") + call pargstr (Memc[cmd]) + } else { + labels = ncmd + call sp_labels (sps, nspec, labels) + } + } + case ULABEL: # Set or show user labels + call gargwrd (Memc[cmd], SZ_LINE) + if (current > 0) { + sp = sps[current] + if (nscan() == 1) { + call printf ("ulabel[%d] %s") + call pargi (SP_INDEX(sp)) + call pargstr (SP_ULABEL(sp)) + } else { + call reset_scan () + call gargwrd (Memc[cmd], SZ_LINE) + call gargstr (Memc[cmd], SZ_LINE) + call strcpy (Memc[cmd], SP_ULABEL(sp), SP_SZULABEL) + if (labels == LABEL_USER) + call strcpy (SP_ULABEL(sp), SP_LABEL(sp), SP_SZLABEL) + } + } else { + if (nscan() == 1) { + call printf ("ulabel:") + do i = 1, nspec { + sp = sps[i] + call printf (" %d=%s") + call pargi (SP_INDEX(sp)) + call pargstr (SP_ULABEL(sp)) + } + } else { + call reset_scan () + call gargwrd (Memc[cmd], SZ_LINE) + call gargstr (Memc[cmd], SZ_LINE) + do i = 1, nspec { + sp = sps[i] + call strcpy (Memc[cmd], SP_ULABEL(sp), SP_SZULABEL) + if (labels == LABEL_USER) + call strcpy (SP_ULABEL(sp), SP_LABEL(sp),SP_SZLABEL) + } + } + } + case PTYPE: # Set or show plotting type + call gargwrd (Memc[cmd], SZ_LINE) + if (current > 0) { + sp = sps[current] + if (nscan() == 1) { + call printf ("ptype[%d] %s") + call pargi (SP_INDEX(sp)) + call pargstr (SP_PTYPE(sp)) + } else { + call strcpy (Memc[cmd], SP_PTYPE(sp), SP_SZPTYPE) + } + } else { + if (nscan() == 1) { + call printf ("ptype:") + do i = 1, nspec { + sp = sps[i] + call printf (" %d=%s") + call pargi (SP_INDEX(sp)) + call pargstr (SP_PTYPE(sp)) + } + } else { + do i = 1, nspec + call strcpy (Memc[cmd], SP_PTYPE(sps[i]), + SP_SZPTYPE) + } + } + case UNITS: # Change plotting units + # Any change of units resets the offset and scale parametes. + call gargstr (Memc[cmd], SZ_LINE) + iferr { + do i = 1, nspec { + if (j > 0 && i != j) + next + sp = sps[i] + sh = SP_SH(sp) + call un_changer (UN(sh), Memc[cmd], Memr[SX(sh)], + SN(sh), YES) + SP_W0(sp) = Memr[SX(sh)] + SP_WPC(sp) = (Memr[SX(sh)+SN(sh)-1] - Memr[SX(sh)]) / + (SN(sh) - 1) + SP_XSCALE(sp) = 1. + SP_XOFFSET(sp) = 0. + call sp_scale (sp, 1, step) + if (i == 1) { + call strcpy (Memc[cmd], units, SZ_FNAME) + call gt_sets (gt, GTXLABEL, UN_LABEL(UN(sh))) + call gt_sets (gt, GTXUNITS, UN_UNITS(UN(sh))) + } + redraw = YES + } + } then + call erract (EA_WARN) + case COLOR: # Set or show color + call gargi (j) + if (current > 0) { + sp = sps[current] + if (nscan() == 1) { + call printf ("color[%d] %d") + call pargi (SP_INDEX(sp)) + call pargi (SP_COLOR(sp)) + } else { + SP_COLOR(sp) = j + } + } else { + if (nscan() == 1) { + call printf ("color:") + do i = 1, nspec { + sp = sps[i] + call printf (" %d=%d") + call pargi (SP_INDEX(sp)) + call pargi (SP_COLOR(sp)) + } + } else { + do i = 1, nspec + SP_COLOR(sps[i]) = j + } + } + default: # Print unknown command + call printf ("Unknown command: %s\007") + call pargstr (cmdstr) + } + +done_ call sfree (stack) +end + + +# SP_GDATA -- Get spectrum and add it to the array of spectrum structures. +# Return an error if the image is not found. If a two or three dimensional +# image enter each line. The spectrum data kept in memory and the image is +# closed. + +procedure sp_gdata (image, units, current, sps, nspec) + +char image[ARB] # Image name +char units[ARB] # Coordinate units +int current # Element to append +pointer sps # Pointer to array of spectra structures +int nspec # Number of spectra + +real scale # Default intensity scale +real offset # Default intensity offset +real xlpos, ylpos # Default position of labels +char ptype[SP_SZPTYPE] # Default plot type + +int i, j, k, l, m, trans +pointer sp, im, mw, sh, stack, aps, bands, str, ptr + +int ctor(), open(), fscan(), nowhite(), clgwrd() +bool rng_elementi(), fp_equalr() +real clgetr(), asumr(), imgetr(), sp_logerr() +pointer immap(), smw_openim(), rng_open() + +errchk immap, smw_openim, open + +extern sp_logerr + +begin + call smark (stack) + call salloc (str, SZ_LINE, TY_CHAR) + + # Map the image and return an error if this fails. + im = immap (image, READ_ONLY, 0) + mw = smw_openim (im) + + # Get parameters. + if (nspec == 0) { + #scale = clgetr ("scale") + #offset = clgetr ("offset") + xlpos = clgetr ("xlpos") + ylpos = clgetr ("ylpos") + call clgstr ("ptype", ptype, SP_SZPTYPE) + trans = clgwrd ("transform", Memc[str], SZ_LINE, TRANSFORMS) + } + + call clgstr ("scale", Memc[str], SZ_LINE) + if (nowhite (Memc[str], Memc[str], SZ_LINE) == 0) + call error (1, "Error in scale parameter") + if (Memc[str] == '@') { + j = open (Memc[str+1], READ_ONLY, TEXT_FILE) + do i = 1, nspec+1 + if (fscan(j) == EOF) + call error (1, "Error reading scale file") + call gargr (scale) + call close (j) + } else if (IS_ALPHA(Memc[str])) { + scale = imgetr (im, Memc[str]) + } else { + i = 1 + if (ctor (Memc[str], i, scale) == 0) + call error (1, "Error in scale parameter") + } + + call clgstr ("offset", Memc[str], SZ_LINE) + if (nowhite (Memc[str], Memc[str], SZ_LINE) == 0) + call error (1, "Error in offset parameter") + if (Memc[str] == '@') { + j = open (Memc[str+1], READ_ONLY, TEXT_FILE) + do i = 1, nspec+1 + if (fscan(j) == EOF) + call error (1, "Error reading offset file") + call gargr (offset) + call close (j) + } else if (IS_ALPHA(Memc[str])) + offset = imgetr (im, Memc[str]) + else { + i = 1 + if (ctor (Memc[str], i, offset) == 0) + call error (1, "Error in offset parameter") + } + + call clgstr ("apertures", Memc[str], SZ_LINE) + iferr (aps = rng_open (Memc[str], INDEF, INDEF, INDEF)) + call error (0, "Bad aperture/record list") + call clgstr ("bands", Memc[str], SZ_LINE) + iferr (bands = rng_open (Memc[str], INDEF, INDEF, INDEF)) + call error (0, "Bad band list") + + # For each line in the image, allocate memory for the spectrum + # structure, get the pixel data, compute the mean and limits, + # set the structure parameters, and add the structure to the + # array of structures. + + do j = 1, SMW_NBANDS(mw) { + if (SMW_FORMAT(mw) != SMW_ND) + if (!rng_elementi (bands, j)) + next + do i = 1, SMW_NSPEC(mw) { + if (SMW_FORMAT(mw) == SMW_ND) { + call smw_mw (mw, i, j, ptr, k, l) + if (!rng_elementi (aps, k) || !rng_elementi (bands, l)) + next + } else { + call shdr_open (im, mw, i, j, INDEFI, SHHDR, sh) + if (!rng_elementi (aps, AP(sh))) + next + } + call shdr_open (im, mw, i, j, INDEFI, SHDATA, sh) + iferr (call shdr_units (sh, units)) + ; + + call sp_alloc (sp, sh) + SP_NPTS(sp) = SN(sh) + SP_W0(sp) = Memr[SX(sh)] + SP_WPC(sp) = (Memr[SX(sh)+SN(sh)-1] - Memr[SX(sh)]) / + (SN(sh) - 1) + switch (trans) { + case TRANS_LOG: + SP_OMIN(sp) = MAX_REAL; SP_OMAX(sp) = -MAX_REAL + ptr = SY(sh); + do m = 1, SP_NPTS(sp) { + if (Memr[ptr] > 0.) { + SP_OMIN(sp) = min (SP_OMIN(sp), Memr[ptr]) + SP_OMAX(sp) = max (SP_OMAX(sp), Memr[ptr]) + } + ptr = ptr + 1 + } + if (SP_OMAX(sp) > 0.) { + call amaxkr (Memr[SY(sh)], SP_OMIN(sp), Memr[SY(sh)], + SN(sh)) + call alogr (Memr[SY(sh)], Memr[SY(sh)], SN(sh), + sp_logerr) + call amovr (Memr[SY(sh)], Memr[SY(SP_SH(sp))], SN(sh)) + } + } + SP_OMEAN(sp) = asumr (Memr[SY(sh)], SN(sh)) / SN(sh) + call alimr (Memr[SY(sh)], SN(sh), SP_OMIN(sp), SP_OMAX(sp)) + + SP_XSCALE(sp) = 1. + SP_XOFFSET(sp) = 0. + SP_SCALE(sp) = scale + SP_OFFSET(sp) = offset + SP_XLPOS(sp) = xlpos + SP_YLPOS(sp) = ylpos + SP_COLOR(sp) = 1 + + call sprintf (SP_IMNAME(sp), SP_SZNAME, "%s%s(%d)") + call pargstr (IMNAME(sh)) + call pargstr (IMSEC(sh)) + call pargi (AP(sh)) + call strcpy (TITLE(sh), SP_IMTITLE(sp), SP_SZTITLE) + call strcpy (ptype, SP_PTYPE(sp), SP_SZPTYPE) + SP_ULABEL(sp) = EOS + + call sp_add (sp, current, sps, nspec) + } + } + + # Close the image. + call shdr_close (sh) + call rng_close (bands) + call rng_close (aps) + call smw_close (mw) + call imunmap (im) + + call sfree (stack) +end + + +# SP_LINEAR -- Reset linear coordinates + +procedure sp_linear (sp) + +pointer sp # SPECPLOT pointer + +int i +pointer x + +begin + x = SX(SP_SH(sp)) + do i = 0, SP_NPTS(sp)-1 + Memr[x+i] = SP_W0(sp) + i * SP_WPC(sp) + SP_XSCALE(sp) = 1. + SP_XOFFSET(sp) = 0. +end + + +# SP_DELETE -- Delete a spectrum from memory. The index numbers are +# decreased to fill the hole. + +procedure sp_delete (current, sps, nspec) + +int current # Element to be deleted +pointer sps # Pointer to array of spectrum structures +int nspec # Number of spectra + +int i + +begin + if (nspec == 0) + return + + for (i = current; i < nspec; i = i + 1) { + Memi[sps+i-1] = Memi[sps+i] + SP_INDEX(Memi[sps+i-1]) = SP_INDEX(Memi[sps+i-1]) - 1 + } + nspec = nspec - 1 +end + + +# SP_ADD -- Add a spectrum structure to the array of structures +# following the specified element. The spectrum index is defined to be +# one higher than the spectrum to be followed and all higher indexed +# spectra are increased by 1. Special cases are when there are no +# spectra in which case the index is set to 1 and when the current +# element to be followed is zero. The current element is set to the +# added spectrum. The array of pointers is expanded in blocks of 100. + +procedure sp_add (sp, current, sps, nspec) + +pointer sp # Spectrum structure to be appended +int current # Element followed (in), added element (out) +pointer sps # Pointer to array of spectrum structures +int nspec # Number of spectra + +int i + +begin + # Reallocate memory for the array of structure pointers in steps of 100. + if (mod (nspec, 100) == 0) + call realloc (sps, nspec + 100, TY_POINTER) + + # Shift higher spectra in the array and increase the index numbers by 1 + # and then add the new spectrum pointer. + for (i = nspec; i > current; i = i - 1) { + Memi[sps+i] = Memi[sps+i-1] + SP_INDEX(Memi[sps+i]) = SP_INDEX(Memi[sps+i]) + 1 + } + Memi[sps+current] = sp + + # Set the new spectrum index. + if (nspec == 0) + SP_INDEX(sp) = 1 + else if (current == 0) + SP_INDEX(sp) = SP_INDEX(Memi[sps+current+1]) - 1 + else + SP_INDEX(sp) = SP_INDEX(Memi[sps+current-1]) + 1 + + # Adjust the current element and number of spectra. + current = current + 1 + nspec = nspec + 1 +end + + +# SP_LABELS -- Set the spectrum labels to the specified type. + +procedure sp_labels (sps, nspec, labels) + +pointer sps[ARB] # Spectrum pointers +int nspec # Number of spectra +int labels # Type of labels + +int i + +begin + for (i = 1; i <= nspec; i = i + 1) { + switch (labels) { + case LABEL_NONE: + SP_LABEL(sps[i]) = EOS + case LABEL_IMNAME: + call strcpy (SP_IMNAME(sps[i]), SP_LABEL(sps[i]), SP_SZLABEL) + case LABEL_IMTITLE: + call strcpy (SP_IMTITLE(sps[i]), SP_LABEL(sps[i]), SP_SZLABEL) + case LABEL_INDEX: + call sprintf (SP_LABEL(sps[i]), SP_SZLABEL, "%-4d") + call pargi (SP_INDEX(sps[i])) + case LABEL_USER: + call strcpy (SP_ULABEL(sps[i]), SP_LABEL(sps[i]), SP_SZULABEL) + } + } +end + + +# SP_ALLOC -- Allocate memory for a spectrum structure with given number of +# data points. The MWCS is not used. + +procedure sp_alloc (sp, sh) + +pointer sp # Spectrum structure pointer to be allocated +pointer sh # Spectrum header pointer + +begin + call calloc (sp, SP_LEN, TY_STRUCT) + call calloc (SP_PX(sp), SN(sh), TY_REAL) + call calloc (SP_PY(sp), SN(sh), TY_REAL) + + call shdr_copy (sh, SP_SH(sp), NO) + MW(SP_SH(sp)) = NULL +end + + +# SP_FREE -- Free a spectrum structure. + +procedure sp_free (sp) + +pointer sp, sh # Spectrum structure pointers + +begin + sh = SP_SH(sp) + call shdr_close (sh) + + call mfree (SP_PX(sp), TY_REAL) + call mfree (SP_PY(sp), TY_REAL) + call mfree (sp, TY_STRUCT) +end + + +# SP_NEAREST -- Find the nearest spectrum to the cursor and return the element. +# Return zero if no spectra are defined. The distance is in NDC. + +int procedure sp_nearest (gp, wx1, wy1, key, cmd, sps, nspec) + +pointer gp # GIO pointer +real wx1, wy1 # Cursor position +int key # Key +char cmd[ARB] # Cursor command +pointer sps[ARB] # Array of structure pointers +int nspec # Number of spectra + +int i, j, k, stridxs() +real wx0, wy0, x0, y0, x1, y1, r2, r2min +pointer sp, px, py + +begin + # Check for explicit specification. + if (key == ':') { + if (stridxs ("[", cmd) > 0) + return (1) + } + + if (IS_INDEFR(wx1)) + wx1 = 0. + if (IS_INDEFR(wy1)) + wy1 = 0. + + # Transform world cursor coordinates to NDC. + call gctran (gp, wx1, wy1, wx0, wy0, 1, 0) + + # Search for nearest point. + k = 0 + r2min = MAX_REAL + do i = 1, nspec { + sp = sps[i] + px = SP_PX(sp) - 1 + py = SP_PY(sp) - 1 + do j = 1, SP_NPTS(sp) { + x1 = Memr[px + j] + y1 = Memr[py + j] + call gctran (gp, x1, y1, x0, y0, 1, 0) + r2 = (x0 - wx0) ** 2 + (y0 - wy0) ** 2 + if (r2 < r2min) { + r2min = r2 + k = i + } + } + } + + return (k) +end + + +# SP_SHOW -- Show parameter information. Clear the screen if the output is +# to the graphics device otherwise append to the specified file. + +procedure sp_show (file, gp, sps, nspec, step) + +char file[ARB] # Optional file +pointer gp # Graphics pointer +pointer sps[ARB] # Spectra data +int nspec # Number of spectra +real step # Separation step + +int i, fd +pointer stack, line, sp + +int open() +errchk open() + +begin + fd = open (file, APPEND, TEXT_FILE) + if (gp != NULL) + call gdeactivate (gp, AW_CLEAR) + + call smark (stack) + call salloc (line, SZ_LINE, TY_CHAR) + call sysid (Memc[line], SZ_LINE) + call fprintf (fd, "%s\n\n") + call pargstr (Memc[line]) + + call fprintf (fd, "Step = %g\n\n") + call pargr (step) + call fprintf (fd, " %16s %7s %7s %7s %7s %s\n") + call pargstr ("Image Name") + call pargstr ("W0") + call pargstr ("WPC") + call pargstr ("Offset") + call pargstr ("Scale") + call pargstr ("Title") + + do i = 1, nspec { + sp = sps[i] + call fprintf (fd, "%2d %16s %7g %7g %7g %7g %s\n") + call pargi (SP_INDEX(sp)) + call pargstr (SP_IMNAME(sp)) + call pargr (SP_W0(sp)*SP_XSCALE(sp)+SP_XOFFSET(sp)) + call pargr (SP_WPC(sp)*SP_XSCALE(sp)) + call pargr (SP_OFFSET(sp)) + call pargr (SP_SCALE(sp)) + call pargstr (SP_IMTITLE(sp)) + } + + call sfree (stack) + + call close (fd) + if (gp != NULL) + call greactivate (gp, AW_PAUSE) +end + + +# SP_VSHOW -- Show verbose parameter information. Clear the screen if the +# output is to the graphics device otherwise append to the specified file. + +procedure sp_vshow (file, gp, sps, nspec, step) + +char file[ARB] # Optional file +pointer gp # Graphics pointer +pointer sps[ARB] # Spectra data +int nspec # Numbeer of spectra +real step # Separation step + +int i, fd +real z, v +pointer stack, line, sp, un + +int open() +errchk open() + +begin + fd = open (file, APPEND, TEXT_FILE) + if (gp != NULL) + call gdeactivate (gp, AW_CLEAR) + + call smark (stack) + call salloc (line, SZ_LINE, TY_CHAR) + call sysid (Memc[line], SZ_LINE) + call fprintf (fd, "%s\n\n") + call pargstr (Memc[line]) + + call fprintf (fd, "Step = %g\n") + call pargr (step) + call fprintf (fd, "\n %16s %7s %7s %7s %7s %s\n") + call pargstr ("Image Name") + call pargstr ("W0") + call pargstr ("WPC") + call pargstr ("Offset") + call pargstr ("Scale") + call pargstr ("Title") + + do i = 1, nspec { + sp = sps[i] + call fprintf (fd, "%2d %16s %7g %7g %7g %7g %s\n") + call pargi (SP_INDEX(sp)) + call pargstr (SP_IMNAME(sp)) + call pargr (SP_W0(sp)*SP_XSCALE(sp)+SP_XOFFSET(sp)) + call pargr (SP_WPC(sp)*SP_XSCALE(sp)) + call pargr (SP_OFFSET(sp)) + call pargr (SP_SCALE(sp)) + call pargstr (SP_IMTITLE(sp)) + } + + call fprintf (fd, "\n %16s %9s %9s %9s %9s\n") + call pargstr ("Image Name") + call pargstr ("Mean") + call pargstr ("DW0") + call pargstr ("Z") + call pargstr ("V(km/s)") + + un = UN(SP_SH(sps[1])) + if (UN_CLASS(un) == UN_VEL) { + do i = 1, nspec { + sp = sps[i] + z = SP_XOFFSET(sp) / UN_SCALE(un) + v = SP_XOFFSET(sp) + call fprintf (fd, "%2d %16s %9g %9g %9g %9g\n") + call pargi (SP_INDEX(sp)) + call pargstr (SP_IMNAME(sp)) + call pargr (SP_OMEAN(sp)) + call pargr (SP_XOFFSET(sp)) + call pargr (z) + call pargr (v) + } + } else if (UN_CLASS(un) == UN_WAVE) { + do i = 1, nspec { + sp = sps[i] + if (SP_XSCALE(sp) != 1.) { + call sprintf (Memc[line], SZ_LINE, "km/s %g %s") + call pargr (SP_W0(sp)) + call pargstr (UN_UNITS(un)) + z = SP_XSCALE(sp) - 1 + v = SP_W0(sp) * SP_XSCALE(sp) + call un_changer (un, Memc[line], v, 1, NO) + } else { + z = 0. + v = 0. + } + call fprintf (fd, "%2d %16s %9g %9g %9g %9g\n") + call pargi (SP_INDEX(sp)) + call pargstr (SP_IMNAME(sp)) + call pargr (SP_OMEAN(sp)) + call pargr (SP_XOFFSET(sp)) + call pargr (z) + call pargr (v) + } + } else if (UN_CLASS(un) == UN_FREQ || UN_CLASS(un) == UN_ENERGY) { + do i = 1, nspec { + sp = sps[i] + if (SP_XSCALE(sp) != 1.) { + call sprintf (Memc[line], SZ_LINE, "km/s %g %s") + call pargr (SP_W0(sp)) + call pargstr (UN_UNITS(un)) + z = 1. / SP_XSCALE(sp) - 1 + v = SP_W0(sp) * SP_XSCALE(sp) + call un_changer (un, Memc[line], v, 1, NO) + } else { + z = 0. + v = 0. + } + call fprintf (fd, "%2d %16s %9g %9g %9g %9g\n") + call pargi (SP_INDEX(sp)) + call pargstr (SP_IMNAME(sp)) + call pargr (SP_OMEAN(sp)) + call pargr (SP_XOFFSET(sp)) + call pargr (z) + call pargr (v) + } + } + + call sfree (stack) + + call close (fd) + if (gp != NULL) + call greactivate (gp, AW_PAUSE) +end + + +# SP_LOGERR -- Value for non-positive values in log function. + +real procedure sp_logerr (x) + +real x + +begin + return (0.) +end -- cgit