Repatch (from linux) of OSX IRAF

1 files changed, 471 insertions, 0 deletions

diff --git a/noao/onedspec/odcombine/src/iccclip.gx b/noao/onedspec/odcombine/src/iccclip.gx
new file mode 100644
index 00000000..609b3448
--- /dev/null
+++ b/noao/onedspec/odcombine/src/iccclip.gx
@@ -0,0 +1,471 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	"../icombine.h"
+
+define	MINCLIP		2	# Mininum number of images for algorithm
+
+$for (sird)
+# IC_ACCDCLIP -- Reject pixels using CCD noise parameters about the average
+
+procedure ic_accdclip$t (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
+$if (datatype == sil)
+real	average[npts]		# Average
+$else
+PIXEL	average[npts]		# Average
+$endif
+
+int	i, j, k, l, jj, n1, n2, nin, nk, maxkeep
+$if (datatype == sil)
+real	d1, low, high, sum, a, s, r, zero
+data	zero /0.0/
+$else
+PIXEL	d1, low, high, sum, a, s, r, zero
+data	zero /0$f/
+$endif
+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 = Mem$t[d[1]+k]
+			do j = 2, n1
+			    sum = sum + Mem$t[d[j]+k]
+			average[i] = sum / n1
+		    }
+		}
+		next
+	    }
+
+	    repeat {
+		if (n1 == 2) {
+		    sum = Mem$t[d[1]+k]
+		    sum = sum + Mem$t[d[2]+k]
+		    a = sum / 2
+		} else {
+		    low = Mem$t[d[1]+k]
+		    high = Mem$t[d[2]+k]
+		    if (low > high) {
+			d1 = low
+			low = high
+			high = d1
+		    }
+		    sum = 0.
+		    do j = 3, n1 {
+		       d1 = Mem$t[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 = Mem$t[dp1]
+			r = (d1 - a) / s
+			if (r < -lsigma || r > hsigma) {
+			    Memr[resid+n1] = abs(r)
+			    if (j < n1) {
+				dp2 = d[n1] + k
+				Mem$t[dp1] = Mem$t[dp2]
+				Mem$t[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 = Mem$t[dp1]
+			r = (d1 - a) / s
+			if (r < -lsigma || r > hsigma) {
+			    Memr[resid+n1] = abs(r)
+			    if (j < n1) {
+				dp2 = d[n1] + k
+				Mem$t[dp1] = Mem$t[dp2]
+				Mem$t[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 = Mem$t[dp1]
+				    Mem$t[dp1] = Mem$t[dp2]
+				    Mem$t[dp2] = d1
+				    mp2 = m[l] + k
+				    s = Memi[mp1]
+				    Memi[mp1] = Memi[mp2]
+				    Memi[mp2] = s
+				}
+			    }
+			}
+			sum = sum + Mem$t[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 = Mem$t[dp1]
+				    Mem$t[dp1] = Mem$t[dp2]
+				    Mem$t[dp2] = d1
+				    if (keepids) {
+					mp1 = m[j] + k
+					mp2 = m[l] + k
+					s = Memi[mp1]
+					Memi[mp1] = Memi[mp2]
+					Memi[mp2] = s
+				    }
+				}
+			    }
+			}
+			sum = sum + Mem$t[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_mccdclip$t (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
+$if (datatype == sil)
+real	median[npts]		# Median
+$else
+PIXEL	median[npts]		# Median
+$endif
+
+int	i, j, k, l, id, n1, n2, n3, nl, nh, nin, maxkeep
+real	r, s
+pointer	sp, resid, mp1, mp2
+$if (datatype == sil)
+real	med, zero
+data	zero /0.0/
+$else
+PIXEL	med, zero
+data	zero /0$f/
+$endif
+
+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 = Mem$t[d[n3-1]+k]
+		    med = (med + Mem$t[d[n3]+k]) / 2.
+		} else
+		    med = Mem$t[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 - Mem$t[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 = (Mem$t[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 - Mem$t[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 = (Mem$t[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 >= 1.)) {
+		j = max (nl, n1 + 1)
+		if (keepids) {
+		    do l = 1, min (n1, nl-1) {
+			Mem$t[d[l]+k] = Mem$t[d[j]+k]
+			if (grow >= 1.) {
+			    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) {
+			Mem$t[d[l]+k] = Mem$t[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
+$endfor