#======================================================================

#                    . / S V D F I T . P L 

#                    doc: Sat Jul 31 22:09:25 1999

#                    dlm: Sat Jul 31 22:45:40 1999

#                    (c) 1999 A.M. Thurnherr

#                    uE-Info: 61 36 NIL 0 0 72 2 2 4 ofnI

#======================================================================

# SVDFIT routine from Numerical Recipes adapted to ANTS

# HISTORY:

#	Jul 31, 1999: - manually converted from c-source

# Notes:

#   - x,y,sig are field numbers for data in $ants_

#   - if sig is a negative number, -sig is used as constant input stddev

#   - @a, @u, @v, @w, &funcs passed as refs

#	- chi square is returned

require "$ANTS/nrutil.pl";

require "$ANTS/svbksb.pl";

require "$ANTS/svdcmp.pl";

{ # BEGIN static scope

my($TOL) = 1.0e-5;

sub svdfit($$$$$$$$)

{

	my($xfnr,$yfnr,$sig,$aR,$uR,$vR,$wR,$funcsR) = @_;

	my($j,$i);										# int

	my($chisq,$wmax,$tmp,$thresh,$sum);					# float

	my(@b,@afunc);									# float[]

	&vector(\@b,1,$#ants_);

	&vector(\@afunc,1,$#{$aR});

	for ($i=0; $i<=$#ants_; $i++) {

		next if ($antsFlagged[$i]);

		&$funcsR($ants_[$i][$xfnr],\@afunc);

		$tmp = 1.0 / (($sig > 0) ? $ants_[$i][$sig] : -$sig);

		for ($j=1; $j<=$#{$aR}; $j++) {

			$uR->[$i][$j] = $afunc[$j]*$tmp;

		}

		$b[$i] = $ants_[$i][$yfnr]*$tmp;

	}

	&svdcmp($uR,$wR,$vR);

	for ($j=1; $j<=$#{$aR}; $j++) {

		$wmax = $wR->[$j] if ($wR->[$j] > $wmax);

	}

	$thresh = $TOL*$wmax;

	for ($j=1; $j<=$#{$aR}; $j++) {

		$wR->[$j] = 0 if ($wR->[$j] < $thresh);

	}

	&svbksb($uR,$wR,$vR,\@b,$aR);

	for ($i=0; $i<=$#ants_; $i++) {

		next if ($antsFlagged[$i]);

		&$funcsR($ants_[$i][$xfnr],\@afunc);

		for ($j=1; $j<=$#{$aR}; $j++) {

			$sum += $aR->[$j]*$afunc[$j];

		}

		$tmp = ($ants_[$i][$yfnr] - $sum) /

				(($sig > 0) ? $ants_[$i][$sig] : -$sig);

		$chisq += $tmp * $tmp;

	}

	return $chisq;

}

} # END static scope