2 #include "AliL3FitUtilities.h"
15 /******************************************************************************/
22 if( !(fp = fopen( "exp_table", "wb" )))
23 printf("build_exp_table: I/O error\n");
25 tab_h.float_size = sizeof( FLOAT_SIZE );
29 tab_h.step_size = (tab_h.dmax - tab_h.dmin)/tab_h.steps;
30 if( !(exp_val = (FLOAT_SIZE *)malloc( tab_h.steps * sizeof( FLOAT_SIZE ) ) ) )
31 printf("build_exp_table: malloc error\n");
32 fwrite((const void *)&tab_h, (size_t)sizeof(exp_header_t), (size_t)1, fp );
33 for( i=0; i<tab_h.steps; ++i ) {
34 exp_val[i] = exp(tab_h.dmin + i*tab_h.step_size);
36 fwrite((const void *)exp_val, (size_t)sizeof(FLOAT_SIZE), (size_t)tab_h.steps, fp );
38 free( (void *)exp_val );
44 /**********************************************************/
46 double i, delta, dmin, dmax, exp_tab();
49 fp = fopen( "exp_test", "w" );
54 for( i=dmin; i<dmax; i+=delta ){
55 fprintf( fp, "%lf\t%lf\t%lf\n", i, exp(i), exp_tab(i) );
60 double exp_tab(double in)
66 /**********************************************************/
69 static exp_header_t tab_h;
70 static FLOAT_SIZE *exp_val=NULL, slope;
75 if( (num_nan /100)*100 == num_nan )
76 fprintf( stderr, "exp_tab: NaN %d\n", num_nan );
85 if( !(fp = fopen( "exp_table", "rb" ))) {
88 if( !(fp = fopen( "exp_table", "rb" )))
89 printf("exp_tab: I/O error\n");
90 fread(&tab_h, (size_t)sizeof(exp_header_t), (size_t)1, fp );
91 if( tab_h.float_size != sizeof( FLOAT_SIZE ) )
94 if( !(exp_val = (FLOAT_SIZE *)malloc( tab_h.steps * sizeof( FLOAT_SIZE ) ) ) )
95 printf("exp_tab: malloc error\n");
96 fread(exp_val, (size_t)sizeof(FLOAT_SIZE), (size_t)tab_h.steps, fp );
97 slope = tab_h.steps / (tab_h.dmax - tab_h.dmin);
99 if( in < tab_h.dmin || in > tab_h.dmax )
102 return exp_val[(int)((in-tab_h.dmin) * slope)];
107 void f2gauss5(double x,double a[],double *y,double dyda[],int na)
110 double fac,fac1,fac2,ex,ex1,ex2,arg1,arg2,u,v;
114 if( index < 0 || index >=FIT_PTS )
116 fprintf( stderr, "ff2gauss: wrong index %ld\n", index );
123 for (i=1;i<=na-1;i+=5)
125 arg1 = (u-a[i+1])/a[i+2];
126 arg2 = (v-a[i+3])/a[i+4];
127 ex1 = exp_tab(-arg1*arg1);
128 ex2 = exp_tab(-arg2*arg2);
135 dyda[i+1] = fac1/a[i+2];
136 dyda[i+2] = fac1*arg1/a[i+2];
137 dyda[i+3] = fac2/a[i+4];
138 dyda[i+4] = fac2*arg2/a[i+4];
142 void nrerror(char error_text[])
143 /* Numerical Recipes standard error handler */
145 printf("%s\n",error_text);
149 void free_vector(double *v, long nl, long nh)
150 /* free a double vector allocated with vector() */
152 free((FREE_ARG) (v+nl-NR_END));
155 void free_ivector(int *v, long nl, long nh)
156 /* free an int vector allocated with ivector() */
158 free((FREE_ARG) (v+nl-NR_END));
161 void free_matrix(double **m, long nrl, long nrh, long ncl, long nch)
162 /* free a double matrix allocated by matrix() */
164 free((FREE_ARG) (m[nrl]+ncl-NR_END));
165 free((FREE_ARG) (m+nrl-NR_END));
168 int *ivector(long nl, long nh)
169 /* allocate an int vector with subscript range v[nl..nh] */
173 v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
174 if (!v) nrerror("allocation failure in ivector()");
178 double *vector(long nl,long nh)
179 /* allocate a double vector with subscript range v[nl..nh] */
183 v=(double *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(double)));
184 if (!v) nrerror("allocation failure in vector()");
188 double **matrix(long nrl,long nrh,long ncl,long nch)
189 /* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
191 long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
194 /* allocate pointers to rows */
195 m=(double **) malloc((unsigned int)((nrow+NR_END)*sizeof(double*)));
196 if (!m) nrerror("allocation failure 1 in matrix()");
200 /* allocate rows and set pointers to them */
201 m[nrl]=(double *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(double)));
202 if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
206 for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
208 /* return pointer to array of pointers to rows */
212 void gaussj(double **a, int n, double **b, int m)
214 int *indxc,*indxr,*ipiv;
215 int i,icol,irow,j,k,l,ll;
216 double big,dum,pivinv,swap;
221 for (j=1;j<=n;j++) ipiv[j]=0;
228 if (fabs(a[j][k]) >= big) {
233 } else if (ipiv[k] > 1) {
234 free_ivector(ipiv,1,n);
235 free_ivector(indxr,1,n);
236 free_ivector(indxc,1,n);
237 nrerror("gaussj: Singular Matrix-1");
242 for (l=1;l<=n;l++) SWAP(a[irow][l],a[icol][l])
243 for (l=1;l<=m;l++) SWAP(b[irow][l],b[icol][l])
247 if (a[icol][icol] == 0.0){
248 free_ivector(ipiv,1,n);
249 free_ivector(indxr,1,n);
250 free_ivector(indxc,1,n);
251 nrerror("gaussj: Singular Matrix-2");
253 if (mabs(a[icol][icol]) < EPSILON )
254 nrerror("gaussj: a[icol][icol] == 0");
255 pivinv=1.0/a[icol][icol];
257 pivinv=1.0/a[icol][icol];
259 for (l=1;l<=n;l++) a[icol][l] *= pivinv;
260 for (l=1;l<=m;l++) b[icol][l] *= pivinv;
261 for (ll=1;ll<=n;ll++)
265 for (l=1;l<=n;l++) a[ll][l] -= a[icol][l]*dum;
266 for (l=1;l<=m;l++) b[ll][l] -= b[icol][l]*dum;
270 if (indxr[l] != indxc[l])
272 SWAP(a[k][indxr[l]],a[k][indxc[l]]);
274 free_ivector(ipiv,1,n);
275 free_ivector(indxr,1,n);
276 free_ivector(indxc,1,n);
279 void covsrt(double **covar, int ma, int ia[], int mfit)
284 for (i=mfit+1;i<=ma;i++)
285 for (j=1;j<=i;j++) covar[i][j]=covar[j][i]=0.0;
287 for (j=ma;j>=1;j--) {
289 for (i=1;i<=ma;i++) SWAP(covar[i][k],covar[i][j])
290 for (i=1;i<=ma;i++) SWAP(covar[k][i],covar[j][i])
296 void mrqcof(double x[], double y[], double sig[], int ndata, double a[], int ia[],
297 int ma, double **alpha, double beta[], double *chisq,
298 void (*funcs)(double, double [], double *, double [], int))
300 int i,j,k,l,m,mfit=0;
301 double ymod,wt,sig2i,dy,*dyda;
307 for (j=1;j<=mfit;j++) {
308 for (k=1;k<=j;k++) alpha[j][k]=0.0;
312 for (i=1;i<=ndata;i++) {
313 (*funcs)(x[i],a,&ymod,dyda,ma);
314 sig2i=1.0/(sig[i]*sig[i]);
316 for (j=0,l=1;l<=ma;l++) {
319 for (j++,k=0,m=1;m<=l;m++)
320 if (ia[m]) alpha[j][++k] += wt*dyda[m];
324 *chisq += dy*dy*sig2i;
326 for (j=2;j<=mfit;j++)
327 for (k=1;k<j;k++) alpha[k][j]=alpha[j][k];
328 free_vector(dyda,1,ma);
331 void mrqmin(double x[], double y[], double sig[], int ndata, double a[], int ia[],
332 int ma, double **covar, double **alpha, double *chisq,
333 void (*funcs)(double, double [], double *, double [], int), double *alamda)
335 void covsrt(double **covar, int ma, int ia[], int mfit);
336 void gaussj(double **a, int n, double **b, int m);
337 void mrqcof(double x[], double y[], double sig[], int ndata, double a[],
338 int ia[], int ma, double **alpha, double beta[], double *chisq,
339 void (*funcs)(double, double [], double *, double [], int));
342 static double ochisq,*atry,*beta,*da,**oneda;
348 for (mfit=0,j=1;j<=ma;j++)
350 oneda=matrix(1,mfit,1,1);
352 mrqcof(x,y,sig,ndata,a,ia,ma,alpha,beta,chisq,funcs);
354 for (j=1;j<=ma;j++) atry[j]=a[j];
356 for (j=0,l=1;l<=ma;l++) {
358 for (j++,k=0,m=1;m<=ma;m++) {
361 covar[j][k]=alpha[j][k];
364 covar[j][j]=alpha[j][j]*(1.0+(*alamda));
368 gaussj(covar,mfit,oneda,1);
369 for (j=1;j<=mfit;j++) da[j]=oneda[j][1];
370 if (*alamda == 0.0) {
371 covsrt(covar,ma,ia,mfit);
372 free_matrix(oneda,1,mfit,1,1);
373 free_vector(da,1,ma);
374 free_vector(beta,1,ma);
375 free_vector(atry,1,ma);
378 for (j=0,l=1;l<=ma;l++)
379 if (ia[l]) atry[l]=a[l]+da[++j];
380 mrqcof(x,y,sig,ndata,atry,ia,ma,covar,da,chisq,funcs);
381 if (*chisq < ochisq) {
384 for (j=0,l=1;l<=ma;l++) {
386 for (j++,k=0,m=1;m<=ma;m++) {
389 alpha[j][k]=covar[j][k];
403 int lev_marq_fit( double x[], double y[], double sig[], int NPT, double a[], int ia[], double dev[], int MA,
404 double *chisq_p, void (*funcs)(double, double [], double *, double [], int) )
407 double alamda,chisq,ochisq,**covar,**alpha;
409 covar=matrix(1,MA,1,MA);
410 alpha=matrix(1,MA,1,MA);
412 if( setjmp(env) == 0 ) {
414 mrqmin(x,y,sig,NPT,a,ia,MA,covar,alpha,&chisq,funcs,&alamda);
420 mrqmin(x,y,sig,NPT,a,ia,MA,covar,alpha,&chisq,funcs,&alamda);
423 else if (fabs(ochisq-chisq) < 0.1)
425 if (itst < 4) continue;
427 mrqmin(x,y,sig,NPT,a,ia,MA,covar,alpha,&chisq,funcs,&alamda);
430 dev[i] = sqrt(covar[i][i]);
433 free_matrix(alpha,1,MA,1,MA);
434 free_matrix(covar,1,MA,1,MA);
438 //if( control_g.print_fit_errors==2 )
439 fprintf( stderr, " runtime error\n" );
441 free_matrix(alpha,1,MA,1,MA);
442 free_matrix(covar,1,MA,1,MA);