# File rvsao/Emvel/emlines.x
# December 4, 1991
# By Doug Mink, Harvard-Smithsonian Center for Astrophysics
# After Bill Wyatt's C and Forth code

# Copyright(c) 1991 Smithsonian Astrophysical Observatory
# You may do anything you like with this file except remove this copyright.
# The Smithsonian Astrophysical Observatory makes no representations about
# the suitability of this software for any purpose.  It is provided "as is"
# without express or implied warranty.

#--- Routines to do simple line finding & fitting

procedure emlines (fdata, fcont, ipix1, ipix2, pts_fit, lsigma, center
		   height, maxlines, nlines, avwidth, avheight, debug)

real	fdata[ARB]		# Spectrum (continuum subtracted)
real	fcont[ARB]		# Continuum fit to spectrum
int	ipix1, ipix2		# First and last pixels to search
int	pts_fit			# Number of points to fit
double	lsigma			# Number of sigma for line acceptance
double	center[ARB]		# Centers for lines
double	height[ARB]		# Heights for lines
int	maxlines		# Maximum number of lines that can be returned
int	nlines			# Number of lines returned
double	avwidth			# Average line width (returned)
double	avheight		# Average line height (returned)
bool	debug 

int     i
int	nl			# number of lines found
double	avw, avh
bool	isline()		# determine whether at center of line
bool	linepar()		# compute line parameters
double	lcenter,lheight,lwidth	# single line center, height, width

begin

	nl = 0
	avw = 0.d0
	avh = 0.d0
#	if (debug) {
#	    call printf("EMLINES: searching pixels %d to %d\n")
#		call pargi (ipix1)
#		call pargi (ipix2)
#	    call printf("    Center    Height   Width   Continuum Slope\n")
#	    }

	i = ipix1
	do i = ipix1, ipix2 {
	    if (isline (i,fdata,fcont,pts_fit,lsigma,debug)) {
		if (linepar (i,fdata,fcont,pts_fit,lcenter,lheight,
			      lwidth,debug)) {
		    if (nl < maxlines) {
			nl = nl + 1
			center[nl] = lcenter
			height[nl] = lheight
			avw = avw + lwidth
			avh = avh + lheight
			}
		    if (debug) {
			call printf(" %9.2f %9.2f %7.2f\n")
			    call pargd (lcenter)
			    call pargd (lheight)
			    call pargd (lwidth)
			}
		    }
		}
	    }
	nlines = nl
	if (nl > 0) {
	    avwidth = avw / double (nl)
	    avheight = avh / double (nl)
	    }
	else {
	    avwidth = 0.d0
	    avheight = 0.d0
	    }
	return
end


#  Given a pixel index, returns 1 if:
#    1) the pixel value is the max of the npts+1 neighbors on each side.
#    2) the sqrt(average counts in the 3 peak pixels) * sigma are <= height.

bool procedure isline (xval, fdata, fcont, npts, lsigma, debug)

int	xval		# index of center pixel
real	fdata[ARB]	# Spectrum with continuum removed
real	fcont[ARB]	# Continuum fit to spectrum
int	npts		# number of pixels to check for peak (/2)
double	lsigma		# number of sigma for line acceptance
double	sigma		# sigma at line
bool	debug

int	i
real	ymax, ycont
double	sqrt()

begin

	ymax = fdata[xval]

#	if (debug) {
#	    call printf ("ISLINE: checking line at pixel %d = %f\n")
#		call pargi (xval)
#		call pargr (fdata[xval])
#	    }

#  If peak is to the left of this pixel, return

	do i = 1, 3 {
	    if (fdata[xval-i] >= ymax) {
		return FALSE
		}
	    }

#  If peak is to the right of this pixel, return

	do i = 1, 3 {
	    if (fdata[xval+i] >= ymax) {
		return FALSE
		}
	    }

	ymax = (fdata[xval-1] + fdata[xval] + fdata[xval + 1])/ 3.0
	ycont = fcont[xval]
        sigma = sqrt (abs (ymax + ycont))

#	if (debug) {
#	    call printf ("ISLINE: line value= %8.3f sigma= %8.3f\n")
#		call pargr (ymax)
#		call pargd (sigma)
#	    }
	
        if (ymax >= sigma * lsigma)
	    return TRUE
	else if (ymax > abs (ycont * 1.e10))
	    return TRUE
	else
	    return FALSE
end


#  from the FORTH program by J. Tonry: 
# 
#  If n = pts_fit (so 2n + 1 points are fit)
#    f(x) = h + a (x-x')**2 
#         = h [ 1 - 2 { (x-x')/w }**2 ]
# 
#  also,  a, b, c are:
#    a = -45/{n(n+1)(2n-1)(2n+1)(2n+3)} * Sum(i=-n, n) {n(n+1)/3 - i**2}yi
#    b = -3/{n(n+1)(2n+1)} * Sum(i=-n,n) {i * yi}
#    c = 1/(2n+1) * Sum(i=-n,n) yi
# 
#  Simplifying, if:
#    u = 3 / { n (n+1) }
#    v =  2n + 1 
#    w = (2n-1)(2n+3) = 4n**2 + 4n - 3
# 
#  Then:
#    a = -15 u / { v (2n-1)(2n+3)} * Sum(i=-n,n) { yi/u - yi * i**2 }
#    a = -u/v * Sum(-n, n){ i * yi }
#    a = Sum(-n, n){ yi } / v
# 
#  Then:
#    x' = center = x0 + b/2a
#    h  = height = c - a [ n(n+1)/3 + (b/2a)**2 ]
#    w  = width  = sqrt(-2h/a)

# returns true for success, false for error

bool procedure linepar (xzero, fdata, fcont, npts, center, height, width, debug)

int	xzero		# Initial guess of line (should be peak!)
real	fdata[ARB]	# Spectrum data to fit
real	fcont[ARB]	# Continuum fit
int	npts		# +/- number pts to fit
double	center		# line center in pixels (returned)
double	height		# line height (returned)
double	width		# line width in pixels (returned)
bool	debug		# true for diagnostic printing

int	i
double	a, b, c
double	u, v, w
double	sumy, sumyi, sumyii
double	y

begin

	if (npts <= 0 || xzero < npts) {
	    if (debug) {
		call printf ("LINEPAR: npts= %d xzero= %d\n")
		    call pargi (npts)
		    call pargi (xzero)
		}
	    return FALSE
	    }

	u = 3.d0 / (npts * npts + npts)
	v = 2.d0 * npts + 1.d0
	w = ((4.d0 * npts) + 4.d0) * npts - 3.d0

	sumy   = 0.d0
	sumyi  = 0.d0
	sumyii = 0.d0
	do i = -npts ,npts {
	    y = fdata[xzero+i]
	    sumy   = sumy + y
	    sumyi  = sumyi + (y * i)
	    sumyii = sumyii + (y * i * i)
	    }

	a = (-15.d0 * u) / (v * w) * (sumy / u - sumyii)
	b = (-u * sumyi) / v
	c = sumy / v

	w = b / ( 2.d0 * a )
	center    = xzero + w
	height    = c - a * (1.d0 / u + (w * w))
	width     = sqrt (abs (-2.d0 * height / a))

#	if (debug) {
#	    call printf ("LINEPAR: line found at pixel %d -> %8.3f\n")
#		call pargi (xzero)
#		call pargd (center)
#	    }

	return TRUE
end

# Dec  4 1991	Use interactive continuum fit to whole spectrum instead
#		of fitting local continuum for each line