include	<mach.h>
include	"refspectra.h"


# REFINTERP -- Assign reference spectra to interpolate between based on sort
# key.  The nearest preceding and following spectra are assigned weights based
# on their distance.  If there is no preceding and following spectrum then
# the nearest spectrum is assigned.

procedure refinterp (input, refs)

pointer	input			# List of input spectra
pointer	refs			# List of reference spectra

bool	ignoreaps		# Ignore apertures?

int	i, i1, i2, nrefs, ap
double	sortval, d, d1, d2
real	wt1, wt2
pointer	sp, image, gval, refimages, refaps, refvals, refgvals

bool	clgetb(), streq(), refginput(), refgref()
int	imtgetim(), imtlen()

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

	# Task parameters
	ignoreaps = clgetb ("ignoreaps")

	# Tabulate reference spectra.  This expands the reference list,
	# checks the spectrum is a reference spectrum of the appropriate
	# aperture.

	call salloc (refimages, imtlen (refs), TY_INT)
	call salloc (refaps, imtlen (refs), TY_INT)
	call salloc (refvals, imtlen (refs), TY_DOUBLE)
	call salloc (refgvals, imtlen (refs), TY_INT)
	nrefs = 0
	while (imtgetim (refs, Memc[image], SZ_FNAME) != EOF) {
	    call refnoextn (Memc[image])
	    if (!refgref (Memc[image], ap, sortval, gval))
		next

	    for (i=0; i<nrefs; i=i+1)
		if (streq (Memc[Memi[refimages+i]], Memc[image]))
		    break
	    if (i == nrefs) {
		call salloc (Memi[refimages+nrefs], SZ_FNAME, TY_CHAR)
		call salloc (Memi[refgvals+nrefs], SZ_FNAME, TY_CHAR)
		call strcpy (Memc[image], Memc[Memi[refimages+i]], SZ_FNAME)
		Memi[refaps+i] = ap
		Memd[refvals+i] = sortval
		call strcpy (Memc[gval], Memc[Memi[refgvals+i]], SZ_FNAME)
		nrefs = i + 1
	    }
	}
	if (nrefs < 1)
	    call error (0, "No reference images specified")


	# Assign following reference spectra to each input spectrum.
	# Skip input spectra which are not of the appropriate aperture
	# or have been assigned previously (unless overriding).

	while (imtgetim (input, Memc[image], SZ_FNAME) != EOF) {
	    call refnoextn (Memc[image])
	    if (!refginput (Memc[image], ap, sortval, gval))
		next

	    i1 = 0
	    i2 = 0
	    d1 = MAX_REAL
	    d2 = -MAX_REAL
	    do i = 1, nrefs {
		if (!streq (Memc[gval], Memc[Memi[refgvals+i-1]]))
		    next
	        if (!ignoreaps  && ap != Memi[refaps+i-1])
		    next
	        d = sortval - Memd[refvals+i-1]
	        if ((d >= 0.) && (d < d1)) {
		    i1 = i
		    d1 = d
	        } else if ((d <= 0.) && (d > d2)) {
		    i2 = i
		    d2 = d
	        }
	    }

	    if (i1 > 0 && i2 > 0) {	# Weight spectra
		if (d1 - d2 == 0.) {
		    wt1 = 0.5
		    wt2 = 0.5
		} else {
		    wt1 = -d2 / (d1 - d2)
		    wt2 =  d1 / (d1 - d2)
		}
		call refspectra (Memc[image], Memc[Memi[refimages+i1-1]], wt1,
		    Memc[Memi[refimages+i2-1]], wt2)
	    } else if (i1 > 0)	# Nearest preceding spectrum
		call refspectra (Memc[image], Memc[Memi[refimages+i1-1]], 1.,
		    Memc[Memi[refimages+i1-1]], 0.)
	    else if (i2 > 0)	# Nearest following spectrum
		call refspectra (Memc[image], Memc[Memi[refimages+i2-1]], 1.,
		    Memc[Memi[refimages+i2-1]], 0.)
	    else {		# No reference spectrum found
		call refprint (STDERR, NO_REFSPEC, Memc[image], "", "", "",
		    ap, 0, "")
		do i = 1, nrefs {
		    if (!streq (Memc[gval], Memc[Memi[refgvals+i-1]])) {
			call refprint (STDERR, REF_GROUP, Memc[image],
			    Memc[Memi[refimages+i-1]], Memc[gval],
			    Memc[Memi[refgvals+i-1]], ap, Memi[refaps+i-1], "")
			next
		    }
		    if (!ignoreaps  && ap != Memi[refaps+i-1])
			call refprint (STDERR, REF_AP, Memc[image],
			    Memc[Memi[refimages+i-1]], Memc[gval],
			    Memc[Memi[refgvals+i-1]], ap, Memi[refaps+i-1], "")
			next
		}
	    }
	}

	call sfree (sp)
end