# File Emsao/gvitrv.x
# February 19, 1997
# Modified by Doug Mink, Harvard-Smithsonian Center for Astrophysics
# Most of this code is from noao.onedspec.splot's gfit.x, voigt.x, and deblend.x

include	<error.h>
include	<mach.h>
include	<gset.h>

define	NSUB	3	# Number of pixel subsamples
define	MC_N	50	# Monte-Carlo samples
define	MC_P	10	# Percent done interval (percent)
define	MC_SIG	68	# Sigma sample point (percent)

# GFIT -- Fit Gaussian

procedure gfit (sh, wx1, wy1, wcs, pix, n, fd1, fd2, xg, yg, sg, lg, pg,ng)

pointer	sh			# SHDR pointer
real	wx1, wy1		# Cursor position
real	wcs[n]			# Spectrum data
real	pix[n]			# Spectrum data
int	n			# Number of points
int	fd1, fd2		# Output file descriptors
pointer	xg, yg, sg, lg, pg	# Pointers to fit parameters
int	ng			# Number of components

int	fit[5], nsub, mc_p, mc_sig, mc_n
int	i, j, i1, npts, nlines, wc, key
long	seed
real	w, dw, wyc, wx, wy, wx2, wy2, v, u
real	slope, peak, flux, cont, gfwhm, lfwhm, eqw, scale, sscale, chisq
real	sigma0, invgain, wyc1, slope1, flux1, cont1, eqw1
bool	fitit
pointer	xg1, yg1, sg1, lg1
pointer	sp, cmd, x, y, s, z, ym, conte, xge, yge, sge, lge, fluxe, eqwe

int	clgeti(), clgcur()
real	clgetr(), model(), gasdev(), asumr()
errchk	dofit, dorefit

define	done_	99

begin
	call smark (sp)
	call salloc (cmd, SZ_FNAME, TY_CHAR)

	# Input cursor is first continuum point now get second continuum point.
	call printf ("k again:")
	if (clgcur ("cursor", wx2, wy2, wc, key, Memc[cmd], SZ_FNAME) == EOF) {
	    call sfree (sp)
	    return
	}

	# Set pixel indices and determine number of points to fit.
	call fixx (sh, wx1, wx2, wy1, wy2, i1, j)
	npts = j - i1 + 1
	if (npts < 3) {
	    call eprintf ("At least 3 points are required\n")
	    call sfree (sp)
	    return
	}

	# Allocate space for the points to be fit.
	call salloc (x, npts, TY_REAL)
	call salloc (y, npts, TY_REAL)
	call salloc (s, npts, TY_REAL)
	call salloc (z, npts, TY_REAL)

	# Scale the data.
	mc_n = clgeti ("nerrsample")
	sigma0 = clgetr ("sigma0")
	invgain = clgetr ("invgain")
	if (IS_INDEF(sigma0) || IS_INDEF(invgain) || sigma0<0. ||
	    invgain<0. || (sigma0 == 0. && invgain == 0.)) {
	    sigma0 = INDEF
	    invgain = INDEF
	}
	scale = 0.
	do i = 1, npts {
	    Memr[x+i-1] = wcs[i1+i-1]
	    Memr[y+i-1] = pix[i1+i-1]
	    if (Memr[y+i-1] <= 0.)
		if (invgain != 0.) {
		    sigma0 = INDEF
		    invgain = INDEF
		    call eprintf (
			"WARNING: Cannot compute errors with non-zero gain")
		    call eprintf (
			" and negative pixel values.\n")
		}
	    scale = max (scale, abs (Memr[y+i-1]))
	}
	if (IS_INDEF(sigma0)) {
	    call amovkr (1., Memr[s], npts)
	    sscale = 1.
	} else {
	    do i = 1, npts
		Memr[s+i-1] = sqrt (sigma0 ** 2 + invgain * Memr[y+i-1])
	    sscale = asumr (Memr[s], npts) / npts
	}
	call adivkr (Memr[y], scale, Memr[y], npts)
	call adivkr (Memr[s], sscale, Memr[s], npts)

	# Allocate memory.
	nlines = 1
	if (ng == 0) {
	    call malloc (xg, nlines, TY_REAL)
	    call malloc (yg, nlines, TY_REAL)
	    call malloc (sg, nlines, TY_REAL)
	    call malloc (lg, nlines, TY_REAL)
	    call malloc (pg, nlines, TY_INT)
	} else if (ng != nlines) {
	    call realloc (xg, nlines, TY_REAL)
	    call realloc (yg, nlines, TY_REAL)
	    call realloc (sg, nlines, TY_REAL)
	    call realloc (lg, nlines, TY_REAL)
	    call realloc (pg, nlines, TY_INT)
	}
	ng = nlines

	# Do fit.
	fit[1] = 1
	fit[2] = 2
	fit[3] = 2
	fit[4] = 2
	fit[5] = 2

	# Setup initial estimates.
	slope = (wy2-wy1) / (wx2-wx1) / scale
	wyc = wy1 / scale - slope * wx1
	wx = 0
	do i = 0, npts-1 {
	    w = Memr[x+i]
	    wy = Memr[y+i] - wyc - slope * w
	    if (abs (wy) > wx) {
		wx = abs (wy)
		j = i
		Memr[xg] = w
		Memr[yg] = wy
	    }
	}

	if (j > 0 && j < npts-1) {
	    w = Memr[x+j-1]
	    wy = min (0.99, max (0.01, abs (Memr[y+j-1] - wyc - slope*w) / wx))
	    gfwhm = 2.355 * sqrt (-0.5 * (w-Memr[xg])**2 / log (wy))
	    w = Memr[x+j+1]
	    wy = min (0.99, max (0.01, abs (Memr[y+j+1] - wyc - slope*w) / wx))
	    gfwhm = (gfwhm + 2.355 * sqrt (-0.5*(w-Memr[xg])**2/log (wy))) / 2
	} else
	    gfwhm = 0.3 * abs (Memr[x+npts-1] - Memr[x])

	switch (key) {
	case 'l':
	    Memr[sg] = 0.
	    Memr[lg] = gfwhm
	    Memi[pg] = 2
	case 'v':
	    Memr[sg] = 0.5 * gfwhm
	    Memr[lg] = 0.5 * gfwhm
	    Memi[pg] = 3
	default:
	    Memr[sg] = gfwhm
	    Memr[lg] = 0.
	    Memi[pg] = 1
	}

	nsub = NSUB
	dw = (wcs[n] - wcs[1]) / (n - 1)
	iferr (call dofit (fit, Memr[x], Memr[y], Memr[s], npts, dw, nsub,
	    wyc, slope, Memr[xg], Memr[yg], Memr[sg], Memr[lg], Memi[pg],
	    ng, chisq)) {
	    fitit = false
	    goto done_
	}

	# Compute Monte-Carlo errors.
	if (mc_n > 9 && !IS_INDEF(sigma0)) {
	    mc_p = nint (mc_n * MC_P / 100.)
	    mc_sig = nint (mc_n * MC_SIG / 100.)

	    call salloc (ym, npts, TY_REAL)
	    call salloc (xg1, ng, TY_REAL)
	    call salloc (yg1, ng, TY_REAL)
	    call salloc (sg1, ng, TY_REAL)
	    call salloc (lg1, ng, TY_REAL)
	    call salloc (conte, mc_n*ng, TY_REAL)
	    call salloc (xge, mc_n*ng, TY_REAL)
	    call salloc (yge, mc_n*ng, TY_REAL)
	    call salloc (sge, mc_n*ng, TY_REAL)
	    call salloc (lge, mc_n*ng, TY_REAL)
	    call salloc (fluxe, mc_n*ng, TY_REAL)
	    call salloc (eqwe, mc_n*ng, TY_REAL)
	    do i = 1, npts {
		w = Memr[x+i-1]
		Memr[ym+i-1] = model (w, dw, nsub, Memr[xg], Memr[yg],
		    Memr[sg], Memr[lg], Memi[pg], ng) + wyc + slope * w
	    }
	    seed = 1
	    do i = 0, mc_n-1 {
		if (i > 0 && mod (i, mc_p) == 0) {
		    call printf ("%2d ")
			call pargi (100 * i / mc_n)
		    call flush (STDOUT)
		}
		do j = 1, npts
		    Memr[y+j-1] = Memr[ym+j-1] +
			sscale / scale * Memr[s+j-1] * gasdev (seed)
		wyc1 = wyc
		slope1 = slope
		call amovr (Memr[xg], Memr[xg1], ng)
		call amovr (Memr[yg], Memr[yg1], ng)
		call amovr (Memr[sg], Memr[sg1], ng)
		call amovr (Memr[lg], Memr[lg1], ng)
		call dorefit (fit, Memr[x], Memr[y], Memr[s], npts,
		    dw, nsub, wyc1, slope1,
		    Memr[xg1], Memr[yg1], Memr[sg1], Memr[lg1], Memi[pg], ng,
		    chisq)

		do j = 0, ng-1 {
		    cont = wyc + slope * Memr[xg+j]
		    cont1 = wyc1 + slope1 * Memr[xg+j]
		    switch (Memi[pg+j]) {
		    case 1:
			flux = 1.064467 * Memr[yg+j] * Memr[sg+j]
			flux1 = 1.064467 * Memr[yg1+j] * Memr[sg1+j]
		    case 2:
			flux = 1.570795 * Memr[yg+j] * Memr[lg+j]
			flux1 = 1.570795 * Memr[yg1+j] * Memr[lg1+j]
		    case 3:
			call voigt (0., 0.832555*Memr[lg+j]/Memr[sg+j], v, u)
			flux = 1.064467 * Memr[yg+j] * Memr[sg+j] / v
			call voigt (0., 0.832555*Memr[lg1+j]/Memr[sg1+j], v, u)
			flux1 = 1.064467 * Memr[yg1+j] * Memr[sg1+j] / v
		    }
		    if (cont > 0. && cont1 > 0.) {
			eqw = -flux / cont
			eqw1 = -flux1 / cont1
		    } else {
			eqw = 0.
			eqw1 = 0.
		    }
		    Memr[conte+j*mc_n+i] = abs (cont1 - cont)
		    Memr[xge+j*mc_n+i] = abs (Memr[xg1+j] - Memr[xg+j])
		    Memr[yge+j*mc_n+i] = abs (Memr[yg1+j] - Memr[yg+j])
		    Memr[sge+j*mc_n+i] = abs (Memr[sg1+j] - Memr[sg+j])
		    Memr[lge+j*mc_n+i] = abs (Memr[lg1+j] - Memr[lg+j])
		    Memr[fluxe+j*mc_n+i] = abs (flux1 - flux)
		    Memr[eqwe+j*mc_n+i] = abs (eqw1 - eqw)
		}
	    }
	    do j = 0, ng-1 {
		call asrtr (Memr[conte+j*mc_n], Memr[conte+j*mc_n], mc_n)
		call asrtr (Memr[xge+j*mc_n], Memr[xge+j*mc_n], mc_n)
		call asrtr (Memr[yge+j*mc_n], Memr[yge+j*mc_n], mc_n)
		call asrtr (Memr[sge+j*mc_n], Memr[sge+j*mc_n], mc_n)
		call asrtr (Memr[lge+j*mc_n], Memr[lge+j*mc_n], mc_n)
		call asrtr (Memr[fluxe+j*mc_n], Memr[fluxe+j*mc_n], mc_n)
		call asrtr (Memr[eqwe+j*mc_n], Memr[eqwe+j*mc_n], mc_n)
	    }
	    call amulkr (Memr[conte], scale, Memr[conte], mc_n*ng)
	    call amulkr (Memr[yge], scale, Memr[yge], mc_n*ng)
	    call amulkr (Memr[fluxe], scale, Memr[fluxe], mc_n*ng)
	}

	call amulkr (Memr[yg], scale, Memr[yg], ng)
	wyc = (wyc + slope * wx1) * scale
	slope = slope * scale

	# Compute model spectrum with continuum and plot.
	fitit = true
	do i = 1, npts {
	    w = wcs[i1+i-1]
	    Memr[z+i-1] = model (w, dw, nsub, Memr[xg], Memr[yg],
		Memr[sg], Memr[lg], Memi[pg], ng) + wyc + slope * (w - wx1)
	}

done_
	# Log computed values
	if (fitit) {
	    do i = 1, nlines {
		w = Memr[xg+i-1]
		cont = wyc + slope * (w - wx1)
		peak = Memr[yg+i-1]
		gfwhm = Memr[sg+i-1]
		lfwhm = Memr[lg+i-1]
		switch (Memi[pg+i-1]) {
		case 1:
		    flux = 1.064467 * peak * gfwhm
		    if (cont > 0.)
			eqw = -flux / cont
		    else
			eqw = INDEF
		    call printf (
	"\n%d: center = %8.6g, flux = %8.4g, eqw = %6.4g, gfwhm = %6.4g")
			call pargi (i)
			call pargr (w)
			call pargr (flux)
			call pargr (eqw)
			call pargr (gfwhm)
		case 2:
		    flux = 1.570795 * peak * lfwhm
		    if (cont > 0.)
			eqw = -flux / cont
		    else
			eqw = INDEF
		    call printf (
	"\n%d: center = %8.6g, flux = %8.4g, eqw = %6.4g, lfwhm = %6.4g")
			call pargi (i)
			call pargr (w)
			call pargr (flux)
			call pargr (eqw)
			call pargr (lfwhm)
		case 3:
		    call voigt (0., 0.832555*lfwhm/gfwhm, v, u)
		    flux = 1.064467 * peak * gfwhm / v
		    if (cont > 0.)
			eqw = -flux / cont
		    else
			eqw = INDEF
		    call printf (
	"\n%d: center = %8.6g, eqw = %6.4g, gfwhm = %6.4g, lfwhm = %6.4g")
			call pargi (i)
			call pargr (w)
			call pargr (eqw)
			call pargr (gfwhm)
			call pargr (lfwhm)
		}
		if (fd1 != NULL) {
		    call fprintf (fd1,
			" %9.7g %9.7g %9.6g %9.4g %9.6g %9.4g %9.4g\n")
			call pargr (w)
			call pargr (cont)
			call pargr (flux)
			call pargr (eqw)
			call pargr (peak)
			call pargr (gfwhm)
			call pargr (lfwhm)
		}
		if (fd2 != NULL) {
		    call fprintf (fd2,
			" %9.7g %9.7g %9.6g %9.4g %9.6g %9.4g %9.4g\n")
			call pargr (w)
			call pargr (cont)
			call pargr (flux)
			call pargr (eqw)
			call pargr (peak)
			call pargr (gfwhm)
			call pargr (lfwhm)
		}
		if (mc_n > 9 && !IS_INDEF(sigma0)) {
		    if (fd1 != NULL) {
			call fprintf (fd1,
		"  (%7.5g) (%7w) (%7.4g) (%7.4g) (%7.4g) (%7.4g) (%7.4g)\n")
			    call pargr (Memr[xge+(i-1)*mc_n+mc_sig])
			    call pargr (Memr[fluxe+(i-1)*mc_n+mc_sig])
			    call pargr (Memr[eqwe+(i-1)*mc_n+mc_sig])
			    call pargr (Memr[yge+(i-1)*mc_n+mc_sig])
			    call pargr (Memr[sge+(i-1)*mc_n+mc_sig])
			    call pargr (Memr[lge+(i-1)*mc_n+mc_sig])
		    }
		    if (fd2 != NULL) {
			call fprintf (fd2,
		"  (%7.5g) (%7w) (%7.4g) (%7.4g) (%7.4g) (%7.4g) (%7.4g)\n")
			    call pargr (Memr[xge+(i-1)*mc_n+mc_sig])
			    call pargr (Memr[fluxe+(i-1)*mc_n+mc_sig])
			    call pargr (Memr[eqwe+(i-1)*mc_n+mc_sig])
			    call pargr (Memr[yge+(i-1)*mc_n+mc_sig])
			    call pargr (Memr[sge+(i-1)*mc_n+mc_sig])
			    call pargr (Memr[lge+(i-1)*mc_n+mc_sig])
		    }
		}
	    }
	} else {
	    call mfree (xg, TY_REAL)
	    call mfree (yg, TY_REAL)
	    call mfree (sg, TY_REAL)
	    call mfree (lg, TY_REAL)
	    call mfree (pg, TY_INT)
	    ng = 0
	}

	call sfree (sp)
end
# VOIGT -- Compute the real (Voigt function) and imaginary parts of the
# complex function w(z)=exp(-z**2)*erfc(-i*z) in the upper half-plane
# z=x+iy.  The maximum relative error of the real part is 2E-6 and the
# imaginary part is 5E-6.
#
# From: Humlicek, J. Quant. Spectrosc. Radiat. Transfer, V21, p309, 1979.

procedure voigt (xarg, yarg, wr, wi)

real	xarg		#I Real part of argument
real	yarg		#I Imaginary part of argument
real	wr		#O Real part of function
real	wi		#O Imaginary part of function

int	i
real	x, y, y1, y2, y3, d, d1, d2, d3, d4, r, r2
real	t[6], c[6], s[6]

data	t/.314240376,.947788391,1.59768264,2.27950708,3.02063703,3.8897249/
data	c/1.01172805,-.75197147,1.2557727e-2,1.00220082e-2,-2.42068135e-4,
	  5.00848061e-7/
data	s/1.393237,.231152406,-.155351466,6.21836624e-3,9.19082986e-5,
	  -6.27525958e-7/

begin
	x = xarg
	y = abs (yarg)
	wr = 0.
	wi = 0.
	y1 = y + 1.5
	y2 = y1 * y1

	# Region II
	if (y < 0.85 && abs(x) > 18.1*y+1.65) {
	    if (abs(x) < 12)
		wr = exp (-x * x)
	    y3 = y + 3
	    do i = 1, 6 {
		r = x - t[i]
		r2 = r * r
		d = 1 / (r2 + y2)
		d1 = y1 * d
		d2 = r * d
		wr = wr + y * (c[i] * (r * d2 - 1.5 * d1) + s[i] * y3 * d2) /
		    (r2 + 2.25)
		r = x + t[i]
		r2 = r * r
		d = 1 / (r2 + y2)
		d3 = y1 * d
		d4 = r * d
		wr = wr + y * (c[i] * (r * d4 - 1.5 * d3) - s[i] * y3 * d4) /
		    (r2 + 2.25)
		wi = wi + c[i] * (d2 + d4) + s[i] * (d1 - d3)
	    }

	# Region I
	} else {
	    do i = 1, 6 {
		r = x - t[i]
		d = 1 / (r * r + y2)
		d1 = y1 * d
		d2 = r * d
		r = x + t[i]
		d = 1 / (r * r + y2)
		d3 = y1 * d
		d4 = r * d
		wr = wr + c[i] * (d1 + d3) - s[i] * (d2 - d4)
		wi = wi + c[i] * (d2 + d4) + s[i] * (d1 - d3)
	    }
	}
end
include	<math.h>
include	<mach.h>

# Profile types.
define	GAUSS		1	# Gaussian profile
define	LORENTZ		2	# Lorentzian profile
define	VOIGT		3	# Voigt profile

# Elements of fit array.
define	BKG		1	# Background
define	POS		2	# Position
define	INT		3	# Intensity
define	GAU		4	# Gaussian FWHM
define	LOR		5	# Lorentzian FWHM

# Type of constraints.
define	FIXED		1	# Fixed parameter
define	SINGLE		2	# Fit a single value for all lines
define	INDEP		3	# Fit independent values for all lines


# DOFIT -- Fit line profiles.  This is an interface to DOFIT1
# which puts parameters into the required form and vice-versa.
# It also implements a constrained approach to the solution.

procedure dofit (fit, x, y, s, npts, dx, nsub, y1, dy,
    xp, yp, gp, lp, tp, np, chisq)

int	fit[5]		# Fit constraints
real	x[npts]		# X data
real	y[npts]		# Y data
real	s[npts]		# Sigma data
int	npts		# Number of points
real	dx		# Pixel size
int	nsub		# Number of subpixels
real	y1		# Continuum offset
real	dy		# Continuum slope
real	xp[np]		# Profile positions
real	yp[np]		# Profile intensities
real	gp[np]		# Profile Gaussian FWHM
real	lp[np]		# Profile Lorentzian FWHM
int	tp[np]		# Profile type	
int	np		# Number of profiles
real	chisq		# Chi squared

int	i, j, fit1[5]
pointer	sp, a, b
errchk	dofit1

begin
	call smark (sp)
	call salloc (a, 8 + 5 * np, TY_REAL)

	# Convert positions and widths relative to first component.
	Memr[a] = dx
	Memr[a+1] = nsub
	Memr[a+2] = y1
	Memr[a+3] = dy
	Memr[a+4] = yp[1]
	Memr[a+5] = xp[1]
	Memr[a+6] = gp[1]
	Memr[a+7] = lp[1]
	do i = 1, np {
	    b = a + 5 * i + 3
	    Memr[b]   = yp[i] / Memr[a+4]
	    Memr[b+1] = xp[i] - Memr[a+5]
	    switch (tp[i]) {
	    case GAUSS:
		if (Memr[a+6] == 0.)
		    Memr[a+6] = gp[i]
		Memr[b+2] = gp[i] / Memr[a+6]
	    case LORENTZ:
		if (Memr[a+7] == 0.)
		    Memr[a+7] = lp[i]
		Memr[b+3] = lp[i] / Memr[a+7]
	    case VOIGT:
		if (Memr[a+6] == 0.)
		    Memr[a+6] = gp[i]
		Memr[b+2] = gp[i] / Memr[a+6]
		if (Memr[a+7] == 0.)
		    Memr[a+7] = lp[i]
		Memr[b+3] = lp[i] / Memr[a+7]
	    }
	    Memr[b+4] = tp[i]
	}

	# Do fit.
	fit1[INT] = fit[INT]
	do i = 1, fit[BKG] {
	    fit1[BKG] = i
	    fit1[GAU] = min (SINGLE, fit[GAU])
	    fit1[LOR] = min (SINGLE, fit[LOR])
	    do j = FIXED, fit[POS] {
		fit1[POS] = j
		if (np > 1 || j != INDEP)
		    call dofit1 (fit1, x, y, s, npts, Memr[a], np, chisq)
	    }
	    if (np > 1 && (fit[GAU] == INDEP || fit[LOR] == INDEP)) {
		fit1[GAU] = fit[GAU]
		fit1[LOR] = fit[LOR]
		call dofit1 (fit1, x, y, s, npts, Memr[a], np, chisq)
	    }
	}

	y1 = Memr[a+2]
	dy = Memr[a+3]
	do i = 1, np {
	    b = a + 5 * i + 3
	    yp[i] = Memr[b] * Memr[a+4]
	    xp[i] = Memr[b+1] + Memr[a+5]
	    switch (tp[i]) {
	    case GAUSS:
		gp[i] = abs (Memr[b+2] * Memr[a+6])
	    case LORENTZ:
		lp[i] = abs (Memr[b+3] * Memr[a+7])
	    case VOIGT:
		gp[i] = abs (Memr[b+2] * Memr[a+6])
		lp[i] = abs (Memr[b+3] * Memr[a+7])
	    }
	}

	call sfree (sp)
end


# DOREFIT -- Refit line profiles.  This assumes the input is very close
# to the final solution and minimizes the number of calls to the
# fitting routines.  This is intended for efficient use in the
# in computing bootstrap error estimates.

procedure dorefit (fit, x, y, s, npts, dx, nsub, y1, dy,
    xp, yp, gp, lp, tp, np, chisq)

int	fit[5]		# Fit constraints
real	x[npts]		# X data
real	y[npts]		# Y data
real	s[npts]		# Sigma data
int	npts		# Number of points
real	dx		# Pixel size
int	nsub		# Number of subpixels
real	y1		# Continuum offset
real	dy		# Continuum slope
real	xp[np]		# Profile positions
real	yp[np]		# Profile intensities
real	gp[np]		# Profile Gaussian FWHM
real	lp[np]		# Profile Lorentzian FWHM
int	tp[np]		# Profile type	
int	np		# Number of profiles
real	chisq		# Chi squared

int	i
pointer	sp, a, b
errchk	dofit1

begin
	call smark (sp)
	call salloc (a, 8 + 5 * np, TY_REAL)

	# Convert positions and widths relative to first component.
	Memr[a] = dx
	Memr[a+1] = nsub
	Memr[a+2] = y1
	Memr[a+3] = dy
	Memr[a+4] = yp[1]
	Memr[a+5] = xp[1]
	Memr[a+6] = gp[1]
	Memr[a+7] = lp[1]
	do i = 1, np {
	    b = a + 5 * i + 3
	    Memr[b]   = yp[i] / Memr[a+4]
	    Memr[b+1] = xp[i] - Memr[a+5]
	    switch (tp[i]) {
	    case GAUSS:
		Memr[b+2] = gp[i] / Memr[a+6]
	    case LORENTZ:
		Memr[b+3] = lp[i] / Memr[a+7]
	    case VOIGT:
		Memr[b+2] = gp[i] / Memr[a+6]
		Memr[b+3] = lp[i] / Memr[a+7]
	    }
	    Memr[b+4] = tp[i]
	}

	# Do fit.
	call dofit1 (fit, x, y, s, npts, Memr[a], np, chisq)

	y1 = Memr[a+2]
	dy = Memr[a+3]
	do i = 1, np {
	    b = a + 5 * i + 3
	    yp[i] = Memr[b] * Memr[a+4]
	    xp[i] = Memr[b+1] + Memr[a+5]
	    switch (tp[i]) {
	    case GAUSS:
		gp[i] = abs (Memr[b+2] * Memr[a+6])
	    case LORENTZ:
		lp[i] = abs (Memr[b+3] * Memr[a+7])
	    case VOIGT:
		gp[i] = abs (Memr[b+2] * Memr[a+6])
		lp[i] = abs (Memr[b+3] * Memr[a+7])
	    }
	}

	call sfree (sp)
end


# MODEL -- Compute model.

real procedure model (x, dx, nsub, xp, yp, gp, lp, tp, np)

real	x		# X value to be evaluated
real	dx		# Pixel width
int	nsub		# Number of subpixels
real	xp[np]		# Profile positions
real	yp[np]		# Profile intensities
real	gp[np]		# Profile Gaussian FWHM
real	lp[np]		# Profile Lorentzian FWHM
int	tp[np]		# Profile type	
int	np		# Number of profiles

int	i, j
real	delta, x1, y, arg1, arg2, v, v0, u

begin
	delta = dx / nsub
	x1 = x - (dx + delta) / 2
	y = 0.
	do j = 1, nsub {
	    x1 = x1 + delta
	    do i = 1, np {
		switch (tp[i]) {
		case GAUSS:
		    arg1 = 1.66511 * abs ((x1 - xp[i]) / gp[i])
		    if (arg1 < 5.)
			y = y + yp[i] * exp (-arg1**2)
		case LORENTZ:
		    arg2 = abs ((x1 - xp[i]) / (lp[i] / 2))
		    y = y + yp[i] / (1 + arg2**2)
		case VOIGT:
		    arg1 = 1.66511 * (x1 - xp[i]) / gp[i]
		    arg2 = 0.832555 * lp[i] / gp[i]
		    call voigt (0., arg2, v0, u)
		    call voigt (arg1, arg2, v, u)
		    y = y + yp[i] * v / v0
		}
	    }
	}
	y = y / nsub
	return (y)
end


# DERIVS -- Compute model and derivatives for MR_SOLVE procedure.
# This could be optimized more for the Voigt profile by reversing
# the do loops since v0 need only be computed once per line.

procedure derivs (x, a, y, dyda, na)

real	x		# X value to be evaluated
real	a[na]		# Parameters
real	y		# Function value
real	dyda[na]	# Derivatives
int	na		# Number of parameters

int	i, j, nsub
real	dx, dx1, delta, x1, wg, wl, arg1, arg2, I0, dI, c, u, v, v0

begin
	dx = a[1]
	nsub = a[2]
	delta = dx / nsub
	dx1 = .1 * delta
	x1 = x - (dx + delta) / 2
	y = 0.
	do i = 1, na
	   dyda[i] = 0.
	do j = 1, nsub {
	    x1 = x1 + delta
	    y = y + a[3] + a[4] * x1
	    dyda[3] = dyda[3] + 1.
	    dyda[4] = dyda[4] + x1
	    do i = 9, na, 5 {
		switch (a[i+4]) {
		case GAUSS:
		    I0 = a[5] * a[i]
		    wg = a[7] * a[i+2]
		    arg1 = 1.66511 * (x1 - a[6] - a[i+1]) / wg
		    if (abs (arg1) < 5.) {
			dI = exp (-arg1**2)
			c = I0 * dI * arg1
			y = y + I0 * dI
			dyda[5] = dyda[5] + a[i] * dI
			dyda[6] = dyda[6] + c / wg
			dyda[7] = dyda[7] + c * arg1 / a[7]
			dyda[i] = dyda[i] + a[5] * dI
			dyda[i+1] = dyda[i+1] + c / wg
			dyda[i+2] = dyda[i+2] + c * arg1 / a[i+2]
		    }
		case LORENTZ:
		    I0 = a[5] * a[i]
		    wl = (a[8] * a[i+3] / 2)
		    arg2 = (x1 - a[6] - a[i+1]) / wl
		    dI = 1 / (1 + arg2**2)
		    c = 2 * I0 * dI * dI * arg2
		    y = y + I0 * dI
		    dyda[5] = dyda[5] + a[i] * dI
		    dyda[6] = dyda[6] + c / wl
		    dyda[8] = dyda[8] + c * arg2 / a[8]
		    dyda[i] = dyda[i] + a[5] * dI
		    dyda[i+1] = dyda[i+1] + c / wl
		    dyda[i+3] = dyda[i+3] + c * arg2 / a[i+3]
		case VOIGT:
		    a[7] = max (dx1, abs(a[7]))
		    a[8] = max (dx1, abs(a[8]))
		    a[i+2] = max (1E-6, abs(a[i+2]))
		    a[i+3] = max (1E-6, abs(a[i+3]))

		    I0 = a[5] * a[i]
		    wg = a[7] * a[i+2]
		    wl = a[8] * a[i+3]
		    arg1 = 1.66511 * (x1 - a[6] - a[i+1]) / wg
		    arg2 = 0.832555 * wl / wg
		    call voigt (0., arg2, v0, u)
		    call voigt (arg1, arg2, v, u)
		    v = v / v0; u = u / v0
		    dI = (1 - v) / (v0 * SQRTOFPI)
		    c = 2 * I0 * arg2
		    y = y + I0 * v
		    dyda[5] = dyda[5] + a[i] * v
		    dyda[6] = dyda[6] + 2 * c * (arg1 * v - arg2 * u) / wl
		    dyda[7] = dyda[7] +
			c * (dI + arg1 * (arg1 / arg2 * v - 2 * u)) / a[7]
		    dyda[8] = dyda[8] + c * (arg1 * u - dI) / a[8]
		    dyda[i] = dyda[i] + a[5] * v
		    dyda[i+1] = dyda[i+1] + 2 * c * (arg1 * v - arg2 * u) / wl
		    dyda[i+2] = dyda[i+2] +
			c * (dI + arg1 * (arg1 / arg2 * v - 2 * u)) / a[i+2]
		    dyda[i+3] = dyda[i+3] + c * (arg1 * u - dI) / a[i+3]
		}
	    }
	}
	y = y / nsub
	do i = 1, na
	   dyda[i] = dyda[i] / nsub
end


# DOFIT1 -- Perform nonlinear iterative fit for the specified parameters.
# This uses the Levenberg-Marquardt method from NUMERICAL RECIPES.

procedure dofit1 (fit, x, y, s, npts, a, nlines, chisq)

int	fit[5]		# Fit constraints
real	x[npts]		# X data
real	y[npts]		# Y data
real	s[npts]		# Sigma data
int	npts		# Number of points
real	a[ARB]		# Fitting parameters
int	nlines		# Number of lines
real	chisq		# Chi squared

int	i, np, nfit
real	mr, chi2
pointer	sp, flags, ptr
errchk	mr_solve

begin
	# Number of terms is 5 for each line plus common background, center,
	# intensity and widths.  Also the pixel size and number of subpixels.

	np = 5 * nlines + 8

	call smark (sp)
	call salloc (flags, np, TY_INT)
	ptr = flags

	# Background.
	switch (fit[BKG]) {
	case SINGLE:
	    Memi[ptr] = 3
	    Memi[ptr+1] = 4
	    ptr = ptr + 2
	}

	# Peaks.
	switch (fit[INT]) {
	case SINGLE:
	    Memi[ptr] = 5
	    ptr = ptr + 1
	case INDEP:
	    do i = 1, nlines {
		Memi[ptr] = 5 * i + 4
		ptr = ptr + 1
	    }
	}

	# Positions.
	switch (fit[POS]) {
	case SINGLE:
	    Memi[ptr] = 6
	    ptr = ptr + 1
	case INDEP:
	    do i = 1, nlines {
		Memi[ptr] = 5 * i + 5
		ptr = ptr + 1
	    }
	}

	# Gaussian FWHM.
	switch (fit[GAU]) {
	case SINGLE:
	    Memi[ptr] = 7
	    ptr = ptr + 1
	case INDEP:
	    do i = 1, nlines {
		Memi[ptr] = 5 * i + 6
		ptr = ptr + 1
	    }
	}

	# Lorentzian FWHM.
	switch (fit[LOR]) {
	case SINGLE:
	    Memi[ptr] = 8
	    ptr = ptr + 1
	case INDEP:
	    do i = 1, nlines {
		Memi[ptr] = 5 * i + 7
		ptr = ptr + 1
	    }
	}

	nfit = ptr - flags
	mr = -1.
	i = 0
	chi2 = MAX_REAL
	repeat {
	    call mr_solve (x, y, s, npts, a, Memi[flags], np, nfit, mr, chisq)
	    if (chi2 - chisq > 0.0001)
		i = 0
	    else
		i = i + 1
	    chi2 = chisq
	} until (i == 5)

	mr = 0.
	call mr_solve (x, y, s, npts, a, Memi[flags], np, nfit, mr, chisq)

	call sfree (sp)
end


# MR_SOLVE -- Levenberg-Marquardt nonlinear chi square minimization.
#
# Use the Levenberg-Marquardt method to minimize the chi squared of a set
# of paraemters.  The parameters being fit are indexed by the flag array.
# To initialize the Marquardt parameter, MR, is less than zero.  After that
# the parameter is adjusted as needed.  To finish set the parameter to zero
# to free memory.  This procedure requires a subroutine, DERIVS, which
# takes the derivatives of the function being fit with respect to the
# parameters.  There is no limitation on the number of parameters or
# data points.  For a description of the method see NUMERICAL RECIPES
# by Press, Flannery, Teukolsky, and Vetterling, p523.

procedure mr_solve (x, y, s, npts, params, flags, np, nfit, mr, chisq)

real	x[npts]			# X data array
real	y[npts]			# Y data array
real	s[npts]			# Sigma data array
int	npts			# Number of data points
real	params[np]		# Parameter array
int	flags[np]		# Flag array indexing parameters to fit
int	np			# Number of parameters
int	nfit			# Number of parameters to fit
real	mr			# MR parameter
real	chisq			# Chi square of fit

int	i
real	chisq1
pointer	new, a1, a2, delta1, delta2

errchk	mr_invert

begin
	# Allocate memory and initialize.
	if (mr < 0.) {
	    call mfree (new, TY_REAL)
	    call mfree (a1, TY_REAL)
	    call mfree (a2, TY_REAL)
	    call mfree (delta1, TY_REAL)
	    call mfree (delta2, TY_REAL)

	    call malloc (new, np, TY_REAL)
	    call malloc (a1, nfit*nfit, TY_REAL)
	    call malloc (a2, nfit*nfit, TY_REAL)
	    call malloc (delta1, nfit, TY_REAL)
	    call malloc (delta2, nfit, TY_REAL)

	    call amovr (params, Memr[new], np)
	    call mr_eval (x, y, s, npts, Memr[new], flags, np, Memr[a2],
	        Memr[delta2], nfit, chisq)
	    mr = 0.001
	}

	# Restore last good fit and apply the Marquardt parameter.
	call amovr (Memr[a2], Memr[a1], nfit * nfit)
	call amovr (Memr[delta2], Memr[delta1], nfit)
	do i = 1, nfit
	    Memr[a1+(i-1)*(nfit+1)] = Memr[a2+(i-1)*(nfit+1)] * (1. + mr)

	# Matrix solution.
	call mr_invert (Memr[a1], Memr[delta1], nfit)

	# Compute the new values and curvature matrix.
	do i = 1, nfit
	    Memr[new+flags[i]-1] = params[flags[i]] + Memr[delta1+i-1]
	call mr_eval (x, y, s, npts, Memr[new], flags, np, Memr[a1],
	    Memr[delta1], nfit, chisq1)

	# Check if chisq has improved.
	if (chisq1 < chisq) {
	    mr = max (EPSILONR, 0.1 * mr)
	    chisq = chisq1
	    call amovr (Memr[a1], Memr[a2], nfit * nfit)
	    call amovr (Memr[delta1], Memr[delta2], nfit)
	    call amovr (Memr[new], params, np)
	} else
	    mr = 10. * mr

	if (mr == 0.) {
	    call mfree (new, TY_REAL)
	    call mfree (a1,  TY_REAL)
	    call mfree (a2,  TY_REAL)
	    call mfree (delta1, TY_REAL)
	    call mfree (delta2, TY_REAL)
	}
end


# MR_EVAL -- Evaluate curvature matrix.  This calls procedure DERIVS.

procedure mr_eval (x, y, s, npts, params, flags, np, a, delta, nfit, chisq)

real	x[npts]			# X data array
real	y[npts]			# Y data array
real	s[npts]			# Sigma data array
int	npts			# Number of data points
real	params[np]		# Parameter array
int	flags[np]		# Flag array indexing parameters to fit
int	np			# Number of parameters
real	a[nfit,nfit]		# Curvature matrix
real	delta[nfit]		# Delta array
int	nfit			# Number of parameters to fit
real	chisq			# Chi square of fit

int	i, j, k
real	ymod, dy, dydpj, dydpk, sig2i
pointer	sp, dydp

begin
	call smark (sp)
	call salloc (dydp, np, TY_REAL)

	do j = 1, nfit {
	   do k = 1, j
	       a[j,k] = 0.
	    delta[j] = 0.
	}

	chisq = 0.
	do i = 1, npts {
	    call derivs (x[i], params, ymod, Memr[dydp], np)
	    if (IS_INDEF(ymod))
		next
	    sig2i = 1. / (s[i] * s[i])
	    dy = y[i] - ymod
	    do j = 1, nfit {
		dydpj = Memr[dydp+flags[j]-1] * sig2i
		delta[j] = delta[j] + dy * dydpj
		do k = 1, j {
		    dydpk = Memr[dydp+flags[k]-1]
		    a[j,k] = a[j,k] + dydpj * dydpk
		}
	    }
	    chisq = chisq + dy * dy * sig2i
	}

	do j = 2, nfit
	    do k = 1, j-1
		a[k,j] = a[j,k]

	call sfree (sp)
end
	    

# MR_INVERT -- Solve a set of linear equations using Householder transforms.

procedure mr_invert (a, b, n)

real	a[n,n]		# Input matrix and returned inverse
real	b[n]		# Input RHS vector and returned solution
int	n		# Dimension of input matrices

int	krank
real	rnorm
pointer	sp, h, g, ip

begin
	call smark (sp)
	call salloc (h, n, TY_REAL)
	call salloc (g, n, TY_REAL)
	call salloc (ip, n, TY_INT)

	call hfti (a, n, n, n, b, n, 1, 1E-10, krank, rnorm,
	    Memr[h], Memr[g], Memi[ip])

	call sfree (sp)
end