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/scombine/README | 17 + noao/onedspec/scombine/generic/icaclip.x | 555 +++++++++++++++++++++++++ noao/onedspec/scombine/generic/icaverage.x | 84 ++++ noao/onedspec/scombine/generic/iccclip.x | 453 +++++++++++++++++++++ noao/onedspec/scombine/generic/icgrow.x | 76 ++++ noao/onedspec/scombine/generic/icmedian.x | 139 +++++++ noao/onedspec/scombine/generic/icmm.x | 152 +++++++ noao/onedspec/scombine/generic/icpclip.x | 224 ++++++++++ noao/onedspec/scombine/generic/icsclip.x | 486 ++++++++++++++++++++++ noao/onedspec/scombine/generic/icsort.x | 275 +++++++++++++ noao/onedspec/scombine/generic/mkpkg | 16 + noao/onedspec/scombine/icgdata.x | 199 +++++++++ noao/onedspec/scombine/iclog.x | 301 ++++++++++++++ noao/onedspec/scombine/icombine.com | 36 ++ noao/onedspec/scombine/icombine.h | 74 ++++ noao/onedspec/scombine/icombine.x | 174 ++++++++ noao/onedspec/scombine/icscale.x | 463 +++++++++++++++++++++ noao/onedspec/scombine/icstat.x | 160 ++++++++ noao/onedspec/scombine/icsum.x | 48 +++ noao/onedspec/scombine/iscombine.key | 23 ++ noao/onedspec/scombine/iscombine.par | 18 + noao/onedspec/scombine/mkpkg | 35 ++ noao/onedspec/scombine/scombine.par | 37 ++ noao/onedspec/scombine/t_scombine.x | 630 +++++++++++++++++++++++++++++ noao/onedspec/scombine/x_scombine.x | 1 + 25 files changed, 4676 insertions(+) create mode 100644 noao/onedspec/scombine/README create mode 100644 noao/onedspec/scombine/generic/icaclip.x create mode 100644 noao/onedspec/scombine/generic/icaverage.x create mode 100644 noao/onedspec/scombine/generic/iccclip.x create mode 100644 noao/onedspec/scombine/generic/icgrow.x create mode 100644 noao/onedspec/scombine/generic/icmedian.x create mode 100644 noao/onedspec/scombine/generic/icmm.x create mode 100644 noao/onedspec/scombine/generic/icpclip.x create mode 100644 noao/onedspec/scombine/generic/icsclip.x create mode 100644 noao/onedspec/scombine/generic/icsort.x create mode 100644 noao/onedspec/scombine/generic/mkpkg create mode 100644 noao/onedspec/scombine/icgdata.x create mode 100644 noao/onedspec/scombine/iclog.x create mode 100644 noao/onedspec/scombine/icombine.com create mode 100644 noao/onedspec/scombine/icombine.h create mode 100644 noao/onedspec/scombine/icombine.x create mode 100644 noao/onedspec/scombine/icscale.x create mode 100644 noao/onedspec/scombine/icstat.x create mode 100644 noao/onedspec/scombine/icsum.x create mode 100644 noao/onedspec/scombine/iscombine.key create mode 100644 noao/onedspec/scombine/iscombine.par create mode 100644 noao/onedspec/scombine/mkpkg create mode 100644 noao/onedspec/scombine/scombine.par create mode 100644 noao/onedspec/scombine/t_scombine.x create mode 100644 noao/onedspec/scombine/x_scombine.x (limited to 'noao/onedspec/scombine') diff --git a/noao/onedspec/scombine/README b/noao/onedspec/scombine/README new file mode 100644 index 00000000..e4031e6a --- /dev/null +++ b/noao/onedspec/scombine/README @@ -0,0 +1,17 @@ +SCOMBINE -- Combine spectra + +This routine is based in large part on IMCOMBINE. The routines in the +generic directory are identical to those in that task except that they +only contain routines for real data. The ic routines in this directory +are similar though modified from IMCOMBINE. + +The iscombine files are for an interactive combine task based on work +by CTIO. Because it is limited currently to linear spectra and is not +organized to take advantage of the IMCOMBINE options it is not installed. +A version of this may someday be added based on the current software. + + +======= + +This version was renamed to OSCOMBINE. It is obsolete and may be removed +at some future time. (4/14/04, Valdes) diff --git a/noao/onedspec/scombine/generic/icaclip.x b/noao/onedspec/scombine/generic/icaclip.x new file mode 100644 index 00000000..41432dd7 --- /dev/null +++ b/noao/onedspec/scombine/generic/icaclip.x @@ -0,0 +1,555 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "../icombine.h" + +define MINCLIP 3 # Minimum number of images for this algorithm + + +# IC_AAVSIGCLIP -- Reject pixels using an average sigma about the average +# The average sigma is normalized by the expected poisson sigma. + +procedure ic_aavsigclipr (d, m, n, scales, zeros, nimages, npts, average) + +pointer d[nimages] # Data pointers +pointer m[nimages] # Image id pointers +int n[npts] # Number of good pixels +real scales[nimages] # Scales +real zeros[nimages] # Zeros +int nimages # Number of images +int npts # Number of output points per line +real average[npts] # Average + +int i, j, k, l, jj, n1, n2, nin, nk, maxkeep +real d1, low, high, sum, a, s, s1, r, one +data one /1.0/ +pointer sp, sums, resid, dp1, dp2, mp1, mp2 + +include "../icombine.com" + +begin + # If there are insufficient pixels go on to the combining. + if (nkeep < 0) + maxkeep = max (0, nimages + nkeep) + else + maxkeep = min (nimages, nkeep) + if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) { + docombine = true + return + } + + call smark (sp) + call salloc (sums, npts, TY_REAL) + call salloc (resid, nimages+1, TY_REAL) + + # Since the unweighted average is computed here possibly skip combining + if (dowts || combine != AVERAGE) + docombine = true + else + docombine = false + + # Compute the unweighted average with the high and low rejected and + # the poisson scaled average sigma. There must be at least three + # pixels at each point to define the average and contributions to + # the mean sigma. Corrections for differences in the image + # scale factors are selected by the doscale1 flag. + + nin = n[1] + s = 0. + n2 = 0 + do i = 1, npts { + k = i - 1 + n1 = n[i] + if (n1 < 3) + next + + # Unweighted average with the high and low rejected + low = Memr[d[1]+k] + high = Memr[d[2]+k] + if (low > high) { + d1 = low + low = high + high = d1 + } + sum = 0. + do j = 3, n1 { + d1 = Memr[d[j]+k] + if (d1 < low) { + sum = sum + low + low = d1 + } else if (d1 > high) { + sum = sum + high + high = d1 + } else + sum = sum + d1 + } + a = sum / (n1 - 2) + sum = sum + low + high + + # Poisson scaled sigma accumulation + if (doscale1) { + do j = 1, n1 { + dp1 = d[j] + k + mp1 = m[j] + k + + d1 = Memr[dp1] + l = Memi[mp1] + s1 = max (one, (a + zeros[l]) / scales[l]) + s = s + (d1 - a) ** 2 / s1 + } + } else { + s1 = max (one, a) + do j = 1, n1 + s = s + (Memr[d[j]+k] - a) ** 2 / s1 + } + n2 = n2 + n1 + + # Save the average and sum for later. + average[i] = a + Memr[sums+k] = sum + } + + # Here is the final sigma. + if (n2 > 1) + s = sqrt (s / (n2 - 1)) + + # Reject pixels and compute the final average (if needed). + # There must be at least three pixels at each point for rejection. + # Iteratively scale the mean sigma and reject pixels + # Compact the data and keep track of the image IDs if needed. + + do i = 1, npts { + k = i - 1 + n1 = n[i] + if (nkeep < 0) + maxkeep = max (0, n1 + nkeep) + else + maxkeep = min (n1, nkeep) + if (n1 <= max (2, maxkeep)) { + if (!docombine) { + if (n1 == 0) + average[i] = blank + else { + sum = Memr[d[1]+k] + do j = 2, n1 + sum = sum + Memr[d[j]+k] + average[i] = sum / n1 + } + } + next + } + + a = average[i] + sum = Memr[sums+k] + + repeat { + n2 = n1 + if (s > 0.) { + if (doscale1) { + for (j=1; j<=n1; j=j+1) { + dp1 = d[j] + k + mp1 = m[j] + k + + d1 = Memr[dp1] + l = Memi[mp1] + s1 = s * sqrt (max (one, (a+zeros[l]) / scales[l])) + r = (d1 - a) / s1 + if (r < -lsigma || r > hsigma) { + Memr[resid+n1] = abs(r) + if (j < n1) { + dp2 = d[n1] + k + Memr[dp1] = Memr[dp2] + Memr[dp2] = d1 + mp2 = m[n1] + k + Memi[mp1] = Memi[mp2] + Memi[mp2] = l + j = j - 1 + } + sum = sum - d1 + n1 = n1 - 1 + } + } + } else { + s1 = s * sqrt (max (one, a)) + for (j=1; j<=n1; j=j+1) { + dp1 = d[j] + k + d1 = Memr[dp1] + r = (d1 - a) / s1 + if (r < -lsigma || r > hsigma) { + Memr[resid+n1] = abs(r) + if (j < n1) { + dp2 = d[n1] + k + Memr[dp1] = Memr[dp2] + Memr[dp2] = d1 + if (keepids) { + mp1 = m[j] + k + mp2 = m[n1] + k + l = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = l + } + j = j - 1 + } + sum = sum - d1 + n1 = n1 - 1 + } + } + } + } + if (n1 > 1) + a = sum / n1 + } until (n1 == n2 || n1 <= max (2, maxkeep)) + + # If too many are rejected add some back in. + # Pixels with equal residuals are added together. + if (n1 < maxkeep) { + nk = maxkeep + if (doscale1) { + for (j=n1+1; j<=nk; j=j+1) { + dp1 = d[j] + k + mp1 = m[j] + k + r = Memr[resid+j] + jj = 0 + do l = j+1, n2 { + s = Memr[resid+l] + if (s < r + TOL) { + if (s > r - TOL) + jj = jj + 1 + else { + jj = 0 + Memr[resid+l] = r + r = s + dp2 = d[l] + k + d1 = Memr[dp1] + Memr[dp1] = Memr[dp2] + Memr[dp2] = d1 + mp2 = m[l] + k + s = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = s + } + } + } + sum = sum + Memr[dp1] + n1 = n1 + 1 + nk = max (nk, j+jj) + } + } else { + for (j=n1+1; j<=nk; j=j+1) { + dp1 = d[j] + k + r = Memr[resid+j] + jj = 0 + do l = j+1, n2 { + s = Memr[resid+l] + if (s < r + TOL) { + if (s > r - TOL) + jj = jj + 1 + else { + jj = 0 + Memr[resid+l] = r + r = s + dp2 = d[l] + k + d1 = Memr[dp1] + Memr[dp1] = Memr[dp2] + Memr[dp2] = d1 + if (keepids) { + mp1 = m[j] + k + mp2 = m[l] + k + s = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = s + } + } + } + } + sum = sum + Memr[dp1] + n1 = n1 + 1 + nk = max (nk, j+jj) + } + } + if (n1 > 1) + a = sum / n1 + } + + # Save the average if needed. + n[i] = n1 + if (!docombine) { + if (n1 > 0) + average[i] = a + else + average[i] = blank + } + } + + # Check if the data flag has to be reset for rejected pixels + if (dflag == D_ALL) { + do i = 1, npts { + if (n[i] != nin) { + dflag = D_MIX + break + } + } + } + + call sfree (sp) +end + + +# IC_MAVSIGCLIP -- Reject pixels using an average sigma about the median +# The average sigma is normalized by the expected poisson sigma. + +procedure ic_mavsigclipr (d, m, n, scales, zeros, nimages, npts, median) + +pointer d[nimages] # Data pointers +pointer m[nimages] # Image id pointers +int n[npts] # Number of good pixels +real scales[nimages] # Scales +real zeros[nimages] # Zeros +int nimages # Number of images +int npts # Number of output points per line +real median[npts] # Median + +int i, j, k, l, id, n1, n2, n3, nl, nh, nin, maxkeep +pointer sp, resid, mp1, mp2 +real med, low, high, r, s, s1, one +data one /1.0/ + +include "../icombine.com" + +begin + # If there are insufficient pixels go on to the combining. + if (nkeep < 0) + maxkeep = max (0, nimages + nkeep) + else + maxkeep = min (nimages, nkeep) + if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) { + docombine = true + return + } + + call smark (sp) + call salloc (resid, nimages+1, TY_REAL) + + # Compute the poisson scaled average sigma about the median. + # There must be at least three pixels at each point to define + # the mean sigma. Corrections for differences in the image + # scale factors are selected by the doscale1 flag. + + s = 0. + n2 = 0 + nin = n[1] + do i = 1, npts { + k = i - 1 + n1 = n[i] + if (n1 < 3) { + if (n1 == 0) + median[i] = blank + else if (n1 == 1) + median[i] = Memr[d[1]+k] + else { + low = Memr[d[1]+k] + high = Memr[d[2]+k] + median[i] = (low + high) / 2. + } + next + } + + # Median + n3 = 1 + n1 / 2 + if (mod (n1, 2) == 0) { + low = Memr[d[n3-1]+k] + high = Memr[d[n3]+k] + med = (low + high) / 2. + } else + med = Memr[d[n3]+k] + + # Poisson scaled sigma accumulation + if (doscale1) { + do j = 1, n1 { + l = Memi[m[j]+k] + s1 = max (one, (med + zeros[l]) / scales[l]) + s = s + (Memr[d[j]+k] - med) ** 2 / s1 + } + } else { + s1 = max (one, med) + do j = 1, n1 + s = s + (Memr[d[j]+k] - med) ** 2 / s1 + } + n2 = n2 + n1 + + # Save the median for later. + median[i] = med + } + + # Here is the final sigma. + if (n2 > 1) + s = sqrt (s / (n2 - 1)) + else + return + + # Compute individual sigmas and iteratively clip. + do i = 1, npts { + k = i - 1 + n1 = n[i] + if (nkeep < 0) + maxkeep = max (0, n1 + nkeep) + else + maxkeep = min (n1, nkeep) + if (n1 < max (3, maxkeep+1)) + next + nl = 1 + nh = n1 + med = median[i] + + repeat { + n2 = n1 + n3 = nl + n1 / 2 + + if (n1 >= max (MINCLIP, maxkeep+1) && s > 0.) { + if (doscale1) { + for (; nl <= n2; nl = nl + 1) { + l = Memi[m[nl]+k] + s1 = s * sqrt (max (one, (med+zeros[l])/scales[l])) + r = (med - Memr[d[nl]+k]) / s1 + if (r <= lsigma) + break + Memr[resid+nl] = r + n1 = n1 - 1 + } + for (; nh >= nl; nh = nh - 1) { + l = Memi[m[nh]+k] + s1 = s * sqrt (max (one, (med+zeros[l])/scales[l])) + r = (Memr[d[nh]+k] - med) / s1 + if (r <= hsigma) + break + Memr[resid+nh] = r + n1 = n1 - 1 + } + } else { + s1 = s * sqrt (max (one, med)) + for (; nl <= n2; nl = nl + 1) { + r = (med - Memr[d[nl]+k]) / s1 + if (r <= lsigma) + break + Memr[resid+nl] = r + n1 = n1 - 1 + } + for (; nh >= nl; nh = nh - 1) { + r = (Memr[d[nh]+k] - med) / s1 + if (r <= hsigma) + break + Memr[resid+nh] = r + n1 = n1 - 1 + } + } + + # Recompute median + if (n1 < n2) { + if (n1 > 0) { + n3 = nl + n1 / 2 + if (mod (n1, 2) == 0) { + low = Memr[d[n3-1]+k] + high = Memr[d[n3]+k] + med = (low + high) / 2. + } else + med = Memr[d[n3]+k] + } else + med = blank + } + } + } until (n1 == n2 || n1 < max (MINCLIP, maxkeep+1)) + + # If too many are rejected add some back in. + # Pixels with equal residuals are added together. + while (n1 < maxkeep) { + if (nl == 1) + nh = nh + 1 + else if (nh == n[i]) + nl = nl - 1 + else { + r = Memr[resid+nl-1] + s = Memr[resid+nh+1] + if (r < s) { + nl = nl - 1 + r = r + TOL + if (s <= r) + nh = nh + 1 + if (nl > 1) { + if (Memr[resid+nl-1] <= r) + nl = nl - 1 + } + } else { + nh = nh + 1 + s = s + TOL + if (r <= s) + nl = nl - 1 + if (nh < n2) { + if (Memr[resid+nh+1] <= s) + nh = nh + 1 + } + } + } + n1 = nh - nl + 1 + + # Recompute median + if (n1 < n2) { + if (n1 > 0) { + n3 = nl + n1 / 2 + if (mod (n1, 2) == 0) { + low = Memr[d[n3-1]+k] + high = Memr[d[n3]+k] + med = (low + high) / 2. + } else + med = Memr[d[n3]+k] + } else + med = blank + } + } + + # Only set median and reorder if needed + n[i] = n1 + if (n1 > 0 && nl > 1 && (combine != MEDIAN || grow > 0)) { + j = max (nl, n1 + 1) + if (keepids) { + do l = 1, min (n1, nl-1) { + Memr[d[l]+k] = Memr[d[j]+k] + if (grow > 0) { + mp1 = m[l] + k + mp2 = m[j] + k + id = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = id + } else + Memi[m[l]+k] = Memi[m[j]+k] + j = j + 1 + } + } else { + do l = 1, min (n1, nl - 1) { + Memr[d[l]+k] = Memr[d[j]+k] + j = j + 1 + } + } + } + + if (combine == MEDIAN) + median[i] = med + } + + # Check if data flag needs to be reset for rejected pixels + if (dflag == D_ALL) { + do i = 1, npts { + if (n[i] != nin) { + dflag = D_MIX + break + } + } + } + + # Flag that the median is computed. + if (combine == MEDIAN) + docombine = false + else + docombine = true + + call sfree (sp) +end + diff --git a/noao/onedspec/scombine/generic/icaverage.x b/noao/onedspec/scombine/generic/icaverage.x new file mode 100644 index 00000000..6c5c870b --- /dev/null +++ b/noao/onedspec/scombine/generic/icaverage.x @@ -0,0 +1,84 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include +include "../icombine.h" + + +# IC_AVERAGE -- Compute the average image line. +# Options include a weight average. + +procedure ic_averager (d, m, n, wts, npts, average) + +pointer d[ARB] # Data pointers +pointer m[ARB] # Image ID pointers +int n[npts] # Number of points +real wts[ARB] # Weights +int npts # Number of output points per line +real average[npts] # Average (returned) + +int i, j, k +real sumwt, wt +real sum + +include "../icombine.com" + +begin + # If no data has been excluded do the average without checking the + # number of points and using the fact that the weights are normalized. + # If all the data has been excluded set the average to the blank value. + + if (dflag == D_ALL) { + if (dowts) { + do i = 1, npts { + k = i - 1 + wt = wts[Memi[m[1]+k]] + sum = Memr[d[1]+k] * wt + do j = 2, n[i] { + wt = wts[Memi[m[j]+k]] + sum = sum + Memr[d[j]+k] * wt + } + average[i] = sum + } + } else { + do i = 1, npts { + k = i - 1 + sum = Memr[d[1]+k] + do j = 2, n[i] + sum = sum + Memr[d[j]+k] + average[i] = sum / n[i] + } + } + } else if (dflag == D_NONE) { + do i = 1, npts + average[i] = blank + } else { + if (dowts) { + do i = 1, npts { + if (n[i] > 0) { + k = i - 1 + wt = wts[Memi[m[1]+k]] + sum = Memr[d[1]+k] * wt + sumwt = wt + do j = 2, n[i] { + wt = wts[Memi[m[j]+k]] + sum = sum + Memr[d[j]+k] * wt + sumwt = sumwt + wt + } + average[i] = sum / sumwt + } else + average[i] = blank + } + } else { + do i = 1, npts { + if (n[i] > 0) { + k = i - 1 + sum = Memr[d[1]+k] + do j = 2, n[i] + sum = sum + Memr[d[j]+k] + average[i] = sum / n[i] + } else + average[i] = blank + } + } + } +end diff --git a/noao/onedspec/scombine/generic/iccclip.x b/noao/onedspec/scombine/generic/iccclip.x new file mode 100644 index 00000000..26b17ba2 --- /dev/null +++ b/noao/onedspec/scombine/generic/iccclip.x @@ -0,0 +1,453 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "../icombine.h" + +define MINCLIP 2 # Mininum number of images for algorithm + + +# IC_ACCDCLIP -- Reject pixels using CCD noise parameters about the average + +procedure ic_accdclipr (d, m, n, scales, zeros, nm, nimages, npts, average) + +pointer d[nimages] # Data pointers +pointer m[nimages] # Image id pointers +int n[npts] # Number of good pixels +real scales[nimages] # Scales +real zeros[nimages] # Zeros +real nm[3,nimages] # Noise model parameters +int nimages # Number of images +int npts # Number of output points per line +real average[npts] # Average + +int i, j, k, l, jj, n1, n2, nin, nk, maxkeep +real d1, low, high, sum, a, s, r, zero +data zero /0.0/ +pointer sp, resid, dp1, dp2, mp1, mp2 + +include "../icombine.com" + +begin + # If there are no pixels go on to the combining. Since the unweighted + # average is computed here possibly skip the combining later. + + # There must be at least max (1, nkeep) pixels. + if (nkeep < 0) + maxkeep = max (0, nimages + nkeep) + else + maxkeep = min (nimages, nkeep) + if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) { + docombine = true + return + } else if (dowts || combine != AVERAGE) + docombine = true + else + docombine = false + + call smark (sp) + call salloc (resid, nimages+1, TY_REAL) + + # There must be at least two pixels for rejection. The initial + # average is the low/high rejected average except in the case of + # just two pixels. The rejections are iterated and the average + # is recomputed. Corrections for scaling may be performed. + # Depending on other flags the image IDs may also need to be adjusted. + + nin = n[1] + do i = 1, npts { + k = i - 1 + n1 = n[i] + if (nkeep < 0) + maxkeep = max (0, n1 + nkeep) + else + maxkeep = min (n1, nkeep) + if (n1 <= max (MINCLIP-1, maxkeep)) { + if (!docombine) { + if (n1 == 0) + average[i] = blank + else { + sum = Memr[d[1]+k] + do j = 2, n1 + sum = sum + Memr[d[j]+k] + average[i] = sum / n1 + } + } + next + } + + repeat { + if (n1 == 2) { + sum = Memr[d[1]+k] + sum = sum + Memr[d[2]+k] + a = sum / 2 + } else { + low = Memr[d[1]+k] + high = Memr[d[2]+k] + if (low > high) { + d1 = low + low = high + high = d1 + } + sum = 0. + do j = 3, n1 { + d1 = Memr[d[j]+k] + if (d1 < low) { + sum = sum + low + low = d1 + } else if (d1 > high) { + sum = sum + high + high = d1 + } else + sum = sum + d1 + } + a = sum / (n1 - 2) + sum = sum + low + high + } + n2 = n1 + if (doscale1) { + for (j=1; j<=n1; j=j+1) { + dp1 = d[j] + k + mp1 = m[j] + k + + l = Memi[mp1] + s = scales[l] + d1 = max (zero, s * (a + zeros[l])) + s = sqrt (nm[1,l] + d1/nm[2,l] + (d1*nm[3,l])**2) / s + + d1 = Memr[dp1] + r = (d1 - a) / s + if (r < -lsigma || r > hsigma) { + Memr[resid+n1] = abs(r) + if (j < n1) { + dp2 = d[n1] + k + Memr[dp1] = Memr[dp2] + Memr[dp2] = d1 + mp2 = m[n1] + k + Memi[mp1] = Memi[mp2] + Memi[mp2] = l + j = j - 1 + } + sum = sum - d1 + n1 = n1 - 1 + } + } + } else { + if (!keepids) { + s = max (zero, a) + s = sqrt (nm[1,1] + s/nm[2,1] + (s*nm[3,1])**2) + } + for (j=1; j<=n1; j=j+1) { + if (keepids) { + l = Memi[m[j]+k] + s = max (zero, a) + s = sqrt (nm[1,l] + s/nm[2,l] + (s*nm[3,l])**2) + } + dp1 = d[j] + k + d1 = Memr[dp1] + r = (d1 - a) / s + if (r < -lsigma || r > hsigma) { + Memr[resid+n1] = abs(r) + if (j < n1) { + dp2 = d[n1] + k + Memr[dp1] = Memr[dp2] + Memr[dp2] = d1 + if (keepids) { + mp1 = m[j] + k + mp2 = m[n1] + k + l = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = l + } + j = j - 1 + } + sum = sum - d1 + n1 = n1 - 1 + } + } + } + } until (n1 == n2 || n1 < max (MINCLIP, maxkeep+1)) + + if (n1 < maxkeep) { + nk = maxkeep + if (doscale1) { + for (j=n1+1; j<=nk; j=j+1) { + dp1 = d[j] + k + mp1 = m[j] + k + r = Memr[resid+j] + jj = 0 + do l = j+1, n2 { + s = Memr[resid+l] + if (s < r + TOL) { + if (s > r - TOL) + jj = jj + 1 + else { + jj = 0 + Memr[resid+l] = r + r = s + dp2 = d[l] + k + d1 = Memr[dp1] + Memr[dp1] = Memr[dp2] + Memr[dp2] = d1 + mp2 = m[l] + k + s = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = s + } + } + } + sum = sum + Memr[dp1] + n1 = n1 + 1 + nk = max (nk, j+jj) + } + } else { + for (j=n1+1; j<=nk; j=j+1) { + dp1 = d[j] + k + r = Memr[resid+j] + jj = 0 + do l = j+1, n2 { + s = Memr[resid+l] + if (s < r + TOL) { + if (s > r - TOL) + jj = jj + 1 + else { + jj = 0 + Memr[resid+l] = r + r = s + dp2 = d[l] + k + d1 = Memr[dp1] + Memr[dp1] = Memr[dp2] + Memr[dp2] = d1 + if (keepids) { + mp1 = m[j] + k + mp2 = m[l] + k + s = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = s + } + } + } + } + sum = sum + Memr[dp1] + n1 = n1 + 1 + nk = max (nk, j+jj) + } + } + } + + n[i] = n1 + if (!docombine) + if (n1 > 0) + average[i] = sum / n1 + else + average[i] = blank + } + + # Check if the data flag has to be reset for rejected pixels + if (dflag == D_ALL) { + do i = 1, npts { + if (n[i] != nin) { + dflag = D_MIX + break + } + } + } + + call sfree (sp) +end + + +# IC_CCDCLIP -- Reject pixels using CCD noise parameters about the median + +procedure ic_mccdclipr (d, m, n, scales, zeros, nm, nimages, npts, median) + +pointer d[nimages] # Data pointers +pointer m[nimages] # Image id pointers +int n[npts] # Number of good pixels +real scales[nimages] # Scales +real zeros[nimages] # Zeros +real nm[3,nimages] # Noise model +int nimages # Number of images +int npts # Number of output points per line +real median[npts] # Median + +int i, j, k, l, id, n1, n2, n3, nl, nh, nin, maxkeep +real r, s +pointer sp, resid, mp1, mp2 +real med, zero +data zero /0.0/ + +include "../icombine.com" + +begin + # There must be at least max (MINCLIP, nkeep+1) pixels. + if (nkeep < 0) + maxkeep = max (0, nimages + nkeep) + else + maxkeep = min (nimages, nkeep) + if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) { + docombine = true + return + } + + call smark (sp) + call salloc (resid, nimages+1, TY_REAL) + + # Compute median and sigma and iteratively clip. + nin = n[1] + do i = 1, npts { + k = i - 1 + n1 = n[i] + if (nkeep < 0) + maxkeep = max (0, n1 + nkeep) + else + maxkeep = min (n1, nkeep) + nl = 1 + nh = n1 + + repeat { + n2 = n1 + n3 = nl + n1 / 2 + + if (n1 == 0) + med = blank + else if (mod (n1, 2) == 0) { + med = Memr[d[n3-1]+k] + med = (med + Memr[d[n3]+k]) / 2. + } else + med = Memr[d[n3]+k] + + if (n1 >= max (MINCLIP, maxkeep+1)) { + if (doscale1) { + for (; nl <= n2; nl = nl + 1) { + l = Memi[m[nl]+k] + s = scales[l] + r = max (zero, s * (med + zeros[l])) + s = sqrt (nm[1,l] + r/nm[2,l] + (r*nm[3,l])**2) / s + r = (med - Memr[d[nl]+k]) / s + if (r <= lsigma) + break + Memr[resid+nl] = r + n1 = n1 - 1 + } + for (; nh >= nl; nh = nh - 1) { + l = Memi[m[nh]+k] + s = scales[l] + r = max (zero, s * (med + zeros[l])) + s = sqrt (nm[1,l] + r/nm[2,l] + (r*nm[3,l])**2) / s + r = (Memr[d[nh]+k] - med) / s + if (r <= hsigma) + break + Memr[resid+nh] = r + n1 = n1 - 1 + } + } else { + if (!keepids) { + s = max (zero, med) + s = sqrt (nm[1,1] + s/nm[2,1] + (s*nm[3,1])**2) + } + for (; nl <= n2; nl = nl + 1) { + if (keepids) { + l = Memi[m[nl]+k] + s = max (zero, med) + s = sqrt (nm[1,l] + s/nm[2,l] + (s*nm[3,l])**2) + } + r = (med - Memr[d[nl]+k]) / s + if (r <= lsigma) + break + Memr[resid+nl] = r + n1 = n1 - 1 + } + for (; nh >= nl; nh = nh - 1) { + if (keepids) { + l = Memi[m[nh]+k] + s = max (zero, med) + s = sqrt (nm[1,l] + s/nm[2,l] + (s*nm[3,l])**2) + } + r = (Memr[d[nh]+k] - med) / s + if (r <= hsigma) + break + Memr[resid+nh] = r + n1 = n1 - 1 + } + } + } + } until (n1 == n2 || n1 < max (MINCLIP, maxkeep+1)) + + while (n1 < maxkeep) { + if (nl == 1) + nh = nh + 1 + else if (nh == n[i]) + nl = nl - 1 + else { + r = Memr[resid+nl-1] + s = Memr[resid+nh+1] + if (r < s) { + nl = nl - 1 + r = r + TOL + if (s <= r) + nh = nh + 1 + if (nl > 1) { + if (Memr[resid+nl-1] <= r) + nl = nl - 1 + } + } else { + nh = nh + 1 + s = s + TOL + if (r <= s) + nl = nl - 1 + if (nh < n2) { + if (Memr[resid+nh+1] <= s) + nh = nh + 1 + } + } + } + n1 = nh - nl + 1 + } + + # Only set median and reorder if needed + n[i] = n1 + if (n1 > 0 && nl > 1 && (combine != MEDIAN || grow > 0)) { + j = max (nl, n1 + 1) + if (keepids) { + do l = 1, min (n1, nl-1) { + Memr[d[l]+k] = Memr[d[j]+k] + if (grow > 0) { + mp1 = m[l] + k + mp2 = m[j] + k + id = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = id + } else + Memi[m[l]+k] = Memi[m[j]+k] + j = j + 1 + } + } else { + do l = 1, min (n1, nl - 1) { + Memr[d[l]+k] = Memr[d[j]+k] + j = j + 1 + } + } + } + + if (combine == MEDIAN) + median[i] = med + } + + # Check if data flag needs to be reset for rejected pixels + if (dflag == D_ALL) { + do i = 1, npts { + if (n[i] != nin) { + dflag = D_MIX + break + } + } + } + + # Flag that the median is computed. + if (combine == MEDIAN) + docombine = false + else + docombine = true + + call sfree (sp) +end + diff --git a/noao/onedspec/scombine/generic/icgrow.x b/noao/onedspec/scombine/generic/icgrow.x new file mode 100644 index 00000000..074bd8c3 --- /dev/null +++ b/noao/onedspec/scombine/generic/icgrow.x @@ -0,0 +1,76 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "../icombine.h" + + +# IC_GROW -- Reject neigbors of rejected pixels. +# The rejected pixels are marked by having nonzero ids beyond the number +# of included pixels. The pixels rejected here are given zero ids +# to avoid growing of the pixels rejected here. The unweighted average +# can be updated but any rejected pixels requires the median to be +# recomputed. When the number of pixels at a grow point reaches nkeep +# no further pixels are rejected. Note that the rejection order is not +# based on the magnitude of the residuals and so a grow from a weakly +# rejected image pixel may take precedence over a grow from a strongly +# rejected image pixel. + +procedure ic_growr (d, m, n, nimages, npts, average) + +pointer d[ARB] # Data pointers +pointer m[ARB] # Image id pointers +int n[npts] # Number of good pixels +int nimages # Number of images +int npts # Number of output points per line +real average[npts] # Average + +int i1, i2, j1, j2, k1, k2, l, is, ie, n2, maxkeep +pointer mp1, mp2 + +include "../icombine.com" + +begin + if (dflag == D_NONE) + return + + do i1 = 1, npts { + k1 = i1 - 1 + is = max (1, i1 - grow) + ie = min (npts, i1 + grow) + do j1 = n[i1]+1, nimages { + l = Memi[m[j1]+k1] + if (l == 0) + next + if (combine == MEDIAN) + docombine = true + + do i2 = is, ie { + if (i2 == i1) + next + k2 = i2 - 1 + n2 = n[i2] + if (nkeep < 0) + maxkeep = max (0, n2 + nkeep) + else + maxkeep = min (n2, nkeep) + if (n2 <= maxkeep) + next + do j2 = 1, n2 { + mp1 = m[j2] + k2 + if (Memi[mp1] == l) { + if (!docombine && n2 > 1) + average[i2] = + (n2*average[i2] - Memr[d[j2]+k2]) / (n2-1) + mp2 = m[n2] + k2 + if (j2 < n2) { + Memr[d[j2]+k2] = Memr[d[n2]+k2] + Memi[mp1] = Memi[mp2] + } + Memi[mp2] = 0 + n[i2] = n2 - 1 + break + } + } + } + } + } +end diff --git a/noao/onedspec/scombine/generic/icmedian.x b/noao/onedspec/scombine/generic/icmedian.x new file mode 100644 index 00000000..e7607340 --- /dev/null +++ b/noao/onedspec/scombine/generic/icmedian.x @@ -0,0 +1,139 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "../icombine.h" + + +# IC_MEDIAN -- Median of lines + +procedure ic_medianr (d, n, npts, median) + +pointer d[ARB] # Input data line pointers +int n[npts] # Number of good pixels +int npts # Number of output points per line +real median[npts] # Median + +int i, j1, j2, j3, k, n1 +bool even +real val1, val2, val3 + +include "../icombine.com" + +begin + if (dflag == D_NONE) { + do i = 1, npts + median[i]= blank + return + } + + # Check for previous sorting + if (mclip) { + if (dflag == D_ALL) { + n1 = n[1] + even = (mod (n1, 2) == 0) + j1 = n1 / 2 + 1 + j2 = n1 / 2 + do i = 1, npts { + k = i - 1 + if (even) { + val1 = Memr[d[j1]+k] + val2 = Memr[d[j2]+k] + median[i] = (val1 + val2) / 2. + } else + median[i] = Memr[d[j1]+k] + } + } else { + do i = 1, npts { + k = i - 1 + n1 = n[i] + if (n1 > 0) { + j1 = n1 / 2 + 1 + if (mod (n1, 2) == 0) { + j2 = n1 / 2 + val1 = Memr[d[j1]+k] + val2 = Memr[d[j2]+k] + median[i] = (val1 + val2) / 2. + } else + median[i] = Memr[d[j1]+k] + } else + median[i] = blank + } + } + return + } + + # Repeatedly exchange the extreme values until there are three + # or fewer pixels. + + do i = 1, npts { + k = i - 1 + n1 = n[i] + while (n1 > 3) { + j1 = 1 + j2 = 1 + val1 = Memr[d[j1]+k] + val2 = val1 + do j3 = 2, n1 { + val3 = Memr[d[j3]+k] + if (val3 > val1) { + j1 = j3 + val1 = val3 + } else if (val3 < val2) { + j2 = j3 + val2 = val3 + } + } + j3 = n1 - 1 + if (j1 < j3 && j2 < j3) { + Memr[d[j1]+k] = val3 + Memr[d[j2]+k] = Memr[d[j3]+k] + Memr[d[j3]+k] = val1 + Memr[d[n1]+k] = val2 + } else if (j1 < j3) { + if (j2 == j3) { + Memr[d[j1]+k] = val3 + Memr[d[n1]+k] = val1 + } else { + Memr[d[j1]+k] = Memr[d[j3]+k] + Memr[d[j3]+k] = val1 + } + } else if (j2 < j3) { + if (j1 == j3) { + Memr[d[j2]+k] = val3 + Memr[d[n1]+k] = val2 + } else { + Memr[d[j2]+k] = Memr[d[j3]+k] + Memr[d[j3]+k] = val2 + } + } + n1 = n1 - 2 + } + + if (n1 == 3) { + val1 = Memr[d[1]+k] + val2 = Memr[d[2]+k] + val3 = Memr[d[3]+k] + if (val1 < val2) { + if (val2 < val3) # abc + median[i] = val2 + else if (val1 < val3) # acb + median[i] = val3 + else # cab + median[i] = val1 + } else { + if (val2 > val3) # cba + median[i] = val2 + else if (val1 > val3) # bca + median[i] = val3 + else # bac + median[i] = val1 + } + } else if (n1 == 2) { + val1 = Memr[d[1]+k] + val2 = Memr[d[2]+k] + median[i] = (val1 + val2) / 2 + } else if (n1 == 1) + median[i] = Memr[d[1]+k] + else + median[i] = blank + } +end diff --git a/noao/onedspec/scombine/generic/icmm.x b/noao/onedspec/scombine/generic/icmm.x new file mode 100644 index 00000000..1c314241 --- /dev/null +++ b/noao/onedspec/scombine/generic/icmm.x @@ -0,0 +1,152 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "../icombine.h" + + +# IC_MM -- Reject a specified number of high and low pixels + +procedure ic_mmr (d, m, n, npts) + +pointer d[ARB] # Data pointers +pointer m[ARB] # Image ID pointers +int n[npts] # Number of good pixels +int npts # Number of output points per line + +int n1, ncombine, npairs, nlow, nhigh, np +int i, i1, j, jmax, jmin +pointer k, kmax, kmin +real d1, d2, dmin, dmax + +include "../icombine.com" + +begin + if (dflag == D_NONE) + return + + if (dflag == D_ALL) { + n1 = n[1] + nlow = flow * n1 + 0.001 + nhigh = fhigh * n1 + 0.001 + ncombine = n1 - nlow - nhigh + npairs = min (nlow, nhigh) + nlow = nlow - npairs + nhigh = nhigh - npairs + } + + do i = 1, npts { + i1 = i - 1 + n1 = n[i] + if (dflag == D_MIX) { + nlow = flow * n1 + 0.001 + nhigh = fhigh * n1 + 0.001 + ncombine = max (ncombine, n1 - nlow - nhigh) + npairs = min (nlow, nhigh) + nlow = nlow - npairs + nhigh = nhigh - npairs + } + + # Reject the npairs low and high points. + do np = 1, npairs { + k = d[1] + i1 + d1 = Memr[k] + dmax = d1; dmin = d1; jmax = 1; jmin = 1; kmax = k; kmin = k + do j = 2, n1 { + d2 = d1 + k = d[j] + i1 + d1 = Memr[k] + if (d1 > dmax) { + dmax = d1; jmax = j; kmax = k + } else if (d1 < dmin) { + dmin = d1; jmin = j; kmin = k + } + } + j = n1 - 1 + if (keepids) { + if (jmax < j) { + if (jmin != j) { + Memr[kmax] = d2 + Memi[m[jmax]+i1] = Memi[m[j]+i1] + } else { + Memr[kmax] = d1 + Memi[m[jmax]+i1] = Memi[m[n1]+i1] + } + } + if (jmin < j) { + if (jmax != n1) { + Memr[kmin] = d1 + Memi[m[jmin]+i1] = Memi[m[n1]+i1] + } else { + Memr[kmin] = d2 + Memi[m[jmin]+i1] = Memi[m[j]+i1] + } + } + } else { + if (jmax < j) { + if (jmin != j) + Memr[kmax] = d2 + else + Memr[kmax] = d1 + } + if (jmin < j) { + if (jmax != n1) + Memr[kmin] = d1 + else + Memr[kmin] = d2 + } + } + n1 = n1 - 2 + } + + # Reject the excess low points. + do np = 1, nlow { + k = d[1] + i1 + d1 = Memr[k] + dmin = d1; jmin = 1; kmin = k + do j = 2, n1 { + k = d[j] + i1 + d1 = Memr[k] + if (d1 < dmin) { + dmin = d1; jmin = j; kmin = k + } + } + if (keepids) { + if (jmin < n1) { + Memr[kmin] = d1 + Memi[m[jmin]+i1] = Memi[m[n1]+i1] + } + } else { + if (jmin < n1) + Memr[kmin] = d1 + } + n1 = n1 - 1 + } + + # Reject the excess high points. + do np = 1, nhigh { + k = d[1] + i1 + d1 = Memr[k] + dmax = d1; jmax = 1; kmax = k + do j = 2, n1 { + k = d[j] + i1 + d1 = Memr[k] + if (d1 > dmax) { + dmax = d1; jmax = j; kmax = k + } + } + if (keepids) { + if (jmax < n1) { + Memr[kmax] = d1 + Memi[m[jmax]+i1] = Memi[m[n1]+i1] + } + } else { + if (jmax < n1) + Memr[kmax] = d1 + } + n1 = n1 - 1 + } + n[i] = n1 + } + + if (dflag == D_ALL && npairs + nlow + nhigh > 0) + dflag = D_MIX +end diff --git a/noao/onedspec/scombine/generic/icpclip.x b/noao/onedspec/scombine/generic/icpclip.x new file mode 100644 index 00000000..d9028a93 --- /dev/null +++ b/noao/onedspec/scombine/generic/icpclip.x @@ -0,0 +1,224 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "../icombine.h" + +define MINCLIP 3 # Minimum number for clipping + + +# IC_PCLIP -- Percentile clip +# +# 1) Find the median +# 2) Find the pixel which is the specified order index away +# 3) Use the data value difference as a sigma and apply clipping +# 4) Since the median is known return it so it does not have to be recomputed + +procedure ic_pclipr (d, m, n, nimages, npts, median) + +pointer d[ARB] # Data pointers +pointer m[ARB] # Image id pointers +int n[npts] # Number of good pixels +int nimages # Number of input images +int npts # Number of output points per line +real median[npts] # Median + +int i, j, k, l, id, n1, n2, n3, n4, n5, nl, nh, nin, maxkeep +bool even, fp_equalr() +real sigma, r, s, t +pointer sp, resid, mp1, mp2 +real med + +include "../icombine.com" + +begin + # There must be at least MINCLIP and more than nkeep pixels. + if (nkeep < 0) + maxkeep = max (0, nimages + nkeep) + else + maxkeep = min (nimages, nkeep) + if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) { + docombine = true + return + } + + call smark (sp) + call salloc (resid, nimages+1, TY_REAL) + + # Set sign of pclip parameter + if (pclip < 0) + t = -1. + else + t = 1. + + # If there are no rejected pixels compute certain parameters once. + if (dflag == D_ALL) { + n1 = n[1] + if (nkeep < 0) + maxkeep = max (0, n1 + nkeep) + else + maxkeep = min (n1, nkeep) + n2 = 1 + n1 / 2 + even = (mod (n1, 2) == 0) + if (pclip < 0.) { + if (even) + n3 = max (1, nint (n2 - 1 + pclip)) + else + n3 = max (1, nint (n2 + pclip)) + } else + n3 = min (n1, nint (n2 + pclip)) + nin = n1 + } + + # Now apply clipping. + do i = 1, npts { + # Compute median. + if (dflag == D_MIX) { + n1 = n[i] + if (nkeep < 0) + maxkeep = max (0, n1 + nkeep) + else + maxkeep = min (n1, nkeep) + if (n1 == 0) { + if (combine == MEDIAN) + median[i] = blank + next + } + n2 = 1 + n1 / 2 + even = (mod (n1, 2) == 0) + if (pclip < 0) { + if (even) + n3 = max (1, nint (n2 - 1 + pclip)) + else + n3 = max (1, nint (n2 + pclip)) + } else + n3 = min (n1, nint (n2 + pclip)) + } + + j = i - 1 + if (even) { + med = Memr[d[n2-1]+j] + med = (med + Memr[d[n2]+j]) / 2. + } else + med = Memr[d[n2]+j] + + if (n1 < max (MINCLIP, maxkeep+1)) { + if (combine == MEDIAN) + median[i] = med + next + } + + # Define sigma for clipping + sigma = t * (Memr[d[n3]+j] - med) + if (fp_equalr (sigma, 0.)) { + if (combine == MEDIAN) + median[i] = med + next + } + + # Reject pixels and save residuals. + # Check if any pixels are clipped. + # If so recompute the median and reset the number of good pixels. + # Only reorder if needed. + + for (nl=1; nl<=n1; nl=nl+1) { + r = (med - Memr[d[nl]+j]) / sigma + if (r < lsigma) + break + Memr[resid+nl] = r + } + for (nh=n1; nh>=1; nh=nh-1) { + r = (Memr[d[nh]+j] - med) / sigma + if (r < hsigma) + break + Memr[resid+nh] = r + } + n4 = nh - nl + 1 + + # If too many pixels are rejected add some back in. + # All pixels with the same residual are added. + while (n4 < maxkeep) { + if (nl == 1) + nh = nh + 1 + else if (nh == n[i]) + nl = nl - 1 + else { + r = Memr[resid+nl-1] + s = Memr[resid+nh+1] + if (r < s) { + nl = nl - 1 + r = r + TOL + if (s <= r) + nh = nh + 1 + if (nl > 1) { + if (Memr[resid+nl-1] <= r) + nl = nl - 1 + } + } else { + nh = nh + 1 + s = s + TOL + if (r <= s) + nl = nl - 1 + if (nh < n2) { + if (Memr[resid+nh+1] <= s) + nh = nh + 1 + } + } + } + n4 = nh - nl + 1 + } + + # If any pixels are rejected recompute the median. + if (nl > 1 || nh < n1) { + n5 = nl + n4 / 2 + if (mod (n4, 2) == 0) { + med = Memr[d[n5-1]+j] + med = (med + Memr[d[n5]+j]) / 2. + } else + med = Memr[d[n5]+j] + n[i] = n4 + } + if (combine == MEDIAN) + median[i] = med + + # Reorder if pixels only if necessary. + if (nl > 1 && (combine != MEDIAN || grow > 0)) { + k = max (nl, n4 + 1) + if (keepids) { + do l = 1, min (n1, nl-1) { + Memr[d[l]+j] = Memr[d[k]+j] + if (grow > 0) { + mp1 = m[l] + j + mp2 = m[k] + j + id = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = id + } else + Memi[m[l]+j] = Memi[m[k]+j] + k = k + 1 + } + } else { + do l = 1, min (n1, nl - 1) { + Memr[d[l]+j] = Memr[d[k]+j] + k = k + 1 + } + } + } + } + + # Check if data flag needs to be reset for rejected pixels. + if (dflag == D_ALL) { + do i = 1, npts { + if (n[i] != nin) { + dflag = D_MIX + break + } + } + } + + # Flag whether the median has been computed. + if (combine == MEDIAN) + docombine = false + else + docombine = true + + call sfree (sp) +end diff --git a/noao/onedspec/scombine/generic/icsclip.x b/noao/onedspec/scombine/generic/icsclip.x new file mode 100644 index 00000000..e38f7935 --- /dev/null +++ b/noao/onedspec/scombine/generic/icsclip.x @@ -0,0 +1,486 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "../icombine.h" + +define MINCLIP 3 # Mininum number of images for algorithm + + +# IC_ASIGCLIP -- Reject pixels using sigma clipping about the average +# The initial average rejects the high and low pixels. A correction for +# different scalings of the images may be made. Weights are not used. + +procedure ic_asigclipr (d, m, n, scales, zeros, nimages, npts, average) + +pointer d[nimages] # Data pointers +pointer m[nimages] # Image id pointers +int n[npts] # Number of good pixels +real scales[nimages] # Scales +real zeros[nimages] # Zeros +int nimages # Number of images +int npts # Number of output points per line +real average[npts] # Average + +int i, j, k, l, jj, n1, n2, nin, nk, maxkeep +real d1, low, high, sum, a, s, r, one +data one /1.0/ +pointer sp, resid, w, wp, dp1, dp2, mp1, mp2 + +include "../icombine.com" + +begin + # If there are insufficient pixels go on to the combining + if (nkeep < 0) + maxkeep = max (0, nimages + nkeep) + else + maxkeep = min (nimages, nkeep) + if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) { + docombine = true + return + } + + # Flag whether returned average needs to be recomputed. + if (dowts || combine != AVERAGE) + docombine = true + else + docombine = false + + # Save the residuals and the sigma scaling corrections if needed. + call smark (sp) + call salloc (resid, nimages+1, TY_REAL) + if (doscale1) + call salloc (w, nimages, TY_REAL) + + # Do sigma clipping. + nin = n[1] + do i = 1, npts { + k = i - 1 + n1 = n[i] + if (nkeep < 0) + maxkeep = max (0, n1 + nkeep) + else + maxkeep = min (n1, nkeep) + + # If there are not enough pixels simply compute the average. + if (n1 < max (3, maxkeep)) { + if (!docombine) { + if (n1 == 0) + average[i] = blank + else { + sum = Memr[d[1]+k] + do j = 2, n1 + sum = sum + Memr[d[j]+k] + average[i] = sum / n1 + } + } + next + } + + # Compute average with the high and low rejected. + low = Memr[d[1]+k] + high = Memr[d[2]+k] + if (low > high) { + d1 = low + low = high + high = d1 + } + sum = 0. + do j = 3, n1 { + d1 = Memr[d[j]+k] + if (d1 < low) { + sum = sum + low + low = d1 + } else if (d1 > high) { + sum = sum + high + high = d1 + } else + sum = sum + d1 + } + a = sum / (n1 - 2) + sum = sum + low + high + + # Iteratively reject pixels and compute the final average if needed. + # Compact the data and keep track of the image IDs if needed. + + repeat { + n2 = n1 + if (doscale1) { + # Compute sigma corrected for scaling. + s = 0. + wp = w - 1 + do j = 1, n1 { + dp1 = d[j] + k + mp1 = m[j] + k + wp = wp + 1 + + d1 = Memr[dp1] + l = Memi[mp1] + r = sqrt (max (one, (a + zeros[l]) / scales[l])) + s = s + ((d1 - a) / r) ** 2 + Memr[wp] = r + } + s = sqrt (s / (n1 - 1)) + + # Reject pixels. Save the residuals and data values. + wp = w - 1 + if (s > 0.) { + for (j=1; j<=n1; j=j+1) { + dp1 = d[j] + k + mp1 = m[j] + k + wp = wp + 1 + + d1 = Memr[dp1] + r = (d1 - a) / (s * Memr[wp]) + if (r < -lsigma || r > hsigma) { + Memr[resid+n1] = abs (r) + if (j < n1) { + dp2 = d[n1] + k + Memr[dp1] = Memr[dp2] + Memr[dp2] = d1 + Memr[wp] = Memr[w+n1-1] + mp2 = m[n1] + k + l = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = l + j = j - 1 + } + sum = sum - d1 + n1 = n1 - 1 + } + } + } + } else { + # Compute the sigma without scale correction. + s = 0. + do j = 1, n1 + s = s + (Memr[d[j]+k] - a) ** 2 + s = sqrt (s / (n1 - 1)) + + # Reject pixels. Save the residuals and data values. + if (s > 0.) { + for (j=1; j<=n1; j=j+1) { + dp1 = d[j] + k + d1 = Memr[dp1] + r = (d1 - a) / s + if (r < -lsigma || r > hsigma) { + Memr[resid+n1] = abs (r) + if (j < n1) { + dp2 = d[n1] + k + Memr[dp1] = Memr[dp2] + Memr[dp2] = d1 + if (keepids) { + mp1 = m[j] + k + mp2 = m[n1] + k + l = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = l + } + j = j - 1 + } + sum = sum - d1 + n1 = n1 - 1 + } + } + } + } + + # Recompute the average. + if (n1 > 1) + a = sum / n1 + } until (n1 == n2 || n1 <= max (2, maxkeep)) + + # If too many pixels are rejected add some back. + # All pixels with equal residuals are added back. + if (n1 < maxkeep) { + nk = maxkeep + if (doscale1) { + for (j=n1+1; j<=nk; j=j+1) { + dp1 = d[j] + k + mp1 = m[j] + k + r = Memr[resid+j] + jj = 0 + do l = j+1, n2 { + s = Memr[resid+l] + if (s < r + TOL) { + if (s > r - TOL) + jj = jj + 1 + else { + jj = 0 + Memr[resid+l] = r + r = s + dp2 = d[l] + k + d1 = Memr[dp1] + Memr[dp1] = Memr[dp2] + Memr[dp2] = d1 + mp2 = m[l] + k + s = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = s + } + } + } + sum = sum + Memr[dp1] + n1 = n1 + 1 + nk = max (nk, j+jj) + } + } else { + for (j=n1+1; j<=nk; j=j+1) { + dp1 = d[j] + k + r = Memr[resid+j] + jj = 0 + do l = j+1, n2 { + s = Memr[resid+l] + if (s < r + TOL) { + if (s > r - TOL) + jj = jj + 1 + else { + jj = 0 + Memr[resid+l] = r + r = s + dp2 = d[l] + k + d1 = Memr[dp1] + Memr[dp1] = Memr[dp2] + Memr[dp2] = d1 + if (keepids) { + mp1 = m[j] + k + mp2 = m[l] + k + s = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = s + } + } + } + } + sum = sum + Memr[dp1] + n1 = n1 + 1 + nk = max (nk, j+jj) + } + } + + # Recompute the average. + if (n1 > 1) + a = sum / n1 + } + + # Save the average if needed. + n[i] = n1 + if (!docombine) { + if (n1 > 0) + average[i] = a + else + average[i] = blank + } + } + + # Check if the data flag has to be reset for rejected pixels + if (dflag == D_ALL) { + do i = 1, npts { + if (n[i] != nin) { + dflag = D_MIX + break + } + } + } + + call sfree (sp) +end + + +# IC_MSIGCLIP -- Reject pixels using sigma clipping about the median + +procedure ic_msigclipr (d, m, n, scales, zeros, nimages, npts, median) + +pointer d[nimages] # Data pointers +pointer m[nimages] # Image id pointers +int n[npts] # Number of good pixels +real scales[nimages] # Scales +real zeros[nimages] # Zeros +int nimages # Number of images +int npts # Number of output points per line +real median[npts] # Median + +int i, j, k, l, id, n1, n2, n3, nl, nh, nin, maxkeep +real r, s +pointer sp, resid, w, mp1, mp2 +real med, one +data one /1.0/ + +include "../icombine.com" + +begin + # If there are insufficient pixels go on to the combining + if (nkeep < 0) + maxkeep = max (0, nimages + nkeep) + else + maxkeep = min (nimages, nkeep) + if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) { + docombine = true + return + } + + # Save the residuals and sigma scaling corrections if needed. + call smark (sp) + call salloc (resid, nimages+1, TY_REAL) + if (doscale1) + call salloc (w, nimages, TY_REAL) + + # Compute median and sigma and iteratively clip. + nin = n[1] + do i = 1, npts { + k = i - 1 + n1 = n[i] + if (nkeep < 0) + maxkeep = max (0, n1 + nkeep) + else + maxkeep = min (n1, nkeep) + nl = 1 + nh = n1 + + repeat { + n2 = n1 + n3 = nl + n1 / 2 + + if (n1 == 0) + med = blank + else if (mod (n1, 2) == 0) + med = (Memr[d[n3-1]+k] + Memr[d[n3]+k]) / 2. + else + med = Memr[d[n3]+k] + + if (n1 >= max (MINCLIP, maxkeep+1)) { + if (doscale1) { + # Compute the sigma with scaling correction. + s = 0. + do j = nl, nh { + l = Memi[m[j]+k] + r = sqrt (max (one, (med + zeros[l]) / scales[l])) + s = s + ((Memr[d[j]+k] - med) / r) ** 2 + Memr[w+j-1] = r + } + s = sqrt (s / (n1 - 1)) + + # Reject pixels and save the residuals. + if (s > 0.) { + for (; nl <= n2; nl = nl + 1) { + r = (med - Memr[d[nl]+k]) / (s * Memr[w+nl-1]) + if (r <= lsigma) + break + Memr[resid+nl] = r + n1 = n1 - 1 + } + for (; nh >= nl; nh = nh - 1) { + r = (Memr[d[nh]+k] - med) / (s * Memr[w+nh-1]) + if (r <= hsigma) + break + Memr[resid+nh] = r + n1 = n1 - 1 + } + } + } else { + # Compute the sigma without scaling correction. + s = 0. + do j = nl, nh + s = s + (Memr[d[j]+k] - med) ** 2 + s = sqrt (s / (n1 - 1)) + + # Reject pixels and save the residuals. + if (s > 0.) { + for (; nl <= n2; nl = nl + 1) { + r = (med - Memr[d[nl]+k]) / s + if (r <= lsigma) + break + Memr[resid+nl] = r + n1 = n1 - 1 + } + for (; nh >= nl; nh = nh - 1) { + r = (Memr[d[nh]+k] - med) / s + if (r <= hsigma) + break + Memr[resid+nh] = r + n1 = n1 - 1 + } + } + } + } + } until (n1 == n2 || n1 < max (MINCLIP, maxkeep+1)) + + # If too many pixels are rejected add some back. + # All pixels with equal residuals are added back. + while (n1 < maxkeep) { + if (nl == 1) + nh = nh + 1 + else if (nh == n[i]) + nl = nl - 1 + else { + r = Memr[resid+nl-1] + s = Memr[resid+nh+1] + if (r < s) { + nl = nl - 1 + r = r + TOL + if (s <= r) + nh = nh + 1 + if (nl > 1) { + if (Memr[resid+nl-1] <= r) + nl = nl - 1 + } + } else { + nh = nh + 1 + s = s + TOL + if (r <= s) + nl = nl - 1 + if (nh < n2) { + if (Memr[resid+nh+1] <= s) + nh = nh + 1 + } + } + } + n1 = nh - nl + 1 + } + + # Only set median and reorder if needed + n[i] = n1 + if (n1 > 0 && nl > 1 && (combine != MEDIAN || grow > 0)) { + j = max (nl, n1 + 1) + if (keepids) { + do l = 1, min (n1, nl-1) { + Memr[d[l]+k] = Memr[d[j]+k] + if (grow > 0) { + mp1 = m[l] + k + mp2 = m[j] + k + id = Memi[mp1] + Memi[mp1] = Memi[mp2] + Memi[mp2] = id + } else + Memi[m[l]+k] = Memi[m[j]+k] + j = j + 1 + } + } else { + do l = 1, min (n1, nl - 1) { + Memr[d[l]+k] = Memr[d[j]+k] + j = j + 1 + } + } + } + + if (combine == MEDIAN) + median[i] = med + } + + # Check if data flag needs to be reset for rejected pixels + if (dflag == D_ALL) { + do i = 1, npts { + if (n[i] != nin) { + dflag = D_MIX + break + } + } + } + + # Flag that the median has been computed. + if (combine == MEDIAN) + docombine = false + else + docombine = true + + call sfree (sp) +end + diff --git a/noao/onedspec/scombine/generic/icsort.x b/noao/onedspec/scombine/generic/icsort.x new file mode 100644 index 00000000..f3d2fb21 --- /dev/null +++ b/noao/onedspec/scombine/generic/icsort.x @@ -0,0 +1,275 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +define LOGPTR 32 # log2(maxpts) (4e9) + + +# IC_SORT -- Quicksort. This is based on the VOPS asrt except that +# the input is an array of pointers to image lines and the sort is done +# across the image lines at each point along the lines. The number of +# valid pixels at each point is allowed to vary. The cases of 1, 2, and 3 +# pixels per point are treated specially. + +procedure ic_sortr (a, b, nvecs, npts) + +pointer a[ARB] # pointer to input vectors +real b[ARB] # work array +int nvecs[npts] # number of vectors +int npts # number of points in vectors + +real pivot, temp, temp3 +int i, j, k, l, p, npix, lv[LOGPTR], uv[LOGPTR] +define swap {temp=$1;$1=$2;$2=temp} +define copy_ 10 + +begin + do l = 0, npts-1 { + npix = nvecs[l+1] + if (npix <= 1) + next + + do i = 1, npix + b[i] = Memr[a[i]+l] + + # Special cases + if (npix <= 3) { + pivot = b[1] + temp = b[2] + if (npix == 2) { + if (temp < pivot) { + b[1] = temp + b[2] = pivot + } else + next + } else { + temp3 = b[3] + if (temp < pivot) { # bac|bca|cba + if (temp < temp3) { # bac|bca + b[1] = temp + if (pivot < temp3) # bac + b[2] = pivot + else { # bca + b[2] = temp3 + b[3] = pivot + } + } else { # cba + b[1] = temp3 + b[3] = pivot + } + } else if (temp3 < temp) { # acb|cab + b[3] = temp + if (pivot < temp3) # acb + b[2] = temp3 + else { # cab + b[1] = temp3 + b[2] = pivot + } + } else + next + } + goto copy_ + } + + # General case + do i = 1, npix + b[i] = Memr[a[i]+l] + + lv[1] = 1 + uv[1] = npix + p = 1 + + while (p > 0) { + if (lv[p] >= uv[p]) # only one elem in this subset + p = p - 1 # pop stack + else { + # Dummy do loop to trigger the Fortran optimizer. + do p = p, ARB { + i = lv[p] - 1 + j = uv[p] + + # Select as the pivot the element at the center of the + # array, to avoid quadratic behavior on an already + # sorted array. + + k = (lv[p] + uv[p]) / 2 + swap (b[j], b[k]) + pivot = b[j] # pivot line + + while (i < j) { + for (i=i+1; b[i] < pivot; i=i+1) + ; + for (j=j-1; j > i; j=j-1) + if (b[j] <= pivot) + break + if (i < j) # out of order pair + swap (b[i], b[j]) # interchange elements + } + + j = uv[p] # move pivot to position i + swap (b[i], b[j]) # interchange elements + + if (i-lv[p] < uv[p] - i) { # stack so shorter done first + lv[p+1] = lv[p] + uv[p+1] = i - 1 + lv[p] = i + 1 + } else { + lv[p+1] = i + 1 + uv[p+1] = uv[p] + uv[p] = i - 1 + } + + break + } + p = p + 1 # push onto stack + } + } + +copy_ + do i = 1, npix + Memr[a[i]+l] = b[i] + } +end + + +# IC_2SORT -- Quicksort. This is based on the VOPS asrt except that +# the input is an array of pointers to image lines and the sort is done +# across the image lines at each point along the lines. The number of +# valid pixels at each point is allowed to vary. The cases of 1, 2, and 3 +# pixels per point are treated specially. A second integer set of +# vectors is sorted. + +procedure ic_2sortr (a, b, c, d, nvecs, npts) + +pointer a[ARB] # pointer to input vectors +real b[ARB] # work array +pointer c[ARB] # pointer to associated integer vectors +int d[ARB] # work array +int nvecs[npts] # number of vectors +int npts # number of points in vectors + +real pivot, temp, temp3 +int i, j, k, l, p, npix, lv[LOGPTR], uv[LOGPTR], itemp +define swap {temp=$1;$1=$2;$2=temp} +define iswap {itemp=$1;$1=$2;$2=itemp} +define copy_ 10 + +begin + do l = 0, npts-1 { + npix = nvecs[l+1] + if (npix <= 1) + next + + do i = 1, npix { + b[i] = Memr[a[i]+l] + d[i] = Memi[c[i]+l] + } + + # Special cases + if (npix <= 3) { + pivot = b[1] + temp = b[2] + if (npix == 2) { + if (temp < pivot) { + b[1] = temp + b[2] = pivot + iswap (d[1], d[2]) + } else + next + } else { + temp3 = b[3] + if (temp < pivot) { # bac|bca|cba + if (temp < temp3) { # bac|bca + b[1] = temp + if (pivot < temp3) { # bac + b[2] = pivot + iswap (d[1], d[2]) + } else { # bca + b[2] = temp3 + b[3] = pivot + itemp = d[2] + d[2] = d[3] + } + } else { # cba + b[1] = temp3 + b[3] = pivot + iswap (d[1], d[3]) + } + } else if (temp3 < temp) { # acb|cab + b[3] = temp + if (pivot < temp3) { # acb + b[2] = temp3 + iswap (d[2], d[3]) + } else { # cab + b[1] = temp3 + b[2] = pivot + itemp = d[2] + d[2] = d[1] + d[1] = d[3] + d[3] = itemp + } + } else + next + } + goto copy_ + } + + # General case + lv[1] = 1 + uv[1] = npix + p = 1 + + while (p > 0) { + if (lv[p] >= uv[p]) # only one elem in this subset + p = p - 1 # pop stack + else { + # Dummy do loop to trigger the Fortran optimizer. + do p = p, ARB { + i = lv[p] - 1 + j = uv[p] + + # Select as the pivot the element at the center of the + # array, to avoid quadratic behavior on an already + # sorted array. + + k = (lv[p] + uv[p]) / 2 + swap (b[j], b[k]); swap (d[j], d[k]) + pivot = b[j] # pivot line + + while (i < j) { + for (i=i+1; b[i] < pivot; i=i+1) + ; + for (j=j-1; j > i; j=j-1) + if (b[j] <= pivot) + break + if (i < j) { # out of order pair + swap (b[i], b[j]) # interchange elements + swap (d[i], d[j]) + } + } + + j = uv[p] # move pivot to position i + swap (b[i], b[j]) # interchange elements + swap (d[i], d[j]) + + if (i-lv[p] < uv[p] - i) { # stack so shorter done first + lv[p+1] = lv[p] + uv[p+1] = i - 1 + lv[p] = i + 1 + } else { + lv[p+1] = i + 1 + uv[p+1] = uv[p] + uv[p] = i - 1 + } + + break + } + p = p + 1 # push onto stack + } + } + +copy_ + do i = 1, npix { + Memr[a[i]+l] = b[i] + Memi[c[i]+l] = d[i] + } + } +end diff --git a/noao/onedspec/scombine/generic/mkpkg b/noao/onedspec/scombine/generic/mkpkg new file mode 100644 index 00000000..4d371363 --- /dev/null +++ b/noao/onedspec/scombine/generic/mkpkg @@ -0,0 +1,16 @@ +$checkout libpkg.a ../ +$update libpkg.a +$checkin libpkg.a ../ +$exit + +libpkg.a: + icaclip.x ../icombine.com ../icombine.h + icaverage.x ../icombine.com ../icombine.h + iccclip.x ../icombine.com ../icombine.h + icgrow.x ../icombine.com ../icombine.h + icmedian.x ../icombine.com ../icombine.h + icmm.x ../icombine.com ../icombine.h + icpclip.x ../icombine.com ../icombine.h + icsclip.x ../icombine.com ../icombine.h + icsort.x + ; diff --git a/noao/onedspec/scombine/icgdata.x b/noao/onedspec/scombine/icgdata.x new file mode 100644 index 00000000..907adc5e --- /dev/null +++ b/noao/onedspec/scombine/icgdata.x @@ -0,0 +1,199 @@ +include +include "icombine.h" + + +# IC_GDATAR - Apply threshold, scaling, and masking + +procedure ic_gdatar (sh, d, id, n, m, lflag, scales, zeros, nimages, npts) + +pointer sh[nimages] # Input spectra structures +pointer d[nimages] # Data pointers +pointer id[nimages] # ID pointers +int n[npts] # Number of good pixels +pointer m[nimages] # Mask pointers +int lflag[nimages] # Empty mask flags +real scales[nimages] # Scale factors +real zeros[nimages] # Zero offset factors +int nimages # Number of spectra +int npts # NUmber of output points + +int i, j, k, l, nused +real a, b +pointer dp, ip, mp + +include "icombine.com" + +begin + # Set data vectors + do i = 1, nimages { + d[i] = SY(sh[i]) + m[i] = SX(sh[i]) + } + + # Apply threshold if needed + if (dothresh) { + do i = 1, nimages { + dp = d[i] + if (lflag[i] == D_ALL) { + do j = 1, npts { + a = Memr[dp] + if (a < lthresh || a > hthresh) { + Memr[m[i]+j-1] = 1 + lflag[i] = D_MIX + dflag = D_MIX + } + dp = dp + 1 + } + } else if (lflag[i] == D_MIX) { + mp = m[i] + do j = 1, npts { + if (Memr[mp] == 0) { + a = Memr[dp] + if (a < lthresh || a > hthresh) { + Memr[m[i]+j-1] = 1 + dflag = D_MIX + } + } + dp = dp + 1 + mp = mp + 1 + } + } + + # Check for completely empty lines + if (lflag[i] == D_MIX) { + lflag[i] = D_NONE + mp = m[i] + do j = 1, npts { + if (Memr[mp] == 0) { + lflag[i] = D_MIX + break + } + mp = mp + 1 + } + } + } + } + + # Apply scaling (avoiding masked pixels which might overflow?) + if (doscale) { + if (dflag == D_ALL) { + do i = 1, nimages { + dp = d[i] + a = scales[i] + b = -zeros[i] + do j = 1, npts { + Memr[dp] = Memr[dp] / a + b + dp = dp + 1 + } + } + } else if (dflag == D_MIX) { + do i = 1, nimages { + dp = d[i] + a = scales[i] + b = -zeros[i] + if (lflag[i] == D_ALL) { + do j = 1, npts { + Memr[dp] = Memr[dp] / a + b + dp = dp + 1 + } + } else if (lflag[i] == D_MIX) { + mp = m[i] + do j = 1, npts { + if (Memr[mp] == 0) + Memr[dp] = Memr[dp] / a + b + dp = dp + 1 + mp = mp + 1 + } + } + } + } + } + + # Sort pointers to exclude unused images. + # Use the lflag array to keep track of the image index. + + if (dflag == D_ALL) + nused = nimages + else { + nused = 0 + do i = 1, nimages + if (lflag[i] != D_NONE) { + nused = nused + 1 + d[nused] = d[i] + m[nused] = m[i] + lflag[nused] = i + } + if (nused == 0) + dflag = D_NONE + } + + # Compact data to remove bad pixels + # Keep track of the image indices if needed + # If growing mark the end of the included image indices with zero + + if (dflag == D_ALL) { + call amovki (nused, n, npts) + if (keepids) + do i = 1, nimages + call amovki (i, Memi[id[i]], npts) + } else if (dflag == D_NONE) + call aclri (n, npts) + else { + call aclri (n, npts) + if (keepids) { + do i = 1, nused { + l = lflag[i] + dp = d[i] + ip = id[i] + mp = m[i] + do j = 1, npts { + if (Memr[mp] == 0) { + n[j] = n[j] + 1 + k = n[j] + if (k < i) { + Memr[d[k]+j-1] = Memr[dp] + Memi[id[k]+j-1] = l + } else + Memi[ip] = l + } + dp = dp + 1 + ip = ip + 1 + mp = mp + 1 + } + } + if (grow > 0) { + do j = 0, npts-1 { + do i = n[i]+1, nimages + Memi[id[i]+j] = 0 + } + } + } else { + do i = 1, nused { + dp = d[i] + mp = m[i] + do j = 1, npts { + if (Memr[mp] == 0) { + n[j] = n[j] + 1 + k = n[j] + if (k < i) + Memr[d[k]+j-1] = Memr[dp] + } + dp = dp + 1 + mp = mp + 1 + } + } + } + } + + # Sort the pixels and IDs if needed + if (mclip) { + call malloc (dp, nimages, TY_REAL) + if (keepids) { + call malloc (ip, nimages, TY_INT) + call ic_2sortr (d, Memr[dp], id, Memi[ip], n, npts) + call mfree (ip, TY_INT) + } else + call ic_sortr (d, Memr[dp], n, npts) + call mfree (dp, TY_REAL) + } +end diff --git a/noao/onedspec/scombine/iclog.x b/noao/onedspec/scombine/iclog.x new file mode 100644 index 00000000..29002c0f --- /dev/null +++ b/noao/onedspec/scombine/iclog.x @@ -0,0 +1,301 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include +include +include "icombine.h" + +# IC_LOG -- Output log information is a log file has been specfied. + +procedure ic_log (sh, shout, ncombine, exptime, sname, zname, wname, + mode, median, mean, scales, zeros, wts, nimages, + dozero, nout, expname, exposure) + +pointer sh[nimages] # Input spectra +pointer shout # Output spectrum +int ncombine[nimages] # Number of previous combined images +real exptime[nimages] # Exposure times +char sname[ARB] # Scale name +char zname[ARB] # Zero name +char wname[ARB] # Weight name +real mode[nimages] # Modes +real median[nimages] # Medians +real mean[nimages] # Means +real scales[nimages] # Scale factors +real zeros[nimages] # Zero or sky levels +real wts[nimages] # Weights +int nimages # Number of images +bool dozero # Zero flag +int nout # Number of images combined in output +char expname[ARB] # Exposure name +real exposure # Output exposure + +int i, j, ctor() +real rval +long clktime() +bool prncombine, prexptime, prmode, prmedian, prmean +bool prrdn, prgain, prsn, prscale, przero, prwts, strne() +pointer sp, fname + +include "icombine.com" + +begin + if (logfd == NULL) + return + + call smark (sp) + call salloc (fname, SZ_LINE, TY_CHAR) + + # Time stamp the log and print parameter information. + + call cnvdate (clktime(0), Memc[fname], SZ_LINE) + call fprintf (logfd, "\n%s: SCOMBINE\n") + call pargstr (Memc[fname]) + switch (combine) { + case AVERAGE: + call fprintf (logfd, " combine = average,") + case MEDIAN: + call fprintf (logfd, " combine = median,") + case SUM: + call fprintf (logfd, " combine = sum\n") + } + if (combine != SUM) { + call fprintf (logfd, " scale = %s, zero = %s, weight = %s\n") + call pargstr (sname) + call pargstr (zname) + call pargstr (wname) + } + + switch (reject) { + case MINMAX: + call fprintf (logfd, " reject = minmax, nlow = %d, nhigh = %d\n") + call pargi (nint (flow * nimages)) + call pargi (nint (fhigh * nimages)) + case CCDCLIP: + call fprintf (logfd, " reject = ccdclip, mclip = %b, nkeep = %d\n") + call pargb (mclip) + call pargi (nkeep) + call fprintf (logfd, + " rdnoise = %s, gain = %s, snoise = %s, sigma = %g, hsigma = %g\n") + call pargstr (Memc[rdnoise]) + call pargstr (Memc[gain]) + call pargstr (Memc[snoise]) + call pargr (lsigma) + call pargr (hsigma) + case CRREJECT: + call fprintf (logfd, + " reject = crreject, mclip = %b, nkeep = %d\n") + call pargb (mclip) + call pargi (nkeep) + call fprintf (logfd, + " rdnoise = %s, gain = %s, snoise = %s, hsigma = %g\n") + call pargstr (Memc[rdnoise]) + call pargstr (Memc[gain]) + call pargstr (Memc[snoise]) + call pargr (hsigma) + case PCLIP: + call fprintf (logfd, " reject = pclip, nkeep = %d\n") + call pargi (nkeep) + call fprintf (logfd, " pclip = %g, lsigma = %g, hsigma = %g\n") + call pargr (pclip) + call pargr (lsigma) + call pargr (hsigma) + case SIGCLIP: + call fprintf (logfd, " reject = sigclip, mclip = %b, nkeep = %d\n") + call pargb (mclip) + call pargi (nkeep) + call fprintf (logfd, " lsigma = %g, hsigma = %g\n") + call pargr (lsigma) + call pargr (hsigma) + case AVSIGCLIP: + call fprintf (logfd, + " reject = avsigclip, mclip = %b, nkeep = %d\n") + call pargb (mclip) + call pargi (nkeep) + call fprintf (logfd, " lsigma = %g, hsigma = %g\n") + call pargr (lsigma) + call pargr (hsigma) + } + if (reject != NONE && grow > 0) { + call fprintf (logfd, " grow = %d\n") + call pargi (grow) + } + if (dothresh) { + if (lthresh > -MAX_REAL && hthresh < MAX_REAL) { + call fprintf (logfd, " lthreshold = %g, hthreshold = %g\n") + call pargr (lthresh) + call pargr (hthresh) + } else if (lthresh > -MAX_REAL) { + call fprintf (logfd, " lthreshold = %g\n") + call pargr (lthresh) + } else { + call fprintf (logfd, " hthreshold = %g\n") + call pargr (hthresh) + } + } + call fprintf (logfd, " blank = %g\n") + call pargr (blank) + call clgstr ("sample", Memc[fname], SZ_LINE) + if (Memc[fname] != EOS) { + call fprintf (logfd, " sample = %s\n") + call pargstr (Memc[fname]) + } + + # Print information pertaining to individual images as a set of + # columns with the image name being the first column. Determine + # what information is relevant and print the appropriate header. + + prncombine = false + prexptime = (expname[1] != EOS) + prscale = (doscale || strne (sname, "none")) + przero = (dozero || strne (zname, "none")) + prwts = (dowts || strne (wname, "none")) + prmode = false + prmedian = false + prmean = false + prrdn = false + prgain = false + prsn = false + do i = 1, nimages { + if (ncombine[i] != ncombine[1]) + prncombine = true + if (exptime[i] != exptime[1]) + prexptime = true + if (mode[i] != mode[1]) + prmode = true + if (median[i] != median[1]) + prmedian = true + if (mean[i] != mean[1]) + prmean = true + if (reject == CCDCLIP || reject == CRREJECT) { + j = 1 + if (ctor (Memc[rdnoise], j, rval) == 0) + prrdn = true + j = 1 + if (ctor (Memc[gain], j, rval) == 0) + prgain = true + j = 1 + if (ctor (Memc[snoise], j, rval) == 0) + prsn = true + } + } + + call fprintf (logfd, " %20s ") + call pargstr ("Images") + if (prncombine) { + call fprintf (logfd, " %6s") + call pargstr ("N") + } + if (prexptime) { + call fprintf (logfd, " %6s") + call pargstr ("Exp") + } + if (prmode) { + call fprintf (logfd, " %7s") + call pargstr ("Mode") + } + if (prmedian) { + call fprintf (logfd, " %7s") + call pargstr ("Median") + } + if (prmean) { + call fprintf (logfd, " %7s") + call pargstr ("Mean") + } + if (prrdn) { + call fprintf (logfd, " %7s") + call pargstr ("Rdnoise") + } + if (prgain) { + call fprintf (logfd, " %6s") + call pargstr ("Gain") + } + if (prsn) { + call fprintf (logfd, " %6s") + call pargstr ("Snoise") + } + if (prscale) { + call fprintf (logfd, " %6s") + call pargstr ("Scale") + } + if (przero) { + call fprintf (logfd, " %7s") + call pargstr ("Zero") + } + if (prwts) { + call fprintf (logfd, " %6s") + call pargstr ("Weight") + } + call fprintf (logfd, "\n") + + do i = 1, nimages { + call fprintf (logfd, " %16s[%3d]") + call pargstr (IMNAME(sh[i])) + call pargi (AP(sh[i])) + if (prncombine) { + call fprintf (logfd, " %6d") + call pargi (ncombine[i]) + } + if (prexptime) { + call fprintf (logfd, " %6.1f") + call pargr (exptime[i]) + } + if (prmode) { + call fprintf (logfd, " %7.5g") + call pargr (mode[i]) + } + if (prmedian) { + call fprintf (logfd, " %7.5g") + call pargr (median[i]) + } + if (prmean) { + call fprintf (logfd, " %7.5g") + call pargr (mean[i]) + } + if (prrdn) { + call fprintf (logfd, " %7g") + call pargr (RA(sh[i])) + } + if (prgain) { + call fprintf (logfd, " %6g") + call pargr (DEC(sh[i])) + } + if (prsn) { + call fprintf (logfd, " %6g") + call pargr (UT(sh[i])) + } + if (prscale) { + call fprintf (logfd, " %6.3f") + call pargr (1./scales[i]) + } + if (przero) { + call fprintf (logfd, " %7.5g") + call pargr (-zeros[i]) + } + if (prwts) { + call fprintf (logfd, " %6.3f") + call pargr (wts[i]) + } + call fprintf (logfd, "\n") + } + + # Log information about the output images. + call fprintf (logfd, "\n Output image = %s, ncombine = %d") + call pargstr (IMNAME(shout)) + call pargi (nout) + if (expname[1] != EOS) { + call fprintf (logfd, ", %s = %g") + call pargstr (expname) + call pargr (exposure) + } + call fprintf (logfd, "\n") + call fprintf (logfd, + " w1 = %g, w2 = %g, dw = %g, nw = %g, dtype = %d\n") + call pargr (W0(shout)) + call pargr (W1(shout)) + call pargr (WP(shout)) + call pargi (SN(shout)) + call pargi (DC(shout)) + + call flush (logfd) + call sfree (sp) +end diff --git a/noao/onedspec/scombine/icombine.com b/noao/onedspec/scombine/icombine.com new file mode 100644 index 00000000..771ada77 --- /dev/null +++ b/noao/onedspec/scombine/icombine.com @@ -0,0 +1,36 @@ +# SCOMBINE Common + +int combine # Combine algorithm +int reject # Rejection algorithm +real blank # Blank value +pointer rdnoise # CCD read noise +pointer gain # CCD gain +pointer snoise # CCD sensitivity noise +real lthresh # Low threshold +real hthresh # High threshold +int nkeep # Minimum to keep +real lsigma # Low sigma cutoff +real hsigma # High sigma cutoff +real pclip # Number or fraction of pixels from median +real flow # Fraction of low pixels to reject +real fhigh # Fraction of high pixels to reject +int grow # Grow radius +bool mclip # Use median in sigma clipping? +real sigscale # Sigma scaling tolerance +int logfd # Log file descriptor + +# These flags allow special conditions to be optimized. + +int dflag # Data flag (D_ALL, D_NONE, D_MIX) +bool doscale # Do the images have to be scaled? +bool doscale1 # Do the sigma calculations have to be scaled? +bool dothresh # Check pixels outside specified thresholds? +bool dowts # Does the final average have to be weighted? +bool keepids # Keep track of the image indices? +bool docombine # Call the combine procedure? +bool sort # Sort data? + +common /scbcom/ combine, reject, blank, rdnoise, gain, snoise, lsigma, hsigma, + lthresh, hthresh, nkeep, pclip, flow, fhigh, grow, logfd, + dflag, sigscale, mclip, doscale, doscale1, + dothresh, dowts, keepids, docombine, sort diff --git a/noao/onedspec/scombine/icombine.h b/noao/onedspec/scombine/icombine.h new file mode 100644 index 00000000..8a45a673 --- /dev/null +++ b/noao/onedspec/scombine/icombine.h @@ -0,0 +1,74 @@ +# SCOMBINE Definitions + +# Grouping options +define GROUP "|all|images|apertures|" +define GRP_ALL 1 +define GRP_IMAGES 2 +define GRP_APERTURES 3 + +# Sorting options +define SORT "|none|increasing|decreasing|" +define SORT_NONE 1 +define SORT_INC 2 +define SORT_DEC 3 + +# Combining modes in interactive mode +define CMB_AGAIN 0 +define CMB_ALL 1 +define CMB_FIRST 2 +define CMB_NEXT 3 +define CMB_SKIP 4 + +# Rejection options: +define REJECT "|none|ccdclip|crreject|minmax|pclip|sigclip|avsigclip|" +define NONE 1 # No rejection algorithm +define CCDCLIP 2 # CCD noise function clipping +define CRREJECT 3 # CCD noise function clipping +define MINMAX 4 # Minmax rejection +define PCLIP 5 # Percentile clip +define SIGCLIP 6 # Sigma clip +define AVSIGCLIP 7 # Sigma clip with average poisson sigma + +# Combine options: +define COMBINE "|average|median|sum|" +define AVERAGE 1 +define MEDIAN 2 +define SUM 3 + +# Scaling options: +define STYPES "|none|mode|median|mean|exposure|" +define ZTYPES "|none|mode|median|mean|" +define WTYPES "|none|mode|median|mean|exposure|" +define S_NONE 1 +define S_MODE 2 +define S_MEDIAN 3 +define S_MEAN 4 +define S_EXPOSURE 5 +define S_FILE 6 +define S_KEYWORD 7 +define S_SECTION "|input|output|overlap|" +define S_INPUT 1 +define S_OUTPUT 2 +define S_OVERLAP 3 + +# Data flag +define D_ALL 0 # All pixels are good +define D_NONE 1 # All pixels are bad or rejected +define D_MIX 2 # Mixture of good and bad pixels + +define TOL 0.001 # Tolerance for equal residuals + +# Spectrum data structure +define NS Memi[$1+$2-1] # Number of spec of given ap +define SH Memi[Memi[$1+$2-1]+$3-1] # Spectrum header structure + +# Combining options +#define COMBINE "|average|sum|" +#define CMB_AVERAGE 1 +#define CMB_SUM 2 + +# Weighting options +#define WT_TYPE "|none|expo|user|" +#define WT_NONE 1 +#define WT_EXPO 2 +#define WT_USER 3 diff --git a/noao/onedspec/scombine/icombine.x b/noao/onedspec/scombine/icombine.x new file mode 100644 index 00000000..5650d3ab --- /dev/null +++ b/noao/onedspec/scombine/icombine.x @@ -0,0 +1,174 @@ +include +include +include "icombine.h" + + +# IC_COMBINE -- Combine images. + +procedure ic_combiner (sh, shout, d, id, n, m, lflag, scales, zeros, wts, + nimages, npts) + +pointer sh[nimages] # Input spectra +pointer shout # Output spectrum +pointer d[nimages] # Data pointers +pointer id[nimages] # Image index ID pointers +int n[npts] # Number of good pixels +pointer m[nimages] # Mask pointers +int lflag[nimages] # Line flags +real scales[nimages] # Scale factors +real zeros[nimages] # Zero offset factors +real wts[nimages] # Combining weights +int nimages # Number of input images +int npts # Number of points per output line + +int i, ctor() +real r +pointer sp, nm +errchk ic_scale + +include "icombine.com" + +begin + call smark (sp) + + # Rebin spectra and set mask arrays + call scb_rebin (sh, shout, lflag, nimages, npts) + + # Set scale and weights and log + call ic_scale (sh, shout, lflag, scales, zeros, wts, nimages) + + # Set combine parameters + switch (combine) { + case AVERAGE: + if (dowts) + keepids = true + else + keepids = false + case MEDIAN: + dowts = false + keepids = false + case SUM: + keepids = false + reject = NONE + grow = 0 + } + docombine = true + + # Set rejection algorithm specific parameters + switch (reject) { + case CCDCLIP, CRREJECT: + call salloc (nm, 3*nimages, TY_REAL) + i = 1 + if (ctor (Memc[rdnoise], i, r) > 0) { + do i = 1, nimages + Memr[nm+3*(i-1)] = r + } else { + do i = 1, nimages + Memr[nm+3*(i-1)] = RA(sh[i]) + } + i = 1 + if (ctor (Memc[gain], i, r) > 0) { + do i = 1, nimages { + Memr[nm+3*(i-1)+1] = r + Memr[nm+3*(i-1)] = (Memr[nm+3*(i-1)] / r) ** 2 + } + } else { + do i = 1, nimages { + r = DEC(sh[i]) + Memr[nm+3*(i-1)+1] = r + Memr[nm+3*(i-1)] = (Memr[nm+3*(i-1)] / r) ** 2 + } + } + i = 1 + if (ctor (Memc[snoise], i, r) > 0) { + do i = 1, nimages + Memr[nm+3*(i-1)+2] = r + } else { + do i = 1, nimages { + r = UT(sh[i]) + Memr[nm+3*(i-1)+2] = r + } + } + if (!keepids) { + if (doscale1 || grow > 0) + keepids = true + else { + do i = 2, nimages { + if (Memr[nm+3*(i-1)] != Memr[nm] || + Memr[nm+3*(i-1)+1] != Memr[nm+1] || + Memr[nm+3*(i-1)+2] != Memr[nm+2]) { + keepids = true + break + } + } + } + } + if (reject == CRREJECT) + lsigma = MAX_REAL + case MINMAX: + mclip = false + if (grow > 0) + keepids = true + case PCLIP: + mclip = true + if (grow > 0) + keepids = true + case AVSIGCLIP, SIGCLIP: + if (doscale1 || grow > 0) + keepids = true + case NONE: + mclip = false + grow = 0 + } + + if (keepids) { + do i = 1, nimages + call salloc (id[i], npts, TY_INT) + } + + call ic_gdatar (sh, d, id, n, m, lflag, scales, zeros, nimages, npts) + + switch (reject) { + case CCDCLIP, CRREJECT: + if (mclip) + call ic_mccdclipr (d, id, n, scales, zeros, Memr[nm], + nimages, npts, Memr[SY(shout)]) + else + call ic_accdclipr (d, id, n, scales, zeros, Memr[nm], + nimages, npts, Memr[SY(shout)]) + case MINMAX: + call ic_mmr (d, id, n, npts) + case PCLIP: + call ic_pclipr (d, id, n, nimages, npts, Memr[SY(shout)]) + case SIGCLIP: + if (mclip) + call ic_msigclipr (d, id, n, scales, zeros, nimages, npts, + Memr[SY(shout)]) + else + call ic_asigclipr (d, id, n, scales, zeros, nimages, npts, + Memr[SY(shout)]) + case AVSIGCLIP: + if (mclip) + call ic_mavsigclipr (d, id, n, scales, zeros, nimages, + npts, Memr[SY(shout)]) + else + call ic_aavsigclipr (d, id, n, scales, zeros, nimages, + npts, Memr[SY(shout)]) + } + + if (grow > 0) + call ic_growr (d, id, n, nimages, npts, Memr[SY(shout)]) + + if (docombine) { + switch (combine) { + case AVERAGE: + call ic_averager (d, id, n, wts, npts, Memr[SY(shout)]) + case MEDIAN: + call ic_medianr (d, n, npts, Memr[SY(shout)]) + case SUM: + call ic_sumr (d, n, npts, Memr[SY(shout)]) + } + } + + call sfree (sp) +end diff --git a/noao/onedspec/scombine/icscale.x b/noao/onedspec/scombine/icscale.x new file mode 100644 index 00000000..009b30c3 --- /dev/null +++ b/noao/onedspec/scombine/icscale.x @@ -0,0 +1,463 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include +include +include +include +include +include "icombine.h" + +# IC_SCALE -- Get the scale factors for the spectra. +# 1. This procedure does CLIO to determine the type of scaling desired. +# 2. The output header parameters for exposure time and NCOMBINE are set. + +procedure ic_scale (sh, shout, lflags, scales, zeros, wts, nimages) + +pointer sh[nimages] # Input spectra +pointer shout # Output spectrum +int lflags[nimages] # Data flags +real scales[nimages] # Scale factors +real zeros[nimages] # Zero or sky levels +real wts[nimages] # Weights +int nimages # Number of images + +int stype, ztype, wtype +int i, j, nout +real mode, median, mean, exposure, zmean +pointer sp, ncombine, exptime, modes, medians, means, expname +pointer str, sname, zname, wname, rg +bool domode, domedian, domean, dozero + +int ic_gscale() +real asumr(), asumi() +pointer ic_wranges() +errchk ic_gscale, ic_statr + +include "icombine.com" + +begin + call smark (sp) + call salloc (ncombine, nimages, TY_INT) + call salloc (exptime, nimages, TY_REAL) + call salloc (modes, nimages, TY_REAL) + call salloc (medians, nimages, TY_REAL) + call salloc (means, nimages, TY_REAL) + call salloc (expname, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (sname, SZ_FNAME, TY_CHAR) + call salloc (zname, SZ_FNAME, TY_CHAR) + call salloc (wname, SZ_FNAME, TY_CHAR) + + # Set the defaults. + call amovki (1, Memi[ncombine], nimages) + call amovkr (0., Memr[exptime], nimages) + call amovkr (INDEF, Memr[modes], nimages) + call amovkr (INDEF, Memr[medians], nimages) + call amovkr (INDEF, Memr[means], nimages) + call amovkr (1., scales, nimages) + call amovkr (0., zeros, nimages) + call amovkr (1., wts, nimages) + + # Set scaling factors. + if (combine == SUM) { + stype = S_NONE + ztype = S_NONE + wtype = S_NONE + do i = 1, nimages + Memr[exptime+i-1] = IT(sh[i]) + } else { + stype = ic_gscale ("scale", Memc[sname], STYPES, sh, Memr[exptime], + scales, nimages) + ztype = ic_gscale ("zero", Memc[zname], ZTYPES, sh, Memr[exptime], + zeros, nimages) + wtype = ic_gscale ("weight", Memc[wname], WTYPES, sh, Memr[exptime], + wts, nimages) + } + + Memc[expname] = EOS + if (combine == SUM || stype == S_EXPOSURE || wtype == S_EXPOSURE) { + call strcpy ("exptime", Memc[expname], SZ_FNAME) + do i = 1, nimages + if (IS_INDEFR(Memr[exptime+i-1])) + Memc[expname] = EOS + } + + # Get image statistics only if needed. + domode = ((stype==S_MODE)||(ztype==S_MODE)||(wtype==S_MODE)) + domedian = ((stype==S_MEDIAN)||(ztype==S_MEDIAN)||(wtype==S_MEDIAN)) + domean = ((stype==S_MEAN)||(ztype==S_MEAN)||(wtype==S_MEAN)) + if (domode || domedian || domean) { + call clgstr ("sample", Memc[str], SZ_LINE) + rg = ic_wranges (Memc[str]) + do i = 1, nimages { + call ic_statr (sh[i], lflags[i], rg, domode, domedian, domean, + mode, median, mean) + if (domode) { + Memr[modes+i-1] = mode + if (stype == S_MODE) + scales[i] = mode + if (ztype == S_MODE) + zeros[i] = mode + if (wtype == S_MODE) + wts[i] = mode + } + if (domedian) { + Memr[medians+i-1] = median + if (stype == S_MEDIAN) + scales[i] = median + if (ztype == S_MEDIAN) + zeros[i] = median + if (wtype == S_MEDIAN) + wts[i] = median + } + if (domean) { + Memr[means+i-1] = mean + if (stype == S_MEAN) + scales[i] = mean + if (ztype == S_MEAN) + zeros[i] = mean + if (wtype == S_MEAN) + wts[i] = mean + } + } + call mfree (rg, TY_REAL) + } + + do i = 1, nimages + if (scales[i] <= 0.) { + call eprintf ("WARNING: Negative scale factors") + call eprintf (" -- ignoring scaling\n") + call amovkr (1., scales, nimages) + break + } + + # Convert to relative factors. + mean = asumr (scales, nimages) / nimages + call adivkr (scales, mean, scales, nimages) + call adivr (zeros, scales, zeros, nimages) + zmean = asumr (zeros, nimages) / nimages + + if (wtype != S_NONE) { + do i = 1, nimages { + if (wts[i] <= 0.) { + call eprintf ("WARNING: Negative weights") + call eprintf (" -- using only NCOMBINE weights\n") + do j = 1, nimages + wts[j] = Memi[ncombine+j-1] + break + } + if (ztype == S_NONE) + wts[i] = Memi[ncombine+i-1] * wts[i] + else { + if (zeros[i] <= 0.) { + call eprintf ("WARNING: Negative zero offsets") + call eprintf (" -- ignoring zero weight adjustments\n") + do j = 1, nimages + wts[j] = Memi[ncombine+j-1] * wts[j] + break + } + wts[i] = Memi[ncombine+i-1] * wts[i] * zmean / zeros[i] + } + } + } + + call asubkr (zeros, zmean, zeros, nimages) + mean = asumr (wts, nimages) + call adivkr (wts, mean, wts, nimages) + + # Because of finite arithmetic it is possible for the zero offsets to + # be nonzero even when they are all equal. Just for the sake of + # a nice log set the zero offsets in this case. + + for (i=2; (i<=nimages)&&(zeros[i]==zeros[1]); i=i+1) + ; + if (i > nimages) + call aclrr (zeros, nimages) + + # Set flags for scaling, zero offsets, sigma scaling, weights. + # Sigma scaling may be suppressed if the scales or zeros are + # different by a specified tolerance. + + doscale = false + dozero = false + doscale1 = false + dowts = false + do i = 2, nimages { + if (scales[i] != scales[1]) + doscale = true + if (zeros[i] != zeros[1]) + dozero = true + if (wts[i] != wts[1]) + dowts = true + } + if (doscale && sigscale != 0.) { + do i = 1, nimages { + if (abs (scales[i] - 1) > sigscale) { + doscale1 = true + break + } + } + if (!doscale1 && zmean > 0.) { + do i = 1, nimages { + if (abs (zeros[i] / zmean) > sigscale) { + doscale1 = true + break + } + } + } + } + + # Set the output header parameters. + nout = asumi (Memi[ncombine], nimages) + call imaddi (IM(shout), "ncombine", nout) + if (Memc[expname] != EOS) { + exposure = 0. + if (combine == SUM) { + do i = 1, nimages + exposure = exposure + Memr[exptime+i-1] + } else { + do i = 1, nimages + exposure = exposure + wts[i] * Memr[exptime+i-1] / scales[i] + } + call imaddr (IM(shout), Memc[expname], exposure) + } else + exposure = INDEF + + # Start the log here since much of the info is only available here. + call ic_log (sh, shout, Memi[ncombine], Memr[exptime], Memc[sname], + Memc[zname], Memc[wname], Memr[modes], Memr[medians], Memr[means], + scales, zeros, wts, nimages, dozero, nout, Memc[expname], exposure) + + doscale = (doscale || dozero) + + call sfree (sp) +end + + +# IC_GSCALE -- Get scale values as directed by CL parameter +# The values can be one of those in the dictionary, from a file specified +# with a @ prefix, or from an image header keyword specified by a ! prefix. + +int procedure ic_gscale (param, name, dic, sh, exptime, values, nimages) + +char param[ARB] #I CL parameter name +char name[SZ_FNAME] #O Parameter value +char dic[ARB] #I Dictionary string +pointer sh[nimages] #I SHDR pointers +real exptime[nimages] #I Exposure times +real values[nimages] #O Values +int nimages #I Number of images + +int type #O Type of value + +int fd, i, nowhite(), open(), fscan(), nscan(), strdic() +real rval +pointer errstr +errchk open + +include "icombine.com" + +begin + call clgstr (param, name, SZ_FNAME) + if (nowhite (name, name, SZ_FNAME) == 0) + type = S_NONE + else if (name[1] == '@') { + type = S_FILE + fd = open (name[2], READ_ONLY, TEXT_FILE) + i = 0 + while (fscan (fd) != EOF) { + call gargr (rval) + if (nscan() != 1) + next + if (i == nimages) { + call eprintf ( + "Warning: Ignoring additional %s values in %s\n") + call pargstr (param) + call pargstr (name[2]) + break + } + i = i + 1 + values[i] = rval + } + call close (fd) + if (i < nimages) { + call salloc (errstr, SZ_LINE, TY_CHAR) + call sprintf (Memc[errstr], SZ_FNAME, + "Insufficient %s values in %s") + call pargstr (param) + call pargstr (name[2]) + call error (1, Memc[errstr]) + } + } else if (name[1] == '!') { + type = S_KEYWORD + do i = 1, nimages { + switch (param[1]) { + case 's': + values[i] = ST(sh[i]) + case 'z': + values[i] = HA(sh[i]) + case 'w': + values[i] = AM(sh[i]) + } + } + } else { + type = strdic (name, name, SZ_FNAME, dic) + if (type == 0) + call error (1, "Unknown scale, zero, or weight type") + if (type==S_EXPOSURE) { + do i = 1, nimages { + if (IS_INDEF(IT(sh[i]))) + call error (1, "Exposure time not defined") + exptime[i] = IT(sh[i]) + values[i] = max (0.001, exptime[i]) + } + } + } + + return (type) +end + + +# IC_WRANGES -- Parse wavelength range string. +# A wavelength range string consists of colon delimited ranges with +# multiple ranges separated by comma and/or whitespace. + +pointer procedure ic_wranges (rstr) + +char rstr[ARB] # Range string +pointer rg # Range pointer + +int i, fd, strlen(), open(), getline() +pointer sp, str, ptr +errchk open, ic_wadd + +begin + call smark (sp) + call salloc (str, max (strlen (rstr), SZ_LINE), TY_CHAR) + call calloc (rg, 1, TY_REAL) + + i = 1 + while (rstr[i] != EOS) { + + # Find beginning and end of a range and copy it to the work string + while (IS_WHITE(rstr[i]) || rstr[i]==',' || rstr[i]=='\n') + i = i + 1 + if (rstr[i] == EOS) + break + + ptr = str + while (!(IS_WHITE(rstr[i]) || rstr[i]==',' || rstr[i]=='\n' || + rstr[i]==EOS)) { + Memc[ptr] = rstr[i] + i = i + 1 + ptr = ptr + 1 + } + Memc[ptr] = EOS + + # Add range(s) + iferr { + if (Memc[str] == '@') { + fd = open (Memc[str+1], READ_ONLY, TEXT_FILE) + while (getline (fd, Memc[str]) != EOF) { + iferr (call ic_wadd (rg, Memc[str])) + call erract (EA_WARN) + } + call close (fd) + } else + call ic_wadd (rg, Memc[str]) + } then + call erract (EA_WARN) + } + + call sfree (sp) + + # Set final structure + i = Memr[rg] + if (i == 0) + call mfree (rg, TY_REAL) + else + call realloc (rg, 1 + 2 * i, TY_REAL) + return (rg) +end + + +# IC_WADD -- Add a range + +procedure ic_wadd (rg, rstr) + +pointer rg # Range descriptor +char rstr[ARB] # Range string + +int i, j, n, strlen(), ctor() +real w1, w2 +pointer sp, str, ptr + +begin + call smark (sp) + call salloc (str, strlen (rstr), TY_CHAR) + + i = 1 + n = Memr[rg] + while (rstr[i] != EOS) { + + # Find beginning and end of a range and copy it to the work string + while (IS_WHITE(rstr[i]) || rstr[i]==',' || rstr[i]=='\n') + i = i + 1 + if (rstr[i] == EOS) + break + + ptr = str + while (!(IS_WHITE(rstr[i]) || rstr[i]==',' || rstr[i]=='\n' || + rstr[i]==EOS)) { + if (rstr[i] == ':') + Memc[ptr] = ' ' + else + Memc[ptr] = rstr[i] + i = i + 1 + ptr = ptr + 1 + } + Memc[ptr] = EOS + + # Parse range + if (Memc[str] == '@') + call error (1, "Cannot nest @files") + else { + # Get range + j = 1 + if (ctor (Memc[str], j, w1) == 0) + call error (1, "Range syntax error") + if (ctor (Memc[str], j, w2) == 0) + call error (1, "Range syntax error") + } + + if (mod (n, 10) == 0) + call realloc (rg, 1+2*(n+10), TY_REAL) + n = n + 1 + Memr[rg+2*n-1] = min (w1, w2) + Memr[rg+2*n] = max (w1, w2) + } + Memr[rg] = n + + call sfree (sp) +end + + +# IC_WISINRANGE -- Is wavelength in range? + +bool procedure ic_wisinrange (rg, w) + +pointer rg # Wavelength range array +real w # Wavelength + +int i, n + +begin + if (rg == NULL) + return (true) + + n = nint (Memr[rg]) + do i = 1, 2*n, 2 + if (w >= Memr[rg+i] && w <= Memr[rg+i+1]) + return (true) + return (false) +end diff --git a/noao/onedspec/scombine/icstat.x b/noao/onedspec/scombine/icstat.x new file mode 100644 index 00000000..3fce4165 --- /dev/null +++ b/noao/onedspec/scombine/icstat.x @@ -0,0 +1,160 @@ +include +include "icombine.h" + + +# IC_STATR -- Compute image statistics within spectrum. + +procedure ic_statr (sh, lflag, rg, domode, domedian, domean, mode, median, mean) + +pointer sh # Spectrum structure +int lflag # Data flag +pointer rg # Wavelength ranges +bool domode, domedian, domean # Statistics to compute +real mode, median, mean # Statistics + +int i, n, npts +real a, w +pointer sp, data, dp, lp, mp +real ic_moder(), asumr() +bool ic_wisinrange() +double shdr_lw() + +include "icombine.com" + +begin + mp = SX(sh) + lp = SY(sh) + npts = SN(sh) + + call smark (sp) + call salloc (data, npts, TY_REAL) + + dp = data + if (lflag == D_ALL && rg == NULL) { + if (dothresh) { + do i = 1, npts { + a = Memr[lp] + if (a >= lthresh && a <= hthresh) { + Memr[dp] = a + dp = dp + 1 + } + lp = lp + 1 + } + } else { + do i = 1, npts { + Memr[dp] = Memr[lp] + dp = dp + 1 + lp = lp + 1 + } + } + } else if (lflag == D_MIX || rg != NULL) { + if (dothresh) { + do i = 1, npts { + if (Memr[mp] == 0) { + a = Memr[lp] + if (a >= lthresh && a <= hthresh) { + w = shdr_lw (sh, double (i)) + if (ic_wisinrange (rg, w)) { + Memr[dp] = a + dp = dp + 1 + } + } + } + mp = mp + 1 + lp = lp + 1 + } + } else { + do i = 1, npts { + if (Memr[mp] == 0) { + w = shdr_lw (sh, double (i)) + if (ic_wisinrange (rg, w)) { + Memr[dp] = Memr[lp] + dp = dp + 1 + } + } + mp = mp + 1 + lp = lp + 1 + } + } + } + + n = dp - data + if (n > 0) { + # Compute only statistics needed. + if (domode || domedian) { + call asrtr (Memr[data], Memr[data], n) + mode = ic_moder (Memr[data], n) + median = Memr[data+n/2-1] + } + if (domean) + mean = asumr (Memr[data], n) / n + } else { + mode = INDEF + median = INDEF + mean = INDEF + } + + call sfree (sp) +end + + +define NMIN 10 # Minimum number of pixels for mode calculation +define ZRANGE 0.8 # Fraction of pixels about median to use +define ZSTEP 0.01 # Step size for search for mode +define ZBIN 0.1 # Bin size for mode. + +# IC_MODE -- Compute mode of an array. The mode is found by binning +# with a bin size based on the data range over a fraction of the +# pixels about the median and a bin step which may be smaller than the +# bin size. If there are too few points the median is returned. +# The input array must be sorted. + +real procedure ic_moder (a, n) + +real a[n] # Data array +int n # Number of points + +int i, j, k, nmax +real z1, z2, zstep, zbin +real mode +bool fp_equalr() + +begin + if (n < NMIN) + return (a[n/2]) + + # Compute the mode. The array must be sorted. Consider a + # range of values about the median point. Use a bin size which + # is ZBIN of the range. Step the bin limits in ZSTEP fraction of + # the bin size. + + i = 1 + n * (1. - ZRANGE) / 2. + j = 1 + n * (1. + ZRANGE) / 2. + z1 = a[i] + z2 = a[j] + if (fp_equalr (z1, z2)) { + mode = z1 + return (mode) + } + + zstep = ZSTEP * (z2 - z1) + zbin = ZBIN * (z2 - z1) + + z1 = z1 - zstep + k = i + nmax = 0 + repeat { + z1 = z1 + zstep + z2 = z1 + zbin + for (; i < j && a[i] < z1; i=i+1) + ; + for (; k < j && a[k] < z2; k=k+1) + ; + if (k - i > nmax) { + nmax = k - i + mode = a[(i+k)/2] + } + } until (k >= j) + + return (mode) +end diff --git a/noao/onedspec/scombine/icsum.x b/noao/onedspec/scombine/icsum.x new file mode 100644 index 00000000..f038b37b --- /dev/null +++ b/noao/onedspec/scombine/icsum.x @@ -0,0 +1,48 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "icombine.h" + + +# IC_SUM -- Compute the summed image line. + +procedure ic_sumr (d, n, npts, sum) + +pointer d[ARB] # Data pointers +int n[npts] # Number of points +int npts # Number of output points per line +real sum[npts] # Average (returned) + +int i, j, k +real s + +include "icombine.com" + +begin + # If no data has been excluded do the sum without checking the + # number of points. If all the data has been excluded set the + # sum to the blank value. + + if (dflag == D_ALL) { + do i = 1, npts { + k = i - 1 + s = Memr[d[1]+k] + do j = 2, n[i] + s = s + Memr[d[j]+k] + sum[i] = s + } + } else if (dflag == D_NONE) { + do i = 1, npts + sum[i] = blank + } else { + do i = 1, npts { + if (n[i] > 0) { + k = i - 1 + s = Memr[d[1]+k] + do j = 2, n[i] + s = s + Memr[d[j]+k] + sum[i] = s + } else + sum[i] = blank + } + } +end diff --git a/noao/onedspec/scombine/iscombine.key b/noao/onedspec/scombine/iscombine.key new file mode 100644 index 00000000..91d7876b --- /dev/null +++ b/noao/onedspec/scombine/iscombine.key @@ -0,0 +1,23 @@ + SCOMBINE CURSOR KEYS + +a - Mark scaling ranges in overlap region for 'v' key +b - Cancel scaling ranges +c - Print cursor position +d - Replace range of pixels by straight cursor line +e - Replace range of pixels by linear interpolation from the endpoint pixels +f - Start over from the first spectrum +j - Fudge a point to vertical cursor value +n - Go on to next spectrum +o - Reset data for current spectrum to initial values +p - Don't include current spectrum in combined image and go on to next spectrum +q - Quit and combine remaining spectra noninteractively (no|yes|YES) +s - Mark accumulation ranges in current spectrum +t - Cancel accumulation ranges +v - Shift overlap average of spectrum vertically to accumulated spectrum +w - Window the graph using gtools commands +x - Shift spectrum horizontally to cursor position +y - Shift spectrum vertically to cursor position +z - Shift spectrum vertically to accumulated spectrum ++ - Set additive scaling for 'v' key +* - Set multiplicative scaling +? - This help page diff --git a/noao/onedspec/scombine/iscombine.par b/noao/onedspec/scombine/iscombine.par new file mode 100644 index 00000000..a9d8846e --- /dev/null +++ b/noao/onedspec/scombine/iscombine.par @@ -0,0 +1,18 @@ +input,s,a,"",,,List of input spectra +output,s,a,"",,,List of output spectra +woutput,s,h,"",,,List of output weight spectra +apertures,s,h,"",,,Apertures to combine +group,s,h,"apertures","all|images|apertures",,Grouping option +combine,s,h,"average","average|sum",,Combining option +scale,s,h,"",,,Header keyword for scaling +weight,s,h,"","",,"Header keyword for weighting +" +w1,r,h,INDEF,,,Starting wavelength of output spectra +w2,r,h,INDEF,,,Ending wavelength of output spectra +dw,r,h,INDEF,,,Wavelength increment of output spectra +nw,i,h,INDEF,,,Length of output spectra +log,b,h,no,,,"Logarithmic increments? +" +interactive,b,h,no,,,Adjust spectra interactively? +sort,s,h,"none","none|increasing|decreasing",,Interactive combining order +cursor,*gcur,h,"",,,Graphics cursor input diff --git a/noao/onedspec/scombine/mkpkg b/noao/onedspec/scombine/mkpkg new file mode 100644 index 00000000..ab60e45b --- /dev/null +++ b/noao/onedspec/scombine/mkpkg @@ -0,0 +1,35 @@ +# SCOMBINE + +$call relink +$exit + +update: + $call relink + $call install + ; + +relink: + $update libpkg.a + $omake x_scombine.x + $link x_scombine.o libpkg.a -lsmw -lxtools -liminterp \ + -o xx_scombine.e + ; + +install: + $move xx_scombine.e noaobin$x_scombine.e + ; + + +libpkg.a: + @generic + + icgdata.x icombine.com icombine.h + iclog.x icombine.com icombine.h + icombine.x icombine.com icombine.h + icscale.x icombine.com icombine.h \ + + icstat.x icombine.com icombine.h + icsum.x icombine.com icombine.h + t_scombine.x icombine.h icombine.com \ + + ; diff --git a/noao/onedspec/scombine/scombine.par b/noao/onedspec/scombine/scombine.par new file mode 100644 index 00000000..932e6e31 --- /dev/null +++ b/noao/onedspec/scombine/scombine.par @@ -0,0 +1,37 @@ +input,s,a,"",,,List of input spectra +output,s,a,"",,,List of output spectra +noutput,s,h,"",,,List of output number combined spectra +logfile,s,h,"STDOUT",,,"Log file +" +apertures,s,h,"",,,Apertures to combine +group,s,h,"apertures","all|images|apertures",,"Grouping option" +combine,s,h,"average","average|median|sum",,Type of combine operation +reject,s,h,"none","none|minmax|ccdclip|crreject|sigclip|avsigclip|pclip",,"Type of rejection +" +first,b,h,no,,,Use first spectrum for dispersion? +w1,r,h,INDEF,,,Starting wavelength of output spectra +w2,r,h,INDEF,,,Ending wavelength of output spectra +dw,r,h,INDEF,,,Wavelength increment of output spectra +nw,i,h,INDEF,,,Length of output spectra +log,b,h,no,,,"Logarithmic increments? +" +scale,s,h,"none",,,Image scaling +zero,s,h,"none",,,Image zero point offset +weight,s,h,"none",,,Image weights +sample,s,h,"",,,"Wavelength sample regions for statistics +" +lthreshold,r,h,INDEF,,,Lower threshold +hthreshold,r,h,INDEF,,,Upper threshold +nlow,i,h,1,0,,minmax: Number of low pixels to reject +nhigh,i,h,1,0,,minmax: Number of high pixels to reject +nkeep,i,h,1,,,Minimum to keep (pos) or maximum to reject (neg) +mclip,b,h,yes,,,Use median in sigma clipping algorithms? +lsigma,r,h,3.,0.,,Lower sigma clipping factor +hsigma,r,h,3.,0.,,Upper sigma clipping factor +rdnoise,s,h,"0.",,,ccdclip: CCD readout noise (electrons) +gain,s,h,"1.",,,ccdclip: CCD gain (electrons/DN) +snoise,s,h,"0.",,,ccdclip: Sensitivity noise (fraction) +sigscale,r,h,0.1,0.,,Tolerance for sigma clipping scaling corrections +pclip,r,h,-0.5,,,pclip: Percentile clipping parameter +grow,i,h,0,,,Radius (pixels) for 1D neighbor rejection +blank,r,h,0.,,,Value if there are no pixels diff --git a/noao/onedspec/scombine/t_scombine.x b/noao/onedspec/scombine/t_scombine.x new file mode 100644 index 00000000..774a5f87 --- /dev/null +++ b/noao/onedspec/scombine/t_scombine.x @@ -0,0 +1,630 @@ +include +include +include +include +include "icombine.h" + + +# T_SCOMBINE - Combine spectra +# The input spectra are combined by medianing, averaging or summing +# with optional rejection, scaling and weighting. The input may be +# grouped by aperture or by image. The combining algorithms are +# similar to those in IMCOMBINE. + +procedure t_scombine() + +int ilist # list of input images +int olist # list of output images +pointer nlist # image name for number combined +pointer aps # aperture ranges +int group # grouping option + +int reject1 +real flow1, fhigh1, pclip1, nkeep1 + +real rval +bool grdn, ggain, gsn +int i, j, k, l, n, naps, npts +pointer im, mw, nout, refim, shin, shout +pointer sp, input, output, noutput, scale, zero, weight, str, logfile, sh, ns +pointer sp1, d, id, nc, m, lflag, scales, zeros, wts + +real clgetr(), imgetr() +bool clgetb(), rng_elementi() +int clgeti(), clgwrd(), ctor() +int imtopenp(), imtgetim(), open(), nowhite() +pointer rng_open(), immap(), smw_openim(), impl2i(), impl2r() +errchk open, immap, smw_openim, shdr_open, imgetr +errchk scb_output, scb_combine, ic_combiner + +include "icombine.com" + +begin + call smark (sp) + call salloc (input, SZ_FNAME, TY_CHAR) + call salloc (output, SZ_FNAME, TY_CHAR) + call salloc (noutput, SZ_FNAME, TY_CHAR) + call salloc (scale, SZ_FNAME, TY_CHAR) + call salloc (zero, SZ_FNAME, TY_CHAR) + call salloc (weight, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (gain, SZ_FNAME, TY_CHAR) + call salloc (snoise, SZ_FNAME, TY_CHAR) + call salloc (rdnoise, SZ_FNAME, TY_CHAR) + call salloc (logfile, SZ_FNAME, TY_CHAR) + + # Get parameters + ilist = imtopenp ("input") + olist = imtopenp ("output") + nlist = imtopenp ("noutput") + call clgstr ("apertures", Memc[str], SZ_LINE) + group = clgwrd ("group", Memc[input], SZ_FNAME, GROUP) + + # IMCOMBINE parameters + call clgstr ("logfile", Memc[logfile], SZ_FNAME) + combine = clgwrd ("combine", Memc[input], SZ_FNAME, COMBINE) + reject1 = clgwrd ("reject", Memc[input], SZ_FNAME, REJECT) + blank = clgetr ("blank") + call clgstr ("scale", Memc[scale], SZ_FNAME) + call clgstr ("zero", Memc[zero], SZ_FNAME) + call clgstr ("weight", Memc[weight], SZ_FNAME) + call clgstr ("gain", Memc[gain], SZ_FNAME) + call clgstr ("rdnoise", Memc[rdnoise], SZ_FNAME) + call clgstr ("snoise", Memc[snoise], SZ_FNAME) + lthresh = clgetr ("lthreshold") + hthresh = clgetr ("hthreshold") + lsigma = clgetr ("lsigma") + hsigma = clgetr ("hsigma") + pclip1 = clgetr ("pclip") + flow1 = clgetr ("nlow") + fhigh1 = clgetr ("nhigh") + nkeep1 = clgeti ("nkeep") + grow = clgeti ("grow") + mclip = clgetb ("mclip") + sigscale = clgetr ("sigscale") + + i = nowhite (Memc[scale], Memc[scale], SZ_FNAME) + i = nowhite (Memc[zero], Memc[zero], SZ_FNAME) + i = nowhite (Memc[weight], Memc[weight], SZ_FNAME) + + # Check parameters, map INDEFs, and set threshold flag + if (combine == SUM) + reject1 = NONE + if (pclip1 == 0. && reject1 == PCLIP) + call error (1, "Pclip parameter may not be zero") + if (IS_INDEFR (blank)) + blank = 0. + if (IS_INDEFR (lsigma)) + lsigma = MAX_REAL + if (IS_INDEFR (hsigma)) + hsigma = MAX_REAL + if (IS_INDEFR (pclip1)) + pclip1 = -0.5 + if (IS_INDEFI (nkeep1)) + nkeep1 = 0 + if (IS_INDEFR (flow1)) + flow1 = 0 + if (IS_INDEFR (fhigh)) + fhigh = 0 + if (IS_INDEFI (grow)) + grow = 0 + if (IS_INDEF (sigscale)) + sigscale = 0. + + if (IS_INDEF(lthresh) && IS_INDEF(hthresh)) + dothresh = false + else { + dothresh = true + if (IS_INDEF(lthresh)) + lthresh = -MAX_REAL + if (IS_INDEF(hthresh)) + hthresh = MAX_REAL + } + + # Get read noise and gain? + grdn = false + ggain = false + gsn = false + if (reject1 == CCDCLIP || reject1 == CRREJECT) { + i = 1 + if (ctor (Memc[rdnoise], i, rval) == 0) + grdn = true + i = 1 + if (ctor (Memc[gain], i, rval) == 0) + ggain = true + i = 1 + if (ctor (Memc[snoise], i, rval) == 0) + gsn = true + } + + # Open the log file. + logfd = NULL + if (Memc[logfile] != EOS) { + iferr (logfd = open (Memc[logfile], APPEND, TEXT_FILE)) { + logfd = NULL + call erract (EA_WARN) + } + } + + iferr (aps = rng_open (Memc[str], INDEF, INDEF, INDEF)) + call error (1, "Error in aperture list") + + # Loop through input images. + while (imtgetim (ilist, Memc[input], SZ_FNAME) != EOF) { + if (imtgetim (olist, Memc[output], SZ_FNAME) == EOF) { + call eprintf ("No output image\n") + break + } + if (imtgetim (nlist, Memc[noutput], SZ_FNAME) == EOF) + Memc[noutput] = EOS + + # Get spectra to combine. + # Because the input images are unmapped we must get all the + # data we need for combining into the spectrum data structures. + # In particular any header keyword parameters that will be + # used. We save the header values in unused elements of + # the spectrum data structure. + + naps = 0 + repeat { + iferr (im = immap (Memc[input], READ_ONLY, 0)) { + if (group == GRP_IMAGES) { + call erract (EA_WARN) + next + } else { + call erract (EA_ERROR) + } + } + mw = smw_openim (im) + shin = NULL + + do i = 1, SMW_NSPEC(mw) { + call shdr_open (im, mw, i, 1, INDEFI, SHDATA, shin) + if (Memc[scale] == '!') + ST(shin) = imgetr (im, Memc[scale+1]) + if (Memc[zero] == '!') + HA(shin) = imgetr (im, Memc[zero+1]) + if (Memc[weight] == '!') + AM(shin) = imgetr (im, Memc[weight+1]) + if (grdn) + RA(shin) = imgetr (im, Memc[rdnoise]) + if (ggain) + DEC(shin) = imgetr (im, Memc[gain]) + if (gsn) + UT(shin) = imgetr (im, Memc[snoise]) + if (!rng_elementi (aps, AP(shin))) + next + if (group == GRP_APERTURES) { + for (j=1; j<=naps; j=j+1) + if (AP(shin) == AP(SH(sh,j,1))) + break + n = 10 + } else { + j = 1 + n = 1 + } + + if (naps == 0) { + call calloc (sh, n, TY_POINTER) + call calloc (ns, n, TY_INT) + } else if (j > naps && mod (naps, n) == 0) { + call realloc (sh, naps+n, TY_POINTER) + call realloc (ns, naps+n, TY_INT) + call aclri (Memi[sh+naps], n) + call aclri (Memi[ns+naps], n) + } + if (j > naps) + naps = naps + 1 + n = NS(ns,j) + if (n == 0) + call malloc (Memi[sh+j-1], 10, TY_POINTER) + else if (mod (n, 10) == 0) + call realloc (Memi[sh+j-1], n+10, TY_POINTER) + + n = n + 1 + SH(sh,j,n) = NULL + NS(ns,j) = n + call shdr_copy (shin, SH(sh,j,n), NO) + } + + call imunmap (IM(shin)) + MW(shin) = NULL + call shdr_close (shin) + + if (group == GRP_IMAGES) + break + } until (imtgetim (ilist, Memc[input], SZ_FNAME) == EOF) + + if (naps < 1) { + call eprintf ("No input spectra to combine\n") + next + } + + # Set the output and combine the spectra. + call scb_output (sh, ns, naps, Memc[output], Memc[noutput], + im, mw, nout, refim) + + do j = 1, naps { + call shdr_open (im, mw, j, 1, INDEFI, SHHDR, shout) + npts = SN(shout) + n = NS(ns,j) + + # Allocate additional memory + call smark (sp1) + call salloc (d, n, TY_POINTER) + call salloc (id, n, TY_POINTER) + call salloc (nc, npts, TY_INT) + call salloc (m, n, TY_POINTER) + call salloc (lflag, n, TY_INT) + call salloc (scales, n, TY_REAL) + call salloc (zeros, n, TY_REAL) + call salloc (wts, n, TY_REAL) + call calloc (SX(shout), npts, TY_REAL) + call calloc (SY(shout), npts, TY_REAL) + call amovki (D_ALL, Memi[lflag], n) + + # Convert the pclip parameter to a number of pixels rather than + # a fraction. This number stays constant even if pixels are + # rejected. The number of low and high pixel rejected, however, + # are converted to a fraction of the valid pixels. + + reject = reject1 + nkeep = nkeep1 + if (nkeep < 0) + nkeep = n + nkeep + if (reject == PCLIP) { + pclip = pclip1 + i = (n - 1) / 2. + if (abs (pclip) < 1.) + pclip = pclip * i + if (pclip < 0.) + pclip = min (-1, max (-i, int (pclip))) + else + pclip = max (1, min (i, int (pclip))) + } + if (reject == MINMAX) { + flow = flow1 + fhigh = fhigh1 + if (flow >= 1) + flow = flow / n + if (fhigh >= 1) + fhigh = fhigh / n + i = flow * n + fhigh * n + if (i == 0) + reject = NONE + else if (i >= n) { + call eprintf ("Bad minmax rejection parameters\n") + call eprintf ("Using no rejection\n") + reject = NONE + } + } + + # Combine spectra + call ic_combiner (SH(sh,j,1), shout, Memi[d], Memi[id], + Memi[nc], Memi[m], Memi[lflag], Memr[scales], Memr[zeros], + Memr[wts], n, npts) + + # Write the results + call amovr (Memr[SY(shout)], Memr[impl2r(im,j)], npts) + if (nout != NULL) + call amovi (Memi[nc], Memi[impl2i(nout,j)], npts) + call sfree (sp1) + } + + # Finish up + call shdr_close (shout) + call smw_close (mw) + call imunmap (im) + call imunmap (refim) + if (nout != NULL) + call imunmap (nout) + + # Find all the distinct SMW pointers and free them. + do j = 1, naps { + do i = 1, NS(ns,j) { + mw = MW(SH(sh,j,i)) + if (mw != NULL) { + do k = 1, naps { + do l = 1, NS(ns,k) { + shin = SH(sh,k,l) + if (MW(shin) == mw) + MW(shin) = NULL + } + } + call smw_close (mw) + } + } + } + do j = 1, naps { + do i = 1, NS(ns,j) + call shdr_close (SH(sh,j,i)) + call mfree (Memi[sh+j-1], TY_POINTER) + } + call mfree (sh, TY_POINTER) + call mfree (ns, TY_INT) + } + + call rng_close (aps) + call imtclose (ilist) + call imtclose (olist) + call imtclose (nlist) + + call sfree (sp) +end + + +# SCB_REBIN - Rebin input spectra to output dispersion +# Use the SX array as mask. If less than 1% of an input +# pixel contributes to an output pixel then flag it as missing data. + +procedure scb_rebin (sh, shout, lflag, ns, npts) + +pointer sh[ns] # Input spectra structures +pointer shout # Output spectrum structure +int lflag[ns] # Empty mask flags +int ns # Number of spectra +int npts # NUmber of output points + +int i, j +real a, b, c +pointer shin +double shdr_wl(), shdr_lw() + +include "icombine.com" + +begin + # Rebin to common dispersion + # Determine overlap with output and set mask arrays + + do i = 1, ns { + shin = sh[i] + c = shdr_wl (shout, shdr_lw (shin, double(0.5))) + b = shdr_wl (shout, shdr_lw (shin, double(SN(shin)+0.5))) + a = max (1, nint (min (b, c) + 0.01)) + b = min (npts, nint (max (b, c) - 0.01)) + j = b - a + 1 + if (j < 1) { + lflag[i] = D_NONE + next + } + else if (j < npts) + lflag[i] = D_MIX + else + lflag[i] = D_ALL + + call shdr_rebin (shin, shout) + call aclrr (Memr[SX(shin)], SN(shin)) + j = a - 1 + if (j > 0) + call amovkr (1.0, Memr[SX(shin)], j) + j = SN(shin) - b + if (j > 0) + call amovkr (1.0, Memr[SX(shin)+SN(shin)-j], j) + } + + dflag = lflag[1] + do i = 2, ns { + if (dflag != lflag[i]) { + dflag = D_MIX + break + } + } +end + + +# SCB_OUTPUT - Set the output spectrum + +procedure scb_output (sh, ns, naps, output, noutput, im, mw, nout, refim) + +pointer sh # spectra structures +int ns # number of spectra +int naps # number of apertures +char output[SZ_FNAME] # output spectrum name +char noutput[SZ_FNAME] # output number combined image name +pointer im # output IMIO pointer +pointer mw # output MWCS pointer +pointer nout # output number combined IMIO pointer +pointer refim # reference image for output image + +int i, ap, beam, dtype, nw, nmax, axis[2] +double w1, dw, z +real aplow[2], aphigh[2] +pointer sp, key, coeff, sh1 +pointer immap(), mw_open(), smw_openim() +errchk immap, smw_openim +data axis/1,2/ + +begin + call smark (sp) + call salloc (key, SZ_FNAME, TY_CHAR) + coeff = NULL + + # Create output image using the first input image as a reference + refim = immap (IMNAME(SH(sh,1,1)), READ_ONLY, 0) + im = immap (output, NEW_COPY, refim) + + # Use smw_openim to clean up old keywords(?). + mw = smw_openim (im) + call smw_close (mw) + + if (naps == 1) + IM_NDIM(im) = 1 + else + IM_NDIM(im) = 2 + call imaddi (im, "SMW_NDIM", IM_NDIM(im)) + IM_LEN(im,2) = naps + if (IM_PIXTYPE(im) != TY_DOUBLE) + IM_PIXTYPE(im) = TY_REAL + + # Set new header. + mw = mw_open (NULL, 2) + call mw_newsystem (mw, "multispec", 2) + call mw_swtype (mw, axis, 2, "multispec", + "label=Wavelength units=Angstroms") + call smw_open (mw, NULL, im) + + nmax = 0 + do i = 1, naps { + sh1 = SH(sh,i,1) + call smw_gwattrs (MW(sh1), APINDEX(sh1), 1, ap, beam, dtype, + w1, dw, nw, z, aplow, aphigh, coeff) + call scb_default (SH(sh,i,1), NS(ns,i), + dtype, w1, dw, nw, z, Memc[coeff]) + call smw_swattrs (mw, i, 1, ap, beam, dtype, + w1, dw, nw, z, aplow, aphigh, Memc[coeff]) + call smw_sapid (mw, i, 1, TITLE(sh1)) + nmax = max (nmax, nw) + } + + IM_LEN(im,1) = nmax + + # Set MWCS header. + call smw_saveim (mw, im) + call smw_close (mw) + mw = smw_openim (im) + + # Create number combined image + if (noutput[1] != EOS) { + nout = immap (noutput, NEW_COPY, im) + IM_PIXTYPE(nout) = TY_INT + call sprintf (IM_TITLE(nout), SZ_LINE, "Number combined for %s") + call pargstr (output) + } + + call mfree (coeff, TY_CHAR) + call sfree (sp) +end + + +# SCB_DEFAULT - Set default values for the starting wavelength, ending +# wavelength, wavelength increment and spectrum length for the output +# spectrum. + +procedure scb_default (shdr, ns, dtype, w1, dw, nw, z, coeff) + +pointer shdr[ARB] # spectra structures +int ns # number of spectra +int dtype # dispersion type +double w1 # starting wavelength +double dw # wavelength increment +int nw # spectrum length +double z # redshift +char coeff[ARB] # nonlinear coefficient array + +bool clgetb() +int i, nwa, clgeti() +double w2, aux, w1a, w2a, dwa, clgetd() +pointer sh + +begin + if (clgetb ("first")) + return + + w1a = clgetd ("w1") + w2a = clgetd ("w2") + dwa = clgetd ("dw") + nwa = clgeti ("nw") + if (clgetb ("log")) + dtype = DCLOG + else + dtype = DCLINEAR + z = 0. + coeff[1] = EOS + + # Dispersion type + if (dtype == DCLINEAR) { + do i = 1, ns { + if (DC(shdr[i]) == DCNO) { + dtype = DCNO + break + } + } + } + + w1 = w1a + w2 = w2a + dw = dwa + nw = nwa + + # Starting wavelength + if (IS_INDEFD (w1)) { + if (IS_INDEFD (dw) || dw > 0.) { + w1 = MAX_REAL + do i = 1, ns { + sh = shdr[i] + if (WP(sh) > 0.) + aux = W0(sh) + else + aux = W1(sh) + if (aux < w1) + w1 = aux + } + } else { + w1 = -MAX_REAL + do i = 1, ns { + sh = shdr[i] + if (WP(sh) > 0.) + aux = W1(sh) + else + aux = W0(sh) + if (aux > w1) + w1 = aux + } + } + } + + # Ending wavelength + if (IS_INDEFD (w2)) { + if (IS_INDEFD (dw) || dw > 0.) { + w2 = -MAX_REAL + do i = 1, ns { + sh = shdr[i] + if (WP(sh) > 0.) + aux = W1(sh) + else + aux = W0(sh) + if (aux > w2) + w2 = aux + } + } else { + w2 = MAX_REAL + do i = 1, ns { + sh = shdr[i] + if (WP(sh) > 0.) + aux = W0(sh) + else + aux = W1(sh) + if (aux < w2) + w2 = aux + } + } + } + + # Wavelength increment + if (IS_INDEFD (dw)) { + dw = MAX_REAL + do i = 1, ns { + aux = abs (WP(shdr[i])) + if (aux < dw) + dw = aux + } + } + if ((w2 - w1) / dw < 0.) + dw = -dw + + # Spectrum length + if (IS_INDEFI (nw)) + nw = int ((w2 - w1) / dw + 0.5) + 1 + + # Adjust the values. + if (IS_INDEFD (dwa)) + dw = (w2 - w1) / (nw - 1) + else if (IS_INDEFD (w2a)) + w2 = w1 + (nw - 1) * dw + else if (IS_INDEFD (w1a)) + w1 = w2 - (nw - 1) * dw + else { + nw = int ((w2 - w1) / dw + 0.5) + 1 + w2 = w1 + (nw - 1) * dw + } +end diff --git a/noao/onedspec/scombine/x_scombine.x b/noao/onedspec/scombine/x_scombine.x new file mode 100644 index 00000000..33943271 --- /dev/null +++ b/noao/onedspec/scombine/x_scombine.x @@ -0,0 +1 @@ +task scombine = t_scombine -- cgit