include	"mkw.h"

# Memory management
define	Line			Memd[line+$1-1]
define	Peak			Memd[peak+$1-1]
define	Sigma			Memd[sigma+$1-1]

#---------------------------------------------------------------------------
.help mkw_mkwave 11Apr95 source
.ih
NAME
.nf
mkw_mkwave -- Make a spectrum.
mkw_near -- Find pixel nearest a specified wavelength.
.fi
.ih
DESCRIPTION
Below is a description of the routine interface
.ls mkw_mkwave (mkw, seed)
Create a spectrum based on the parameters found in the MKW object.
.ls ARGUMENTS
.ls mkw (i: pointer)
Pointer to the MKW object, created by mkw_alloc.
.le
.ls seed (i: long)
The seed for the random number generator
.le
.le
.ls ERRORS
Errors generated by the routines malloc.
.le
.le
.ls pointer mkw_near (warray, w, n)
Find the pixel nearest a specified wavelength in an array.  Use linear
interpolation to find sub-pixel positions.
.ls ARGUMENTS
.ls warray (i: double[n])
Array containing wavelengths.  The wavelengths must be monotonic,
though whether it is decreasing or increasing is not important.
.le
.ls w (i: double)
Wavelength to find in the array.
.le
.ls n (i: int)
Length of the warray array.
.le
.le
.ls RETURNS
The pixel position of the specified wavelength in the array.
.le
.ls ERRORS
Error #1 - wavelength is outside of array
.le
.le
.endhelp
#---------------------------------------------------------------------------
procedure mkw_mkwave (mkw, seed)

pointer	mkw			# I:  MKW object.
long	seed			# I:  Random seed.

# Declarations.
double	d			# Local dispersion.
double	fnu			# Exponential conversion.
int	ip			# Integer pixel.
int	l			# Current line.
pointer line			# Relevant set of lines.
double	mkw_near()		# Find pixel position of a wavelength.
int	nl			# Number of lines in spectrum.
double	p			# Pixel position of center of line.
pointer	peak			# Peak of line.
double	s, c			# Gaussian parameters.
pointer	sigma			# Width of lines.
real	tseed			# Seed holder.
real	urand()			# Get random number.
double	w			# Wavelength.
double	x, x1, x2		# Generic.

errchk	malloc

begin
	tseed = seed

	# Redefine the wavelengths, strengths, and sigmas necessary
	# for creating lines.
	call malloc (line, MKW_NLINES(mkw), TY_DOUBLE)
	call malloc (peak, MKW_NLINES(mkw), TY_DOUBLE)
	call malloc (sigma, MKW_NLINES(mkw), TY_DOUBLE)
	nl = 0
	do l = 1, MKW_NLINES(mkw) {

	    # If line is not in the domain, forget it.
	    w = MKW_LINES(mkw,l) * MKW_Z(mkw)
	    if (w < min (MKW_WAVE(mkw,1), MKW_WAVE(mkw,MKW_NPIX(mkw))) ||
		w > max (MKW_WAVE(mkw,1), MKW_WAVE(mkw,MKW_NPIX(mkw))))
		next

	    # Add to list.
	    nl = nl + 1
	    Line(nl) = w

	    # Get Strength.
	    if (MKW_INT_PTR(mkw) != NULL)
		Peak(nl) = MKW_INT(mkw,l)
	    else
		Peak(nl) = MKW_PEAK(mkw) * urand (tseed)
	    Peak(nl) = Peak(nl) / MKW_Z(mkw) * MKW_SUBSAMPLE(mkw)

	    # Get Sigma
	    if (MKW_FWHM_PTR(mkw) != NULL)
		Sigma(nl) = MKW_FWHM(mkw,l)
	    else
		Sigma(nl) = MKW_SIGMA(mkw)
	    Sigma(nl) = Sigma(nl) * MKW_Z(mkw)
	}
	
	# Make the lines in the list.
	do l = 1, nl {

	    # Find pixel position of line and linear dispersion.
	    iferr (p = mkw_near (MKW_WAVE(mkw,1), Line(l), MKW_NPIX(mkw)))
		next
	    ip = p + 0.5d0
	    if (ip == MKW_NPIX(mkw) || ip == 1)
		next

	    # Create the line.
	    s = Sigma(l) / (MKW_WAVE(mkw,ip+1) - MKW_WAVE(mkw,ip))
	    c = -0.5 / s**2
	    x1 = max (1.0d0, p - MKW_NSIGMA(mkw) * s)
	    x2 = min (double (MKW_NPIX(mkw)), p + MKW_NSIGMA(mkw) * s)
	    for (x = x1; x <= x2; x = x+MKW_SUBSAMPLE(mkw)) {
		ip = x + 0.5
		MKW_SPEC(mkw,ip) = MKW_SPEC(mkw,ip) + Peak(l) *
				   exp (c * (x-p)**2)
	    }
	}
	
	# Now add the overall continuum.
	fnu = 5
	if (MKW_FNU(mkw) == YES)
	    fnu = 3

	if (MKW_TEMP(mkw) > 0.d0) {
	    w = MKW_WAVE(mkw,1) * 1.0d-8
	    x1 = exp (1.4388 / (w * MKW_TEMP(mkw)))
	    x2 = w**fnu * (x1-1.d0)
	}

	do ip = 1, MKW_NPIX(mkw) {
	    w = MKW_WAVE(mkw,ip) / MKW_Z(mkw)
	    if (ip == 1)
		d = (MKW_WAVE(mkw,2) - MKW_WAVE(mkw,1))
	    else
		d = (MKW_WAVE(mkw,ip) - MKW_WAVE(mkw,1)) / (ip-1)
	    x = MKW_CONT(mkw) + MKW_SLOPE(mkw) / d * (w - MKW_WAVE(mkw,1))

	    if (MKW_TEMP(mkw) > 0.d0) {
		w = w * 1.0d-8
		x1 = exp (1.4388 / (w * MKW_TEMP(mkw)))
		x = x * (x2 / w**fnu / (x1 - 1.d0))
	    }

	    if (x > 0.d0)
		MKW_SPEC(mkw,ip) = max (0.d0, x * (1.0 + MKW_SPEC(mkw,ip)))
	}

	# That's all folks.
	call mfree (sigma, TY_DOUBLE)
	call mfree (peak, TY_DOUBLE)
	call mfree (line, TY_DOUBLE)
end
#---------------------------------------------------------------------------
# End of mkw_mkwave
#---------------------------------------------------------------------------
double procedure mkw_near (warray, w, n)

double	warray[n]		# I:  Wavelength array.
double	w			# I:  Wavelength to find.
int	n			# I:  Size of array.

# Declarations.
int	high			# Pixel of greater wavelength.
int	low			# Pixel of lesser wavelength.
int	mid			# Mid pixel of current search range.
int	np			# Pixels in the search range.
double	p			# Closest pixel to specified wavelength.

begin
	# Now find the two pixels the specified wavelength is between.
	np = n
	if (warray[1] < warray[np]) {
	    if (w < warray[1] || w > warray[np])
		call error (1, "wavelength is outside of array")
	    low = 1
	    high = np
	} else {
	    if (w > warray[1] || w < warray[np])
		call error (1, "wavelength is outside of array")
	    low = np
	    high = 1
	}
	while (np > 2) {
	    mid = (low + high) / 2
	    if (warray[mid] < w)
		low = mid
	    else
		high = mid
	    np = abs (high - low) + 1
	}

	# Linearly interpolate for the subpixel.
	p = low + ((w - warray[low]) / (warray[high] - warray[low]))

	# That's all folks.
	return (p)
end
#---------------------------------------------------------------------------
# End of mkw_near
#---------------------------------------------------------------------------