|
1 #====================================================================== |
|
2 # . / S V D F I T . P L |
|
3 # doc: Sat Jul 31 22:09:25 1999 |
|
4 # dlm: Sat Jul 31 22:45:40 1999 |
|
5 # (c) 1999 A.M. Thurnherr |
|
6 # uE-Info: 61 36 NIL 0 0 72 2 2 4 ofnI |
|
7 #====================================================================== |
|
8 |
|
9 # SVDFIT routine from Numerical Recipes adapted to ANTS |
|
10 |
|
11 # HISTORY: |
|
12 # Jul 31, 1999: - manually converted from c-source |
|
13 |
|
14 # Notes: |
|
15 # - x,y,sig are field numbers for data in $ants_ |
|
16 # - if sig is a negative number, -sig is used as constant input stddev |
|
17 # - @a, @u, @v, @w, &funcs passed as refs |
|
18 # - chi square is returned |
|
19 |
|
20 require "$ANTS/nrutil.pl"; |
|
21 require "$ANTS/svbksb.pl"; |
|
22 require "$ANTS/svdcmp.pl"; |
|
23 |
|
24 { # BEGIN static scope |
|
25 |
|
26 my($TOL) = 1.0e-5; |
|
27 |
|
28 sub svdfit($$$$$$$$) |
|
29 { |
|
30 my($xfnr,$yfnr,$sig,$aR,$uR,$vR,$wR,$funcsR) = @_; |
|
31 my($j,$i); # int |
|
32 my($chisq,$wmax,$tmp,$thresh,$sum); # float |
|
33 my(@b,@afunc); # float[] |
|
34 |
|
35 &vector(\@b,1,$#ants_); |
|
36 &vector(\@afunc,1,$#{$aR}); |
|
37 for ($i=0; $i<=$#ants_; $i++) { |
|
38 next if ($antsFlagged[$i]); |
|
39 &$funcsR($ants_[$i][$xfnr],\@afunc); |
|
40 $tmp = 1.0 / (($sig > 0) ? $ants_[$i][$sig] : -$sig); |
|
41 for ($j=1; $j<=$#{$aR}; $j++) { |
|
42 $uR->[$i][$j] = $afunc[$j]*$tmp; |
|
43 } |
|
44 $b[$i] = $ants_[$i][$yfnr]*$tmp; |
|
45 } |
|
46 &svdcmp($uR,$wR,$vR); |
|
47 for ($j=1; $j<=$#{$aR}; $j++) { |
|
48 $wmax = $wR->[$j] if ($wR->[$j] > $wmax); |
|
49 } |
|
50 $thresh = $TOL*$wmax; |
|
51 for ($j=1; $j<=$#{$aR}; $j++) { |
|
52 $wR->[$j] = 0 if ($wR->[$j] < $thresh); |
|
53 } |
|
54 &svbksb($uR,$wR,$vR,\@b,$aR); |
|
55 for ($i=0; $i<=$#ants_; $i++) { |
|
56 next if ($antsFlagged[$i]); |
|
57 &$funcsR($ants_[$i][$xfnr],\@afunc); |
|
58 for ($j=1; $j<=$#{$aR}; $j++) { |
|
59 $sum += $aR->[$j]*$afunc[$j]; |
|
60 } |
|
61 $tmp = ($ants_[$i][$yfnr] - $sum) / |
|
62 (($sig > 0) ? $ants_[$i][$sig] : -$sig); |
|
63 $chisq += $tmp * $tmp; |
|
64 } |
|
65 return $chisq; |
|
66 } |
|
67 |
|
68 } # END static scope |