1 /**
2 @file
3 @ingroup nr
4 */
5 #include <bpm/bpm_messages.h>
6 #include <bpm/bpm_nr.h>
7
8
9 /**
10 Fit data to a straight line. Nicked from numerical recipes, C15.2, p665
11 See: http://www.library.cornell.edu/nr/cbookcpdf.html
12
13 @param x array with x values
14 @param y array with corresponding y values
15 @param ndata number of datapoints
16 @param sig array with errors on y datapoints
17 @param mwt used weighted (so including errors on datapoints ?)
18 @param a fitted slope
19 @param b fitted intercept
20 @param siga error on fitted slope
21 @param sigb error on fitted intercept
22 @param chi2 chi2 of fit
23 @param q quality factor of fit
24
25 @return BPM_FAILURE upon failure, BPM_SUCCESS upon success
26 */
nr_fit(double * x,double y[],int ndata,double sig[],int mwt,double * a,double * b,double * siga,double * sigb,double * chi2,double * q)27 int nr_fit ( double *x, double y[], int ndata, double sig[], int mwt,
28 double *a, double *b, double *siga, double *sigb,
29 double *chi2, double *q ) {
30 int i;
31 double wt, t, sxoss, sx=0.0, sy=0.0, st2=0.0, ss, sigdat;
32
33 if( ! ( x && y ) ) {
34 bpm_error( "Invalid arguments in nr_fit(...)",
35 __FILE__, __LINE__ );
36 return BPM_FAILURE;
37 }
38
39 if( ( mwt ) && ! sig ) {
40 bpm_error( "Invalid arguments using sig[] in nr_fit(...)",
41 __FILE__, __LINE__ );
42 return BPM_FAILURE;
43 }
44
45 if( ndata < 3 ) {
46 bpm_error( "Number of datapoints to small (<3) in nr_fit(...)",
47 __FILE__, __LINE__ );
48 return BPM_FAILURE;
49 }
50
51 *b=0.0;
52 if ( mwt ) {
53 ss=0.0;
54
55 for ( i=0; i<ndata;i++ ) {
56 if( ABS( sig[i] ) > 0. ) {
57 wt=1.0/SQR(sig[i]);
58 } else {
59 bpm_error( "sig[] contains 0 values in nr_fit(...)",
60 __FILE__, __LINE__ );
61 return BPM_FAILURE;
62 }
63 ss += wt;
64 sx += x[i]*wt;
65 sy += y[i]*wt;
66 }
67
68 } else {
69
70 for (i=0; i<ndata;i++) {
71 sx += x[i];
72 sy += y[i];
73 }
74
75 ss=ndata;
76 } /* if ( mwt ) */
77
78 if( ABS(ss) > 0. ) {
79 sxoss=sx/ss;
80 } else {
81 bpm_error( "ss is zero in nr_fit(...)", __FILE__, __LINE__ );
82 return BPM_FAILURE;
83 }
84
85 if ( mwt ) {
86
87 for (i=0; i<ndata;i++) {
88 t=(x[i]-sxoss)/sig[i];
89 st2 += t*t;
90 *b += t*y[i]/sig[i];
91 }
92
93 } else {
94 for (i=0; i<ndata; i++) {
95 t=x[i]-sxoss;
96 st2 += t*t;
97 *b += t*y[i];
98 }
99 } /* if ( mwt ) */
100
101 if( ABS(st2) > 0. ) {
102 *b /= st2;
103 *a = (sy-sx*(*b))/ss;
104 *siga=sqrt((1.0+sx*sx/(ss*st2))/ss);
105 *sigb=sqrt(1.0/st2);
106 } else {
107 bpm_error( "st2 is zero in nr_fit(...)", __FILE__, __LINE__ );
108 return BPM_FAILURE;
109 }
110
111 *chi2=0.0;
112 *q=1.0;
113
114 if ( mwt == 0 ) {
115 for (i=0;i<ndata;i++)
116 *chi2 +=SQR(y[i]-(*a)-(*b)*x[i]);
117
118 sigdat=sqrt((*chi2)/(ndata-2));
119 *siga *= sigdat;
120 *sigb *= sigdat;
121
122 } else {
123
124 for (i=0; i<ndata;i++) *chi2 += SQR((y[i]-(*a)-(*b)*x[i])/sig[i]);
125
126 *q = nr_gammq(0.5*(ndata-2), 0.5*(*chi2));
127
128 } /* if ( mwt == 0 ) */
129
130
131 return BPM_SUCCESS;
132 }
133
134