// Filename: fitfunc.cpp // Author: Michael P. Schubmehl // Date(s): July 2002 // Purpose: This file contains the implementaion for the fittingFunction // class described in fitfunc.h. // // Notes: Updated version uses class inheritance. See fitfunc.h for details. #include // Streams #include // Output formatting #include // File streams #include // String manipulation routines #include // For exit() and NULL #include // Absolute value, exp, etc. #include "fitfunc.h" // Fitting functions #ifndef PI // If there is no numerical #define PI 3.1415926535 // value for pi yet, define #endif // it to be 3.1415926535. #define R 0 // Name indices for easy reference #define ALPHA 1 // later. For example, the #define SIGMA 2 // first parameter in each fit // function is r, which can be // accessed by // _currentParams[R] or // _currentParams[0]. // Table of Contents *********************************************************** // class fittingFunction; // The abstract fitting function // // class deltaFunction; // Derived classes for several // class gaussianFunction; // common types of functions. // class logonormalFunction; double min(double x, double y); // Returns the minimum of x and y double max(double x, double y); // Returns the maximum of x and y // END Table of Contents ******************************************************* // class fittingFunction ******************************************************* fittingFunction::~fittingFunction() { // Input: None. // Output: Frees up space allocated for this fitting function by deleting // the arrays used for the parameters. _numParams = 0; delete[] _minimumParams; delete[] _maximumParams; delete[] _currentParams; } const int fittingFunction::numParams() { // Input: None. // Output: Returns a constant integer equal to the number of parameters // present in this particular fitting function. The user cannot // modify this value directly. return _numParams; } double& fittingFunction::minimumParams(int i) { // Input: An integer i. // Output: Returns by reference the array element _minimumParams[i], or // prints an out-of-range error and exits with error code 2. // The return by reference allows direct access to the elements // of the array through syntax like f->minimumParams(1) = 2; if((i < _numParams) && (i >= 0)) { return _minimumParams[i]; } else { cerr << "Error: Attempt to access parameter " << i << " in a list with "; cerr << _numParams << " elements." << endl; exit(2); } } double& fittingFunction::maximumParams(int i) { // Input: An integer i. // Output: Returns by reference the array element _maximumParams[i], or // prints an out-of-range error and exits with error code 2. if((i < _numParams) && (i >= 0)) { return _maximumParams[i]; } else { cerr << "Error: Attempt to access parameter " << i << " in a list with "; cerr << _numParams << " elements." << endl; exit(2); } } double& fittingFunction::currentParams(int i) { // Input: An integer i. // Output: Returns by reference the array element _currentParams[i], or // prints an out-of-range error and exits with error code 2. if((i < _numParams) && (i >= 0)) { return _currentParams[i]; } else { cerr << "Error: Attempt to access parameter " << i << " in a list with "; cerr << _numParams << " elements." << endl; exit(2); } } bool& fittingFunction::fitting(int i) { // Input: An integer i. // Output: Returns by reference the boolean _fitting[i]. if((i < _numParams) && (i >= 0)) { return _fitting[i]; } else { cerr << "Error: Attempt to access parameter " << i << " in a list with "; cerr << _numParams << " elements." << endl; exit(2); } } ostream& operator<<(ostream& stream, fittingFunction& source) { // Input: A fitting function and a stream to write it to. // Output: Writes the fitting function to the stream and returns the stream. source.output(stream); // Output source func to stream stream << setprecision(6); // Reset to default precision return stream; // Return stream for next << } // END class fittingFunction *************************************************** // class deltaFunction ********************************************************* deltaFunction::deltaFunction(double r, double alpha) { // Input: Values for r and alpha to take on. // Output: Constructor. Initializes the delta function so that parameter // ranges are set to [r, r] and [alpha, alpha], and neither parameter // is being fit. _numParams = 2; _currentParams = new double[2]; _currentParams[R] = r; _currentParams[ALPHA] = alpha; _maximumParams = new double[2]; _maximumParams[R] = r; _maximumParams[ALPHA] = alpha; _minimumParams = new double[2]; _minimumParams[R] = r; _minimumParams[ALPHA] = alpha; _fitting = new bool[2]; _fitting[R] = false; _fitting[ALPHA] = false; _inMicronsOfRadius = false; } deltaFunction::deltaFunction() { // Input: None. // Output: Constructor. Initializes the delta function to have parameter // ranges r in [0, 0] and alpha in [1, 1], and neither parameter // is being fit. Useful for setting up an array, because this is // the default constructor. _numParams = 2; _currentParams = new double[2]; _currentParams[R] = 0; _currentParams[ALPHA] = 1; _maximumParams = new double[2]; _maximumParams[R] = 0; _maximumParams[ALPHA] = 1; _minimumParams = new double[2]; _minimumParams[R] = 0; _minimumParams[ALPHA] = 1; _fitting = new bool[2]; _fitting[R] = false; _fitting[ALPHA] = false; _inMicronsOfRadius = false; } void deltaFunction::readFromFile(ifstream& file) { // Input: An ifstream (input file stream), positioned so that the next item // to be read is the name of one of the delta function parameters. // Output: Reads two lines from the file, determines which corresponds to // which parameter, if possible, and stores the values. If this // process fails, an error is printed and execution ends. char* buffer = new char[20]; // For storing item names char* garbage = new char[20]; // For discarding = signs int currentParam = -1; bool doneR = false; // True once r is read bool doneAlpha = false; // True once alpha is read for(int i = 0; i < 2; i++) { file >> buffer; // Read item name and try to if(strncmp(buffer, "r", 1) == 0) { // determine what it is. currentParam = R; // Keep track of the index of doneR = true; // the item to be read. } else if(strncmp(buffer, "alpha", 5) == 0) { currentParam = ALPHA; doneAlpha = true; } else { // Parameter wasn't r or alpha cerr << "Error: Unknown parameter '" << buffer << "'." << endl; cerr << "Valid Delta parameters are r and alpha." << endl; exit(2); } file >> garbage; // Read = sign file >> _minimumParams[currentParam]; // Read minimum and maximum into file >> _maximumParams[currentParam]; // the appropriate arrays. if(_minimumParams[currentParam] > _maximumParams[currentParam]) { cerr << "Error: Delta: Minimum value for " << buffer << " is greater"; cerr << "than maximum value." << endl; exit(2); } if(_minimumParams[currentParam] == _maximumParams[currentParam]) { _fitting[currentParam] = false; // No need to fit if the min } else { // and max parameter values _fitting[currentParam] = true; // are the same. } } if(!(doneR && doneAlpha)) { // Some parameters missing. cerr << "Error: Delta needs parameters r and alpha." << endl; exit(2); } _inMicronsOfRadius = true; // File is in microns of radius delete[] buffer; // Free up memory allocated delete[] garbage; // for reading in strings. } void deltaFunction::toSizeParams(double n, double lambda) { // Input: The index of refraction of the scatterers and the wavelength // of the incident light (nm). // Output: Converts the parameter arrays from microns of radius into // the dimensionless size parameter // // n pi (2r) // x = ----------- . // lambda // // The conversions for each parameter are determined by using the // existing function f(x) to compute f(n pi (2r)/lambda) and then // solve for the values that put f into an equivalent form. In this // case, the proper transformations are given by // // 2 pi n R // alpha |--> alpha R |--> ---------- . // lambda lambda *= 0.001; // Convert lambda from nm to um. if(_inMicronsOfRadius) { // Only perform conversion if double factor = (2.0*PI*n)/lambda; // the function is not already // in size parameters. _minimumParams[R] *= factor; _currentParams[R] *= factor; _maximumParams[R] *= factor; _inMicronsOfRadius = false; // Conversion completed } } void deltaFunction::toMicronsR(double n, double lambda) { // Input: The index of refraction of the scatterers and the wavelength // of the incident light (nm). // Output: Converts the parameter arrays from dimensionless size parameter // // n pi (2r) // x = ----------- // lambda // // to microns of radius. The conversions for each parameter are // determined by inverting the transformations from the previous // conversion function. // // lambda R // alpha |--> alpha R |--> ---------- . // 2 pi n lambda *= 0.001; // Convert lambda to microns if(!(_inMicronsOfRadius)) { // Only perform conversion if double factor = lambda/(2.0*PI*n); // not already in microns. _minimumParams[R] *= factor; _currentParams[R] *= factor; _maximumParams[R] *= factor; _inMicronsOfRadius = true; // Conversion completed } } bool deltaFunction::inRange(double& minX, double& maxX) { // Input: A size parameter range over which to check the numerical // integrals of the functions. // Output: Returns true iff the R values of the delta function fall // within the specified range of size parameters. Since lambda // and n are unknown, it is impossible to convert microns of // radius to size parameters, so performing this check while the // function is expressed in microns generates an error message. if(_inMicronsOfRadius) { cerr << "Error: Unable to perform range check on delta function " << "when parameters are expressed in microns." << endl; exit(2); } bool isOkay = (_minimumParams[R] >= minX) && (_maximumParams[R] <= maxX); if(!isOkay) { // If there is a range error minX = _minimumParams[R]; // reset minX and maxX to be maxX = _maximumParams[R]; // the range of values needed. } return isOkay; // Return false if not okay } double deltaFunction::evaluateWith(double* params, double x) { // Input: A length 2 array [r, alpha] ane a size parameter x. // Output: Returns alpha if x == r (approximately) and 0 otherwise. if(_inMicronsOfRadius) { // Can't compare x and r cerr << "Error: Can't evaluate delta function at size " << x << " when parameters are in microns." << endl; exit(2); } if(fabs(x-params[R]) < 0.000001) { // Exact equality is too much return params[ALPHA]; // to ask, due to roundoff } // errors, so allow some else { // tolerance. return 0; } } double deltaFunction::evaluate(double x) { // Input: A size parameter x. // Output: Returns alpha if x == r (approximately) and 0 otherwise. return evaluateWith(_currentParams, x); } void deltaFunction::output(ostream& stream) { // Input: A deltaFunction to be written to the output stream. // Output: Writes the deltaFunction to that stream in the form // // Delta: r = 1.00 alpha = 3.5e-007 stream.setf(ios::fixed); // Set for fixed-point mode stream << "Delta: r = " << setw(5) // Set field width to 5 and number << setprecision(2) // of decimal places to 2, then << _currentParams[R]; // output r. stream.unsetf(ios::fixed); // Remove fixed mode specifier stream.setf(ios::floatfield); // and set to a float with stream << " alpha = " << setw(8) // exponent and width 8, then << _currentParams[ALPHA]; // output alpha. stream.unsetf(ios::floatfield); // Remove floatfield specifier. } // END class deltaFunction ***************************************************** // class gaussianFunction ****************************************************** gaussianFunction::gaussianFunction() { // Input: None. // Output: Default constructor. Initializes the Gaussian function to have // parameter ranges r in [0, 0], alpha in [1, 1] and sigma in [1, 1], // with no parameters being fit. _numParams = 3; _currentParams = new double[3]; _currentParams[R] = 0; _currentParams[ALPHA] = 1; _currentParams[SIGMA] = 1; _maximumParams = new double[3]; _maximumParams[R] = 0; _maximumParams[ALPHA] = 1; _maximumParams[SIGMA] = 1; _minimumParams = new double[3]; _minimumParams[R] = 0; _minimumParams[ALPHA] = 1; _minimumParams[SIGMA] = 1; _fitting = new bool[3]; _fitting[R] = false; _fitting[ALPHA] = false; _fitting[SIGMA] = false; _inMicronsOfRadius = false; } gaussianFunction::gaussianFunction(double r, double alpha, double sigma) { // Input: Values for r, alpha, and sigma. // Output: Constructor. Sets to r in [r, r] alpha in [alpha, alpha] and // sigma in [sigma, sigma] with no fitting. _numParams = 3; _currentParams = new double[3]; _currentParams[R] = r; _currentParams[ALPHA] = alpha; _currentParams[SIGMA] = sigma; _maximumParams = new double[3]; _maximumParams[R] = r; _maximumParams[ALPHA] = alpha; _maximumParams[SIGMA] = sigma; _minimumParams = new double[3]; _minimumParams[R] = r; _minimumParams[ALPHA] = alpha; _minimumParams[SIGMA] = sigma; _fitting = new bool[3]; _fitting[R] = false; _fitting[ALPHA] = false; _fitting[SIGMA] = false; _inMicronsOfRadius = false; } void gaussianFunction::readFromFile(ifstream& file) { // Input: An ifstream, positioned so that the next item to be read is the // name of one of the Gaussian function parameters. // Output: Reads three lines of the file and tries to sort out parameters. // If successful, it stores the parameters into this Gaussian. // Otherwise, an error is printed and the program exits. char* buffer = new char[20]; // For reading parameter names char* garbage = new char[20]; // For discarding the = sign int currentParam = -1; bool doneR = false; // Flags to be set to true when bool doneAlpha = false; // the corresponding parameter bool doneSigma = false; // has been found in the file. for(int i = 0; i < 3; i++) { // Read three lines and determine file >> buffer; // which parameter is on each. if(strncmp(buffer, "r", 1) == 0) { currentParam = R; doneR = true; } else if(strncmp(buffer, "alpha", 5) == 0) { currentParam = ALPHA; doneAlpha = true; } else if(strncmp(buffer, "sigma", 5) == 0) { currentParam = SIGMA; doneSigma = true; } else { // Not r, alpha, or sigma cerr << "Error: Unknown parameter '" << buffer << "'." << endl; cerr << "Valid Gaussian parameters are r, alpha, and sigma." << endl; exit(2); } file >> garbage; // Read = sign file >> _minimumParams[currentParam]; // Store parameter range in the file >> _maximumParams[currentParam]; // appropriate arrays. if(_minimumParams[currentParam] > _maximumParams[currentParam]) { cerr << "Error: Gaussian: Minimum value for " << buffer << " is greater"; cerr << "than maximum value." << endl; exit(2); } if(_minimumParams[currentParam] == _maximumParams[currentParam]) { _fitting[currentParam] = false; // No need to fit } else { _fitting[currentParam] = true; } } if(!(doneR && doneAlpha && doneSigma)) { // Missing some parameters cerr << "Error: Gaussian needs parameters r, alpha, and sigma." << endl; exit(2); } _inMicronsOfRadius = true; // Files are in microns delete[] buffer; // Free up memory for strings delete[] garbage; } void gaussianFunction::toSizeParams(double n, double lambda) { // Input: The index of refraction of the scatterers and the wavelength // of the incident light (nm). // Output: Converts the parameter arrays from microns of radius into // the dimensionless size parameter // // n pi (2r) // x = ----------- . // lambda // // The conversions for each parameter are determined by using the // existing function f(x) to compute f(n pi (2r)/lambda) and then // solve for the values that put f into an equivalent form. In this // case, the proper transformations are given by // // 2 pi n alpha 2 pi n R 2 pi n sigma // alpha |--> -------------- R |--> ---------- sigma |--> -------------- // lambda lambda lambda lambda *= 0.001; // Convert lambda to microns if(_inMicronsOfRadius) { // Perform conversion to size double factor = (2.0*PI*n)/lambda; // parameters only if the // function is currently in _minimumParams[ALPHA] *= factor; // microns of radius. _currentParams[ALPHA] *= factor; _maximumParams[ALPHA] *= factor; _minimumParams[R] *= factor; _currentParams[R] *= factor; _maximumParams[R] *= factor; _minimumParams[SIGMA] *= factor; _currentParams[SIGMA] *= factor; _maximumParams[SIGMA] *= factor; _inMicronsOfRadius = false; // Conversion completed } } void gaussianFunction::toMicronsR(double n, double lambda) { // Input: The index of refraction of the scatterers and the wavelength // of the incident light (nm). // Output: Converts the parameter arrays from dimensionless size parameter // // n pi (2r) // x = ----------- // lambda // // to microns of radius. The conversions for each parameter are // determined by inverting the transformations from the opposite // conversion function. // // lambda alpha lambda R lambda sigma // alpha |--> -------------- R |--> ---------- sigma |--> -------------- // 2 pi n 2 pi n 2 pi n lambda *= 0.001; // Convert lambda to microns if(!(_inMicronsOfRadius)) { // Switch to microns only if double factor = lambda/(2.0*PI*n); // not already in microns. _minimumParams[ALPHA] *= factor; _currentParams[ALPHA] *= factor; _maximumParams[ALPHA] *= factor; _minimumParams[R] *= factor; _currentParams[R] *= factor; _maximumParams[R] *= factor; _minimumParams[SIGMA] *= factor; _currentParams[SIGMA] *= factor; _maximumParams[SIGMA] *= factor; _inMicronsOfRadius = true; // Conversion completed } } bool gaussianFunction::inRange(double& minX, double& maxX) { // Input: Size parameter range. // Output: Returns true iff the mean +/- 1.5 sigma of the Gaussian is // within the specified range. This corresponds to about 95% of the // area of the Gaussian falling within the specified range. if(_inMicronsOfRadius) { cerr << "Error: Unable to perform range check on gaussian function " << "when parameters are expressed in microns." << endl; exit(2); } bool isOkay = ((_minimumParams[R]-1.5*_maximumParams[SIGMA]) >= minX) && ((_maximumParams[R]+1.5*_maximumParams[SIGMA]) <= maxX); if(!isOkay) { // If there is a range problem, minX = (_minimumParams[R] // set minX and maxX to hold -1.5*_maximumParams[SIGMA]); // the range of size parameters maxX = (_maximumParams[R] // that would be needed. +1.5*_maximumParams[SIGMA]); } return isOkay; } double gaussianFunction::evaluateWith(double* params, double x) { // Input: An array of length 3 containing [r, alpha, sigma] and a size // parameter x at which to evaluate the Gaussian. // Output: Returns the value of the Gaussian with mean r, area alpha, and // standard deviation sigma, evaluated at x: // // (x - r)² // alpha - ---------- // ----------------- e 2 sigma² // sigma sqrt(2pi) return (params[ALPHA]/(params[SIGMA]*sqrt(2.0*PI)))* exp(-pow(x-params[R],2.0)/(2.0*pow(params[SIGMA],2.0))); } double gaussianFunction::evaluate(double x) { // Input: A size parameter x. // Output: Returns the value of the Gaussian with mean r, area alpha, and // standard deviation sigma, given by _currentParams, evaluated at x: // // (x - r)² // alpha - ---------- // ----------------- e 2 sigma² // sigma sqrt(2pi) return evaluateWith(_currentParams, x); } void gaussianFunction::output(ostream& stream) { // Input: A gaussianFunction to be written to the output stream. // Output: Writes the gaussianFunction to that stream in the form // // Gaussian: r = 1.34 alpha = 3.5e-007 sigma = 0.49 stream.setf(ios::fixed); // Use fixed-decimal output stream << "Gaussian: r = " << setw(5) // Set width to 5 and precision << setprecision(2) // to 2, then output r and << _currentParams[R]; // clear flags. stream.unsetf(ios::fixed); stream.setf(ios::floatfield); // Set to floating point output stream << " alpha = " << setw(8) // with an exponent, and field << _currentParams[ALPHA]; // width to 8, then output alpha stream.unsetf(ios::floatfield); // and unset floatfield. stream.setf(ios::fixed); // Fixed output of sigma stream << " sigma = " << setw(6) << _currentParams[SIGMA]; stream.unsetf(ios::fixed); } // END class gaussianFunction ************************************************** // class logonormalFunction **************************************************** logonormalFunction::logonormalFunction() { // Input: None. // Output: Default constructor. Initializes the logonormal function to have // parameter ranges r in [0, 0], alpha in [1, 1] and sigma in [1, 1], // with no parameters being fit. _numParams = 3; _currentParams = new double[3]; _currentParams[R] = 0; _currentParams[ALPHA] = 1; _currentParams[SIGMA] = 1; _maximumParams = new double[3]; _maximumParams[R] = 0; _maximumParams[ALPHA] = 1; _maximumParams[SIGMA] = 1; _minimumParams = new double[3]; _minimumParams[R] = 0; _minimumParams[ALPHA] = 1; _minimumParams[SIGMA] = 1; _fitting = new bool[3]; _fitting[R] = false; _fitting[ALPHA] = false; _fitting[SIGMA] = false; _inMicronsOfRadius = false; } logonormalFunction::logonormalFunction(double r, double alpha, double sigma) { // Input: Values for r, alpha, and sigma. // Output: Constructor. Sets to r in [r, r] alpha in [alpha, alpha] and // sigma in [sigma, sigma] with no fitting. _numParams = 3; _numParams = 3; _currentParams = new double[3]; _currentParams[R] = r; _currentParams[ALPHA] = alpha; _currentParams[SIGMA] = sigma; _maximumParams = new double[3]; _maximumParams[R] = r; _maximumParams[ALPHA] = alpha; _maximumParams[SIGMA] = sigma; _minimumParams = new double[3]; _minimumParams[R] = r; _minimumParams[ALPHA] = alpha; _minimumParams[SIGMA] = sigma; _fitting = new bool[3]; _fitting[R] = false; _fitting[ALPHA] = false; _fitting[SIGMA] = false; _inMicronsOfRadius = false; } void logonormalFunction::readFromFile(ifstream& file) { // Input: An ifstream, positioned so that the next item to be read is the // name of one of the logonormal function parameters. // Output: Reads three lines of the file and tries to sort out parameters. // If successful, it stores the parameters into this logonormal. // Otherwise, an error is printed and the program exits. char* buffer = new char[20]; // For reading parameter names char* garbage = new char[20]; // For discarding the = sign int currentParam = -1; bool doneR = false; // Flags to be set to true when bool doneAlpha = false; // the corresponding parameter bool doneSigma = false; // has been found in the file. for(int i = 0; i < 3; i++) { // Read lines and determine which file >> buffer; // corresponds to which param. if(strncmp(buffer, "r", 1) == 0) { currentParam = R; doneR = true; } else if(strncmp(buffer, "alpha", 5) == 0) { currentParam = ALPHA; doneAlpha = true; } else if(strncmp(buffer, "sigma", 5) == 0) { currentParam = SIGMA; doneSigma = true; } else { cerr << "Error: Unknown parameter '" << buffer << "'." << endl; cerr << "Valid logonormal parameters are r, alpha, and sigma." << endl; exit(2); } file >> garbage; // Read = sign file >> _minimumParams[currentParam]; // Read parameter ranges and store file >> _maximumParams[currentParam]; // into the appropriate places. if(_minimumParams[currentParam] > _maximumParams[currentParam]) { cerr << "Error: Logonormal: Minimum value for " << buffer << " is greater"; cerr << "than maximum value." << endl; exit(2); } if(_minimumParams[currentParam] == _maximumParams[currentParam]) { _fitting[currentParam] = false; // No need to fit } else { _fitting[currentParam] = true; } } if(!(doneR && doneAlpha && doneSigma)) { cerr << "Error: Logonormal needs parameters r, alpha, and sigma." << endl; exit(2); } _inMicronsOfRadius = true; // File is in microns delete[] buffer; // Free up memory for strings delete[] garbage; } void logonormalFunction::toSizeParams(double n, double lambda) { // Input: The index of refraction of the scatterers and the wavelength // of the incident light (nm). // Output: Converts the parameter arrays from microns of radius into // the dimensionless size parameter // // n pi (2r) // x = ----------- . // lambda // // The conversions for each parameter are determined by using the // existing function f(x) to compute f(n pi (2r)/lambda) and then // solve for the values that put f into an equivalent form. In this // case, the proper transformations are given by // // 2 pi n alpha 2 pi n R // alpha |--> -------------- R |--> ---------- sigma |--> sigma // lambda lambda lambda *= 0.001; // Convert lambda to microns if(_inMicronsOfRadius) { // Only switch to size params if double factor = (2.0*PI*n)/lambda; // currently in microns. _minimumParams[ALPHA] *= factor; _currentParams[ALPHA] *= factor; _maximumParams[ALPHA] *= factor; _minimumParams[R] *= factor; _currentParams[R] *= factor; _maximumParams[R] *= factor; _inMicronsOfRadius = false; // Conversion completed } } void logonormalFunction::toMicronsR(double n, double lambda) { // Input: The index of refraction of the scatterers and the wavelength // of the incident light (nm). // Output: Converts the parameter arrays from dimensionless size parameter // // n pi (2r) // x = ----------- // lambda // // to microns of radius. The conversions for each parameter are // determined by inverting the transformations from the opposite // conversion function. // // lambda alpha lambda R // alpha |--> -------------- R |--> ---------- sigma |--> sigma // 2 pi n 2 pi n lambda *= 0.001; // Convert lambda to microns if(!(_inMicronsOfRadius)) { // Only switch to microns if double factor = lambda/(2.0*PI*n); // not already there. _minimumParams[ALPHA] *= factor; _currentParams[ALPHA] *= factor; _maximumParams[ALPHA] *= factor; _minimumParams[R] *= factor; _currentParams[R] *= factor; _maximumParams[R] *= factor; _inMicronsOfRadius = true; // Conversion completed } } double min(double x, double y) { // Input: Doubles x and y. // Output: Their minimum. if(x < y) { return x; } else { return y; } } double max(double x, double y) { // Input: Doubles x and y. // Output: Their maximum. if(x > y) { return x; } else { return y; } } bool logonormalFunction::inRange(double& minX, double& maxX) { // Input: A range of size parameters to be tested. // Output: True iff this logonormal function has 95% of its area in the // specified range. The bounds on 95% of the area are given by // // 1/2(2 M ± 2.7718 S + S²) // e , // // where M is the natural log of the R parameter, and S is the // SIGMA parameter. Uses the min and max functions above. Note // that, since S^2 - 2.7718 S has a local minimum, exp of this // quantity does as well, and this must be dealt with carefully. // Finally, since lambda and n are unknown here, this check can // only be performed while the function is in size parameters. if(_inMicronsOfRadius) { cerr << "Error: Unable to perform range check on gaussian function " << "when parameters are expressed in microns." << endl; exit(2); } double lowerEnd = min(exp(1.0/2.0*(2.0*log(_minimumParams[R]) + pow(_minimumParams[SIGMA],2.0) - 2.7718*_minimumParams[SIGMA])) , exp(1.0/2.0*(2.0*log(_minimumParams[R]) + pow(_maximumParams[SIGMA],2.0) - 2.7718*_maximumParams[SIGMA])) ); if((_minimumParams[SIGMA] < 1.386) && (_maximumParams[SIGMA] > 1.386)) { lowerEnd = min(lowerEnd, exp(1.0/2.0*(2.0*log(_minimumParams[R]) + pow(1.386, 2.0)- 2.7718*1.386)) ); } double upperEnd = exp(1.0/2.0*(2.0*log(_maximumParams[R]) + pow(_maximumParams[SIGMA],2.0) + 2.7718*_maximumParams[SIGMA])); if((lowerEnd >= minX) && (upperEnd <= maxX)) { return true; } else { minX = lowerEnd; // Reset minX and maxX to hold maxX = upperEnd; // the range required. return false; } } double logonormalFunction::evaluateWith(double* params, double x) { // Input: An array of length 3 containing [r, alpha, sigma] and a size // parameter x. // Output: Returns the value of the logonormal with mean r, area alpha, and // standard deviation sigma, evaluated at x: // // -ln(r) - sigma²/4 (ln(x) - ln(r))² // alpha e - ------------------ // -------------------------- e sigma² . // sigma sqrt(pi) // return (params[ALPHA] * exp(-log(params[R]) - pow(params[SIGMA],2.0)/4.0)/ (params[SIGMA] * sqrt(PI)) * exp(-pow(log(x) - log(params[R]), 2.0)/pow(params[SIGMA],2.0)) ); } double logonormalFunction::evaluate(double x) { // Input: A size parameter x. // Output: Returns the value of the logonormal with mean r, area alpha, and // standard deviation sigma given by _currentParams, evaluated at x. return evaluateWith(_currentParams, x); } void logonormalFunction::output(ostream& stream) { // Input: A logonormalFunction to be written to the output stream. // Output: Writes the logonormalFunction to that stream in the form // // Logonormal: r = 1.34 alpha = 3.5e-007 sigma = 0.49 stream.setf(ios::fixed); // Fixed decimal for r stream << "Logonormal: r = " << setw(5) << setprecision(2) << _currentParams[R]; stream.unsetf(ios::fixed); stream.setf(ios::floatfield); // Exponent for alpha stream << " alpha = " << setw(8) << _currentParams[ALPHA]; stream.unsetf(ios::floatfield); stream.setf(ios::fixed); // Fixed decimal sigma stream << " sigma = " << setw(6) << _currentParams[SIGMA]; stream.unsetf(ios::fixed); } // END class logonormalFunction ************************************************