// 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 <iostream.h>                         // Streams
#include <iomanip.h>                          // Output formatting
#include <fstream.h>                          // File streams
#include <string.h>                           // String manipulation routines
#include <stdlib.h>                           // For exit() and NULL
#include <math.h>                             // 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 ************************************************
