include <error.h>
include	<foc.h>
include "../dbfit.h"

# T_RFIXX -- fixes holes in reseau entry by evaluating distortion fit
# at reference reseau on the holes (INDEF positions), and replaces holes
# by the results - this leaves the non-hole (i.e. good) reseau marks untouched.
# An entry name in the database must match an entry name in the input reseau
# table.  If no match is found, the output reseau table will be deleted.
#
# Dave Giaretta			subroutine created
# Phil Hodge,  8-Dec-89		call to db_readr corrected:  add type & axis
# Phil Hodge,  1-Mar-90		fix several bugs
# Phil Hodge, 18-Mar-91		move defines to ../dbfit.h
# Phil Hodge, 20-Mar-92		delete output table if no entry found in db

procedure t_rfixx ()

char	input1[SZ_FNAME]		# input reseau with holes
char	entry[SZ_R_PATTERN]		# transform entries to use
char	input2[SZ_FNAME]		# fit file 
char	input3[SZ_FNAME]		# reference reseau file for fit
char	refentry[SZ_R_PATTERN]		# reference reseau entry for fit
char	output[SZ_FNAME]		# output reseau file
#--
pointer	rs_open(), dtmap()
pointer	respt			# ptr for output positions (copied or evaluated)
pointer rentp			# ptr for input positions (to be fixed)
pointer refpt			# ptr for reference positions
pointer	rp, rpref, rpout
pointer rs_template(), rlist, db_template(), dlist
pointer fit, dt
int	direction		# direction of fit:  backwards
int	axis			# has to do with long-slit mode; ignored
bool	not_found		# true if no reseau entry found in database
int	nowhite(), i
int	rs_readr(), rs_list(), db_list(), db_readr(), npts
bool	streq()

begin
	# This is for the call to db_readr.  It isn't really important
	# whether forward or backward is specified.
	direction = BACKWARD

	# get infile1 and check for blanks
	call clgstr( "input1", input1, SZ_FNAME)
	if (nowhite (input1, input1, SZ_FNAME) == 0) 
	    call error( 0, "input1 reseau file must be given")

	call clgstr ("entry", entry, SZ_R_PATTERN)

	call clgstr ("input2", input2, SZ_FNAME)
	if (nowhite (input2,  input2, SZ_FNAME) == 0)
	    call error (0, "input2 fit file must be given")
	else
	    # open database for fit coefficients
	    dt = dtmap (input2, READ_ONLY)
	
	# if input3 is blank then use input1
	call clgstr( "input3", input3, SZ_FNAME)
	if ( nowhite( input3,  input3, SZ_FNAME) == 0) 
	    call strcpy( input1, input3, SZ_FNAME)

	call clgstr( "refentry", refentry, SZ_R_PATTERN)

	# open the reseau table(s)
	rp = rs_open (input1, READ_ONLY, NULL) 

	# input2 replaced by input3, PEH 1-Mar-90
	if (streq (input1, input3))
	    rpref = rp
	else
	    rpref = rs_open (input3, READ_ONLY, NULL)

	# open output reseau file
	call clgstr ("output", output, SZ_FNAME)
	rpout = rs_open( output, NEW_COPY, rp)

	# allocate space for reseau data 
	call rs_alloc( rp   , TY_REAL, rentp)
	call rs_alloc( rpref, TY_REAL, refpt)
	call rs_alloc( rpout, TY_REAL, respt)
	# space for fit pointers
	call salloc( fit, LEN_FIT, TY_LONG)

	# read reference
	rlist = rs_template( rpref, refentry)
	if ( rs_readr ( rpref, refpt, EOS, rs_list( rlist ) ) == 
							RES_F_NONEXTENTRY )
	    call error ( 0, " cannot read reference entry ")

	npts = RES_NROWS(rpref)*RES_NCOLS(rpref)

	not_found = true			# initial value

	# expand entry template
	rlist = rs_template( rp, entry)
	# process until we reach the end of list
	while ( rs_readr( rp, rentp, EOS, rs_list( rlist) ) != 
							RES_F_NONEXTENTRY ) {

	    # expand db template with this entry name -
	    # process data if we find a fit
	    dlist = db_template( dt, RES_ENT_ENTRY(rentp) )
	    if (db_readr (dt, direction, EOS, db_list( dlist), axis, fit) != 
				RES_F_NONEXTENTRY) {

		# We found at least one entry in the database.
		not_found = false
	        iferr { 
		    # Memb[RES_FPT(respt)] was replaced by Memb[RES_FPT(rentp)]
		    # by PEH 1-Mar-90.
	            call rfixx (fit, npts,
				Memr[RES_XPT(refpt)], Memr[RES_YPT(refpt)],
				Memb[RES_FPT(rentp)],
				Memr[RES_XPT(respt)], Memr[RES_YPT(respt)])
	    	} then {
		    call eprintf ("Cannot process coordinate list: %s\n")
		    	call pargstr (RES_ENT_ENTRY(refpt))
		    call erract (EA_WARN)
		    next
	    	}

		# fix up the holes in the input1 entry
		do i = 1, npts {
		    if (Memb[RES_FPT(rentp)+i-1] ) {
			Memr[RES_XPT(rentp)+i-1] = Memr[RES_XPT(respt)+i-1]
			Memr[RES_YPT(rentp)+i-1] = Memr[RES_YPT(respt)+i-1]
		    }
		}

		# Now we write out the fixed up entry
	    	call rs_write( rpout, rentp, 0)
	    }
	}

	# close up
	call rs_close (rp)
	call dtunmap  (dt)
	call rs_close (rpref)
	call rs_close (rpout)
	if (not_found) {
	    call eprintf (
"Warning:  entry not found in input reseau table; output deleted.\n")
	    call tbtdel (output)
	}
end

# rfixx -- evaluate the coordinate transformations

procedure rfixx (fit, npts, x, y, flag, xout, yout)

pointer	fit			# i: pointer to fit parameters
int	npts			# i: number of points
real	x[ARB]			# i: x coords to eval fit
real	y[ARB]			# i: y coords to eval fit
bool	flag[ARB]		# i: flag values
real	xout[ARB]		# o: evaluated x positions
real	yout[ARB]		# o: evaluated y positions
#--
pointer	sx1, sy1, sx2, sy2
int	i
real	gseval()
errchk	gseval

begin
	sx1 = FIT_SX1(fit)
	sy1 = FIT_SY1(fit)
	sx2 = FIT_SX2(fit)
	sy2 = FIT_SY2(fit)

	# PEH 1-Mar-90:  reverse the logic here; it used to be that if flag[i]
	# was true then xout[i] & yout[i] were set to INDEFR.
	do i = 1, npts {
	    if ( flag[i] ) {
		xout[i] = gseval (sx1, x[i], y[i])
		if ( sx2 != NULL)
		    xout[i] = xout[i] + gseval (sx2, x[i], y[i])
		yout[i] = gseval (sy1, x[i], y[i])
		if ( sy2 != NULL)
		    yout[i] = yout[i] + gseval (sy2, x[i], y[i])
	    } else {
		xout[i] = x[i]
		yout[i] = y[i]
	    }
	}
end