diff options
Diffstat (limited to 'noao/onedspec/ecidentify/ecffit')
| -rw-r--r-- | noao/onedspec/ecidentify/ecffit/ecfcolon.x | 102 | ||||
| -rw-r--r-- | noao/onedspec/ecidentify/ecffit/ecfeval.x | 68 | ||||
| -rw-r--r-- | noao/onedspec/ecidentify/ecffit/ecffit.com | 23 | ||||
| -rw-r--r-- | noao/onedspec/ecidentify/ecffit/ecffit.h | 20 | ||||
| -rw-r--r-- | noao/onedspec/ecidentify/ecffit/ecffit.key | 53 | ||||
| -rw-r--r-- | noao/onedspec/ecidentify/ecffit/ecffit.x | 193 | ||||
| -rw-r--r-- | noao/onedspec/ecidentify/ecffit/ecfgdata.x | 37 | ||||
| -rw-r--r-- | noao/onedspec/ecidentify/ecffit/ecfget.x | 84 | ||||
| -rw-r--r-- | noao/onedspec/ecidentify/ecffit/ecfgraph.x | 50 | ||||
| -rw-r--r-- | noao/onedspec/ecidentify/ecffit/ecfnearest.x | 85 | ||||
| -rw-r--r-- | noao/onedspec/ecidentify/ecffit/ecfreject.x | 53 | ||||
| -rw-r--r-- | noao/onedspec/ecidentify/ecffit/ecfrms.x | 26 | ||||
| -rw-r--r-- | noao/onedspec/ecidentify/ecffit/ecfset.x | 92 | ||||
| -rw-r--r-- | noao/onedspec/ecidentify/ecffit/ecfshift.x | 55 | ||||
| -rw-r--r-- | noao/onedspec/ecidentify/ecffit/ecfsolve.x | 196 | ||||
| -rw-r--r-- | noao/onedspec/ecidentify/ecffit/ecftitle.x | 48 | ||||
| -rw-r--r-- | noao/onedspec/ecidentify/ecffit/mkpkg | 21 |
17 files changed, 1206 insertions, 0 deletions
diff --git a/noao/onedspec/ecidentify/ecffit/ecfcolon.x b/noao/onedspec/ecidentify/ecffit/ecfcolon.x new file mode 100644 index 00000000..4307335b --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecfcolon.x @@ -0,0 +1,102 @@ +include <error.h> +include <gset.h> + +# List of colon commands +define CMDS "|show|function|xorder|yorder|niterate|lowreject|highreject|" + +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 NITERATE 5 # Set or show rejection iterations +define LOW 6 # Set or show low rejection threshold +define HIGH 7 # Set or show high rejection threshold + +# ECF_COLON -- Processes colon commands. + +procedure ecf_colon (cmdstr, gp) + +char cmdstr[ARB] # Command string +pointer gp # GIO pointer + +double dval +int ncmd, ival +int nscan(), strdic() +include "ecffit.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 (ecfstr, SZ_LINE) + ncmd = strdic (ecfstr, ecfstr, SZ_LINE, CMDS) + + switch (ncmd) { + case SHOW: # :show - Show the values of the fitting parameters. + call gdeactivate (gp, AW_CLEAR) + call printf ("function %s\nxorder %d\nyorder %d\n") + call pargstr (function) + call pargi (xorder) + call pargi (yorder) + call printf ("niterate %d\nlowreject %g\nhighreject\nnreject %d\n") + call pargi (niterate) + call pargd (low) + call pargd (high) + call pargi (nreject) + call printf ("slope %d\noffset %d\nshift %g\n") + call pargi (slope) + call pargi (offset) + call pargd (shift) + call printf ("rms %g\n") + call pargd (rms) + call greactivate (gp, AW_PAUSE) + case FUNCTION: # :function - List or set the fitting function. + call gargwrd (ecfstr, SZ_LINE) + if (nscan() == 1) { + call printf ("function = %s\n") + call pargstr (function) + } else { + iferr (call ecf_sets ("function", ecfstr)) + call erract (EA_WARN) + } + case XORDER: # xorder: List or set the function order. + call gargi (ival) + if (nscan() == 1) { + call printf ("xorder %d\n") + call pargi (xorder) + } 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 + yorder = ival + case NITERATE: # niterate: List or set rejection iterations. + call gargi (ival) + if (nscan() == 1) { + call printf ("niterate %d\n") + call pargi (niterate) + } else + niterate = ival + case LOW: # low: List or set low rejection threshold. + call gargd (dval) + if (nscan() == 1) { + call printf ("lowreject %g\n") + call pargd (low) + } else + low = dval + case HIGH: # highreject: List or set high rejection threshold. + call gargd (dval) + if (nscan() == 1) { + call printf ("highreject %g\n") + call pargd (high) + } else + high = dval + default: + call printf ("Unrecognized or ambiguous command\007") + } +end diff --git a/noao/onedspec/ecidentify/ecffit/ecfeval.x b/noao/onedspec/ecidentify/ecffit/ecfeval.x new file mode 100644 index 00000000..1901522f --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecfeval.x @@ -0,0 +1,68 @@ +# ECF_EVAL -- Evaluate wavelength at a given order and pixel position. + +double procedure ecf_eval (ecf, order, x) + +pointer ecf # GSURFIT pointer +int order # Order +double x # X point + +int ecf_oeval() +double y, dgseval() +include "ecffit.com" + +begin + y = ecf_oeval (ecf, order) + if (ecf == NULL) + return (x + shift / y) + else + return ((dgseval (ecf, x, y) + shift) / y) +end + + +# ECF_VECTOR -- Evaluate echelle dispersion function for a vector of points of +# the same order. + +procedure ecf_vector (ecf, order, x, fit, npts) + +pointer ecf # GSURFIT pointer +int order # Order +double x[npts] # X points +double fit[npts] # Fitted points +int npts # Number of points + +double yval +pointer sp, y +int ecf_oeval() +include "ecffit.com" + +begin + call smark (sp) + call salloc (y, npts, TY_DOUBLE) + + yval = ecf_oeval (ecf, order) + if (ecf == NULL) + call amovd (x, fit, npts) + else { + call amovkd (yval, Memd[y], npts) + call dgsvector (ecf, x, Memd[y], fit, npts) + call adivkd (fit, yval, fit, npts) + } + if (shift != 0.) + call aaddkd (fit, shift / yval, fit, npts) + + call sfree (sp) +end + + +# ECF_OEVAL -- Evaluate the fit order. + +int procedure ecf_oeval (ecf, order) + +pointer ecf # GSURFIT pointer +int order # User order + +include "ecffit.com" + +begin + return (slope * order + offset) +end diff --git a/noao/onedspec/ecidentify/ecffit/ecffit.com b/noao/onedspec/ecidentify/ecffit/ecffit.com new file mode 100644 index 00000000..61f3104a --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecffit.com @@ -0,0 +1,23 @@ +# Common parameters. + +char function[SZ_FNAME] # Fitting function +char ecfstr[SZ_LINE] # Working char string +int gstype # Surface function type +int xorder # X order of surface function +int yorder # Y order of surface function +int niterate # Number of rejection iterations +int nreject # Number of rejected points +int xtype # X axis type +int ytype # Y axis type +int slope # Slope of order +int offset # Order offset of fit +double low, high # Low and high rejection thresholds +double xmin, xmax # X range +double ymin, ymax # Y range +double shift # First order shift +double rms # RMS of fit + + +common /ecfcom/ low, high, xmin, xmax, ymin, ymax, shift, rms, gstype, + xorder, yorder, niterate, nreject, xtype, ytype, slope, offset, + function, ecfstr diff --git a/noao/onedspec/ecidentify/ecffit/ecffit.h b/noao/onedspec/ecidentify/ecffit/ecffit.h new file mode 100644 index 00000000..20825c71 --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecffit.h @@ -0,0 +1,20 @@ +define IGSPARAMS 7 + +define FEATURE 1 +define X 2 +define Y 3 +define Z 4 +define W 5 +define S 6 +define R 7 + +define IGS_FUNCTION 1 +define IGS_XORDER 2 +define IGS_YORDER 3 +define IGS_XMIN 4 +define IGS_XMAX 5 +define IGS_YMIN 6 +define IGS_YMAX 7 +define IGS_OFFSET 8 + +define SFTYPES "|chebyshev|legendre|" # Surface types diff --git a/noao/onedspec/ecidentify/ecffit/ecffit.key b/noao/onedspec/ecidentify/ecffit/ecffit.key new file mode 100644 index 00000000..f24407b9 --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecffit.key @@ -0,0 +1,53 @@ + ECHELLE DISPERSION FUNCTION FITTING KEYS + + +CURSOR KEY SUMMARY + +? Help c Print coordinates d Delete point +f Fit dispersion o Fit with fixed order offset q Quit +r Redraw graph u Undelete point w Window graph +x Set ordinate y Set abscissa I Interrupt + + +COLON COMMAND SUMMARY + +:show :function [value] :highreject [value] :lowreject [value] +:niterate [value] :xorder [value] :yorder [value] + + +CURSOR KEYS + +? Print this list of cursor keys +c Print cursor coordinates +d Delete the nearest undeleted point to the cursor +f Fit dispersion function including determining the order offset +o Fit dispersion function with the order offset fixed +q Quit and return to the spectrum display +r Redraw the graph +u Undelete the nearest deleted point to the cursor (may be outside the window) +w Window the graph (type ? to the window prompt for more help) +x Set the quantity plotted along the ordinate (x axis) +y Set the quantity plotted along the abscissa (y axis) +I Interrupt the task immediately + + +COLON COMMANDS + +:show Print current function and orders +:function [value] Print or set the function type (chebyshev|legendre) +:highreject [value] Print or set high rejection limit +:lowreject [value] Print or set high rejection limit +:niterate [value] Print or set number of rejection iterations +:xorder [value] Print or set the order for the dispersion dependence +:yorder [value] Print or set the order for the echelle order dependence + + +The dispersion function fitted is given by a two dimensional function +(either chebyshev or legendre) of the pixel position along the +dispersion of an order (called x) and the order number (called y). The +order number is determined from the aperture number by an offset and +direction of increasing order number. The basic order dependence is +separated from the surface function as given below. + + y = offset +/- aperture + wavelength = f (x, y) / y diff --git a/noao/onedspec/ecidentify/ecffit/ecffit.x b/noao/onedspec/ecidentify/ecffit/ecffit.x new file mode 100644 index 00000000..408a1b77 --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecffit.x @@ -0,0 +1,193 @@ +include <error.h> +include <pkg/gtools.h> + +define HELP "noao$onedspec/ecidentify/ecffit/ecffit.key" +define PROMPT "fitcoords surface fitting options" + +# EC_FIT -- Echelle dispersion fitting. +# +# X - Pixel coordinates along dispersion +# Y - Relative order number +# Z - Wavelength + +procedure ecf_fit (ecf, gp, gt, xd, yd, zd, wd, npts, fixedorder) + +pointer ecf # GSURFIT pointer +pointer gp # GIO pointer +pointer gt # GTOOLS pointer +double xd[npts] # Pixel coordinates along dispersion +double yd[npts] # Order number +double zd[npts] # Wavelength +double wd[npts] # Weights +int npts # Number of points +int fixedorder # Fixed order? + +real wx, wy +int wcs, key +int i, newgraph +pointer sp, wd1, rd, xr, yr +char cmd[SZ_LINE] + +int ecf_nearest() +int clgcur(), scan(), nscan() +errchk ecf_solve() +include "ecffit.com" + +begin + # Allocate residuals and weights with rejected points arrays + call smark (sp) + call salloc (wd1, npts, TY_DOUBLE) + call salloc (rd, npts, TY_DOUBLE) + call amovd (wd, Memd[wd1], npts) + + # Compute a solution and return if not interactive. + if (gp == NULL) { + call ecf_solve (ecf, xd, yd, zd, Memd[wd1], Memd[rd], npts, + fixedorder) + call ecf_reject (ecf, xd, yd, zd, Memd[wd1], Memd[rd], npts, + fixedorder) + do i = 1, npts + if (Memd[wd1+i-1] != wd[i]) + wd[i] = -1. + call sfree (sp) + return + } + + # Allocate real graph vectors. + call salloc (xr, npts, TY_REAL) + call salloc (yr, npts, TY_REAL) + + # Read cursor commands. + key = 'f' + repeat { + switch (key) { + case 'o': + call printf ("Order offset (%d): ") + call pargi (offset) + call flush (STDOUT) + if (scan() != EOF) { + call gargi (i) + if (nscan() == 1) + offset = i + call amovd (wd, Memd[wd1], npts) + call ecf_solve (ecf, xd, yd, zd, Memd[wd1], Memd[rd], npts, + YES) + call ecf_reject (ecf, xd, yd, zd, Memd[wd1], Memd[rd], npts, + YES) + call ecf_gdata (ecf, xtype, xd, yd, zd, Memd[rd], + Memr[xr], npts) + call ecf_gdata (ecf, ytype, xd, yd, zd, Memd[rd], + Memr[yr], npts) + call ecf_title (gt) + newgraph = YES + } + + 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 ecf_colon (cmd, gp) + + case 'x': # Set ordinate + call printf ("Ordinate - ") + call printf ( + "(p)ixel, (o)rder, (w)avelength, (r)esidual, (v)elocity: ") + if (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF) + break + + if (key != xtype) { + if (key=='p'||key=='o'||key=='w'||key=='r'||key=='v') { + xtype = key + call gt_setr (gt, GTXMIN, INDEF) + call gt_setr (gt, GTXMAX, INDEF) + call ecf_gdata (ecf, xtype, xd, yd, zd, Memd[rd], + Memr[xr], npts) + call ecf_title (gt) + newgraph = YES + } else + call printf ("\007") + } + + case 'y': # Set abscissa + call printf ("Abscissa - ") + call printf ( + "(p)ixel, (o)rder, (w)avelength, (r)esidual, (v)elocity: ") + if(clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF) + break + + if (key != ytype) { + if (key=='p'||key=='o'||key=='w'||key=='r'||key=='v') { + ytype = key + call gt_setr (gt, GTYMIN, INDEF) + call gt_setr (gt, GTYMAX, INDEF) + call ecf_gdata (ecf, ytype, xd, yd, zd, Memd[rd], + Memr[yr], npts) + call ecf_title (gt) + newgraph = YES + } else + call printf ("\007") + } + + case 'r': # Redraw + newgraph = YES + + case 'c': # Cursor coordinates + i = ecf_nearest (gp, gt, wx, wy, wcs, key, Memr[xr], Memr[yr], + wd, npts) + call printf ("%10.2g %d %10.8g\n") + call pargd (xd[i]) + call pargd (yd[i]) + call pargd (zd[i]) + + case 'd': # Delete + i = ecf_nearest (gp, gt, wx, wy, wcs, key, Memr[xr], Memr[yr], + wd, npts) + if (i > 0) + Memd[wd1+i-1] = wd[i] + + case 'u': # Undelete + i = ecf_nearest (gp, gt, wx, wy, wcs, key, Memr[xr], Memr[yr], + wd, npts) + if (i > 0) + Memd[wd1+i-1] = wd[i] + + case 'f': # Fit + call amovd (wd, Memd[wd1], npts) + call ecf_solve (ecf, xd, yd, zd, Memd[wd1], Memd[rd], npts, + fixedorder) + call ecf_reject (ecf, xd, yd, zd, Memd[wd1], Memd[rd], npts, + fixedorder) + call ecf_gdata (ecf, xtype, xd, yd, zd, Memd[rd], + Memr[xr], npts) + call ecf_gdata (ecf, ytype, xd, yd, zd, Memd[rd], + Memr[yr], npts) + call ecf_title (gt) + newgraph = YES + + case 'w': # Window graph + call gt_window (gt, gp, "cursor", newgraph) + + case 'q': # Quit + break + + case 'I': # Interrupt + call fatal (0, "Interrupt") + + default: # Ring the bell. + call printf ("\07\n") + } + + if (newgraph == YES) { + call ecf_graph (gp, gt, Memr[xr], Memr[yr], wd, Memd[wd1], npts) + newgraph = NO + } + } until (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF) + + do i = 1, npts + if (Memd[wd1+i-1] != wd[i]) + wd[i] = -1. + call sfree (sp) +end diff --git a/noao/onedspec/ecidentify/ecffit/ecfgdata.x b/noao/onedspec/ecidentify/ecffit/ecfgdata.x new file mode 100644 index 00000000..eebb34d6 --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecfgdata.x @@ -0,0 +1,37 @@ +include <pkg/gtools.h> + +# ECF_GDATA -- Get graph data for the specified axis type from the fitting data. + +procedure ecf_gdata (ecf, type, x, y, z, r, data, npts) + +pointer ecf # GSURFIT pointer +int type # Axis type +double x[npts] # X fit data +double y[npts] # Y fit data +double z[npts] # Z fit data +double r[npts] # Residuals +real data[npts] # Graph data +int npts # Number of points + +pointer sp, v +include "ecffit.com" + +begin + switch (type) { + case 'p': + call achtdr (x, data, npts) + case 'o': + call achtdr (y, data, npts) + case 'w': + call achtdr (z, data, npts) + case 'r': + call achtdr (r, data, npts) + case 'v': + call smark (sp) + call salloc (v, npts, TY_DOUBLE) + call adivd (r, z, Memd[v], npts) + call amulkd (Memd[v], 300000.D0, Memd[v], npts) + call achtdr (Memd[v], data, npts) + call sfree (sp) + } +end diff --git a/noao/onedspec/ecidentify/ecffit/ecfget.x b/noao/onedspec/ecidentify/ecffit/ecfget.x new file mode 100644 index 00000000..025059df --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecfget.x @@ -0,0 +1,84 @@ +# ECF_GETI -- Get the value of an integer parameter. + +int procedure ecf_geti (param) + +char param[ARB] # ECF parameter + +int i, strdic() +include "ecffit.com" + +begin + i = strdic (param, ecfstr, SZ_LINE, + "|slope|offset|xorder|yorder|niterate|") + switch (i) { + case 1: + return (slope) + case 2: + return (offset) + case 3: + return (xorder) + case 4: + return (yorder) + case 5: + return (niterate) + default: + call error (0, "ecf_geti: Unknown parameter") + } +end + + +# ECF_GETS -- Get the value of a string parameter. + +procedure ecf_gets (param, str, maxchar) + +char param[ARB] # ECF parameter +char str[maxchar] # String +int maxchar # Maximum number of characters + +int i, strdic() +include "ecffit.com" + +begin + i = strdic (param, ecfstr, SZ_LINE, "|function|") + switch (i) { + case 1: + call strcpy (function, str, maxchar) + default: + call error (0, "ecf_gets: Unknown parameter") + } +end + + +# ECF_GETD -- Get the values of double valued fitting parameters. + +double procedure ecf_getd (param) + +char param[ARB] # ECF parameter + +int i, strdic() +include "ecffit.com" + +begin + i = strdic (param, ecfstr, SZ_LINE, + "|xmin|xmax|ymin|ymax|shift|rms|low|high|") + switch (i) { + case 1: + return (xmin) + case 2: + return (xmax) + case 3: + return (ymin) + case 4: + return (ymax) + case 5: + return (shift) + case 6: + return (rms) + case 7: + return (low) + case 8: + return (high) + default: + call error (0, "ecf_gets: Unknown parameter") + } +end diff --git a/noao/onedspec/ecidentify/ecffit/ecfgraph.x b/noao/onedspec/ecidentify/ecffit/ecfgraph.x new file mode 100644 index 00000000..22749527 --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecfgraph.x @@ -0,0 +1,50 @@ +include <gset.h> +include <mach.h> +include <pkg/gtools.h> + +# ECF_GRAPH -- Graph the fitted data. + +procedure ecf_graph (gp, gt, x, y, w, rej, npts) + +pointer gp # GIO pointer +pointer gt # GTOOLS pointer +real x[npts] # X data +real y[npts] # Y data +double w[npts] # Weights +double rej[npts] # Rejected points +int npts # Number of pts points + +int i +real xsize, ysize, ymin, ymax, gt_getr() + +begin + xsize = gt_getr (gt, GTXSIZE) + ysize = gt_getr (gt, GTYSIZE) + + call gclear (gp) + + ymin = MAX_REAL + ymax = -MAX_REAL + do i = 1, npts + if (w[i] > 0.) { + ymin = min (ymin, y[i]) + ymax = max (ymax, y[i]) + } + + call gascale (gp, x, npts, 1) + call gswind (gp, INDEF, INDEF, ymin, ymax) + call gt_swind (gp, gt) + call gt_labax (gp, gt) + + do i = 1, npts { + if (rej[i] == 0.) { + if (y[i] >= ymin && y[i] <= ymax) { + if (w[i] == 0.) + call gmark (gp, x[i], y[i], GM_CROSS, xsize, ysize) + else + call gmark (gp, x[i], y[i], GM_DIAMOND, xsize, ysize) + } + } else + call gmark (gp, x[i], y[i], GM_PLUS, xsize, ysize) + } +end diff --git a/noao/onedspec/ecidentify/ecffit/ecfnearest.x b/noao/onedspec/ecidentify/ecffit/ecfnearest.x new file mode 100644 index 00000000..af1b1f78 --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecfnearest.x @@ -0,0 +1,85 @@ +include <mach.h> +include <gset.h> +include <pkg/gtools.h> + +# ECF_NEAREST -- Find nearest point to the cursor. + +int procedure ecf_nearest (gp, gt, wx, wy, wcs, key, x, y, w, npts) + +pointer gp # GIO pointer +pointer gt # GTOOLS pointer +real wx, wy # Cursor coordinates +int wcs # WCS +int key # Nearest key +real x[npts] # Data points +real y[npts] # Data points +double w[npts] # Weight +int npts # Number of data points + +int i, j +real r2, r2min, x0, y0, xsize, ysize, gt_getr() + +begin + call gctran (gp, wx, wy, wx, wy, wcs, 0) + r2min = MAX_REAL + j = 0 + + switch (key) { + case 'c': + do i = 1, npts { + call gctran (gp, x[i], y[i], x0, y0, wcs, 0) + r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2 + if (r2 < r2min) { + r2min = r2 + j = i + } + } + call gscur (gp, x[j], y[j]) + case 'd': + do i = 1, npts { + if (w[i] == 0.) + next + call gctran (gp, x[i], y[i], x0, y0, wcs, 0) + r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2 + if (r2 < r2min) { + r2min = r2 + j = i + } + } + if (j > 0) { + xsize = gt_getr (gt, GTXSIZE) + ysize = gt_getr (gt, GTYSIZE) + + call gseti (gp, G_PMLTYPE, 0) + call gmark (gp, x[j], y[j], GM_PLUS, xsize, ysize) + call gseti (gp, G_PMLTYPE, 1) + call gmark (gp, x[j], y[j], GM_CROSS, xsize, ysize) + w[j] = 0. + call gscur (gp, x[j], y[j]) + } + case 'u': + do i = 1, npts { + if (w[i] != 0.) + next + call gctran (gp, x[i], y[i], x0, y0, wcs, 0) + r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2 + if (r2 < r2min) { + r2min = r2 + j = i + } + } + if (j > 0) { + xsize = gt_getr (gt, GTXSIZE) + ysize = gt_getr (gt, GTYSIZE) + + call gseti (gp, G_PMLTYPE, 0) + call gmark (gp, x[j], y[j], GM_CROSS, xsize, ysize) + call gseti (gp, G_PMLTYPE, 1) + call gmark (gp, x[j], y[j], GM_PLUS, xsize, ysize) + w[j] = 1. + call gscur (gp, x[j], y[j]) + } + } + + return (j) +end diff --git a/noao/onedspec/ecidentify/ecffit/ecfreject.x b/noao/onedspec/ecidentify/ecffit/ecfreject.x new file mode 100644 index 00000000..a772069e --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecfreject.x @@ -0,0 +1,53 @@ +include <mach.h> + +# ECF_REJECT -- Reject points with large residuals from the fit. + +procedure ecf_reject (ecf, x, y, z, w, r, npts, fixedorder) + +pointer ecf # GSURFIT pointer +double x[npts] # X points +double y[npts] # Y points +double z[npts] # Z points +double w[npts] # Weights +double r[npts] # Residuals +int npts # Number of points +int fixedorder # Fixed order? + +int i, j, newreject +double low_cut, high_cut +include "ecffit.com" + +begin + # Return if rejection is not desired. + nreject = 0 + if (niterate == 0 || (low == 0. && high == 0.)) + return + + # Reject points. + do i = 1, niterate { + if (low > 0.) + low_cut = -low * rms + else + low_cut = -MAX_REAL + if (high > 0.) + high_cut = high * rms + else + high_cut = MAX_REAL + + newreject = 0 + do j = 1, npts { + if (w[j] == 0.) + next + if ((r[j] > high_cut) || (r[j] < low_cut)) { + w[j] = 0. + newreject = newreject + 1 + } + } + + if (newreject == 0) + break + + call ecf_solve (ecf, x, y, z, w, r, npts, fixedorder) + nreject = nreject + newreject + } +end diff --git a/noao/onedspec/ecidentify/ecffit/ecfrms.x b/noao/onedspec/ecidentify/ecffit/ecfrms.x new file mode 100644 index 00000000..1140dc29 --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecfrms.x @@ -0,0 +1,26 @@ +# ECF_RMS -- Compute the rms with deleted points ignored. + +double procedure ecf_rms (r, w, npts) + +double r[npts] # Residuals +double w[npts] # Weights +int npts # Number of points + +int i, n +double rms + +begin + n = 0 + rms = 0. + do i = 1, npts { + if (w[i] == 0.) + next + n = n + 1 + rms = rms + r[i] * r[i] + } + if (n > 0) + rms = sqrt (rms / n) + else + rms = INDEFD + return (rms) +end diff --git a/noao/onedspec/ecidentify/ecffit/ecfset.x b/noao/onedspec/ecidentify/ecffit/ecfset.x new file mode 100644 index 00000000..4b6402b1 --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecfset.x @@ -0,0 +1,92 @@ +# ECF_SETS -- Set the values of string valued fitting parameters. + +procedure ecf_sets (param, str) + +char param[ARB] # Parameter to be set +char str[ARB] # String value + +char temp[10] +int i, strdic() +include "ecffit.com" + +begin + i = strdic (param, temp, 10, "|function|") + switch (i) { + case 1: + i = strdic (str, str, SZ_FNAME, "|chebyshev|legendre|") + if (i == 0) + call error (0, "Unknown function type") + call strcpy (str, function, SZ_LINE) + gstype = i + default: + call error (0, "ecf_sets: Unknown parameter") + } +end + + +# ECF_SETI -- Set the values of integer valued fitting parameters. + +procedure ecf_seti (param, ival) + +char param[ARB] # Parameter to be set +int ival # Integer value + +int i, strdic() +include "ecffit.com" + +begin + i = strdic (param, ecfstr, SZ_LINE, + "|slope|offset|xorder|yorder|xtype|ytype|niterate|") + switch (i) { + case 1: + slope = ival + case 2: + offset = ival + case 3: + xorder = ival + case 4: + yorder = ival + case 5: + xtype = ival + case 6: + ytype = ival + case 7: + niterate = max (0, ival) + default: + call error (0, "ecf_seti: Unknown parameter") + } +end + + +# ECF_SETD -- Set the values of double valued fitting parameters. + +procedure ecf_setd (param, dval) + +char param[ARB] # Parameter to be set +double dval # Double value + +int i, strdic() +include "ecffit.com" + +begin + i = strdic (param, ecfstr, SZ_LINE, + "|xmin|xmax|ymin|ymax|shift|low|high|") + switch (i) { + case 1: + xmin = dval + case 2: + xmax = dval + case 3: + ymin = dval + case 4: + ymax = dval + case 5: + shift = dval + case 6: + low = max (0.D0, dval) + case 7: + high = max (0.D0, dval) + default: + call error (0, "ecf_setd: Unknown parameter") + } +end diff --git a/noao/onedspec/ecidentify/ecffit/ecfshift.x b/noao/onedspec/ecidentify/ecffit/ecfshift.x new file mode 100644 index 00000000..75655703 --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecfshift.x @@ -0,0 +1,55 @@ +# ECF_GSHIFT -- Return the shift for the given order. + +double procedure ecf_gshift (ecf, order) + +pointer ecf # GSURFIT pointer +int order # User order + +include "ecffit.com" + +begin + return (shift / (slope * order + offset)) +end + + +# ECF_PSHIFT -- Put the shift for the given order. + +procedure ecf_pshift (ecf, order, shft) + +pointer ecf # GSURFIT pointer +int order # User order +double shft # Shift at given order + +include "ecffit.com" + +begin + shift = shft * (slope * order + offset) +end + + +procedure ecf_vector (ecf, order, x, fit, npts) + +pointer ecf # GSURFIT pointer +int order # Order +double x[npts] # X points +double fit[npts] # Fitted points +int npts # Number of points + +double yval +pointer sp, y + +include "ecffit.com" + +begin + call smark (sp) + call salloc (y, npts, TY_DOUBLE) + + yval = slope * order + offset + call amovkd (yval, Memd[y], npts) + call dgsvector (ecf, x, Memd[y], fit, npts) + call adivkd (fit, yval, fit, npts) + if (shift != 0.) + call aaddkd (fit, shift / yval, fit, npts) + + call sfree (sp) +end diff --git a/noao/onedspec/ecidentify/ecffit/ecfsolve.x b/noao/onedspec/ecidentify/ecffit/ecfsolve.x new file mode 100644 index 00000000..1c844e76 --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecfsolve.x @@ -0,0 +1,196 @@ +include <mach.h> +include <math/gsurfit.h> + +define ECFTYPES "|chebyshev|legendre|" # Fit types + + +# ECF_SOLVE -- Fit +# +# f(x, slope*y+offset) = (y+slope*offset)*z +# +# with offset minimizing the RMS. + +procedure ecf_solve (ecf, x, y, z, w, r, npts, fixedorder) + +pointer ecf # GSURFIT pointer +double x[npts] # X points +double y[npts] # Y points +double z[npts] # Z points +double w[npts] # Weights +double r[npts] # Residuals +int npts # Number of points +int fixedorder # Fixed order? + +int i, j, k, err +double ya, yb, newrms, ecf_rms() +pointer sp, y1, ecf1 +errchk ecf_solve1 +include "ecffit.com" +define fit_ 99 + +begin + if (fixedorder == YES) { + call ecf_solve1 (ecf, x, y, z, w, r, npts) + return + } + + call smark (sp) + call salloc (y1, npts, TY_DOUBLE) + + # Determine if the orders are reversed. + j = 1 + k = 1 + do i = 1, npts { + if (z[i] < z[j]) + j = i + if (z[i] > z[k]) + k = i + } + if (y[j] >= y[k]) { + slope = 1 + offset = max (offset, int(1. - ymin)) + } else { + slope = -1 + offset = max (offset, int(1. + ymax)) + } + + call dgsfree (ecf) + shift = 0. + + rms = MAX_DOUBLE + j = 1 + k = 0 + + for (i=offset;;i=i+j) { + if (slope == 1) { + ya = i + ymin + yb = i + ymax + } else { + ya = i - ymax + yb = i - ymin + } + if (ya < 1.) + break + + call altmd (y, Memd[y1], npts, double(slope), double(i)) + call amuld (Memd[y1], z, r, npts) + +fit_ call dgsinit (ecf1, gstype, xorder, yorder, YES, xmin, xmax, ya, yb) + call dgsfit (ecf1, x, Memd[y1], r, w, npts, WTS_USER, err) + + if (err != OK) { + if (xorder > 2 || yorder > 2) { + call dgsfree (ecf) + xorder = max (2, xorder - 1) + yorder = max (2, yorder - 1) + goto fit_ + } + + switch (err) { + case SINGULAR: + call dgsfree (ecf) + ecf = ecf1 + call eprintf ("Singular solution\n") + case NO_DEG_FREEDOM: + call sfree (sp) + call error (0, "No degrees of freedom") + } + } + + call dgsvector (ecf1, x, Memd[y1], r, npts) + call adivd (r, Memd[y1], r, npts) + call asubd (z, r, r, npts) + + newrms = ecf_rms (r, w, npts) + k = k + 1 + + if (newrms / rms < 0.999) { + call dgsfree (ecf) + ecf = ecf1 + offset = i + rms = newrms + } else { + call dgsfree (ecf1) + if (k > 2) + break + i = offset + j = -j + } + } + + call altmd (y, Memd[y1], npts, double(slope), double(offset)) + call dgsvector (ecf, x, Memd[y1], r, npts) + call adivd (r, Memd[y1], r, npts) + call asubd (z, r, r, npts) + + call sfree (sp) + +end + + +# ECF_SOLVE1 -- Fit f(x, y+offset) = (y+offset)*z with offset fixed. + +procedure ecf_solve1 (ecf, x, y, z, w, r, npts) + +pointer ecf # GSURFIT pointer +double x[npts] # X points +double y[npts] # Y points +double z[npts] # Z points +double w[npts] # Weights +double r[npts] # Residuals +int npts # Number of points + +int err +pointer sp, y1 +double ya, yb, ecf_rms() +include "ecffit.com" +define fit_ 99 + +begin + call smark (sp) + call salloc (y1, npts, TY_DOUBLE) + + call dgsfree (ecf) + shift = 0. + + if (slope == 1) { + offset = max (offset, int(1. - ymin)) + ya = offset + ymin + yb = offset + ymax + } else { + offset = max (offset, int(1. + ymax)) + ya = offset - ymax + yb = offset - ymin + } + + call altmd (y, Memd[y1], npts, double (slope), double (offset)) + call amuld (Memd[y1], z, r, npts) + +fit_ call dgsinit (ecf, gstype, xorder, yorder, YES, xmin, xmax, + min (ya, yb), max (ya, yb)) + call dgsfit (ecf, x, Memd[y1], r, w, npts, WTS_USER, err) + + if (err != OK) { + if (xorder > 2 || yorder > 2) { + call dgsfree (ecf) + xorder = max (2, xorder - 1) + yorder = max (2, yorder - 1) + goto fit_ + } + + switch (err) { + case SINGULAR: + call eprintf ("Singular solution\n") + case NO_DEG_FREEDOM: + call sfree (sp) + call error (0, "No degrees of freedom") + } + } + + call dgsvector (ecf, x, Memd[y1], r, npts) + call adivd (r, Memd[y1], r, npts) + call asubd (z, r, r, npts) + rms = ecf_rms (r, w, npts) + + call sfree (sp) +end diff --git a/noao/onedspec/ecidentify/ecffit/ecftitle.x b/noao/onedspec/ecidentify/ecffit/ecftitle.x new file mode 100644 index 00000000..3b754f31 --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecftitle.x @@ -0,0 +1,48 @@ +include <pkg/gtools.h> + +# ECF_TITLE -- Set the GTOOLS parameter string. + +procedure ecf_title (gt) + +pointer gt # GTOOLS pointer + +include "ecffit.com" + +begin + call sprintf (ecfstr, SZ_LINE, + "Function=%s, xorder=%d, yorder=%d, slope=%d, offset=%d, rms=%6g") + call pargstr (function) + call pargi (xorder) + call pargi (yorder) + call pargi (slope) + call pargi (offset) + call pargd (rms) + call gt_sets (gt, GTPARAMS, ecfstr) + call gt_sets (gt, GTTITLE, "Echelle Dispersion Function Fitting") + + switch (xtype) { + case 'p': + call gt_sets (gt, GTXLABEL, "Pixel") + case 'o': + call gt_sets (gt, GTXLABEL, "Order") + case 'w': + call gt_sets (gt, GTXLABEL, "Wavelength") + case 'r': + call gt_sets (gt, GTXLABEL, "Residual") + case 'v': + call gt_sets (gt, GTXLABEL, "Velocity") + } + + switch (ytype) { + case 'p': + call gt_sets (gt, GTYLABEL, "Pixel") + case 'o': + call gt_sets (gt, GTYLABEL, "Order") + case 'w': + call gt_sets (gt, GTYLABEL, "Wavelength") + case 'r': + call gt_sets (gt, GTYLABEL, "Residual") + case 'v': + call gt_sets (gt, GTYLABEL, "Velocity") + } +end diff --git a/noao/onedspec/ecidentify/ecffit/mkpkg b/noao/onedspec/ecidentify/ecffit/mkpkg new file mode 100644 index 00000000..40324cb8 --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/mkpkg @@ -0,0 +1,21 @@ +# Echelle Dispersion Fitting Package + +$checkout libpkg.a ../../ +$update libpkg.a +$checkin libpkg.a ../../ +$exit + +libpkg.a: + ecfcolon.x ecffit.com <error.h> <gset.h> + ecfeval.x ecffit.com + ecffit.x ecffit.com <error.h> <pkg/gtools.h> + ecfgdata.x ecffit.com <pkg/gtools.h> + ecfget.x ecffit.com + ecfgraph.x <gset.h> <mach.h> <pkg/gtools.h> + ecfnearest.x <gset.h> <mach.h> <pkg/gtools.h> + ecfreject.x ecffit.com <mach.h> + ecfrms.x + ecfset.x ecffit.com + ecfsolve.x ecffit.com <mach.h> <math/gsurfit.h> + ecftitle.x ecffit.com <pkg/gtools.h> + ; |