include	<evexpr.h>
include	<imhdr.h>

define	SZ_UNITS	8

# COMPSPEC -- Create a one dimensional image from an input function
#
# B.Simon	01-Jun-87	Original

procedure compspec()

#--
include	"compspec.com"

char	spectrum[SZ_FNAME]		# Name of the file to be created
char	expr[SZ_LINE]			# Algebraic expression to be evaluated
int	npoints				# Number of points in output spectrun
real	fstwave				# Wavelength of first point in spectrum
real	lstwave				# Wavelength of last point in spectrum
char	units[SZ_UNITS]			# Units of wavelength

int	ipoint
long	aeon
pointer	im, buffer, op
real	p, value, minpix, maxpix, delta

data	aeon	/ 0 /

int	clgeti(), locpr(), imaccf()
long	clktime()
pointer	immap(), impl1r(), evexpr()
real	clgetr()

extern	get_specpt()

begin
	# Read the parameter file

	call clgstr ("spectrum",spectrum,SZ_FNAME)
	call clgstr ("expr",expr,SZ_LINE)
	npoints = clgeti ("npoints")
	fstwave = clgetr ("fstwave")
	lstwave = clgetr ("lstwave")
	call clgstr ("units",units,SZ_UNITS)

	# Case conversion of input parameters for case insensitivity

	call strlwr (expr)
	call strupr (units)

	# Open the output image

	im = immap (spectrum,NEW_IMAGE,NULL)
	IM_NDIM(im) = 1
	IM_LEN(im,1) = npoints

	buffer = impl1r (im)

	# Initialize flags passed thru common block

	first_point = true
	first_var = true

	do ipoint = 0, npoints-1 {

	    # Compute the wavelength at current point in spectrum

	    p = real(ipoint) / real(npoints-1)
	    wave = fstwave * (1.0 - p) + lstwave * p

	    # Evaluate the function at this point

	    op = evexpr (expr, locpr(get_specpt), NULL)

	    # Retrieve the result from the structure

	    switch (O_TYPE(op)) {
	    case TY_BOOL :
		call error(1,"Boolean expressions not allowed")
	    case TY_CHAR :
		call error(1,"Character expressions not allowed")
	    case TY_INT :
		value = O_VALI(op)
	    case TY_REAL :
		value = O_VALR(op)
	    }

	    Memr[buffer+ipoint] = value

	    # Compute minimum and maximum of spectrum

	    if (first_point) {
		minpix = value
		maxpix = value
	    } else {
		minpix = min(value,minpix)
		maxpix = max(value,maxpix)
	    }

	    first_point = false
	}

	# Copy minimum and maximum to spectrum

	IM_MIN(im) = minpix
	IM_MAX(im) = maxpix
	IM_LIMTIME(im) = clktime (aeon)

	# Compute the world coordinates and copy to spectrum

	if (imaccf (im, "CRPIX1") == YES)
	    call imputr (im,"CRPIX1",1.0)
	if (imaccf (im, "CRVAL1") == YES)
	    call imputr (im,"CRVAL1",fstwave)
	if (imaccf (im, "CTYPE1") == YES)
	    call impstr (im,"CTYPE1",units)

	delta = (lstwave - fstwave) / real(npoints - 1)
	if (imaccf (im, "CD1_1") == YES)
	   call imputr (im,"CD1_1",delta)

	call imunmap (im)
end