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/twodspec/apextract/apfit.x | 737 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 737 insertions(+) create mode 100644 noao/twodspec/apextract/apfit.x (limited to 'noao/twodspec/apextract/apfit.x') diff --git a/noao/twodspec/apextract/apfit.x b/noao/twodspec/apextract/apfit.x new file mode 100644 index 00000000..67bf149d --- /dev/null +++ b/noao/twodspec/apextract/apfit.x @@ -0,0 +1,737 @@ +include +include +include +include "apertures.h" + + +# AP_FITSPEC -- Fit a spectrum by a smoothing function. + +procedure ap_fitspec (ap, in, spec, ny) + +pointer ap # Aperture (used for labels) +pointer in # Input image (used for labels) +real spec[ny] # spectrum +int ny # Number of points in spectra + +int i, fd, apaxis, clgeti() +real clgetr() +pointer sp, str, x, wts, cv, gp, gt, ic, ic1, gt_init() +bool ap_answer() +data ic1 /NULL/ +errchk icg_fit, ic_fit + +common /apn_com/ ic, gt + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (x, ny, TY_REAL) + call salloc (wts, ny, TY_REAL) + + do i = 1, ny { + Memr[x+i-1] = i + Memr[wts+i-1] = 1 + } + + if (ic == NULL || ic1 == NULL) { + call ic_open (ic) + ic1 = ic + call clgstr ("function", Memc[str], SZ_LINE) + call ic_pstr (ic, "function", Memc[str]) + call ic_puti (ic, "order", clgeti ("order")) + call clgstr ("sample", Memc[str], SZ_LINE) + call ic_pstr (ic, "sample", Memc[str]) + call ic_puti (ic, "naverage", clgeti ("naverage")) + call ic_puti (ic, "niterate", clgeti ("niterate")) + call ic_putr (ic, "low", clgetr ("low_reject")) + call ic_putr (ic, "high", clgetr ("high_reject")) + call ic_putr (ic, "grow", clgetr ("grow")) + call ic_pstr (ic, "ylabel", "") + + gt = gt_init() + } + + call ic_putr (ic, "xmin", 1.) + call ic_putr (ic, "xmax", real (ny)) + apaxis = AP_AXIS(ap) + switch (apaxis) { + case 1: + call ic_pstr (ic, "xlabel", "Line") + case 2: + call ic_pstr (ic, "xlabel", "Column") + } + call gt_sets (gt, GTTYPE, "line") + + # Fit spectrum by a smoothing function. + call sprintf (Memc[str], SZ_LINE, + "%s: %s - Aperture %s") + call pargstr (IM_HDRFILE(in)) + call pargstr (IM_TITLE(in)) + call pargi (AP_ID(ap)) + call gt_sets (gt, GTTITLE, Memc[str]) + + # Query the user to fit the spectrum interactively. + call sprintf (Memc[str], SZ_LINE, + "Fit spectrum for aperture %d for %s interactively?") + call pargi (AP_ID(ap)) + call pargstr (IM_HDRFILE(in)) + if (ap_answer ("ansfitspec1", Memc[str])) { + call ap_gopen (gp) + call icg_fit (ic, gp, "gcur", gt, cv, Memr[x], spec, + Memr[wts], ny) + call amovkr (1., Memr[wts], ny) + } else + call ic_fit (ic, cv, Memr[x], spec, Memr[wts], ny, + YES, YES, YES, YES) + + # Make a graph to the plot log. + call ap_popen (gp, fd, "fitspec") + if (gp != NULL) { + call icg_graphr (ic, gp, gt, cv, Memr[x], spec, Memr[wts], ny) + call ap_pclose (gp, fd) + } + + call cvvector (cv, Memr[x], spec, ny) + call cvfree (cv) +end + + +procedure ap_fitfree () + +pointer ic, gt +common /apn_com/ ic, gt + +begin + call ic_closer (ic) + call gt_free (gt) +end + + +# AP_LNORM -- Normalize the input line apertures by the norm spectra. + +procedure ap_lnorm (ap, out, gain, dbuf, nc, nl, c1, l1, spec, ny, ys, init) + +pointer ap # Aperture structure +pointer out # Output IMIO pointer +real gain # Gain +pointer dbuf # Data buffer +int nc, nl # Size of data buffer +int c1, l1 # Origin of data buffer +real spec[ny] # Normalization spectrum +int ny # Size of profile array +int ys # Start of spectrum in image +int init # Fill between apertures with 1? + +bool clgetb() # Center normalize? +real threshold, clgetr() # Division threshold + +int i, ncols, nlines, ix1, ix2, iy, nsum +real cen, low, high, s, x1, x2, sum, ap_cveval(), asumr() +pointer cv, datain, dataout, imps2r(), impl2r() + +begin + threshold = clgetr ("threshold") + + cen = AP_CEN(ap,1) + low = AP_CEN(ap,1) + AP_LOW(ap,1) + high = AP_CEN(ap,1) + AP_HIGH(ap,1) + cv = AP_CV(ap) + ncols = IM_LEN(out, 1) + nlines = IM_LEN(out, 2) + + # Normalize by the aperture width and apply threshold. + call adivkr (spec, high - low, spec, ny) + if (clgetb ("cennorm")) { + sum = 0. + nsum = 0 + do i = 1, nlines { + iy = i - ys + 1 + if (iy < 1 || iy > ny) + next + s = cen + ap_cveval (cv, real (i)) + ix1 = max (1, int (s)) + ix2 = min (ncols, int (s + 1)) + if (ix1 > ix2) + next + datain = dbuf + (i - l1) * nc + ix1 - c1 + if (ix1 == ix2) + sum = sum + Memr[datain] + else + sum = sum + (ix2-s)*Memr[datain] + (s-ix1)*Memr[datain+1] + nsum = nsum + 1 + } + if (nsum > 0) { + sum = (asumr (spec, ny) / ny) / (sum / nsum / gain) + call adivkr (spec, sum, spec, ny) + } + } + if (!IS_INDEF (threshold)) + call arltr (spec, ny, threshold, threshold) + + do i = 1, nlines { + if (init == YES) { + dataout = impl2r (out, i) + call amovkr (1., Memr[dataout], ncols) + } + + iy = i - ys + 1 + if (iy < 1 || iy > ny) + next + s = ap_cveval (cv, real (i)) + x1 = max (0.5, low + s) + x2 = min (ncols + 0.49, high + s) + if (x1 > x2) + next + + ix1 = nint (x1) + ix2 = nint (x2) + + datain = dbuf + (i - l1) * nc + ix1 - c1 + if (init == YES) + dataout = dataout + ix1 - 1 + else + dataout = imps2r (out, ix1, ix2, i, i) + call adivkr (Memr[datain], spec[iy] * gain, Memr[dataout], + ix2-ix1+1) + } + + call imaddr (out, "CCDMEAN", 1.) +end + + +# AP_CNORM -- Normalize the input column apertures by the norm spectra. + +procedure ap_cnorm (ap, out, gain, dbuf, nc, nl, c1, l1, spec, ny, ys, init) + +pointer ap # Aperture structure +pointer out # Output IMIO pointer +real gain # Gain +pointer dbuf # Data buffer +int nc, nl # Size of data buffer +int c1, l1 # Origin of data buffer +real spec[ny] # Normalization spectrum +int ny # Size of profile array +int ys # Start of spectrum in image +int init # Fill between apertures with 1? + +bool clgetb() # Center normalize? +real threshold, clgetr() # Division threshold + +int ncols, nlines, ix, iy, ix1, ix2, iy1, iy2, nsum +real cen, low, high, s, sum, ap_cveval(), asumr() +pointer sp, y1, y2, cv, datain, dataout, buf, imps2r(), impl2r() + +begin + threshold = clgetr ("threshold") + + call smark (sp) + call salloc (y1, 2 * ny, TY_INT) + y1 = y1 - ys + y2 = y1 + ny + + cen = AP_CEN(ap,2) + low = AP_CEN(ap,2) + AP_LOW(ap,2) + high = AP_CEN(ap,2) + AP_HIGH(ap,2) + cv = AP_CV(ap) + ncols = IM_LEN(out, 1) + nlines = IM_LEN(out, 2) + + # Normalize by the aperture width and apply threshold. + call adivkr (spec, high - low, spec, ny) + if (clgetb ("cennorm")) { + sum = 0. + nsum = 0 + do ix = ys, ys+ny-1 { + s = cen + ap_cveval (cv, real (ix)) + iy1 = max (1, int (s)) + iy2 = min (nlines, int (s + 1)) + if (iy1 > iy2) + next + datain = dbuf + (ix - l1) * nc + iy1 - c1 + if (iy1 == iy2) + sum = sum + Memr[datain] + else + sum = sum + (iy2-s)*Memr[datain] + (s-iy1)*Memr[datain+1] + nsum = nsum + 1 + } + if (nsum > 0) { + sum = (asumr (spec, ny) / ny) / (sum / nsum / gain) + call adivkr (spec, sum, spec, ny) + } + } + if (!IS_INDEF (threshold)) + call arltr (spec, ny, threshold, threshold) + + do ix = ys, ys+ny-1 { + s = ap_cveval (cv, real (ix)) + Memi[y1+ix] = nint (low + s) + Memi[y2+ix] = nint (high + s) + } + call alimi (Memi[y1+ys], 2 * ny, iy1, iy2) + + do iy = 1, nlines { + if (init == YES) { + buf = impl2r (out, iy) + call amovkr (1., Memr[buf], ncols) + } + + if (iy < iy1 || iy > iy2) + next + + for (ix1=ys; ix1<=ys+ny-1; ix1=ix1+1) { + if (iy < Memi[y1+ix1] || iy > Memi[y2+ix1]) + next + for (ix2=ix1+1; ix2<=ys+ny-1; ix2=ix2+1) + if (iy < Memi[y1+ix2] || iy > Memi[y2+ix2]) + break + ix2 = ix2 - 1 + + datain = dbuf + (ix1 - l1) * nc + iy - c1 + if (init == YES) + dataout = buf + ix1 - 1 + else + dataout = imps2r (out, ix1, ix2, iy, iy) + do ix = ix1, ix2 { + Memr[dataout] = Memr[datain] / spec[ix-ys+1] / gain + datain = datain + nc + dataout = dataout + 1 + } + ix1 = ix2 + } + } + + call imaddr (out, "CCDMEAN", 1.) + + call sfree (sp) +end + + +# AP_LFLAT -- Flatten the input line apertures by the norm spectra. + +procedure ap_lflat (ap, out, dbuf, nc, nl, c1, l1, spec, sbuf, profile, nx, ny, + xs, ys, init) + +pointer ap # Aperture structure +pointer out # Output IMIO pointer +pointer dbuf # Data buffer +int nc, nl # Size of data buffer +int c1, l1 # Origin of data buffer +real spec[ny] # Normalization spectrum +pointer sbuf # Sky buffer +real profile[ny,nx] # Profile +int nx, ny # Size of profile array +int xs[ny], ys # Start of spectrum in image +int init # Fill between apertures with 1? + +real threshold, clgetr() # Division threshold + +int i, ncols, nlines, ix, iy, ix1, ix2 +real low, high, s, x1, x2, ap_cveval() +pointer cv, datain, dataout, sky, imps2r(), impl2r() + +begin + threshold = clgetr ("threshold") + if (IS_INDEF(threshold)) + threshold = 0. + threshold = max (0., threshold) + + low = AP_CEN(ap,1) + AP_LOW(ap,1) + high = AP_CEN(ap,1) + AP_HIGH(ap,1) + cv = AP_CV(ap) + ncols = IM_LEN(out, 1) + nlines = IM_LEN(out, 2) + + do i = 1, nlines { + if (init == YES) { + dataout = impl2r (out, i) + call amovkr (1., Memr[dataout], ncols) + } + + iy = i - ys + 1 + if (iy < 1 || iy > ny) + next + s = ap_cveval (cv, real (i)) + x1 = max (0.5, low + s) + x2 = min (ncols + 0.49, high + s) + if (x1 > x2) + next + + ix1 = nint (x1) + ix2 = nint (x2) + + datain = dbuf + (i - l1) * nc + ix1 - c1 + if (init == YES) + dataout = dataout + ix1 - 1 + else + dataout = imps2r (out, ix1, ix2, i, i) + if (sbuf != NULL) + sky = sbuf + (iy - 1) * nx - xs[iy] + do ix = ix1, ix2 { + s = spec[iy] * profile[iy, ix-xs[iy]+1] + if (sbuf != NULL) + s = s + Memr[sky+ix] + if (s > threshold) + Memr[dataout] = Memr[datain] / s + else + Memr[dataout] = 1. + datain = datain + 1 + dataout = dataout + 1 + } + } + + call imaddr (out, "CCDMEAN", 1.) +end + + +# AP_CFLAT -- Flatten the input column apertures by the norm spectra. + +procedure ap_cflat (ap, out, dbuf, nc, nl, c1, l1, spec, sbuf, profile, nx, ny, + xs, ys, init) + +pointer ap # Aperture structure +pointer out # Output IMIO pointer +pointer dbuf # Data buffer +int nc, nl # Size of data buffer +int c1, l1 # Origin of data buffer +real spec[ny] # Normalization spectrum +pointer sbuf # Sky buffer +real profile[ny,nx] # Profile +int nx, ny # Size of profile array +int xs[ny], ys # Start of spectrum in image +int init # Fill between apertures with 1? + +real threshold, clgetr() # Division threshold + +int ncols, nlines, ix, iy, ix1, ix2, iy1, iy2 +real low, high, s, ap_cveval() +pointer sp, y1, y2, cv, datain, dataout, sky, buf, imps2r(), impl2r() + +begin + threshold = clgetr ("threshold") + if (IS_INDEF(threshold)) + threshold = 0. + threshold = max (0., threshold) + + call smark (sp) + call salloc (y1, 2 * ny, TY_INT) + y1 = y1 - ys + y2 = y1 + ny + + low = AP_CEN(ap,2) + AP_LOW(ap,2) + high = AP_CEN(ap,2) + AP_HIGH(ap,2) + cv = AP_CV(ap) + ncols = IM_LEN(out, 1) + nlines = IM_LEN(out, 2) + + do ix = ys, ys+ny-1 { + s = ap_cveval (cv, real (ix)) + Memi[y1+ix] = nint (low + s) + Memi[y2+ix] = nint (high + s) + } + call alimi (Memi[y1+ys], 2 * ny, iy1, iy2) + + do iy = 1, nlines { + if (init == YES) { + buf = impl2r (out, iy) + call amovkr (1., Memr[buf], ncols) + } + + if (iy < iy1 || iy > iy2) + next + + for (ix1=ys; ix1<=ys+ny-1; ix1=ix1+1) { + if (iy < Memi[y1+ix1] || iy > Memi[y2+ix1]) + next + for (ix2=ix1+1; ix2<=ys+ny-1; ix2=ix2+1) + if (iy < Memi[y1+ix2] || iy > Memi[y2+ix2]) + break + ix2 = ix2 - 1 + + datain = dbuf + (ix1 - l1) * nc + iy - c1 + if (init == YES) + dataout = buf + ix1 - 1 + else + dataout = imps2r (out, ix1, ix2, iy, iy) + if (sbuf != NULL) + sky = sbuf - ys * nx + iy - xs[iy] + do ix = ix1, ix2 { + s = spec[ix-ys+1] * profile[ix-ys+1, iy-xs[ix-ys+1]+1] + if (sbuf != NULL) + s = s + Memr[sky+ix*nx] + if (s > threshold) + Memr[dataout] = Memr[datain] / s + else + Memr[dataout] = 1. + datain = datain + nc + dataout = dataout + 1 + } + ix1 = ix2 + } + } + + call imaddr (out, "CCDMEAN", 1.) + + call sfree (sp) +end + + +# AP_LDIFF -- Model residuals. + +procedure ap_ldiff (ap, out, gain, dbuf, nc, nl, c1, l1, spec, profile, nx, ny, + xs, ys, init) + +pointer ap # Aperture structure +pointer out # Output IMIO pointer +real gain # Gain +pointer dbuf # Data buffer +int nc, nl # Size of data buffer +int c1, l1 # Origin of data buffer +real spec[ny] # Normalization spectrum +real profile[ny,nx] # Profile +int nx, ny # Size of profile array +int xs[ny], ys # Start of spectrum in image +int init # Fill between apertures with 1? + +int i, ncols, nlines, ix, iy, ix1, ix2 +real low, high, s, x1, x2, ap_cveval() +pointer cv, datain, dataout, imps2r(), impl2r() + +begin + low = AP_CEN(ap,1) + AP_LOW(ap,1) + high = AP_CEN(ap,1) + AP_HIGH(ap,1) + cv = AP_CV(ap) + ncols = IM_LEN(out, 1) + nlines = IM_LEN(out, 2) + + do i = 1, nlines { + if (init == YES) { + dataout = impl2r (out, i) + call aclrr (Memr[dataout], ncols) + } + + iy = i - ys + 1 + if (iy < 1 || iy > ny) + next + s = ap_cveval (cv, real (i)) + x1 = max (0.5, low + s) + x2 = min (ncols + 0.49, high + s) + if (x1 > x2) + next + + ix1 = nint (x1) + ix2 = nint (x2) + + datain = dbuf + (i - l1) * nc + ix1 - c1 + if (init == YES) + dataout = dataout + ix1 - 1 + else + dataout = imps2r (out, ix1, ix2, i, i) + do ix = ix1, ix2 { + s = spec[iy] * profile[iy, ix-xs[iy]+1] + Memr[dataout] = (Memr[datain] - s) / gain + datain = datain + 1 + dataout = dataout + 1 + } + } +end + + +# AP_CDIFF -- Model residuals + +procedure ap_cdiff (ap, out, gain, dbuf, nc, nl, c1, l1, spec, profile, nx, ny, + xs, ys, init) + +pointer ap # Aperture structure +pointer out # Output IMIO pointer +real gain # Gain +pointer dbuf # Data buffer +int nc, nl # Size of data buffer +int c1, l1 # Origin of data buffer +real spec[ny] # Normalization spectrum +real profile[ny,nx] # Profile +int nx, ny # Size of profile array +int xs[ny], ys # Start of spectrum in image +int init # Fill between apertures with 1? + +int ncols, nlines, ix, iy, ix1, ix2, iy1, iy2 +real low, high, s, ap_cveval() +pointer sp, y1, y2, cv, datain, dataout, buf, imps2r(), impl2r() + +begin + call smark (sp) + call salloc (y1, 2 * ny, TY_INT) + y1 = y1 - ys + y2 = y1 + ny + + low = AP_CEN(ap,2) + AP_LOW(ap,2) + high = AP_CEN(ap,2) + AP_HIGH(ap,2) + cv = AP_CV(ap) + ncols = IM_LEN(out, 1) + nlines = IM_LEN(out, 2) + + do ix = ys, ys+ny-1 { + s = ap_cveval (cv, real (ix)) + Memi[y1+ix] = nint (low + s) + Memi[y2+ix] = nint (high + s) + } + call alimi (Memi[y1+ys], 2 * ny, iy1, iy2) + + do iy = 1, nlines { + if (init == YES) { + buf = impl2r (out, iy) + call aclrr (Memr[buf], ncols) + } + + if (iy < iy1 || iy > iy2) + next + + for (ix1=ys; ix1<=ys+ny-1; ix1=ix1+1) { + if (iy < Memi[y1+ix1] || iy > Memi[y2+ix1]) + next + for (ix2=ix1+1; ix2<=ys+ny-1; ix2=ix2+1) + if (iy < Memi[y1+ix2] || iy > Memi[y2+ix2]) + break + ix2 = ix2 - 1 + + datain = dbuf + (ix1 - l1) * nc + iy - c1 + if (init == YES) + dataout = buf + ix1 - 1 + else + dataout = imps2r (out, ix1, ix2, iy, iy) + do ix = ix1, ix2 { + s = spec[ix-ys+1] * profile[ix-ys+1, iy-xs[ix-ys+1]+1] + Memr[dataout] = (Memr[datain] - s) / gain + datain = datain + nc + dataout = dataout + 1 + } + ix1 = ix2 + } + } + + call sfree (sp) +end + + +# AP_LFIT -- Model fit + +procedure ap_lfit (ap, out, gain, spec, profile, nx, ny, xs, ys, init) + +pointer ap # Aperture structure +pointer out # Output IMIO pointer +real gain # Gain +real spec[ny] # Normalization spectrum +real profile[ny,nx] # Profile +int nx, ny # Size of profile array +int xs[ny], ys # Start of spectrum in image +int init # Fill between apertures with 1? + +int i, ncols, nlines, ix, iy, ix1, ix2 +real low, high, s, x1, x2, ap_cveval() +pointer cv, dataout, imps2r(), impl2r() + +begin + low = AP_CEN(ap,1) + AP_LOW(ap,1) + high = AP_CEN(ap,1) + AP_HIGH(ap,1) + cv = AP_CV(ap) + ncols = IM_LEN(out, 1) + nlines = IM_LEN(out, 2) + + do i = 1, nlines { + if (init == YES) { + dataout = impl2r (out, i) + call aclrr (Memr[dataout], ncols) + } + + iy = i - ys + 1 + if (iy < 1 || iy > ny) + next + s = ap_cveval (cv, real (i)) + x1 = max (0.5, low + s) + x2 = min (ncols + 0.49, high + s) + if (x1 > x2) + next + + ix1 = nint (x1) + ix2 = nint (x2) + + if (init == YES) + dataout = dataout + ix1 - 1 + else + dataout = imps2r (out, ix1, ix2, i, i) + do ix = ix1, ix2 { + s = spec[iy] * profile[iy, ix-xs[iy]+1] + Memr[dataout] = s / gain + dataout = dataout + 1 + } + } +end + + +# AP_CFIT -- Model fit + +procedure ap_cfit (ap, out, gain, spec, profile, nx, ny, xs, ys, init) + +pointer ap # Aperture structure +pointer out # Output IMIO pointer +real gain # Gain +real spec[ny] # Normalization spectrum +real profile[ny,nx] # Profile +int nx, ny # Size of profile array +int xs[ny], ys # Start of spectrum in image +int init # Fill between apertures with 1? + +int ncols, nlines, ix, iy, ix1, ix2, iy1, iy2 +real low, high, s, ap_cveval() +pointer sp, y1, y2, cv, dataout, buf, imps2r(), impl2r() + +begin + call smark (sp) + call salloc (y1, 2 * ny, TY_INT) + y1 = y1 - ys + y2 = y1 + ny + + low = AP_CEN(ap,2) + AP_LOW(ap,2) + high = AP_CEN(ap,2) + AP_HIGH(ap,2) + cv = AP_CV(ap) + ncols = IM_LEN(out, 1) + nlines = IM_LEN(out, 2) + + do ix = ys, ys+ny-1 { + s = ap_cveval (cv, real (ix)) + Memi[y1+ix] = nint (low + s) + Memi[y2+ix] = nint (high + s) + } + call alimi (Memi[y1+ys], 2 * ny, iy1, iy2) + + do iy = 1, nlines { + if (init == YES) { + buf = impl2r (out, iy) + call aclrr (Memr[buf], ncols) + } + + if (iy < iy1 || iy > iy2) + next + + for (ix1=ys; ix1<=ys+ny-1; ix1=ix1+1) { + if (iy < Memi[y1+ix1] || iy > Memi[y2+ix1]) + next + for (ix2=ix1+1; ix2<=ys+ny-1; ix2=ix2+1) + if (iy < Memi[y1+ix2] || iy > Memi[y2+ix2]) + break + ix2 = ix2 - 1 + + if (init == YES) + dataout = buf + ix1 - 1 + else + dataout = imps2r (out, ix1, ix2, iy, iy) + do ix = ix1, ix2 { + s = spec[ix-ys+1] * profile[ix-ys+1, iy-xs[ix-ys+1]+1] + Memr[dataout] = s / gain + dataout = dataout + 1 + } + ix1 = ix2 + } + } + + call sfree (sp) +end -- cgit