Initial commit

1 files changed, 199 insertions, 0 deletions

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	<smw.h>
+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