1 files changed, 57 insertions, 0 deletions
diff --git a/pkg/xtools/icfit/icrejectd.x b/pkg/xtools/icfit/icrejectd.x
new file mode 100644
index 00000000..36985923
--- /dev/null
+++ b/pkg/xtools/icfit/icrejectd.x
@@ -0,0 +1,57 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<math/curfit.h>
+include	"names.h"
+
+# IC_REJECT -- Reject points with large residuals from the fit.
+#
+# The sigma of the fit residuals is calculated.  The rejection thresholds
+# are set at low_reject*sigma and high_reject*sigma.  Points outside the
+# rejection threshold are rejected from the fit and flagged in the rejpts
+# array.  Finally, the remaining points are refit.
+
+procedure ic_rejectd (cv, x, y, w, rejpts, npts, low_reject, high_reject,
+    niterate, grow, nreject)
+
+pointer	cv				# Curve descriptor
+double	x[npts]				# Input ordinates
+double	y[npts]				# Input data values
+double	w[npts]				# Weights
+int	rejpts[npts]			# Points rejected
+int	npts				# Number of input points
+real	low_reject, high_reject		# Rejection threshold
+int	niterate			# Number of rejection iterations
+real	grow				# Rejection radius
+int	nreject				# Number of points rejected
+
+int	i, ierr, nit, newreject
+errchk	ic_deviantd
+
+begin
+	# Initialize rejection.
+	nreject = 0
+	call amovki (NO, rejpts, npts)
+
+	if (niterate <= 0)
+	    return
+
+	# Find deviant points.  If an error occurs reduce the number of
+	# iterations and start again.
+	iferr {
+	    nit = 0
+	    do i = 1, niterate {
+		call ic_deviantd (cv, x, y, w, rejpts, npts, low_reject,
+		    high_reject, grow, YES, nreject, newreject)
+		nit = nit + 1
+		if (newreject == 0)
+		    break
+	    }
+	} then {
+	    call dcvfit (cv, x, y, w, npts, WTS_USER, ierr)
+	    nreject = 0
+	    call amovki (NO, rejpts, npts)
+	    do i = 1, nit
+		call ic_deviantd (cv, x, y, w, rejpts, npts, low_reject,
+		    high_reject, grow, YES, nreject, newreject)
+	}
+end