#======================================================================
#                    . L S F I T . B I L I N 
#                    doc: Wed Feb 24 09:40:06 1999
#                    dlm: Fri Jul 28 13:36:36 2006
#                    (c) 1999 A.M. Thurnherr
#                    uE-Info: 32 41 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================

# What you need to provide if you wanna fit a different
# linear model function to your data:
#	- a number of global variables to be set during loading
#	- a number of subs to perform admin tasks (usage, init, ...)
#	- a sub to evaluate the basis funs at a given x value; each
#	  y value must be stored in @A (beginning with
#	  A[1]!!!).
#	- the interface is documented between +++++++ lines

# fit bi-linear function to data, i.e. y = A + B*x1 + C*x2

# HISTORY:
#	Aug 01, 1999: - adapted from [.lsfit.poly]
#	Aug 02, 1999: - added &antsDescription()
#	Sep 26, 1999: - cosmetics
#				  - added vars & covars
#	Sep 27, 1999: - changed from covar to sigmas
#	Oct 01, 1999: - cosmetics
#	Oct 06, 1999: - added -l
#   Mar 17, 2001: - param->arg
#	May  2, 2001: - updated doc
#	Nov 17, 2005: - commented out antsDescription()
#				  - updated stats on -p
#	Jul 28, 2006: - Version 3.3 [HISTORY]

# NOTES:
#	- could be easily extended to multidimensional linear fit but
#	  what's the use?
#	- -p zeroes param C for bilinear spice method (T = A + B*sig + C*neph)
#	  so that residual field becomes spice anomaly (linearly correlated
#	  with neph)

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
#
#	$modelOpts			string of allowed options
#	$modelOptsUsage		usage information string for options
#	$modelMinArgs		min # of arguments of model
#	$modelArgsUsage		usage information string for arguments
#
# The following variables may be set later but not after &modelInit()
#
#	$modelNFit			number of params to fit in model
#	@nameA				symbolic names of model parameters
#
# You should call &antsDescription() for the -ct options here
#
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

$modelOpts = "pl:";
$modelOptsUsage = "[-p)artial f/r] [-l)imit r_BC <min,max>]";
$modelMinArgs = 1;
$modelArgsUsage = "<2nd x-field>";
$modelNFit = 3;
$nameA[1] = "A"; $nameA[2] = "B"; $nameA[3] = "C";
$A[1] = nan; $A[2] = nan; $A[3] = nan;
#&antsDescription("c","bilin_$nameA[1]",
#				 "c","bilin_$nameA[2]",
#				 "c","bilin_$nameA[3]");
#&antsDescription("t","bilin_sigma_$nameA[1]",
#				 "t","bilin_sigma_$nameA[2]",
#				 "t","bilin_sigma_$nameA[3]",
#				 "t","bilin_ccc_$nameA[1]_$nameA[2]",
#				 "t","bilin_ccc_$nameA[1]_$nameA[3]",
#				 "t","bilin_ccc_$nameA[2]_$nameA[3]");

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
# &modelUsage()		mangle parameters; NB: there may be `infinite' # of
#					filenames after model arguments; this usually sets
#					@A (the model parameters) but these can later be
#					calculated heuristically during &modelInit()
#
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

sub modelUsage()
{
	if (defined($opt_l)) {
		($minBC,$maxBC) = split(',',$opt_l);
		&antsUsageError("\n>>> error with -l")
			unless (defined($maxBC) && $maxBC > $minBC);
	}
	$x2fnr = &antsFieldArg();
	&antsUsageError() unless ($#ARGV < 0 || -r $ARGV[0]);
}

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
# &modelInit()		initializes model after reading of data
#
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

sub modelInit() {}

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#
# &modelEvaluate(idx,xfnr,vals)	evaluate basis funs at x (NB: x1, x2)
#		idx				       	current index in @ants_
#		xfnr					field number of x field
#		vals					reference to return values (1-relative!)
#
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

sub modelEvaluate($$$)
{
	my($idx,$xfnr,$valsR) = @_;
	my($i);

	$valsR->[1] = 1;
	$valsR->[2] = $ants_[$idx][$xfnr];
	$valsR->[3] = $ants_[$idx][$x2fnr];
}

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# &modelCleanup()	cleans up after fitting but before output
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

sub modelCleanup() 
{
	if (defined($opt_l)) {
		my($ccc) = $covar[2][3]/(sqrt($covar[2][2])*sqrt($covar[3][3]));
		if ($ccc < $minBC || $ccc > $maxBC) {
			&antsInfo("CCC B/C = %.3g out of range, fit discarded",$ccc);
			$suppressFit = 1;
		}
	}
	if ($opt_p) {
		&antsInfo("parameter $nameA[3] = %.3g discarded",$A[3]);
		$A[3] = 'discarded';
		$covar[3][3] = $RMS = $sig = nan;
	}
}