Initial commit

1 files changed, 240 insertions, 0 deletions

diff --git a/noao/onedspec/splot/eqwidthcp.x b/noao/onedspec/splot/eqwidthcp.x
new file mode 100644
index 00000000..4fc4fd3d
--- /dev/null
+++ b/noao/onedspec/splot/eqwidthcp.x
@@ -0,0 +1,240 @@
+include	<gset.h>
+
+
+# EQWIDTH_CP -- Equivalent width following algorithm provided by
+#               Caty Pilachowski. This assumes a Gaussian line profile
+#               and fits to the specified core level, the width at the
+#		specified flux level above the core, and the specified
+#		continuum.  The line position is found by searching
+#               near the vertical cursor for the nearest minimum.
+
+define	LEFT	1	# Fit to left edge
+define	RIGHT	2	# Fit to right edge
+define	BOTH	3	# Fit to both edges
+
+procedure eqwidth_cp (sh, gfd, center, cont, ylevel, y, n, key, fd1, fd2,
+	xg, yg, sg, lg, pg, ng)
+
+pointer	sh
+int	gfd
+real	center, cont, ylevel
+real	y[n]
+int	n
+int	key
+int	fd1, fd2
+pointer	xg, yg, sg, lg, pg	# Pointers to fit parameters
+int	ng			# Number of components
+
+int	i, i1, i2, isrch, icore, edge
+double	xleft, xright, rcore, rinter, yl, gfwhm, lfwhm, flux, eqw, w, w1, w2
+double	xpara[3], ypara[3], coefs[3], xcore, ycore
+double	shdr_lw(), shdr_wl()
+
+# Initialize reasonable values
+#	isrch -- nr of pixels on either side of cursor to search for min
+
+data	isrch /3/
+
+begin
+	# Check continuum.
+	if (cont <= 0.) {
+	    call eprintf ("Continuum cannot be less than zero.\n")
+	    return
+	}
+
+	# Determine which edges of the line to use.
+	switch (key) {
+	case 'a', 'l':
+	    edge = LEFT
+	case 'b', 'r':
+	    edge = RIGHT
+	default:
+	    edge = BOTH
+	}
+
+	# Search for local minimum or maximum
+	icore = max (1, min (n, nint (shdr_wl (sh, double(center)))))
+	i1 = max (1, icore-isrch)
+	i2 = min (n, icore+isrch)
+
+	# If half lines is selected, restrict the search
+	if (edge == LEFT)
+	    i2 = max (i2-2, icore+1)
+	if (edge == RIGHT)
+	    i1 = min (i1+2, icore-1)
+
+	# Search for core.
+	# Someday it may be desirable to use parabolic interpolation
+	# to locate an estimated minimum or maximum for the region
+	do i = i1, i2 {
+	    if (abs (y[i] - cont) > abs (y[icore] - cont))
+		icore = i
+	}
+
+	# Fit parabola to three points around minimum pixel
+	xpara[1] = icore - 1
+	xpara[2] = icore
+	xpara[3] = icore + 1
+	ypara[1] = y[icore-1]
+	ypara[2] = y[icore]
+	ypara[3] = y[icore+1]
+
+	call para (xpara, ypara, coefs)
+
+	# Compute pixel value at minimum
+	xcore = -coefs[2] / 2.0 / coefs[3]
+	ycore = coefs[1] + coefs[2] * xcore + coefs[3] * xcore**2
+
+	# Locate left and right line edges.  If the ylevel is INDEF then use
+	# the half flux point.
+	if (IS_INDEF (ylevel))
+	    yl = (cont + ycore) / 2.
+	else
+	    yl = ylevel
+
+	rcore = abs (ycore - cont)
+	rinter = abs (yl - cont)
+
+	if (rcore <= rinter) {
+	    call eprintf (
+		"Y cursor must be between the continuum and the line core\n")
+	    return
+	}
+
+	# Bound flux level of interest
+	if ((edge == LEFT) || (edge == BOTH)) {
+	    for (i=icore; i >= 1; i=i-1)
+		if (abs (y[i] - cont) < rinter)
+		    break
+
+	    if (i < 1) {
+		call eprintf ("Can't find left edge of line\n")
+		return
+	    }
+
+	    xleft = float (i) + (yl - y[i]) / (y[i+1] - y[i])
+	    if (edge == LEFT)
+	        xright = xcore + (xcore - xleft)
+	}
+
+	# Now bound the right side
+	if ((edge == RIGHT) || (edge == BOTH)) {
+	    for (i=icore; i <= n; i=i+1)
+		if (abs (y[i] - cont) < rinter)
+		    break
+
+	    if (i > n) {
+		call eprintf ("Can't find right edge of line\n")
+		return
+	    }
+
+	    xright = float (i) - (yl - y[i]) / (y[i-1] - y[i])
+	    if (edge == RIGHT)
+	        xleft  = xcore - (xright - xcore)
+	}
+
+	# Compute in wavelength
+	w = shdr_lw (sh, double(xcore))
+	w1 = shdr_lw (sh, double(xleft))
+	w2 = shdr_lw (sh, double(xright))
+
+	# Apply Gaussian model
+	gfwhm = 1.665109 * abs (w2 - w1) / 2. / sqrt (log (rcore/rinter))
+	lfwhm = 0.
+	rcore = ycore - cont
+	flux = 1.064467 * rcore * gfwhm
+	eqw = -flux / cont
+
+	call printf (
+	    "center = %9.7g, eqw = %9.4g, gfwhm = %9.4g\n")
+	    call pargd (w)
+	    call pargd (eqw)
+	    call pargd (gfwhm)
+
+	if (fd1 != NULL) {
+	    call fprintf (fd1, " %9.7g %9.7g %9.6g %9.4g %9.6g %9.4g %9.4g\n")
+		call pargd (w)
+		call pargr (cont)
+		call pargd (flux)
+		call pargd (eqw)
+		call pargd (ycore - cont)
+		call pargd (gfwhm)
+		call pargd (lfwhm)
+	}
+	if (fd2 != NULL) {
+	    call fprintf (fd2, " %9.7g %9.7g %9.6g %9.4g %9.6g %9.4g %9.4g\n")
+		call pargd (w)
+		call pargr (cont)
+		call pargd (flux)
+		call pargd (eqw)
+		call pargd (ycore - cont)
+		call pargd (gfwhm)
+		call pargd (lfwhm)
+	}
+
+	# Mark line computed
+	call gline (gfd, real(w), cont, real(w), real(ycore))
+	call gline (gfd, real(w1), real(yl), real(w2), real(yl))
+
+	w1 = w - 2 * gfwhm
+	w2 = cont + rcore * exp (-(1.665109*(w1-w)/gfwhm)**2)
+	call gseti (gfd, G_PLTYPE, 2)
+	call gseti (gfd, G_PLCOLOR, 2)
+	call gamove (gfd, real(w1), real(w2))
+	for (; w1 <= w+2*gfwhm; w1=w1+0.05*gfwhm) {
+	    w2 = cont + rcore * exp (-(1.665109*(w1-w)/gfwhm)**2)
+	    call gadraw (gfd, real(w1), real(w2))
+	}
+	call gseti (gfd, G_PLTYPE, 1)
+	call gseti (gfd, G_PLCOLOR, 1)
+
+	# Save fit parameters
+	if (ng == 0) {
+	    call malloc (xg, 1, TY_REAL)
+	    call malloc (yg, 1, TY_REAL)
+	    call malloc (sg, 1, TY_REAL)
+	    call malloc (lg, 1, TY_REAL)
+	    call malloc (pg, 1, TY_INT)
+	} else if (ng != 1) {
+	    call realloc (xg, 1, TY_REAL)
+	    call realloc (yg, 1, TY_REAL)
+	    call realloc (sg, 1, TY_REAL)
+	    call realloc (lg, 1, TY_REAL)
+	    call realloc (pg, 1, TY_INT)
+	}
+	Memr[xg] = w
+	Memr[yg] = rcore
+	Memr[sg] = gfwhm
+	Memr[lg] = lfwhm
+	Memi[pg] = 1
+	ng = 1
+end
+
+# PARA -- Fit a parabola to three points
+
+procedure para (x, y, c)
+
+double	x[3], y[3], c[3]
+double  x12, x13, x23, x213, x223, y13, y23
+
+begin
+        x12 = x[1] - x[2]
+        x13 = x[1] - x[3]
+        x23 = x[2] - x[3]
+
+	if (x12 == 0. || x13 == 0. || x23 == 0.)
+	    call error (1, "X points are not distinct")
+
+	# Compute relative to an origin at x[3]
+        x213 = x13 * x13
+        x223 = x23 * x23
+        y13 = y[1] - y[3]
+        y23 = y[2] - y[3]
+        c[3] = (y13 - y23 * x13 / x23) / (x213 - x223 * x13 / x23)
+        c[2] = (y23 - c[3] * x223) / x23
+        c[1] = y[3]
+
+	# Compute relative to an origin at 0.
+	c[1] = c[1] - x[3] * (c[2] - c[3] * x[3])
+	c[2] = c[2] - 2 * c[3] * x[3]
+end