// Routine for fitting to an exponential (for loading curves) // Using Numerical Recipes in C, 15.5 Nonlinear Models Example #include #include #include #include // NUMERICAL RECIPES ROUTINES - CHECK PATHNAMES BEFORE USING #include "/home/gensemer/NumRec/Other/NRUTIL.C" #include "/home/gensemer/NumRec/Other/nr.h" #include "/home/gensemer/NumRec/ANSI/MRQCOF.C" #include "/home/gensemer/NumRec/ANSI/MRQMIN.C" #include "/home/gensemer/NumRec/ANSI/GAUSSJ.C" #include "/home/gensemer/NumRec/ANSI/COVSRT.C" // USER-SUPPLIED FUNCTION HERE #include "/home/gensemer/datanalysis/analysis/betafit/beta.cpp" void betafit (char *filename, char *outfile, int verbose) { // OPEN INPUT AND OUTPUT FILES--------------------------------- ifstream file (filename); if (!file) cerr << "Cannot open infile for input"; ofstream file2 (outfile); if (!file2) cerr << "Cannot open outfile for output"; int i, xmax, j; float *y,*x; y = vector (1, 25000); x = vector (1, 25000); // LOAD DATA FROM FILE---------------------------------------- i = 0; while (file.good () == 1 && i <= 25000) { i++; file >> y[i]; } const int ndata = i-1; if (i == 25000) cout << "Data Overflow (file contains > 25,000 points)\n"; // DECLARE VARIABLES AND ARRAYS--------------------------------- const int one = 1; const int ma = 3; // number of fitting params in func float *yfit, *sig, *a, *dyda; yfit = vector (one, ndata); sig = vector (one, ndata); a = vector (one, ma); dyda = vector (one, ma); int *ia; ia = ivector (one, ma); float **covar, **alpha; covar = matrix (one, ma, one, ma); alpha = matrix (one, ma, one, ma); void (*funcs) (); float chi, oldchi, lamda, oldlamda, *chisq, *alamda; alamda = &lamda; chisq = χ // SET x[], sig[], AND ia[] VALUES--------------------------------- for (i = 1; i <= ndata; i++) x[i] = i/10; // x axis values (optional) for (i = 1; i <= ndata; i++) sig[i] = 0.02; // estimated experimental uncertainty for (i = 1; i <= ma; i++) ia[i] = 1; // weighting // GUESS FIT PARAMETERS--------------------------------- double v1,v2; v1 = 0; v2 = 0; a[1] = y[1]; // Initial value for (i=1;i<=5;i++) v1 += y[4*ndata/5-i]/5; for (i=1;i<=5;i++) v2 += y[ndata-i]/5; a[2] = log(v1/v2)/(ndata/50); // Gamma (exponential time constant) a[3] = a[2]*(1/v2 - exp(a[2]*ndata/10)/a[1])/(exp(a[2]*ndata/10) - 1); //Beta (non-exponential correction) // INITIALIZE ROUTINE--------------------------------- *alamda = -1; if (verbose == 1) { cout << "i alamda chisq amatrix\n"; cout << -1 << " " << *alamda << " " << *chisq << " "; for (j = 1; j <= ma; j++) cout << a[j] << " "; cout << "\n"; } i = 0; mrqmin (x, y, sig, ndata, a, ia, ma, covar, alpha, chisq, beta, alamda); if (verbose == 1) { cout << i << " " << *alamda << " " << *chisq << " "; for (j = 1; j <= ma; j++) cout << a[j] << " "; cout << "\n"; } i++; // FIT WHILE CONDITION--------------------------------- oldchi = chi * 2; while (((oldchi - chi) / oldchi > .001 || chi > oldchi) && i < 40) { oldchi = chi; mrqmin (x, y, sig, ndata, a, ia, ma, covar, alpha, chisq, beta, alamda); if (verbose == 1) { cout << i << " " << *alamda << " " << *chisq << " "; for (j = 1; j <= ma; j++) cout << a[j] << " "; cout << "\n"; } i++; } if (i == 40) cout << "No Convergence\n"; // FINAL CALL TO DETERMINE CONVERGED VALUES---------------------------- *alamda = 0; mrqmin (x, y, sig, ndata, a, ia, ma, covar, alpha, chisq, beta, alamda); // OUTPUT RESULTS TO STDOUT (pipe to file if needed) if (verbose == 1) file2 << "Results:" << filename << "\n i a[i] error[i]\n"; for (i = 1; i <= ma; i++) file2 << i << " " << a[i] << " " << sqrt (3.53 * covar[i][i]) << "\n"; }