include	<pkg/dttext.h>
include	<foc.h>
include "../dbfit.h"

define MAXITER 40       # Maximum iterations
define OKERROR 0.005    # Acceptable error in Y position

# This file contains db_template, db_list, db_readr,
# which are used for obtaining polynomial fit information from
# text database files.
#
# Dave Giaretta wrote the original.
# Richard Hook, ST-ECF, Jan 88  Add 1-D and backwards options.
# Phil Hodge,  3-Mar-1991  Remove print statements; extract from revalfitx.
# Phil Hodge, 20-Mar-1992  db_template:  allocate one more element, and
#		store numrows in correct location if nothing matches template.

# DB_TEMPLATE -- expands template for database entries

pointer procedure db_template (dt, entry) 

pointer	dt		# i: database pointer
char	entry[ARB]	# i: template for transform
#--
pointer	sel
int	icount, nsel, patsize, patmake(), patmatch()
char	template[SZ_R_PATTERN], name[SZ_R_PATTERN]
int	numrows

begin
	numrows = DT_NRECS(dt)
	# Allocate one for each entry in the database, plus LEN_RES_TINDEX
	# (=1) for the count of entries which have been returned by db_list
	# and one for the RES_F_NONEXTENTRY list terminator.
	call salloc( sel, numrows + LEN_RES_TINDEX + 1, TY_INT)

	# This is the index to the last item in the list that was read.
	nsel = 1
	RES_TEMPL_COUNT(sel) = nsel

	# make pattern stuff
	patsize = patmake( entry, template, SZ_R_PATTERN)

	# Include the entry number for every entry whose name matches
	# the template.
	for ( icount = 1; icount <= numrows; icount=icount+1) {
	    call strcpy( DT_NAME(dt, icount), name, SZ_R_ENTRY )
	    if ( patmatch( name, template ) != 0  ) {
		Memi[ sel + nsel] = icount 
	        nsel = nsel + 1
	    }
	}

	# If there was no match between the input entry (template) and
	# any entry in the database, include the last database entry.
	# It's not obvious to me that it should work this way, but it is
	# useful to be able to fix or evaluate all reseau table entries
	# with one database entry.
	# This used to be Memi[sel] = numrows, which caused db_list to
	# return garbage; fixed 3/20/92 by PEH.
	if (nsel == 1 && numrows > 0) {
	    Memi[sel+nsel] = numrows
	    nsel = nsel + 1
	}

	# This terminates the list of entry numbers that match the template.
	Memi[sel + nsel] = RES_F_NONEXTENTRY
	
	return (sel)
end

# DB_LIST -- give next entry number to read

int procedure db_list (dlist)

pointer	dlist		# i:   list pointer
#--
int	ip		# use count

begin
	ip = RES_TEMPL_COUNT(dlist)
 	ip = ip + 1
	if ( ip <= 0 )
	    call error(0, "illegal index in DB_LIST")

	RES_TEMPL_COUNT(dlist) = ip
	return ( Memi[dlist + ip - 1] )
end

# DB_READR -- reads specified database entry
#
# This routine modified for the longslit case, Richard Hook, ST-ECF
# August 1988. The transform information is headed 'surface' and
# represents lambda as a function of x,y.
# If a value is given for TRANSFORM it will be used; otherwise,
# the entry is gotten from REC.
#
# Further modifications to handle 1d shifts either in X or Y added - 
# Richard Hook, ST-ECF, Jan 89

int procedure db_readr (dt, direction, transform, rec, axis, fit)

pointer	dt		# i: database pointer
int     direction	# i: The direction of the transform
char	transform[ARB]	# i: name of transform wanted
int	rec		# i: record number of transform
int 	axis		# o: axis number (1=X,2=Y, 1d case only)
pointer	fit		# o: pointer to trans. coefficients
#--
pointer	sx1, sy1	# pointer to linear part of surface fit
pointer	sx2, sy2	# pointer to higher order surface

int	i, ncoeff, junk
pointer	xcoeff, ycoeff, coeffs, newsx1, newsy1  # coeffs added for longslit

int	dtlocate(), dtgeti(), dtscan()
errchk  dtgeti

begin
	# if rec is NONEXTENTRY then just return
	if ( rec == RES_F_NONEXTENTRY)
	    return( RES_F_NONEXTENTRY)

	# if name given then we must use it
	if ( transform[1] != EOS ) 
	    rec = dtlocate (dt, transform)

	call strcpy( DT_NAME(dt, rec), FIT_NAME(fit), SZ_R_ENTRY)

	# Check dimensionality; "axis" keyword will only be found
	# (and set to 1 or 2) in the 1d case.
	# Create the full 'type' code from the axis number and dimensionality.
	iferr (axis = dtgeti (dt, rec, "axis"))
	    FIT_TYPE(fit) = direction * 2
	else
	    FIT_TYPE(fit) = direction * 2 - 1

	if (FIT_TYPE(fit) == FOR2D || FIT_TYPE(fit) == BACK2D) {
	    # Simple case as in original program

	    # get linear part of fit
	    ncoeff = dtgeti (dt, rec, "surface1")
	    call malloc (xcoeff, ncoeff, TY_REAL)
	    call malloc (ycoeff, ncoeff, TY_REAL)
	    do i = 1, ncoeff {
		junk = dtscan (dt)
		call gargr (Memr[xcoeff+i-1])
		call gargr (Memr[ycoeff+i-1])
	    }
	    call gsrestore (sx1, Memr[xcoeff])
	    call gsrestore (sy1, Memr[ycoeff])

	    # adjust linear part of fit
	    call gscopy (sx1, newsx1)
	    call gscopy (sy1, newsy1)

	    # get higher order part of fit
	    ncoeff = dtgeti (dt, rec, "surface2")
	    if (ncoeff > 0 ) {

		# get distortion coefficients
		call realloc (xcoeff, ncoeff, TY_REAL)
		call realloc (ycoeff, ncoeff, TY_REAL)
		do i = 1, ncoeff {
		    junk = dtscan(dt)
		    call gargr (Memr[xcoeff+i-1])
		    call gargr (Memr[ycoeff+i-1])
		}
		call gsrestore (sx2, Memr[xcoeff])
		call gsrestore (sy2, Memr[ycoeff])

	    } else {

		# set higher order surface pointers to null
		sx2 = NULL
		sy2 = NULL
	    }

	    # redefine the surfaces
	    call gsfree (sx1)
	    call gscopy (newsx1, sx1)
	    call gsfree (newsx1)
	    call gsfree (sy1)
	    call gscopy (newsy1, sy1)
	    call gsfree (newsy1)

	    call mfree (xcoeff, TY_REAL)
	    call mfree (ycoeff, TY_REAL)
	    
	    # insert info into structure
	    FIT_SX1(fit) = sx1
	    FIT_SY1(fit) = sy1
	    FIT_SX2(fit) = sx2
	    FIT_SY2(fit) = sy2
	} else {
	    # 1d case - just one set of coeffs
	    # first find out which direction the shift is in

	    axis = dtgeti (dt, rec, "axis")

	    ncoeff = dtgeti (dt, rec, "surface")
	    call malloc (coeffs, ncoeff, TY_REAL)
	    call dtgar (dt, rec, "surface", Memr[coeffs], ncoeff, ncoeff)
	    call gsrestore (sx1, Memr[coeffs])
	    FIT_SX1(fit) = sx1
	    call mfree (coeffs, TY_REAL)
        }

	return ( rec )
end