include "ffit.h"

define	LEN_XY_INCR	20	# size of increment in X, Y arrays

# f_rd_sharp -- read file of "sharpness" data
# This routine reads the output from either fradius or fsquares and saves
# the data (x = refoc position, y = blurriness or sharpness) in the fnl
# struct.  The F_XPT and F_YPT pointers are allocated by this routine.

procedure f_rd_sharp (fnl, infile)

pointer fnl		# i: pointer to fitting struct
char	infile[ARB]	# i: name of ascii file
#--
pointer sp
pointer buf		# scratch for input buffer
real	x, y		# values read from file
real	xmin, xmax	# min & max values of refoc position
real	ymin, ymax	# min & max values of dependent variable
int	fd		# fd for input file
int	npts		# local counter for number of values read
int	maxpts		# current size of X, Y arrays
int	nvpar		# number of variable parameters
int	ip		# for ctor
int	k
int	open(), getline(), ctor()

begin
	call smark (sp)
	call salloc (buf, SZ_LINE, TY_CHAR)

	fd = open (infile, READ_ONLY, TEXT_FILE)

	# Set initial size and allocate arrays for X & Y values.
	maxpts = LEN_XY_INCR
	call malloc (F_XPT(fnl), maxpts, TY_REAL)
	call malloc (F_YPT(fnl), maxpts, TY_REAL)

	npts = 0

	while (getline (fd, Memc[buf]) != EOF) {

	    if (Memc[buf] == '#')
		next

	    ip = 1
	    if (ctor (Memc[buf], ip, x) < 1)
		next

	    if (ctor (Memc[buf], ip, y) < 1) {
		call close (fd)
		call error (1, "fewer than 2 values on line in input file")
	    }
	    npts = npts + 1
	    if (npts > maxpts) {
		maxpts = maxpts + LEN_XY_INCR
		call realloc (F_XPT(fnl), maxpts, TY_REAL)
		call realloc (F_YPT(fnl), maxpts, TY_REAL)
	    }
	    F_X(fnl,npts) = x
	    F_Y(fnl,npts) = y
	}
	call close (fd)

	if (npts == 0)
	    call error (1, "no data in input file")

	F_NPTS(fnl) = npts
	call realloc (F_XPT(fnl), npts, TY_REAL)
	call realloc (F_YPT(fnl), npts, TY_REAL)

	# Count the number of variable parameters.
	nvpar = 0
	do k = 1, F_NPAR(fnl) {
	    if (F_VAR(fnl,k))
		nvpar = nvpar + 1
	}
	if (F_NPTS(fnl) < nvpar+1)
	    call error (1, "not enough data to fit that many parameters")

	# Find min & max of x & y.
	xmin = F_X(fnl,1)		# initial values
	xmax = xmin
	ymin = F_Y(fnl,1)
	ymax = ymin
	do k = 1, F_NPTS(fnl) {
	    if (F_X(fnl,k) < xmin)
		xmin = F_X(fnl,k)
	    if (F_X(fnl,k) > xmax)
		xmax = F_X(fnl,k)
	    if (F_Y(fnl,k) < ymin)
		ymin = F_Y(fnl,k)
	    if (F_Y(fnl,k) > ymax)
		ymax = F_Y(fnl,k)
	}

	F_XMIN(fnl) = xmin
	F_XMAX(fnl) = xmax
	F_YMIN(fnl) = ymin
	F_YMAX(fnl) = ymax
	if (xmax == xmin)
	    call error (1, "range of X values is zero")
	if (ymax == ymin)
	    call error (1, "range of Y values is zero")

	call sfree (sp)
end