include "ffit.h"

define	FROM_PAR_FILE	1	# get initial estimates from par file?
define	FROM_DATA	2	# get initial estimates from data?

# fcsinit -- initial values
# This file contains routines for getting the initial parameter values,
# either from the par file or from the data, and routines for getting
# an estimate of the range of each parameter value.
# Routines in this file:
# f_4g_init	initial values for Gaussian
# f_4g_range	range of variable parameters for Gaussian
# f_3h_init	initial values for hyperbolic function
# f_3h_range	range of variable parameters for hyperbolic function

# f_4g_init -- get initial parameter values for Gaussian

procedure f_4g_init (fnl, iinit)

pointer fnl		# i: pointer to fit struct
int	iinit		# i: get init est of param values from par file or data
#--
real	locn		# X location of maximum Y value
int	npts		# number of data values
int	k		# loop index
real	clgetr()

begin
	if (iinit == FROM_PAR_FILE) {
	    F_PAR(fnl,1) = clgetr ("locn")
	    F_PAR(fnl,2) = clgetr ("ampl")
	    F_PAR(fnl,3) = clgetr ("sigma")
	    F_PAR(fnl,4) = clgetr ("baseline")

	} else {			# estimate from data
	    # Look for the location of the maximum value in such a way that
	    # we check the middle point last, in case of duplicate values.
	    npts = F_NPTS(fnl)
	    locn = F_X(fnl,1)
	    do k = 1, npts/2
		if (F_Y(fnl,k) == F_YMAX(fnl))
		    locn = F_X(fnl,k)
	    do k = npts, npts/2 + 1, -1
		if (F_Y(fnl,k) == F_YMAX(fnl))
		    locn = F_X(fnl,k)

	    if (F_VAR(fnl,1))
		F_PAR(fnl,1) = locn		# location of center
	    else
		F_PAR(fnl,1) = clgetr ("locn")

	    if (F_VAR(fnl,2))
		F_PAR(fnl,2) = F_YMAX(fnl) - F_YMIN(fnl)
	    else
		F_PAR(fnl,2) = clgetr ("ampl")

	    if (F_VAR(fnl,4))
		F_PAR(fnl,4) = F_YMIN(fnl)
	    else
		F_PAR(fnl,4) = clgetr ("baseline")

	    # This is a pretty crude estimate of sigma!
	    if (F_VAR(fnl,3))
		F_PAR(fnl,3) = (F_XMAX(fnl) - F_XMIN(fnl)) / 4.
	    else
		F_PAR(fnl,3) = clgetr ("sigma")
	}
end

# f_4g_range -- get an estimate of the range of each parameter
# This is used to get initial positions for the simplex vertices.

procedure f_4g_range (fnl, range)

pointer fnl		# i: pointer to fit struct
real	range[ARB]	# o: the range (roughly) of each variable parameter
#--
int	k

begin
	k = 0
	if (F_VAR(fnl,1)) {
	    k = k + 1
	    range[k] = F_XMAX(fnl) - F_XMIN(fnl)	# location
	}
	if (F_VAR(fnl,2)) {
	    k = k + 1
	    range[k] = F_YMAX(fnl) - F_YMIN(fnl)	# amplitude
	}
	if (F_VAR(fnl,3)) {
	    k = k + 1
	    range[k] = F_XMAX(fnl) - F_XMIN(fnl)	# sigma
	}
	if (F_VAR(fnl,4)) {
	    k = k + 1
	    range[k] = F_YMAX(fnl) - F_YMIN(fnl)	# baseline
	}
end


# f_3h_init -- get initial parameter values
# This routine is for a hyperbola.  The four parameters are p[1] = location
# of center, p[2] = a, p[3] = b as in (y/a)**2 - (x/b)**2 = 1, and p[4]
# is the baseline.

procedure f_3h_init (fnl, iinit)

pointer fnl		# i: pointer to fit struct
int	iinit		# i: get init est of param values from par file or data
#--
real	slope		# a / b
real	locmin		# X location of minimum Y value
real	locmax		# X location of maximum Y value
int	npts		# number of data values
int	k		# loop index
real	clgetr()

begin
	if (iinit == FROM_PAR_FILE) {
	    F_PAR(fnl,1) = clgetr ("locn")
	    F_PAR(fnl,2) = clgetr ("yscale")
	    F_PAR(fnl,3) = clgetr ("xscale")
	    F_PAR(fnl,4) = clgetr ("baseline")

	} else {			# estimate from data
	    # Look for the location of the minimum value in such a way that
	    # we check the middle point last, in case of duplicate values.
	    npts = F_NPTS(fnl)
	    locmin = F_X(fnl,1)			# initial value
	    do k = 1, npts/2
		if (F_Y(fnl,k) == F_YMIN(fnl))
		    locmin = F_X(fnl,k)
	    do k = npts, npts/2 + 1, -1
		if (F_Y(fnl,k) == F_YMIN(fnl))
		    locmin = F_X(fnl,k)

	    locmax = F_X(fnl,1)			# initial value
	    do k = 1, npts {
		if (F_Y(fnl,k) == F_YMAX(fnl))
		    locmax = F_X(fnl,k)
	    }

	    if (F_VAR(fnl,1))
		F_PAR(fnl,1) = locmin		# location of center
	    else
		F_PAR(fnl,1) = clgetr ("locn")

	    if (F_VAR(fnl,2))
		F_PAR(fnl,2) = abs (F_YMIN(fnl))
	    else
		F_PAR(fnl,2) = clgetr ("yscale")

	    if (F_VAR(fnl,3)) {
		slope = abs ((F_YMAX(fnl) - F_YMIN(fnl)) / (locmax - locmin))
		F_PAR(fnl,3) = F_PAR(fnl,2) / slope
	    } else {
		F_PAR(fnl,3) = clgetr ("xscale")
	    }

	    if (F_VAR(fnl,4))
		F_PAR(fnl,4) = 0.
	    else
		F_PAR(fnl,4) = clgetr ("baseline")
	}
end

# f_3h_range -- get an estimate of the range of each parameter
# This is used to get initial positions for the simplex vertices.

procedure f_3h_range (fnl, range)

pointer fnl		# i: pointer to fit struct
real	range[ARB]	# o: the range (roughly) of each variable parameter
#--
int	k

begin
	k = 0
	if (F_VAR(fnl,1)) {
	    k = k + 1
	    range[k] = F_XMAX(fnl) - F_XMIN(fnl)	# location
	}
	if (F_VAR(fnl,2)) {
	    k = k + 1
	    range[k] = F_YMAX(fnl) - F_YMIN(fnl)	# y-scale (= a)
	}
	if (F_VAR(fnl,3)) {
	    k = k + 1
	    range[k] = F_XMAX(fnl) - F_XMIN(fnl)	# x-scale (= b)
	}
	if (F_VAR(fnl,4)) {
	    k = k + 1
	    range[k] = (F_YMAX(fnl) - F_YMIN(fnl)) / 3.	# baseline
	}
end