include "ffit.h"

define	BIGVALUE	1.e10	# return this value if parameter out of range

# fcsfunk -- functions
# This file contains the following functions:
#
# f_gaus_funk	function to be minimized for a Gaussian
# fcs_gauss	Gaussian function
# f_hyp_funk	function to be minimized for a hyperbola
# fcs_hyper	hyperbolic function

# f_gaus_funk -- sum of squares for Gaussian
# This function evaluates the Gaussian function defined by the parameters SP
# and returns either the sum of the squared deviations or the sum of the
# absolute values of the deviations from the data.

real procedure f_gaus_funk (fnl, vp)

pointer fnl		# pointer to struct for fitting parameters
real	vp[MPAR]	# array of parameter values (only the variable ones)
#--
real	p[MPAR]		# array of parameter values (F_NPAR of them)
real	sum		# for accumulating sum of squares
real	g		# value of Gaussian function
int	j, k
real	fcs_gauss()

begin
	# Copy parameters to local array.  We do this so we can have some
	# variable parameters and some fixed ones.  The SP array contains
	# only the variable parameters, and we need them all.
	j = 0
	do k = 1, F_NPAR(fnl) {
	    if (F_VAR(fnl,k)) {
		j = j + 1
		p[k] = vp[j]
	    } else {
		p[k] = F_PAR(fnl,k)
	    }
	}

	# Check for unreasonable parameters.
	if (p[1] < 2 * F_XMIN(fnl) - F_XMAX(fnl))	# x too small?
	    return (BIGVALUE)
	if (p[1] > 2 * F_XMAX(fnl) - F_XMIN(fnl))	# x too big?
	    return (BIGVALUE)
	if (p[2] > 10. * (F_YMAX(fnl) - F_YMIN(fnl)))	# ampl too big?
	    return (BIGVALUE)
	if (p[3] < (F_XMAX(fnl) - F_XMIN(fnl)) / 100.)	# sigma too small?
	    return (BIGVALUE)

	# Accumulate the sum of squared differences (or the sum of the
	# absolute values) from the Gaussian.
	sum = 0.
	do k = 1, F_NPTS(fnl) {

	    g = fcs_gauss (F_X(fnl,k), p)	# Gaussian function

	    if (F_OPTION(fnl) == F_LEAST_SQUARES)
		sum = sum + (g - F_Y(fnl,k)) ** 2
	    else if (F_OPTION(fnl) == F_MINSUM)
		sum = sum + abs (g - F_Y(fnl,k))
	}

	return (sum)
end

# fcs_gauss -- evaluate a Gaussian function
# This routine evaluates a Gaussian using parameters passed in the
# calling sequence.

real procedure fcs_gauss (x, p)

real	x		# i: value at which function is to be evaluated
real	p[4]		# i: parameters for Gaussian
#--
real	arg		# argument for Gaussian function

begin
	arg = (x - p[1]) / p[3]
	arg = arg * arg / 2.
	return (p[2] * exp (-arg) + p[4])
end



# f_hyp_funk -- sum of squares for hyperbolic function
# This function evaluates the hyperbolic function defined by the parameters
# SP and returns either the sum of the squared deviations or the sum of the
# absolute values of the deviations from the data.
# The function is (y/b)**2 - (x/a)**2 = 1, where in this equation x is
# referred to an origin at p[1], and p[2] = a, p[3] = b.  There is also
# an offset (baseline) for y.

real procedure f_hyp_funk (fnl, vp)

pointer fnl		# pointer to struct for fitting parameters
real	vp[MPAR]	# array of parameter values (only the variable ones)
#--
real	p[MPAR]		# array of parameter values (F_NPAR of them)
real	sum		# for accumulating sum of squares
real	h		# value of hyperbolic function
int	j, k
real	fcs_hyper()

begin
	# Copy parameters to local array.  We do this so we can have some
	# variable parameters and some fixed ones.  The SP array contains
	# only the variable parameters, and we need them all.
	j = 0
	do k = 1, F_NPAR(fnl) {
	    if (F_VAR(fnl,k)) {
		j = j + 1
		p[k] = vp[j]
	    } else {
		p[k] = F_PAR(fnl,k)
	    }
	}

	# Check for unreasonable parameters.
	if (p[1] < 2 * F_XMIN(fnl) - F_XMAX(fnl))	# x too small?
	    return (BIGVALUE)
	if (p[1] > 2 * F_XMAX(fnl) - F_XMIN(fnl))	# x too big?
	    return (BIGVALUE)
	if (p[2] <= 0.)					# a too small?
	    return (BIGVALUE)
	if (p[3] <= 0.)					# b too small?
	    return (BIGVALUE)

	# Accumulate the sum of squared differences (or the sum of the
	# absolute values) from the hyperbolic function.
	sum = 0.
	do k = 1, F_NPTS(fnl) {

	    h = fcs_hyper (F_X(fnl,k), p)		# hyperbolic function

	    if (F_OPTION(fnl) == F_LEAST_SQUARES)
		sum = sum + (h - F_Y(fnl,k)) ** 2
	    else if (F_OPTION(fnl) == F_MINSUM)
		sum = sum + abs (h - F_Y(fnl,k))
	}

	return (sum)
end

# fcs_hyper -- evaluate a hyperbolic function
# This routine evaluates a hyperbolic function  using parameters passed in
# the calling sequence.

real procedure fcs_hyper (x, p)

real	x		# i: value at which function is to be evaluated
real	p[4]		# i: parameters for hyperbolic function
#--
real	arg		# argument for hyperbolic function

begin
	arg = (x - p[1]) / p[3]
	arg = arg * arg + 1.
	return (p[2] * sqrt (arg) + p[4])
end