include	<error.h>
include	<imhdr.h>
include	<imset.h>
include	<pkg/dttext.h>
include	<smw.h>
include	<units.h>
include	"dispcor.h"

# Symbol table structure for the dispersion solutions.
define	LEN_DC		11		# Length of dispersion solution struct.
define	DC_FORMAT	Memi[$1]	# Type of dispersion
define	DC_PAPS		Memi[$1+1]	# Pointer to aperture numbers
define	DC_PAPCEN	Memi[$1+2]	# Pointer to aperture centers
define	DC_PUN		Memi[$1+3]	# Pointer to units
define	DC_PSHIFT	Memi[$1+4]	# Pointer to shifts
define	DC_PCOEFF	Memi[$1+5]	# Pointer to coefficients
define	DC_NAPS		Memi[$1+6]	# Number of apertures
define  DC_OFFSET       Memi[$1+7]      # Aperture to order offset
define  DC_SLOPE        Memi[$1+8]      # Aperture to order slope
define  DC_COEFFS       Memi[$1+9]      # Dispersion coefficients
define  DC_SHIFT        Memr[$1+10]      # Dispersion function shift


# DC_OPEN    -- Initialize the dispersion data structures
# DC_CLOSE   -- Close the dispersion data structures
# DC_GMS     -- Get a multispec spectrum
# DC_GMSDB   -- Get a multispec dispersion database entry
# DC_REFSHFT -- Get a reference shift
# DC_GEC     -- Get an echelle spectrum
# DC_GECDB   -- Get an echelle dispersion database entry
# DC_ECMS    -- Convert echelle dispersion coeffs to multispec coeffs


# DC_OPEN -- Initialize the dispersion routines.  This consists
# of opening a symbol table for the dispersion solution functions.  A
# symbol table is used since the same dispersion reference (arc image)
# may be be used multiple times and the database access is slow.

procedure dc_open (stp, db)

pointer	stp			# Symbol table pointer
char	db[SZ_FNAME]		# Database name

pointer	sym, stopen(), stenter(), stpstr()

begin
	stp = stopen ("disp", 10, 10, 10*SZ_FNAME)
	sym = stenter (stp, "database", 1)
	Memi[sym] = stpstr (stp, db, 0)
end


# DC_CLOSE -- Close the dispersion data structures.

procedure dc_close (stp)

int	i
pointer stp, sym, sthead, stnext

begin
	# Close each dispersion function and then the symbol table.
	for (sym = sthead (stp); sym != NULL; sym = stnext (stp, sym)) {
	    if (DC_FORMAT(sym) == 1) {
		do i = 1, DC_NAPS(sym) {
		    call un_close (Memi[DC_PUN(sym)+i-1])
		    call mfree (Memi[DC_PCOEFF(sym)+i-1], TY_DOUBLE)
		}
		call mfree (DC_PAPS(sym), TY_INT)
		call mfree (DC_PAPCEN(sym), TY_REAL)
		call mfree (DC_PUN(sym), TY_POINTER)
		call mfree (DC_PSHIFT(sym), TY_DOUBLE)
		call mfree (DC_PCOEFF(sym), TY_POINTER)
	    } else if (DC_FORMAT(sym) == 2) {
		call un_close (DC_PUN(sym))
		call mfree (DC_COEFFS(sym), TY_DOUBLE)
	    }
	}
	call stclose (stp)
end


# DC_GMS -- Get a multispec spectrum.  This consists of mapping the image
# and setting a MWCS coordinate transformation.  If not dispersion corrected
# the dispersion function is found in the database for the reference
# spectra and set in the SMW.

procedure dc_gms (spec, im, smw, stp, ignoreaps, ap, fd1, fd2)

char	spec[ARB]	#I Spectrum name
pointer	im		#I IMIO pointer
pointer	smw		#I SMW pointer
pointer	stp		#I Dispersion symbol table
int	ignoreaps	#I Ignore aperture numbers?
pointer	ap		#O Aperture data structure
int	fd1		#I Logfile descriptor
int	fd2		#I Logfile descriptor

double	wt1, wt2, dval
int	i, j, k, k1, k2, l, dc, sfd, axis[2], naps, naps1, naps2, ncoeffs
pointer	sp, str1, str2, papcen, pshift, coeffs, ct1, ct2, un, un1, un2
pointer	paps1, paps2, punits1, punits2, pshift1, pshift2, pcoeff1, pcoeff2

bool	un_compare()
double	smw_c1trand()
int	imaccf(), nscan(), stropen()
pointer	smw_sctran(), un_open()
errchk	dc_gmsdb, dc_refshft, imgstr, smw_sctran, un_open

data	axis/1,2/
define	done_	90

begin
	call smark (sp)
	call salloc (str1, SZ_LINE, TY_CHAR)
	call salloc (str2, SZ_LINE, TY_CHAR)

	# Set WCS attributes
	naps = IM_LEN(im,2)
	call calloc (ap, LEN_AP(naps), TY_STRUCT)
	do i = 1, naps {
	    DC_PL(ap,i) = i
	    DC_CF(ap,i) = NULL
	    call smw_gwattrs (smw, DC_PL(ap,i), 1, DC_AP(ap,i), DC_BM(ap,i),
		DC_DT(ap,i), DC_W1(ap,i), DC_DW(ap,i), DC_NW(ap,i), DC_Z(ap,i),
		DC_LW(ap,i), DC_UP(ap,i), DC_CF(ap,i))
	    if (i == 1) {
		iferr (call mw_gwattrs (SMW_MW(smw,0), 1, "units", Memc[str1],
		    SZ_LINE))
		    Memc[str1] = EOS
		DC_UN(ap,i) = un_open (Memc[str1])
	    }
	    dc = DC_DT(ap,i)
	}

	# Check if the spectra have been dispersion corrected
	# by an earlier version of DISPCOR.  If so then don't allow
	# another database dispersion correction.  This assumes all
	# spectra have the same dispersion type.  Check for a
	# reference spectrum.

	if ((imaccf (im, "REFSPEC1") == NO) ||
	    (dc > -1 && imaccf (im, "DCLOG1") == NO)) {
	    if (fd1 != NULL) {
		call fprintf (fd1,
		    "%s: Resampling using current coordinate system\n")
		    call pargstr (spec)
	    }
	    if (fd2 != NULL) {
		call fprintf (fd2,
		    "%s: Resampling using current coordinate system\n")
		    call pargstr (spec)
	    }
	    goto done_
	}
	    
	# Get the reference spectra dispersion function from the database
	# and determine a reference shift.

	iferr {
	    call imgstr (im, "REFSPEC1", Memc[str1], SZ_LINE)
	    call sscan (Memc[str1])
	    call gargwrd (Memc[str1], SZ_LINE)
	    call gargd (wt1)
	    if (nscan() == 1)
		wt1 = 1.
	} then {
	    call strcpy (spec, Memc[str1], SZ_FNAME)
	    wt1 = 1.
	}
	iferr (call dc_gmsdb (Memc[str1], stp, paps1, papcen, punits1, pshift,
	    pcoeff1, naps1)) {
	    call sfree (sp)
	    call erract (EA_ERROR)
	}
	call salloc (pshift1, naps1, TY_DOUBLE)
	call amovd (Memd[pshift], Memd[pshift1], naps1)
	if (fd1 != NULL) {
	    call fprintf (fd1, "%s: REFSPEC1 = '%s %.8g'\n")
		call pargstr (spec)
		call pargstr (Memc[str1])
		call pargd (wt1)
	}
	if (fd2 != NULL) {
	    call fprintf (fd2, "%s: REFSPEC1 = '%s %.8g'\n")
		call pargstr (spec)
		call pargstr (Memc[str1])
		call pargd (wt1)
	}

	iferr (call  dc_refshft (spec, stp, Memc[str1], "REFSHFT1", im,
	    Memi[paps1], Memr[papcen], Memd[pshift1], naps1, fd1, fd2))
	    ;

	iferr {
            call imgstr (im, "REFSPEC2", Memc[str1], SZ_LINE)
            call sscan (Memc[str1])
            call gargwrd (Memc[str1], SZ_LINE)
            call gargd (wt2)
            if (nscan() == 1)
                wt2 = 1.
            call dc_gmsdb (Memc[str1], stp, paps2, papcen, punits2, pshift,
		pcoeff2, naps2)
            call salloc (pshift2, naps2, TY_DOUBLE)
            call amovd (Memd[pshift], Memd[pshift2], naps2)
	    if (fd1 != NULL) {
		call fprintf (fd1, "%s: REFSPEC2 = '%s %.8g'\n")
		    call pargstr (spec)
		    call pargstr (Memc[str1])
		    call pargd (wt2)
	    }
	    if (fd2 != NULL) {
		call fprintf (fd2, "%s: REFSPEC2 = '%s %.8g'\n")
		    call pargstr (spec)
		    call pargstr (Memc[str1])
		    call pargd (wt2)
	    }
	    iferr (call  dc_refshft (spec, stp, Memc[str1],
		"REFSHFT2", im, Memi[paps2], Memr[papcen], Memd[pshift2],
		naps2, fd1, fd2))
	        ;
	} then
	    wt2 = 0.

	# Adjust weights to unit sum.
	dval = wt1 + wt2
	wt1 = wt1 / dval
	wt2 = wt2 / dval

	# Enter dispersion function in the MWCS.
	do i = 1, naps {
	    j = DC_AP(ap,i)
	    for (k1=0; k1<naps1 && Memi[paps1+k1]!=j; k1=k1+1)
		;
	    if (k1 == naps1)
		for (k1=0; k1<naps1 && !IS_INDEFI(Memi[paps1+k1]); k1=k1+1)
		    ;
	    if (k1 == naps1) {
		if (ignoreaps == YES)
		    k1 = 0
		else {
		    call sprintf (Memc[str1], SZ_LINE,
			"%s - Missing reference for aperture %d")
			call pargstr (spec)
			call pargi (j)
		    call fatal (1, Memc[str1])
		}
	    }
	    un1 = Memi[punits1+k1]

	    # The following assumes some knowledge of the data structure in
	    # order to shortten the the attribute string.
	    coeffs = Memi[pcoeff1+k1]
	    if (coeffs == NULL) {
		if (DC_DT(ap,i) == 2) {
		    sfd = NULL
		    if (wt2 <= 0.)
			call sshift1 (Memd[pshift1+k1], DC_CF(ap,i))
		} else {
		    ncoeffs = 6
		    l = 20 * (ncoeffs + 2)
		    if (wt2 > 0.)
			l = 2 * l
		    call realloc (DC_CF(ap,i), l, TY_CHAR)
		    call aclrc (Memc[DC_CF(ap,i)], l)
		    sfd = stropen (Memc[DC_CF(ap,i)], l, NEW_FILE)
		    call fprintf (sfd, "%.8g %g")
			call pargd (wt1)
			call pargd (Memd[pshift1+k1])
		    dval = DC_DW(ap,i) * (DC_NW(ap,i) - 1) / 2.
		    call fprintf (sfd, " 1 2 1 %d %g %g")
			call pargi (DC_NW(ap,i))
			call pargd (DC_W1(ap,i) + dval)
			call pargd (dval)
		}
	    } else {
		ncoeffs = nint (Memd[coeffs])
		l = 20 * (ncoeffs + 2)
		if (wt2 > 0.)
		    l = 2 * l
		call realloc (DC_CF(ap,i), l, TY_CHAR)
		call aclrc (Memc[DC_CF(ap,i)], l)
		sfd = stropen (Memc[DC_CF(ap,i)], l, NEW_FILE)
		call fprintf (sfd, "%.8g %g %d %d")
		    call pargd (wt1)
		    call pargd (Memd[pshift1+k1])
		    call pargi (nint (Memd[coeffs+1]))
		    call pargi (nint (Memd[coeffs+2]))
		do k = 3, ncoeffs {
		    call fprintf (sfd, " %.15g")
			call pargd (Memd[coeffs+k])
		}
	    }

	    if (wt2 > 0.) {
		for (k2=0; k2<naps2 && Memi[paps2+k2]!=j; k2=k2+1)
		    ;
		if (k2 == naps2)
		    for (k2=0; k2<naps2 && !IS_INDEFI(Memi[paps2+k2]); k2=k2+1)
			;
		if (k2 == naps2) {
		    if (ignoreaps == YES)
			k2 = 0
		    else {
			call sprintf (Memc[str1], SZ_LINE,
			    "%s - Missing reference for aperture %d")
			    call pargstr (spec)
			    call pargi (j)
			if (sfd != NULL)
			    call strclose (sfd)
			call sfree (sp)
			call fatal (1, Memc[str1])
		    }
		}
		un2 = Memi[punits2+k2]
		if (!un_compare (un1, un2)) {
		    call sfree (sp)
		    call error (2,
			"Can't combine references with different units")
		}
		if (DC_DT(ap,i)==2 && !(coeffs==NULL&&Memi[pcoeff2+k2]==NULL)) {
		    call sfree (sp)
		    call error (2,
			"Can't combine references with non-linear dispersions")
		}
		coeffs = Memi[pcoeff2+k2]
		if (coeffs == NULL) {
		    if (DC_DT(ap,i) == 2) {
			dval = (wt1*Memd[pshift1+k1] + wt2*Memd[pshift2+k2]) /
			    (wt1 + wt2)
			call sshift1 (dval, DC_CF(ap,i))
		    } else {
			call fprintf (sfd, " %.8g %g")
			    call pargd (wt2)
			    call pargd (Memd[pshift2+k2])
			dval = DC_DW(ap,i) * (DC_NW(ap,i) - 1) / 2.
			call fprintf (sfd, " 1 2 1 %d %g %g")
			    call pargi (DC_NW(ap,i))
			    call pargd (DC_W1(ap,i) + dval)
			    call pargd (dval)
		    }
		} else {
		    call fprintf (sfd, " %.8g %g %d %d")
			call pargd (wt2)
			call pargd (Memd[pshift2+k2])
			call pargi (nint (Memd[coeffs+1]))
			call pargi (nint (Memd[coeffs+2]))
		    ncoeffs = nint (Memd[coeffs])
		    do k = 3, ncoeffs {
			call fprintf (sfd, " %.15g")
			    call pargd (Memd[coeffs+k])
		    }
		}
	    }

	    if (i == 1) {
		un = un1
		if (UN_LABEL(un) != EOS)
		    call mw_swattrs (SMW_MW(smw,0), 1, "label", UN_LABEL(un))
		if (UN_UNITS(un) != EOS)
		    call mw_swattrs (SMW_MW(smw,0), 1, "units", UN_UNITS(un))
		call un_close (DC_UN(ap,i))
		DC_UN(ap,i) = un
	    } else if (!un_compare (un, un1)) {
		call sfree (sp)
		call error (3, "Units must be the same for all apertures")
	    }
	    DC_DT(ap,i) = 2
	    call smw_swattrs (smw, DC_PL(ap,i), 1, DC_AP(ap,i), DC_BM(ap,i),
		DC_DT(ap,i), DC_W1(ap,i), DC_DW(ap,i), DC_NW(ap,i), DC_Z(ap,i),
		DC_LW(ap,i), DC_UP(ap,i), Memc[DC_CF(ap,i)])
	    if (sfd != NULL)
		call strclose (sfd)
	}

	# Update the linear part of WCS.
	ct1 = smw_sctran (smw, "logical", "physical", 2)
	ct2 = smw_sctran (smw, "physical", "world", 3)
	do i = 1, naps {
	    call smw_gwattrs (smw, DC_PL(ap,i), 1, DC_AP(ap,i), DC_BM(ap,i),
		DC_DT(ap,i), DC_W1(ap,i), DC_DW(ap,i), DC_NW(ap,i), DC_Z(ap,i),
		DC_LW(ap,i), DC_UP(ap,i), DC_CF(ap,i))
	    wt1 = nint (smw_c1trand (ct1, double(i)))
	    call smw_c2trand (ct2, double(DC_NW(ap,i)), wt1, DC_W2(ap,i), wt2)
	    DC_DW(ap,i) = (DC_W2(ap,i) - DC_W1(ap,i)) / (DC_NW(ap,i) - 1)
	    call smw_swattrs (smw, DC_PL(ap,i), 1, DC_AP(ap,i), DC_BM(ap,i),
		DC_DT(ap,i), DC_W1(ap,i), DC_DW(ap,i), DC_NW(ap,i), DC_Z(ap,i),
		DC_LW(ap,i), DC_UP(ap,i), Memc[DC_CF(ap,i)])
	}
	call smw_ctfree (ct1)
	call smw_ctfree (ct2)

done_	# Set aperture parameters in terms of logical image.
	ct1 = smw_sctran (smw, "physical", "logical", 1)
	j = nint (smw_c1trand (ct1, 1D0))
	do i = 1, naps {
	    k = nint (smw_c1trand (ct1, double(DC_NW(ap,i))))
	    DC_NW(ap,i) = min (IM_LEN(im,1), max (j, k))
	}
	call smw_ctfree (ct1)

	ct1 = smw_sctran (smw, "logical", "world", 3)
	do i = 1, naps {
	    wt1 = i
	    call smw_c2trand (ct1, 1D0, wt1, DC_W1(ap,i), wt2)
	    call smw_c2trand (ct1, double(DC_NW(ap,i)), wt1, DC_W2(ap,i), wt2)
	    DC_DW(ap,i) = (DC_W2(ap,i) - DC_W1(ap,i)) / (DC_NW(ap,i) - 1)
	}
	call smw_ctfree (ct1)

	do i = 1, naps
	    call mfree (DC_CF(ap,i), TY_CHAR)
	call sfree (sp)
end


# DC_GMSDB -- Get a dispersion database entry.
# The database entry is read only once from the database and stored in a
# symbol table keyed by the spectrum name.  Subsequent requests for the
# reference spectrum returns the data from the symbol table.

procedure dc_gmsdb (spec, stp, paps, papcen, punits, pshift, pcoeff, naps)

char	spec[ARB]	# Spectrum image name
pointer	stp		# Symbol table pointer
pointer	paps		# Pointer to aperture numbers
pointer	papcen		# Pointer to aperture centers
pointer	punits		# Pointer to units
pointer	pshift		# Pointer to shifts
pointer	pcoeff		# Pointer to coefficients
int	naps		# Number of apertures

double	dval
int	i, n, dtgeti(), getline(), ctod()
real	low, high, dtgetr()
pointer	sp, str, coeffs, sym, db, dt, dt1
pointer	stfind(), stenter(), strefsbuf(), dtmap1(), un_open()
errchk	dtmap1, dtgeti, dtgad, un_open

begin
	# Check if dispersion solution is in the symbol table from a previous
	# call.  If not in the symbol table get it from the database and
	# store it in the symbol table.

	sym = stfind (stp, spec)
	if (sym == NULL) {
	    call smark (sp)
	    call salloc (str, SZ_LINE, TY_CHAR)
	    call strcpy ("id", Memc[str], SZ_LINE)
	    call imgcluster (spec, Memc[str+2], SZ_LINE-2)
	    call xt_imroot (Memc[str+2], Memc[str+2], SZ_LINE-2)
	    db = strefsbuf (stp, Memi[stfind (stp, "database")])
	    dt = dtmap1 (Memc[db], Memc[str], READ_ONLY)
	    call strcpy ("ec", Memc[str], SZ_LINE)
	    call imgcluster (spec, Memc[str+2], SZ_LINE-2)
	    call xt_imroot (Memc[str+2], Memc[str+2], SZ_LINE-2)
	    ifnoerr (dt1 = dtmap1 (Memc[db], Memc[str], READ_ONLY)) {
		call sprintf (Memc[str], SZ_LINE,
		    "Ambiguous database files: %s/%s and %s/%s")
		    call pargstr (DT_DNAME(dt))
		    call pargstr (DT_FNAME(dt))
		    call pargstr (DT_DNAME(dt1))
		    call pargstr (DT_FNAME(dt1))
		call dtunmap (dt)
		call dtunmap (dt1)
		call fatal (3, Memc[str])
	    }

	    naps = max (1, DT_NRECS(dt))
	    call calloc (paps, naps, TY_INT)
	    call calloc (papcen, naps, TY_REAL)
	    call calloc (punits, naps, TY_POINTER)
	    call calloc (pshift, naps, TY_DOUBLE)
	    call calloc (pcoeff, naps, TY_POINTER)
	    if (DT_NRECS(dt) > 0) {
		for (i = 1; i <= naps; i = i + 1) {
		    iferr (Memi[paps+i-1] = dtgeti (dt, i, "aperture"))
			Memi[paps+i-1] = INDEFI
		    iferr (low = dtgetr (dt, i, "aplow"))
			low = INDEF
		    iferr (high = dtgetr (dt, i, "aphigh"))
			high = INDEF
		    if (IS_INDEF(low) || IS_INDEF(high))
			Memr[papcen+i-1] = 0.
		    else
			Memr[papcen+i-1] = (low + high) / 2.
		    iferr (call dtgstr (dt, i, "units", Memc[str], SZ_LINE))
			call strcpy ("Angstroms", Memc[str], SZ_LINE)
		    Memi[punits+i-1] = un_open (Memc[str])
		    iferr (Memd[pshift+i-1] = dtgetr (dt, i, "shift"))
			Memd[pshift+i-1] = 0.
		    iferr {
			n = dtgeti (dt, i, "coefficients")
			call malloc (coeffs, 1+n, TY_DOUBLE)
			Memd[coeffs] = n
			call dtgad (dt, i, "coefficients", Memd[coeffs+1], n, n)
			Memi[pcoeff+i-1] = coeffs
		    } then
			Memi[pcoeff+i-1] = NULL
		}
	    } else {
		Memi[paps] = INDEFI
		Memr[papcen] = INDEFR
		Memi[punits] = un_open ("")
		Memd[pshift] = 0.
		call malloc (coeffs, 100, TY_DOUBLE)
		n = 3
		call seek (Memi[dt], BOF)
		while (getline (Memi[dt], Memc[str]) != EOF) {
		    i = 1
		    if (ctod (Memc[str], i, dval) == 0)
			next
		    if (mod (n, 100) == 0)
			call realloc (coeffs, n+100, TY_DOUBLE)
		    Memd[coeffs+n] = dval
		    n = n + 1
		}
		Memd[coeffs] = n - 1
		Memd[coeffs+1] = 5
		Memd[coeffs+2] = n - 3
		Memi[pcoeff] = coeffs
	    }

	    call dtunmap (dt)
	    call sfree (sp)

	    sym = stenter (stp, spec, LEN_DC)
	    DC_FORMAT(sym) = 1
	    DC_PAPS(sym) = paps
	    DC_PAPCEN(sym) = papcen
	    DC_PUN(sym) = punits
	    DC_PSHIFT(sym) = pshift
	    DC_PCOEFF(sym) = pcoeff
	    DC_NAPS(sym) = naps
	} else {
	    if (DC_FORMAT(sym) != 1)
		call error (1, "Not a multispec dispersion function")
	    paps = DC_PAPS(sym)
	    papcen = DC_PAPCEN(sym)
	    punits = DC_PUN(sym)
	    pshift = DC_PSHIFT(sym)
	    pcoeff = DC_PCOEFF(sym)
	    naps = DC_NAPS(sym)
	}
end


# DC_REFSHFT -- Compute dispersion shift.

procedure dc_refshft (spec, stp, refspec, keywrd, im, aps, apcens, shifts,
	naps, fd1, fd2)

char	spec[ARB]		# Spectrum to be corrected
pointer	stp			# Symbol table pointer
char	refspec[ARB]		# Reference spectrum
char	keywrd[ARB]		# Header keyword (for log only)
pointer	im			# IMIO pointer to spectrum to be corrected
int	aps[naps]		# Reference apertures
real	apcens[naps]		# Reference aperture centers
double	shifts[naps]		# Reference aperture shifts (to be modified)
int	naps			# Number of refernce apertures
int	fd1			# Logfile descriptor
int	fd2			# Logfile descriptor

int	i, j, k, pnaps
double	apcen, shift, sumx, sumy, sumxx, sumyy, sumxy, a, b
pointer	sp, refshft, option, paps, papcen, punits, pshift, pcoeff
bool	streq()
errchk	imgstr, dc_gmsdb

begin
	call smark (sp)
	call salloc (refshft, SZ_FNAME, TY_CHAR)
	call salloc (option, SZ_FNAME, TY_CHAR)

	# Parse header parameter.
	call imgstr (im, keywrd, Memc[refshft], SZ_FNAME)
	call sscan (Memc[refshft])
	call gargwrd (Memc[refshft], SZ_FNAME)
	if (streq (Memc[refshft], refspec)) {
	    call sfree (sp)
	    return
	}
	call gargwrd (Memc[option], SZ_FNAME)

	# Get reference shift apertures.
	call dc_gmsdb (Memc[refshft], stp, paps, papcen, punits, pshift,
	    pcoeff, pnaps)
	if (pnaps == 0) {
	    call sfree (sp)
	    return
	}
	
	# Compute mean shift and RMS.
	sumy = 0.
	sumyy = 0.
	do i = 1, pnaps {
	    sumy = sumy + Memd[pshift+i-1]
	    sumyy = sumyy + Memd[pshift+i-1] ** 2
	}
	sumy = sumy / pnaps
	sumyy = sqrt (max (0.D0, sumyy / pnaps - sumy ** 2))

	# Print.
	if (fd1 != NULL) {
	    call fprintf (fd1, "%s: %s = '%s %s', shift = %.6g, rms = %.6g\n")
		call pargstr (spec)
		call pargstr (keywrd)
		call pargstr (Memc[refshft])
		call pargstr (Memc[option])
		call pargd (sumy)
		call pargd (sumyy)
	}
	if (fd2 != NULL) {
	    call fprintf (fd2, "%s: %s = '%s %s', shift = %.6g, rms = %.6g\n")
	        call pargstr (spec)
		call pargstr (keywrd)
	        call pargstr (Memc[refshft])
	        call pargstr (Memc[option])
	        call pargd (sumy)
	        call pargd (sumyy)
	}

	if (streq (Memc[option], "interp")) {
	    if (pnaps > 1) {
	        sumx = 0.
	        sumy = 0.
	        sumxx = 0.
	        sumyy = 0.
	        sumxy = 0.
	        do i = 0, pnaps-1 {
		    apcen = Memr[papcen+i]
		    shift = Memd[pshift+i]
		    sumx = sumx + apcen
		    sumy = sumy + shift
		    sumxx = sumxx + apcen * apcen
		    sumyy = sumyy + shift * shift
		    sumxy = sumxy + apcen * shift
	        }
	        b = pnaps * sumxx - sumx * sumx
	        a = (sumy * sumxx - sumx * sumxy) / b
	        b = (pnaps * sumxy - sumx * sumy) / b
	    } else {
	        a = sumy
	        b = 0.
	    }
	    do i = 1, naps
	        shifts[i] = shifts[i] + a + b * apcens[i]
	    if (fd1 != NULL) {
	        call fprintf (fd1, "\tintercept = %.6g, slope = %.6g\n")
		    call pargd (a)
		    call pargd (b)
	    }
	    if (fd2 != NULL) {
	        call fprintf (fd2, "\tintercept = %.6g, slope = %.6g\n")
		    call pargd (a)
		    call pargd (b)
	    }
	} else if (streq (Memc[option], "nearest")) {
	    do i = 1, naps {
	        k = 0
	        sumy = abs (apcens[i] - Memr[papcen])
	        for (j = 1; j < pnaps; j = j + 1)
		    if (abs (apcens[i] - Memr[papcen+j]) < sumy) {
		        k = j
		        sumy = abs (apcens[i] - Memr[papcen+k])
		    }
	        shifts[i] = shifts[i] + Memd[pshift+k]
	        if (fd1 != NULL) {
		    call fprintf (fd1, "\t%4d %7.2f %4d %7.2f %.6g\n")
		        call pargi (aps[i])
		        call pargr (apcens[i])
		        call pargi (Memi[paps+k])
		        call pargr (Memr[papcen+k])
		        call pargd (Memd[pshift+k])
	        }
	        if (fd2 != NULL) {
		    call fprintf (fd2, "\t%4d %7.2f %4d %7.2f %.6g\n")
		        call pargi (aps[i])
		        call pargr (apcens[i])
		        call pargi (Memi[paps+k])
		        call pargr (Memr[papcen+k])
		        call pargd (Memd[pshift+k])
	        }
	    }
	} else
	    call aaddkd (shifts, sumy, shifts, naps)

	call sfree (sp)
end


# DC_GEC -- Get an echelle spectrum.  This consists of mapping the image
# and setting a MWCS coordinate transformation.  If not dispersion corrected
# the dispersion function is found in the database for the reference
# spectra and set in the SMW.

procedure dc_gec (spec, im, smw, stp, ap, fd1, fd2)

char	spec[ARB]	#I Spectrum name
pointer	im		#I IMIO pointer
pointer	smw		#I SMW pointers
pointer	stp		#I Symbol table
pointer	ap		#O Aperture data structure
int	fd1		#I Logfile descriptor
int	fd2		#I Logfile descriptor

double	wt1, wt2, dval
int	i, j, k, l, dc, sfd, axis[2], naps, ncoeffs, offset, slope
pointer	sp, str1, str2, coeff, coeffs, ct1, ct2, un1, un2, un3
pointer	pshift1, pshift2, pshift3, pcoeff1, pcoeff2, pcoeff3

bool	un_compare()
double	smw_c1trand()
int	imaccf(), nscan(), stropen()
pointer	smw_sctran(), un_open()
errchk	dc_gecdb, imgstr, smw_sctran, un_open

data	axis/1,2/
define	done_	90

begin
	call smark (sp)
	call salloc (str1, SZ_LINE, TY_CHAR)
	call salloc (str2, SZ_LINE, TY_CHAR)
	coeff = NULL

	# Set WCS attributes
	naps = IM_LEN(im,2)
	call calloc (ap, LEN_AP(naps), TY_STRUCT)
	do i = 1, naps {
	    DC_PL(ap,i) = i
	    call smw_gwattrs (smw, DC_PL(ap,i), 1, DC_AP(ap,i), DC_BM(ap,i),
		DC_DT(ap,i), DC_W1(ap,i), DC_DW(ap,i), DC_NW(ap,i), DC_Z(ap,i),
		DC_LW(ap,i), DC_UP(ap,i), coeff)
	    if (i == 1) {
		iferr (call mw_gwattrs (SMW_MW(smw,0), 1, "units", Memc[str1],
		    SZ_LINE))
		    Memc[str1] = EOS
		DC_UN(ap,i) = un_open (Memc[str1])
	    }
	    dc = DC_DT(ap,i)
	}

	# Check if the spectra have been dispersion corrected
	# by an earlier version of DISPCOR.  If so then don't allow
	# another database dispersion correction.  This assumes all
	# spectra have the same dispersion type.  Check for a
	# reference spectrum.

	if ((imaccf (im, "REFSPEC1") == NO) ||
	    (dc > -1 && imaccf (im, "DCLOG1") == NO)) {
	    if (fd1 != NULL) {
		call fprintf (fd1,
		    "%s: Resampling using current coordinate system\n")
		    call pargstr (spec)
	    }
	    if (fd2 != NULL) {
		call fprintf (fd2,
		    "%s: Resampling using current coordinate system\n")
		    call pargstr (spec)
	    }
	    goto done_
	}

	# Get the reference spectra dispersion function from the database
	# and determine a reference shift.

	iferr {
	    call imgstr (im, "REFSPEC1", Memc[str1], SZ_LINE)
	    call sscan (Memc[str1])
	    call gargwrd (Memc[str1], SZ_LINE)
	    call gargd (wt1)
	    if (nscan() == 1)
		wt1 = 1.
	} then {
	    call strcpy (spec, Memc[str1], SZ_LINE)
	    wt1 = 1.
	}
	call salloc (pshift1, naps, TY_DOUBLE)
	call salloc (pcoeff1, naps, TY_POINTER)
	slope = 0
	iferr (call dc_gecdb (Memc[str1], stp, ap, un1, Memd[pshift1],
	    Memi[pcoeff1], naps, offset, slope)) {
	    call sfree (sp)
	    call erract (EA_ERROR)
	}
	if (fd1 != NULL) {
	    call fprintf (fd1, "%s: REFSPEC1 = '%s %.8g'\n")
		call pargstr (spec)
		call pargstr (Memc[str1])
		call pargd (wt1)
	}
	if (fd2 != NULL) {
	    call fprintf (fd2, "%s: REFSPEC1 = '%s %.8g'\n")
		call pargstr (spec)
		call pargstr (Memc[str1])
		call pargd (wt1)
	}

        iferr {
            call imgstr (im, "refshft1", Memc[str1], SZ_LINE)
	    call salloc (pshift3, naps, TY_DOUBLE)
	    call salloc (pcoeff3, naps, TY_POINTER)
	    call dc_gecdb (Memc[str1], stp, ap, un3, Memd[pshift3],
		Memi[pcoeff3], naps, offset, slope)
            if (fd1 != NULL) {
                call fprintf (fd1, "%s: REFSHFT1 = '%s', shift = %.6g\n")
                    call pargstr (spec)
                    call pargstr (Memc[str1])
                    call pargd (Memd[pshift3])
            }
            if (fd2 != NULL) {
                call fprintf (fd2, "%s: REFSHFT1 = '%s', shift = %.6g\n")
                    call pargstr (spec)
                    call pargstr (Memc[str1])
                    call pargd (Memd[pshift3])
            }
	    call aaddd (Memd[pshift1], Memd[pshift3], Memd[pshift1], naps)
        } then
            ;

	iferr {
            call imgstr (im, "REFSPEC2", Memc[str1], SZ_LINE)
            call sscan (Memc[str1])
            call gargwrd (Memc[str1], SZ_LINE)
            call gargd (wt2)
            if (nscan() == 1)
                wt2 = 1.
	    call salloc (pshift2, naps, TY_DOUBLE)
	    call salloc (pcoeff2, naps, TY_POINTER)
	    call dc_gecdb (Memc[str1], stp, ap, un2, Memd[pshift2],
		Memi[pcoeff2], naps, offset, slope)
	    if (fd1 != NULL) {
		call fprintf (fd1, "%s: REFSPEC2 = '%s %.8g'\n")
		    call pargstr (spec)
		    call pargstr (Memc[str1])
		    call pargd (wt2)
	    }
	    if (fd2 != NULL) {
		call fprintf (fd2, "%s: REFSPEC2 = '%s %.8g'\n")
		    call pargstr (spec)
		    call pargstr (Memc[str1])
		    call pargd (wt2)
	    }

	    iferr {
		call imgstr (im, "refshft2", Memc[str1], SZ_LINE)
		call salloc (pshift3, naps, TY_DOUBLE)
		call salloc (pcoeff3, naps, TY_POINTER)
		call dc_gecdb (Memc[str1], stp, ap, un3, Memd[pshift3],
		    Memi[pcoeff3], naps, offset, slope)
		if (fd1 != NULL) {
		    call fprintf (fd1, "%s: REFSHFT2 = '%s', shift = %.6g\n")
			call pargstr (spec)
                    call pargstr (Memc[str1])
                    call pargd (Memd[pshift3])
		}
		if (fd2 != NULL) {
		    call fprintf (fd2, "%s: REFSHFT2 = '%s', shift = %.6g\n")
			call pargstr (spec)
                    call pargstr (Memc[str1])
                    call pargd (Memd[pshift3])
		}
		call aaddd (Memd[pshift1], Memd[pshift3], Memd[pshift1], naps)
	    } then
		;
	} then
	    wt2 = 0.

	# Adjust weights to unit sum.
	dval = wt1 + wt2
	wt1 = wt1 / dval
	wt2 = wt2 / dval

	# Enter dispersion function in the MWCS.
	do i = 1, naps {
	    coeffs = Memi[pcoeff1+i-1]
	    ncoeffs = nint (Memd[coeffs])
	    l = 20 * (ncoeffs + 2)
	    if (wt2 > 0.)
		l = 2 * l
	    call realloc (coeff, l, TY_CHAR)
	    call aclrc (Memc[coeff], l)
	    sfd = stropen (Memc[coeff], l, NEW_FILE)
	    call fprintf (sfd, "%.8g %g")
		call pargd (wt1)
		call pargd (Memd[pshift1+i-1])

	    # The following assumes some knowledge of the data structure in
	    # order to shortten the the attribute string.

	    call fprintf (sfd, " %d %d %.8g %.8g")
		call pargi (nint (Memd[coeffs+1]))
		call pargi (nint (Memd[coeffs+2]))
		call pargd (Memd[coeffs+3])
		call pargd (Memd[coeffs+4])
	    do j = 5, ncoeffs {
		call fprintf (sfd, " %.15g")
		    call pargd (Memd[coeffs+j])
	    }

	    if (wt2 > 0.) {
		coeffs = Memi[pcoeff2+i-1]
		ncoeffs = nint (Memd[coeffs])
		call fprintf (sfd, "%.8g %g")
		    call pargd (wt2)
		    call pargd (Memd[pshift2+i-1])
		call fprintf (sfd, " %d %d %.8g %.8g")
		    call pargi (nint (Memd[coeffs+1]))
		    call pargi (nint (Memd[coeffs+2]))
		    call pargd (Memd[coeffs+3])
		    call pargd (Memd[coeffs+4])
		do j = 5, ncoeffs {
		    call fprintf (sfd, " %.15g")
			call pargd (Memd[coeffs+j])
		}
		if (!un_compare (un1, un2)) {
		    call sfree (sp)
		    call error (2,
			"Can't combine references with different units")
		}
	    }

	    if (i == 1) {
		if (UN_LABEL(un1) != EOS)
		    call mw_swattrs (SMW_MW(smw,0), 1, "label", UN_LABEL(un1))
		if (UN_UNITS(un1) != EOS)
		    call mw_swattrs (SMW_MW(smw,0), 1, "units", UN_UNITS(un1))
		call un_close (DC_UN(ap,i))
		DC_UN(ap,i) = un1
	    }
	    DC_DT(ap,i) = 2
	    call smw_swattrs (smw, DC_PL(ap,i), 1, DC_AP(ap,i), DC_BM(ap,i),
		DC_DT(ap,i), DC_W1(ap,i), DC_DW(ap,i), DC_NW(ap,i), DC_Z(ap,i),
		DC_LW(ap,i), DC_UP(ap,i), Memc[coeff])
	    call strclose (sfd)
	}

	# Update the linear part of WCS.
	ct1 = smw_sctran (smw, "logical", "physical", 2)
	ct2 = smw_sctran (smw, "physical", "world", 3)
	do i = 1, naps {
	    call smw_gwattrs (smw, DC_PL(ap,i), 1, DC_AP(ap,i), DC_BM(ap,i),
		DC_DT(ap,i), DC_W1(ap,i), DC_DW(ap,i), DC_NW(ap,i), DC_Z(ap,i),
		DC_LW(ap,i), DC_UP(ap,i), coeff)
	    wt1 = nint (smw_c1trand (ct1, double(i)))
	    call smw_c2trand (ct2, 1D0, wt1, DC_W1(ap,i), wt2)
	    call smw_c2trand (ct2, double(DC_NW(ap,i)), wt1, DC_W2(ap,i), wt2)
	    DC_DW(ap,i) = (DC_W2(ap,i) - DC_W1(ap,i)) / (DC_NW(ap,i) - 1)
	    call smw_swattrs (smw, DC_PL(ap,i), 1, DC_AP(ap,i), DC_BM(ap,i),
		DC_DT(ap,i), DC_W1(ap,i), DC_DW(ap,i), DC_NW(ap,i), DC_Z(ap,i),
		DC_LW(ap,i), DC_UP(ap,i), Memc[coeff])
	}
	call smw_ctfree (ct1)
	call smw_ctfree (ct2)

done_	# Set aperture parameters in terms of logical image.
	ct1 = smw_sctran (smw, "physical", "logical", 1)
	j = nint (smw_c1trand (ct1, 1D0))
	do i = 1, naps {
	    k = nint (smw_c1trand (ct1, double(DC_NW(ap,i))))
	    DC_NW(ap,i) = min (IM_LEN(im,1), max (j, k))
	}
	call smw_ctfree (ct1)

	ct1 = smw_sctran (smw, "logical", "world", 3)
	do i = 1, naps {
	    wt1 = i
	    call smw_c2trand (ct1, 1D0, wt1, DC_W1(ap,i), wt2)
	    call smw_c2trand (ct1, double(DC_NW(ap,i)), wt1, DC_W2(ap,i), wt2)
	    DC_DW(ap,i) = (DC_W2(ap,i) - DC_W1(ap,i)) / (DC_NW(ap,i) - 1)
	}
	call smw_ctfree (ct1)

	call mfree (coeff, TY_CHAR)
	call sfree (sp)
end


# DC_GECDB -- Get a dispersion database entry.
# The database entry is read only once from the database and stored in a
# symbol table keyed by the spectrum name.  Subsequent requests for the
# reference spectrum returns the data from the symbol table.

procedure dc_gecdb (spec, stp, ap, un, shifts, pcoeff, naps, offset, slope)

char	spec[ARB]	# Spectrum image name
pointer	stp		# Symbol table pointer
pointer	ap		# Aperture data structure
pointer	un		# Units
double	shifts[naps]	# Shifts
pointer	pcoeff[naps]	# Pointer to coefficients
int	naps		# Number of apertures
int     offset          # Aperture to order offset
int     slope           # Aperture to order slope

double  shift
real    dtgetr()
int     i, rec, offst, slpe, n, dtlocate(), dtgeti()
pointer sp, str, coeffs, sym, db, dt
pointer	stfind(), stenter(), strefsbuf(), dtmap1(), un_open()
errchk  dtmap1, dtlocate, dtgeti, dtgad, un_open

begin
	# Check if dispersion solution is in the symbol table from a previous
	# call.  If not in the symbol table get it from the database and
	# store it in the symbol table.

	sym = stfind (stp, spec)
	if (sym == NULL) {
	    call smark (sp)
	    call salloc (str, SZ_LINE, TY_CHAR)
	    call strcpy ("ec", Memc[str], SZ_LINE)
	    call imgcluster (spec, Memc[str+2], SZ_LINE-2)
	    call xt_imroot (Memc[str+2], Memc[str+2], SZ_LINE-2)
	    db = strefsbuf (stp, Memi[stfind (stp, "database")])
	    dt = dtmap1 (Memc[db], Memc[str], READ_ONLY)

            call sprintf (Memc[str], SZ_LINE, "ecidentify %s")
                call pargstr (spec)
            iferr (rec = dtlocate (dt, Memc[str])) {
		call sprintf (Memc[str], SZ_LINE,
		    "DISPCOR: Echelle dispersion function not found (%s/%s)")
		    call pargstr (DT_DNAME(dt))
		    call pargstr (DT_FNAME(dt))
                call fatal (0, Memc[str])
	    }

	    iferr (call dtgstr (dt, rec, "units", Memc[str], SZ_LINE))
		call strcpy ("Angstroms", Memc[str], SZ_LINE)
	    un = un_open (Memc[str])
            iferr (offst = dtgeti (dt, rec, "offset"))
                offst = 0
            iferr (slpe = dtgeti (dt, rec, "slope"))
                slpe = 1
            iferr (shift = dtgetr (dt, rec, "shift"))
                shift = 0.
            n = dtgeti (dt, rec, "coefficients")
            call malloc (coeffs, n, TY_DOUBLE)
            call dtgad (dt, rec, "coefficients", Memd[coeffs], n, n)

            sym = stenter (stp, spec, LEN_DC)
	    DC_FORMAT(sym) = 2
	    DC_PUN(sym) = un
            DC_OFFSET(sym) = offst
            DC_SLOPE(sym) = slpe
            DC_SHIFT(sym) = shift
            DC_COEFFS(sym) = coeffs

	    call dtunmap (dt)
	    call sfree (sp)
	} else {
	    if (DC_FORMAT(sym) != 2)
		call error (1, "Not an echelle dispersion function")
	    un = DC_PUN(sym)
            offst = DC_OFFSET(sym)
            slpe = DC_SLOPE(sym)
            coeffs = DC_COEFFS(sym)
            shift = DC_SHIFT(sym)
	}

	# Check aperture to order parameters.
	if (slope == 0) {
	    offset = offst
	    slope = slpe
	} else if (offset != offst || slope != slpe) {
	    call eprintf (
		"WARNING: Echelle order offsets/slopes are not the same.\n")
	}

	# Convert to multispec coefficients
	do i = 1, naps {
	    DC_BM(ap,i) = offset + slope * DC_AP(ap,i)
	    call dc_ecms (DC_BM(ap,i), Memd[coeffs], pcoeff[i])
	    shifts[i] = shift / DC_BM(ap,i)
	}
end


# DC_ECMS -- Convert echelle dispersion coefficients to multispec coefficients

procedure dc_ecms (order, eccoeff, mscoeff)

int	order			# Echelle order
double	eccoeff[ARB]		# Echelle dispersion coefficients
pointer	mscoeff			# Pointer to multispec coefficients

int	i, j, k, type, xorder, yorder
double	xmin, xmax, ymin, ymax, ymaxmin, yrange, y, coeff, a, b, c

begin
	type = nint (eccoeff[1])
	xorder = nint (eccoeff[2])
	yorder = nint (eccoeff[3])
	xmin = eccoeff[5]
	xmax = eccoeff[6]
	ymin = eccoeff[7]
	ymax = eccoeff[8]

	yrange = 2. / (ymax - ymin)
	ymaxmin = (ymax + ymin) / 2
	y = (order - ymaxmin) * yrange

	call malloc (mscoeff, 5+xorder, TY_DOUBLE)
	Memd[mscoeff] = 4+xorder
	Memd[mscoeff+1] = type
	Memd[mscoeff+2] = xorder
	Memd[mscoeff+3] = xmin
	Memd[mscoeff+4] = xmax

	switch (type) {
	case 1:
	    do k = 1, xorder {
		j = 9 + k - 1
		coeff = eccoeff[j]
		if (yorder > 1) {
		    j = j + xorder
		    coeff = coeff + eccoeff[j] * y
		}
		if (yorder > 2) {
		    a = 1
		    b = y
		    do i = 3, yorder {
			c = 2 * y * b - a
			j = j + xorder
			coeff = coeff + eccoeff[j] * c
			a = b
			b = c
		    }
		}
		Memd[mscoeff+4+k] = coeff / order
	    }
	case 2:
	    do k = 1, xorder {
		j = 9 + k - 1
		coeff = eccoeff[j]
		if (yorder > 1) {
		    j = j + xorder
		    coeff = coeff + eccoeff[j] * y
		}
		if (yorder > 2) {
		    a = 1
		    b = y
		    do i = 3, yorder {
			c = ((2 * i - 3) * y * b - (i - 2) * a) / (i - 1)
			j = j + xorder
			coeff = coeff + eccoeff[j] * c
			a = b
			b = c
		    }
		}
		Memd[mscoeff+4+k] = coeff / order
	    }
	}
end