00001
00005 #include <bpm/bpm_messages.h>
00006 #include <bpm/bpm_nr.h>
00007
00008
00027 int nr_fit ( double *x, double y[], int ndata, double sig[], int mwt,
00028 double *a, double *b, double *siga, double *sigb,
00029 double *chi2, double *q ) {
00030 int i;
00031 double wt, t, sxoss, sx=0.0, sy=0.0, st2=0.0, ss, sigdat;
00032
00033 if( ! ( x && y ) ) {
00034 bpm_error( "Invalid arguments in nr_fit(...)",
00035 __FILE__, __LINE__ );
00036 return BPM_FAILURE;
00037 }
00038
00039 if( ( mwt ) && ! sig ) {
00040 bpm_error( "Invalid arguments using sig[] in nr_fit(...)",
00041 __FILE__, __LINE__ );
00042 return BPM_FAILURE;
00043 }
00044
00045 if( ndata < 3 ) {
00046 bpm_error( "Number of datapoints to small (<3) in nr_fit(...)",
00047 __FILE__, __LINE__ );
00048 return BPM_FAILURE;
00049 }
00050
00051 *b=0.0;
00052 if ( mwt ) {
00053 ss=0.0;
00054
00055 for ( i=0; i<ndata;i++ ) {
00056 if( ABS( sig[i] ) > 0. ) {
00057 wt=1.0/SQR(sig[i]);
00058 } else {
00059 bpm_error( "sig[] contains 0 values in nr_fit(...)",
00060 __FILE__, __LINE__ );
00061 return BPM_FAILURE;
00062 }
00063 ss += wt;
00064 sx += x[i]*wt;
00065 sy += y[i]*wt;
00066 }
00067
00068 } else {
00069
00070 for (i=0; i<ndata;i++) {
00071 sx += x[i];
00072 sy += y[i];
00073 }
00074
00075 ss=ndata;
00076 }
00077
00078 if( ABS(ss) > 0. ) {
00079 sxoss=sx/ss;
00080 } else {
00081 bpm_error( "ss is zero in nr_fit(...)", __FILE__, __LINE__ );
00082 return BPM_FAILURE;
00083 }
00084
00085 if ( mwt ) {
00086
00087 for (i=0; i<ndata;i++) {
00088 t=(x[i]-sxoss)/sig[i];
00089 st2 += t*t;
00090 *b += t*y[i]/sig[i];
00091 }
00092
00093 } else {
00094 for (i=0; i<ndata; i++) {
00095 t=x[i]-sxoss;
00096 st2 += t*t;
00097 *b += t*y[i];
00098 }
00099 }
00100
00101 if( ABS(st2) > 0. ) {
00102 *b /= st2;
00103 *a = (sy-sx*(*b))/ss;
00104 *siga=sqrt((1.0+sx*sx/(ss*st2))/ss);
00105 *sigb=sqrt(1.0/st2);
00106 } else {
00107 bpm_error( "st2 is zero in nr_fit(...)", __FILE__, __LINE__ );
00108 return BPM_FAILURE;
00109 }
00110
00111 *chi2=0.0;
00112 *q=1.0;
00113
00114 if ( mwt == 0 ) {
00115 for (i=0;i<ndata;i++)
00116 *chi2 +=SQR(y[i]-(*a)-(*b)*x[i]);
00117
00118 sigdat=sqrt((*chi2)/(ndata-2));
00119 *siga *= sigdat;
00120 *sigb *= sigdat;
00121
00122 } else {
00123
00124 for (i=0; i<ndata;i++) *chi2 += SQR((y[i]-(*a)-(*b)*x[i])/sig[i]);
00125
00126 *q = nr_gammq(0.5*(ndata-2), 0.5*(*chi2));
00127
00128 }
00129
00130
00131 return BPM_SUCCESS;
00132 }
00133