// Filename: fitdist.cpp
   #define   VERSION 2.2
// Author :  Mike Schubmehl
// Date(s):  July 2002
// Purpose:  Fits data from Mie scattering experiments to a distribution
//           of sizes. The program takes as input a data file that contains
//           the polarization, incident light wavelength, and relative
//           index of refraction for experimental data, as well as
//           intensity data as a function of angle, with uncertainty.
//           For example, a data file might look like:
//
//                  // Some description of the data here.
//                  // This is from a 2.4 MHz board.
//
//                  polarization	 = unpolarized
//                  wavelength  	 = 488
//                  ambientIndex	 = 1.00
//                  scattererIndex = 1.33
//                  relativeIndex  = 1.33
//
//                  2.    0.8934722517541993       0.04462386510041677
//                  3.    0.7374139804418689       0.047095193908696334
//                  4.    0.49943374858437145      0.035329931847326135
//                  5.    0.3472889498970487       0.016730686710008463
//                  6.    0.23342999171499584      0.011860555740763163
//
//           The program also requires theoretical tables of intensity as a
//           function of angle at varous size parameters, in the format used
//           by the mietable program. That is,
//
//        *   1.000
//          0.00      0.52613690E-01      0.52613690E-01      0.52613690E-01
//          1.00      0.52610300E-01      0.52595302E-01      0.52602801E-01
//          2.00      0.52600145E-01      0.52540164E-01      0.52570157E-01
//            .             .                   .                   .
//            .             .                   .                   .
//            .             .                   .                   .
//        *   2.000 (size parameter)
//          0.00      0.39379461E+01      0.39379461E+01      0.39379461E+01
//          1.00      0.39368415E+01      0.39360471E+01      0.39364443E+01
//            .             .                   .                   .
//            .             .                   .                   .
//            .             .                   .                   .
//          (angle)   (perpendicular)       (parallel)        (unpolarized)
//
//           If a table is not explicitly provided, the program will use its
//           own table (fitdist.dat). In addition to data and theoretical Mie
//           tables, the program requires a fit file specifying the type of
//           fit and the ranges for all relevant parameters. This file
//           has the format
//
//                    [Gaussian]
//                    alpha = 0   1
//                    sigma = 0.5 3
//                    r     = 1   10
//
//                    [Delta]
//                    alpha = 0 1
//                    r     = 3 5
//
//                    [Delta]
//                    alpha = 0 1
//                    r     = 5 7
//
//                    [Logonormal]
//                    alpha = 0   1
//                    sigma = 0.5 3
//                    r     = 1   10
//
//                    .
//                    .
//                    .
//           
//           This particular file would fit a distribution consisting of one
//           Gaussian peak, two discrete delta peaks, and one logonormal
//           function, using the specified parameter ranges. Formulas for the
//           various functions are given in the accompanying documentation,
//           and in the evaluate() function comments in fitfunc.cpp.
//
//           In practice, no more than one or two Gaussians, or two or three
//           delta functions should actually be used, as there are too many
//           degrees of freedom to perform a good fit.
//
//           As it proceeds with the fitting process, the program outputs its
//           current status to cerr, on a single line. After obtaining a fit,
//           a summary of the results is displayed, and the user is asked
//           whether the output should be written to file, and then whether
//           a modified fit is desired. This method saves time, because the
//           data and mie tables do not need to be read in from file again. If
//           the data are written to file, the output is
//
//              Distributed Fit with fitdist #VERSION
//
//              Data:
//              (entire contents of data file)
//
//              Parameters:
//              (entire contents of fit file)
//
//              Mie Tables:
//              Size parameters range from #### to #### in steps of ####.
//
//              Summary of Results:
//                Chi^2:     ####
//                FitFunc1:  r = ####  alpha = #### sigma = ####
//                FitFunc2:  r = ####  alpha = ####
//
//              Size Distribution:
//                radius (um)   relative abundance
//              ####            ####
//              ####            ####
//                   .                  .
//                   .                  .
//                   .                  .
//              (relative abundance as a function of radius in microns)
//
//              Theoretical Intensity Profile:
//                Angle         Intensity
//              ##              #####
//              ##              #####
//                  .              .
//                  .              .
//                  .              .
//              (intensity as a function of angle for above distribution)
//
//           If another fit is requested, the user is given a chance to read
//           a new fit file. The program terminates if the user does not
//           wish to perform another fit.
//
// Notes:    Requires stdlib.h, iostream.h, fstream.h, string.h.

#include <stdlib.h>                           // For exit() and NULL
#include <iostream.h>                         // Input and output streams
#include <iomanip.h>                          // For output formatting
#include <fstream.h>                          // File streams
#include <string.h>                           // String manipulation
#include <ctype.h>                            // For isdigit() etc
#include <math.h>                             // Mathematical functions.
#include "fitfunc.h"                          // For fittingFunction class

#define EXTENDED true                         // Flags for detail level of the
#define SIMPLE   false                        //   usage message.

#define SIMPLE_FIT 0                          // Fitting methods
#define LEVMAR_FIT 1

#define PERPENDICULAR 1                       // Perpendicular to scatter plane
#define PARALLEL      2                       // Parallel to scatter plane
#define UNPOLARIZED   3                       // Unpolarized

#define EPSILON 0.0001                        // A "small" number
#define N_AIR 1.00                            // The index of refraction of air
#ifndef PI                                    
#define PI 3.1415926535                       // If no pi defined, use this
#endif

struct experimentalData {                     // Struct to hold data
  int    polarization;                        // Polarization of incident light
  double wavelength;                          // Wavelength of incident light
  double relativeIndex;                       // (scatter index)/(ambient index),
                                              //   so water in air is 1.33.
  double scattererIndex;                      // Index of sphere
  double ambientIndex;                        // Index of surrounding medium

  double* angles;                             // Array of angles present in data
  double* intensities;                        // Corresponding intensity and
  double* uncertainties;                      //   uncertainty information.

  int numPoints;                              // Number of data points.
} data;                                       // Call the actual variable data

struct theoreticalMieTables {                 // Struct to hold tables
  double* angles;                             // Angles contained in the table
  int     numAngles;                          // Number of angles present
  double* sizeParameters;                     // Size parameters in table
  int     numSizeParameters;                  // Number of size parameters

  double  mieTablesMaxAngle;                  // Record table properties for
  double  mieTablesAngleStep;                 //   output to final results.

  double** intensities;                       // Scattered intensity, indexed
} theory;                                     //   by [sizeparam#][angle#].

struct fittingSetup {                         // Struct to hold fit parameters
  char* dataFilename;                         // Filename for data file
  char* fitFilename;                          // Filename for fit setup file
  char* tablesFilename;                       // Filename for Mie tables

  char* outputFilename;                       // Filename for output

  int fitMethod;                              // SIMPLE_FIT or LEVMAR_FIT
  int pointsPerParameter;                     // Points per parameter if
                                              //   doing a simple fit.

  double bestChiSquared;                      // Best chi^2 obtained in fit
  double overallAlpha;                        // Overall scaling factor

  int fittingFunctionCount;                   // Number of fitting functions
  fittingFunction** fittingFunctions;         // Array of pointers to fitting
} setup;                                      //   functions. Done with pointers
                                              //   so non-default constructor
                                              //   can be called with new.

// Table of Contents ***********************************************************
void printUsage(bool extended);               // Print usage message, with
                                              //   details if extended is true
void processArguments(int argc, char* argv[]);// Read command line arguments and
                                              //   transfer data into setup.
void openFile(char* filename,                 // Open specified file and make
              ifstream& file);                //   sure it can be read.
void readData(char* dataFilename);            // Read data from specified file
                                              //   and check for validity.
void readNecessaryMieTables(                  // Read from file only those parts
              char* tablesFilename);          //   of the Mie tables needed to
                                              //   compute chi^2 values. Need to
                                              //   examine data for this.
void readFit(char* fitFilename);              // Read fit from file, storing
                                              //   fit info in setup.
void performFit();                            // Perform the actual of theory to
                                              //   data, as specified by the
                                              //   fitting configuration.
void performLevenbergMarquardtFit();          // Fit using the Levenberg
                                              //   Marquardt method of steepest
                                              //   descent/inverse Hessian.
void performSimpleFit();                      // Fit by breaking parameter
                                              //   space into a grid and
                                              //   checking all possibilities.
void outputSummary();                         // Write brief summary of fit.
bool fullOutputRequested();                   // Prompt user. True if full fit
                                              //   results are to be output.
char* constructDefaultOutputName();           // Create a sensible default name
                                              //   for the output.
char* getFilename();                          // Prompt user for a filename.
char* getFilename(char* defaultFilename);     // Prompt with default.
void outputFullResults();                     // Write full fitting results
                                              //   to a file.
bool additionalFitRequested();                // Prompt user. True if another
                                              //   fit is to be performed.

double chiSquared();                          // Compute chi^2 per DOF
void computeLevenbergMarquardtCoefficients(double* params, int numParams,
                                           double** alpha, double* beta,
                                           double& chiSquared);
                                              // Compute coefficients needed for
                                              //   fitting by LevMar algorithm.
double evaluateTheory(double* params,         // Takes fitting functions and the
                      int angleNum);          //   array of parameters, stores
                                              //   parameters to the current
                                              //   params of the functions,
                                              //   and evaluates at angle.
// END Table of Contents *******************************************************



// Main Program ****************************************************************
int main(int argc, char* argv[]) {
  if(argc < 2) {                              // If there are fewer than 2
    printUsage(EXTENDED);                     //   arugments, counting the full
    exit(1);                                  //   path to run the program,
  }                                           //   print a usage message.

  processArguments(argc, argv);               // Process command arguments
  readData(setup.dataFilename);               // Read in experimental data
  readNecessaryMieTables(setup.tablesFilename);//  and the needed Mie tables.

  bool doneFitting = false;                   // Flag to indicate completion
  while(!doneFitting) {
    readFit(setup.fitFilename);               // Read from current fit file
    performFit();                             // Do the fit by the selected
    outputSummary();                          //   method and output results.
    if(fullOutputRequested()) {               // Output to file if desired
      setup.outputFilename = constructDefaultOutputName();
      setup.outputFilename = getFilename(setup.outputFilename);
      outputFullResults();
    }
    if(additionalFitRequested()) {            // Shortcut to perform another fit
      setup.fitFilename = getFilename(setup.fitFilename);
    } else {
      doneFitting = true;
    }
  }

  return 0;                                   // Exit without error status
}
// END Main Program ************************************************************











// Mathematical Functions ******************************************************
double chiSquared() {
  // Input:  The fitting functions, with their current parameter values,
  //         Mie tables containing theoretical scattered intensity, and
  //         experimental data points with uncertainty.
  // Output: The value of chi^2, computed between the theoretical size
  //         distribution's Mie amplitudes and the data, and then divided
  //         by the number of degrees of freedom, given by
  //
  //                DOF = (# data points) - (# parameters in fit).


  int i = 0;
  double sumSquaredResiduals = 0;
  double degreesOfFreedom = data.numPoints;   // Compute degrees of freedom
  for(i = 0; i < setup.fittingFunctionCount; i++) {
    degreesOfFreedom -= setup.fittingFunctions[i]->numParams();
  }
      
  double* theoryDistribution = new double[theory.numAngles];

  for(i = 0; i < theory.numAngles; i++) {     // Declare and initialize Mie
    theoryDistribution[i] = 0;                //   theory distribution
  }                                           //   (intensity vs. angle).

  for(i = 0; i < theory.numSizeParameters; i++) {
    double currentCoefficient = 0;            // Multiplier at this size param
    int j = 0;
    for(j = 0; j < setup.fittingFunctionCount; j++) {
      currentCoefficient += setup.fittingFunctions[j]->
                               evaluate(theory.sizeParameters[i]);
    }                                         // Sum up all functions' values

    for(j = 0; j < theory.numAngles; j++) {   // Use this multiplier with tables
      theoryDistribution[j] += currentCoefficient*theory.intensities[i][j];
    }
  }

  for(i = 0; i < theory.numAngles; i++) {     // Compute sumSquaredResiduals
    sumSquaredResiduals += pow((data.intensities[i]-theoryDistribution[i])/
                                data.uncertainties[i],2.0);
  }

  delete[] theoryDistribution;                // Free memory for theory

  return sumSquaredResiduals /                // Divide by degrees of freedom
         degreesOfFreedom;
}


void performFit() {
  // Input:  None.
  // Output: Performs the fit using selected algorithm.

  switch(setup.fitMethod) {
  case SIMPLE_FIT:
    performSimpleFit();
    break;
  case LEVMAR_FIT:
    performLevenbergMarquardtFit();
    break;
  default:
    cerr << "Error: Unknown fitting method selected." << endl;
    exit(1);
  }
}

// Levenberg-Marquardt Fitting Routine -----------------------------------------
inline void swap(double& a, double& b) {
  // Input:  Doubles a and b, passed by reference.
  // Output: Swaps the value of a and b.

  double temp = a;
  a = b;
  b = temp;
}

void solveGaussJordan(double** A, double* b, int n) {
  // Input:  An n × n matrix a and an n-item vector b.
  // Output: Solves the system
  //
  //                                 A x = b
  //
  //         by Gaussian elimination and returns the solution vector x,
  //         storing it in the vector b.
  // Notes:  Adapted from Numerical Recipes in C, second edition.

  double largest;                             // Largest element found so far
  int    pivotRow = 0;                        // Row and column on which to
  int    pivotCol = 0;                        //   pivot, and the number of
  int*   pivot = new int[n];                  //   times the row has pivoted.
                                              
  int i = 0;
  for(i = 0; i < n; i++) {
    pivot[i] = 0;
  }

  int j = 0;
  int l = 0;
  for(i = 0; i < n; i++) {
    largest = 0;
    for(j = 0; j < n; j++) {
      if(pivot[j] != 1) {
        for(int k = 0; k < n; k++) {
          if(pivot[k] == 0) {
            if(fabs(A[j][k]) >= largest) {    // Find largest element, and
              largest = fabs(A[j][k]);        //   mark new location for pivot.
              pivotRow = j;
              pivotCol = k;
            }
          } else if(pivot[k] > 1) {           // Already pivoted here. Error.
            cerr << "Error: Singular matrix A in the system Ax = b." << endl;
            exit(2);
          }
        }
      }
    }
    (pivot[pivotCol])++;                      

    if(pivotRow != pivotCol) {                // Have pivot element. Need to
      for(l = 0; l < n; l++) {                //   interchange rows and columns
        swap(A[pivotRow][l], A[pivotCol][l]); //   to put it on the diagonal.
      }
      swap(b[pivotRow], b[pivotCol]);
    }

    if(A[pivotCol][pivotCol] == 0.0) {        // Numerical recipes claims that
      return;                                 //   this is also an error. It
    }                                         //   looks to me like you're just
                                              //   done with the reduction.

    double pivotInverse = 1.0/A[pivotCol][pivotCol];
    for(l = 0; l < n; l++) {              // Divide through by the inverse
      A[pivotCol][l] *= pivotInverse;         //   of the pivot element.
    }
    b[pivotCol] *= pivotInverse;              // Divide vector entries as well

    for(l = 0; l < n; l++) {                  // Now subtract off the row
      if(l != pivotCol) {                     //   to reduce all rows below it.
        double temp = A[l][pivotCol];
        for(int m = 0; m < n; m++) {
          A[l][m] -= A[pivotCol][m] * temp;
        }
        b[l] -= b[pivotCol] * temp;
      }
    }
  }

  delete[] pivot;                             // Free memory for pivot array
}

double evaluateTheory(double* params, int angleNum) {
  // Input:  Fitting setup, Mie tables, an array of parameters, and an angle.
  // Output: Substitutes params for the current params of the fitting
  //         functions, then evaluates the theory at the given angle.

  fittingFunction* f;                         // Shorthand for fitting functions

  setup.overallAlpha = 1;                     // Set overall scaling factor

  int i = 0;
  int j = 0;
  int k = 0;
  for(i = 0; i < setup.fittingFunctionCount; i++) {
    f = setup.fittingFunctions[i];            // Distribute parameters in the
    for(j = 0; j < f->numParams(); j++) {     //   array back into parameters
      if(f->fitting(j)) {                     //   of individual fit functions.
        f->currentParams(j) = params[k];
        k++;
      }
    }
  }

  double result = 0;                          // Theory value at this angle
  for(i = 0; i < theory.numSizeParameters; i++) {
    double currentCoefficient = 0;            // Multiplier for this size param
    for(j = 0; j < setup.fittingFunctionCount; j++) {
      currentCoefficient += setup.fittingFunctions[j]->
                            evaluate(theory.sizeParameters[i]);
    }
    result += setup.overallAlpha*currentCoefficient*
              theory.intensities[i][angleNum];
  }

  return result;
}
                   
void computeLevenbergMarquardtCoefficients(double* params, int numParams,
                                           double** alpha, double* beta,
                                           double& chiSquared) {
  // Input:  The fitting setup, to be used in evaluating the fitting functions,
  //         the array of parameters to be fit, the size of this array,
  //         the experimental data to be fit to, and the alpha, beta, and
  //         chiSquared coefficients (to be updated).
  // Output: Computes alpha, beta, and chiSquared, as defined and discussed in
  //         Numerical Recipes in C, Second Ed. (c) 1992.

  chiSquared = 0;
  int i = 0;
  int j = 0;
  for(i = 0; i < numParams; i++) {
    for(j = 0; j <= i; j++) {
      alpha[i][j] = 0;
    }
    beta[i] = 0;
  }

  double* partialDs = new double[numParams];
  double* tweakedParams = new double[numParams];

  for(i = 0; i < data.numPoints; i++) {
    for(j = 0; j < numParams; j++) {
      for(int k = 0; k < numParams; k++) {    // Numerical approximation to
        if(k == j) {                          //   partial derivative.
          tweakedParams[k] = params[k]*(1.0+EPSILON);
        } else {
          tweakedParams[k] = params[k];
        }
      }
      partialDs[j] = (evaluateTheory(tweakedParams, i) - 
                      evaluateTheory(params, i))/
                     (EPSILON*params[j]);
    }

    double sigmaSquared = (data.uncertainties[i] * data.uncertainties[i]);
    double residual = evaluateTheory(params, i) - data.intensities[i];

    for(j = 0; j < numParams; j++) {
      double weight = partialDs[j]/sigmaSquared;
      for(int k = 0; k <= j; k++) {
        alpha[j][k] += weight * partialDs[k];
      }
      beta[j] += residual * weight;
    }
    chiSquared += residual * residual / sigmaSquared;
  }

  chiSquared = chiSquared /                   // Compute chi squared per DOF
               (data.numPoints - numParams);

  delete[] partialDs;                         // Free memory
  delete[] tweakedParams;

  for(i = 1; i < numParams; i++) {            // Fill in symmetric half of
    for(j = 0; j < i; j++) {                  //   the alpha matrix.
      alpha[j][i] = alpha[i][j];
    }
  }
}

void performLevenbergMarquardtFit() {
  // Input:  Fitting functions with parameter ranges, theoretical Mie tables,
  //         and experimental data.
  // Output: Performs the fit and stores the new values of the parameters in
  //         the fitting functions.


  fittingFunction* f;

  //                  STEP 1: ESTABLISHING PARAMETERS

  int parameterCount = 0;
  int i = 0;
  int j = 0;
  for(i = 0; i < setup.fittingFunctionCount; i++) {
    f = setup.fittingFunctions[i];            // Count up parameters being fit
    for(j = 0; j < f->numParams(); j++) {
      if(f->fitting(j)) {
        parameterCount++;
      }
    }
  }

  double*    params = new double[parameterCount];
  double* minParams = new double[parameterCount];
  double* maxParams = new double[parameterCount];

  int k = 0;
  for(i = 0; i < setup.fittingFunctionCount; i++) {
    f = setup.fittingFunctions[i];
    for(j = 0; j < f->numParams(); j++) {
      f->currentParams(j) = (f->maximumParams(j) + f->minimumParams(j))/2.0;
      if(f->fitting(j)) {                     // Store values in the params
        params[k] = f->currentParams(j);      //   arrays for use by the fitting
        minParams[k] = f->minimumParams(j);   //   routine.
        maxParams[k] = f->maximumParams(j);
        k++;
      }
    }
  }

  //                       STEP 2: FITTING SETUP

  double*  testParams = new double[parameterCount];
  double* dTestParams = new double[parameterCount];
  double** covariance = new double*[parameterCount];
  for(i = 0; i < parameterCount; i++) {
    covariance[i] = new double[parameterCount];
    testParams[i] = params[i];
  }
  double  chiSquared    =  0;                 // Current chi^2 and old chi^2 for
  double  oldChiSquared = -1;                 //   comparison. Start at -1 so
                                              //   initially chi^2 > old chi^2.

  double lambda = 0.001;                      // Lambda is the parameter that
                                              //   chooses a point on the
                                              //   continuum from steepest-
                                              //   descent to inverse Hessian
                                              //   fitting methods. This value
                                              //   is a standard initial value.

  //                        STEP 3: ITERATIVE FIT

  bool* hitBoundary = new bool[parameterCount];
  bool  hitAllBoundaries = false;             // Keep track of which boundaries
  for(i = 0; i < parameterCount; i++) {       //   have been hit.
    hitBoundary[i] = false;
  }

  int iterationCount = 0;
  double testChiSquared;
  while((chiSquared > oldChiSquared) ||
        (fabs(1.0 - chiSquared/oldChiSquared) > EPSILON * EPSILON)) {

    computeLevenbergMarquardtCoefficients(params, parameterCount,
                                          covariance, dTestParams, chiSquared);

    for(i = 0; i < parameterCount; i++) {
      covariance[i][i] *= (1.0+lambda);       // Augment diagonal
    }

    solveGaussJordan(covariance, dTestParams, // Solve for the change in test
                     parameterCount);         //   parameters needed.

    hitAllBoundaries = true;                  // Not really true. Needed to
                                              //   && together the array.
    for(i = 0; i < parameterCount; i++) {     // Solution gives the needed
      testParams[i] = params[i] -             //   adjustment to the test
                      dTestParams[i];         //   parameters, so adjust the
      if(testParams[i] < minParams[i]) {      //   test parameters accordingly.
        testParams[i] = minParams[i];
        hitBoundary[i] = true;
      }
      if(testParams[i] > maxParams[i]) {
        testParams[i] = maxParams[i];
        hitBoundary[i] = true;
      }

      hitAllBoundaries = hitAllBoundaries && hitBoundary[i];
    }

    computeLevenbergMarquardtCoefficients(testParams, parameterCount,
                                          covariance, dTestParams,
                                          testChiSquared);

    if(testChiSquared < chiSquared) {         // Computed new chi^2. If it's
      lambda *= 0.1;                          //   better than the old one,
      oldChiSquared = chiSquared;             //   adopt the test parameters,
      chiSquared = testChiSquared;            //   store chi^2, and decrease
      for(i = 0; i < parameterCount; i++) {   //   lambda (move closer to the
        params[i] = testParams[i];            //   inverse-Hessian method).
      }
    } else {                                  // If the adjustment made chi^2
      lambda *= 10.0;                         //   worse, reject test params
    }                                         //   and increase lambda.

    iterationCount++;                         // Keep track of how many
    cerr << "Fitting: ";                      //   iterations have taken place
    for(i = 0; i < 3; i++) {                  //   and use this to generate a
      if(i <= (iterationCount % 4)) {         //   progress indicator of sorts.
        cerr << ".";
      } else {
        cerr << " ";
      }
    }
    cerr << " " << flush;
    cerr << "  Best chi^2:  " << chiSquared;   // This is the real indicator
    cerr << "            ";
    cerr << "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b";
    cerr << "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b" << flush;
  }


  //                       STEP 4: REPORT RESULTS

  cerr << "Fitting: Done  Best chi^2: " << chiSquared << "         " << flush;

  setup.bestChiSquared = chiSquared;          // Store best chi^2 for output

  delete[] hitBoundary;                       // Free memory
  delete[] params;
  delete[] testParams;
  delete[] dTestParams;
  for(i = 0; i < parameterCount; i++) {
    delete[] covariance[i];
  }
  delete[] covariance;
}
// END Levenberg-Marquardt Fitting Routine -------------------------------------


// Simple Chi^2 Exploration Fitting Routine ------------------------------------
void performSimpleFit() {
  // Input:  Fitting functions with parameter ranges, theoretical Mie tables,
  //         and experimental data.
  // Output: Performs the fit and stores the new values of the parameters in
  //         the fitting functions.

  fittingFunction* f;                         // For shorthand

  int parameterCount = 0;                     // Count parameters to be fit
  int i = 0;
  int j = 0;
  for(i = 0; i < setup.fittingFunctionCount; i++) {
    f = setup.fittingFunctions[i];
    f->toSizeParams(data.ambientIndex, data.wavelength);
    for(j = 0; j < f->numParams(); j++) {
      if(f->fitting(j)) {
        parameterCount++;
      }
    }
  }

  const int pointsPerParameter = setup.pointsPerParameter;
  double currentChiSquared;                   // Number of points into which
  double bestChiSquared  = -1;                //   each parameter range is to
  double percentComplete = 0;                 //   be divided.
  bool   doneAdding      = false;             // True when next lattice point
                                              //   has been computed.

  int totalPoints = (int) pow(pointsPerParameter, parameterCount);
  double sizeParamSpacing = (theory.sizeParameters[theory.numSizeParameters-1] -
                             theory.sizeParameters[0])/
                             (theory.numSizeParameters-1);
  double initialSizeParam = theory.sizeParameters[0];

  if(parameterCount < 1) {                    // Nothing to fit
    cerr << "Warning: No parameters specified for fitting." << flush;
    for(i = 0; i < setup.fittingFunctionCount; i++) {
      f = setup.fittingFunctions[i];
      for(j = 0; j < f->numParams(); j++) {
        f->currentParams(j) = (f->minimumParams(j) + f->maximumParams(j))/2;
        if((f->numParams()==2) && (j==0)) {   // Check for delta function and
                                              //   make sure it falls on a size
                                              //   that's present in Mie tables.
          f->currentParams(j) = theory.sizeParameters[0] +
                                sizeParamSpacing * (int)
                                floor((f->currentParams(j)-initialSizeParam)/
                                sizeParamSpacing + 0.5);
        }
      }
    }
    setup.bestChiSquared = chiSquared();      // Compute chi^2 with these params
  } else {                                    // In this case, some parameters
                                              //   were specified for fitting,
                                              //   so explore parameter space.
    int* latticePoint = new int[parameterCount];
    for(i = 0; i < parameterCount; i++) {     // Location in parameter space
      latticePoint[i] = 0;
    }

    int* bestLatticePoint = new int[parameterCount];
                                              // Location of best fit so far


    for(i = 0; i < totalPoints; i++) {        // Check all of parameter space
      int m = 0;
      for(j = 0; j < setup.fittingFunctionCount; j++) {
        f = setup.fittingFunctions[j];
        for(int k = 0; k < f->numParams(); k++) {
          if(f->fitting(k)) {
            f->currentParams(k) = f->minimumParams(k) +
                                  (f->maximumParams(k) - f->minimumParams(k))/
                                         (pointsPerParameter-1)
                                   *latticePoint[m];
            m++;
          } else {
            f->currentParams(k) = (f->minimumParams(k) + f->maximumParams(k))/
                                                      2.0;
          }
          if((f->numParams()==2) && (k==0)) { // Check for delta functions and
                                              //   align them to the Mie tables.
            f->currentParams(k) = theory.sizeParameters[0] +
                                  sizeParamSpacing * (int)
                                  floor((f->currentParams(k)-initialSizeParam)/
                                  sizeParamSpacing + 0.5);
          }
        }
      }
      currentChiSquared = chiSquared();       // Compute chi^2 at test point and
                                              //   save results if best yet.
      if((currentChiSquared < bestChiSquared) || (bestChiSquared == -1)) {
        for(int m = 0; m < parameterCount; m++) {
          bestLatticePoint[m] = latticePoint[m];
        }
        bestChiSquared = currentChiSquared;
      }

      percentComplete = (int) (100 * (double) (i+1)/(double) totalPoints);
      cerr << "Fitting: ";                    // Progress report
      if(percentComplete < 100) cerr << " ";
      if(percentComplete < 10) cerr << "0";
      cerr << percentComplete << "%  Best chi^2: " << bestChiSquared;
      cerr << "            ";
      cerr << "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b";
      cerr << "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b";
      cerr << flush;

      doneAdding = false;                     // Now, add 1 to the lattice point
      int digit = parameterCount - 1;         //   and carry as needed to
      while(!doneAdding && (digit > 0)) {     //   generate the newest point.
        if(latticePoint[digit] >= (pointsPerParameter-1)) {
          latticePoint[digit] = 0;
        } else {
          latticePoint[digit]++;
          doneAdding = true;
        }
        digit--;
      }
      if(!doneAdding) {
        latticePoint[0]++;
      }
    }
                                              // Update the actual parameters
                                              //   of the functions to reflect
    int m = 0;                                //   the best fit found.
    for(j = 0; j < setup.fittingFunctionCount; j++) {
      for(int k = 0; k < setup.fittingFunctions[j]->numParams(); k++) {
        f = setup.fittingFunctions[j];
        if(f->fitting(k)) {
          f->currentParams(k) = f->minimumParams(k) +
                                (f->maximumParams(k) - f->minimumParams(k))/(pointsPerParameter-1)*
                                bestLatticePoint[m];
          m++;
        }
        if((f->numParams()==2) && (k==0)) {
          f->currentParams(k) = theory.sizeParameters[0] +
                                sizeParamSpacing * (int)
                                floor((f->currentParams(k)-initialSizeParam)/
                                sizeParamSpacing + 0.5);
        }
      }
    }
    setup.bestChiSquared = bestChiSquared;    // Record best chi^2 for output

    delete[] latticePoint;
    delete[] bestLatticePoint;
  } // end if fitting any params  
}
// END Simple Chi^2 Exploration Fitting Routine --------------------------------
// END Mathematical Functions **************************************************














// Support Functions ***********************************************************
void printUsage(bool extended) {
  // Input:  A boolean, telling whether a detailed usage message is desired.
  // Output: If extended is true, a detailed usage message is printed to cerr.
  //         If its value is false, a shorter message is used.

  cerr << "Usage: fitdist datafile fitfile [-tables tablefile]" 
       << " [-simple #]" << endl;

  if(extended) {
    cerr << "   datafile            - Experimental data file." << endl;
    cerr << "   fitfile             - Fit type and parameter ranges. " << endl;
    cerr << "   [-tables tablefile] - Use specified tables instead "
                                  << "of fitdist.dat." << endl;
    cerr << "   [-simple #]         - Disable steepest descent and fit "
                                  << "by sampling the" << endl;
    cerr << "                         parameter ranges at # possible points "
                                  << "and checking" << endl;
    cerr << "                         all possible combinations." << endl;
    cerr << "File formats are described in readme.txt" << endl;
  }
}

void processArguments(int argc, char* argv[]) {
  // Input:  The argc and argv from the command line.
  // Output: Sets the filenames in setup to the appropriate values, and
  //         decects a request for a simple fit, if it is present.

  setup.pointsPerParameter = 15;              // Default values
  setup.fittingFunctionCount = 0;
  setup.fittingFunctions = NULL;
  setup.fitMethod = LEVMAR_FIT;               // Do steepest descent by default

  setup.dataFilename = argv[1];               // If we've reached this point,
  setup.fitFilename = argv[2];                //   there are at least three
                                              //   arguments (counting argv[0],
                                              //   which is the path to the
                                              //   program on the command line.
  switch(argc) {
  case 3:                                     // In this case, use default Mie
    setup.tablesFilename = new char[12];      //   tables and default fit.
    strncpy(setup.tablesFilename, "fitdist.dat", 12);
    break;
  case 5:                                     // Either -tables or -simple (and
    if(strncmp(argv[3], "-tables", 7) == 0) { //   the corresponding argument)
      setup.tablesFilename = argv[4];         //   was given. Check which it is.
    } else if(strncmp(argv[3], "-simple", 7) == 0) {
      setup.fitMethod = SIMPLE_FIT;           // Switch was -simple. Get # pts.
      setup.pointsPerParameter = atoi(argv[4]);
      if(setup.pointsPerParameter < 3) {      // At least get ends and average
        cerr << "Error: Specify a number of points per paremeter, "
             << "3 or greater." << endl;
        exit(1);
      }
      setup.tablesFilename = new char[12];
      strncpy(setup.tablesFilename, "fitdist.dat", 12);
    } else {                                  // Unknown switch
      printUsage(EXTENDED);
      exit(1);
    }
    break;
  case 7:                                     // Both -tables and -simple were
    if(strncmp(argv[3], "-tables", 7) == 0) { //   provided, so figure out which
      setup.tablesFilename = argv[4];         //   is first and act accordingly.
    } else if(strncmp(argv[3], "-simple", 7) == 0) {
      setup.fitMethod = SIMPLE_FIT;
      setup.pointsPerParameter = atoi(argv[4]);
      if(setup.pointsPerParameter < 3) {
        cerr << "Error: Specify a number of points per paremeter, "
             << "3 or greater." << endl;
        exit(1);
      }
    } else {
      printUsage(EXTENDED);
      exit(1);
    }
    if(strncmp(argv[5], "-tables", 7) == 0) {
      setup.tablesFilename = argv[6];
    } else if(strncmp(argv[5], "-simple", 7) == 0) {
      setup.fitMethod = SIMPLE_FIT;
      setup.pointsPerParameter = atoi(argv[6]);
      if(setup.pointsPerParameter < 3) {
        cerr << "Error: Specify a number of points per paremeter, "
             << "3 or greater." << endl;
        exit(1);
      }
    } else {                                  // Unknown switch
      printUsage(EXTENDED);
      exit(1);
    }
    break;
  default:                                    // Improper number of arguments
    printUsage(EXTENDED);
    exit(1);
  }
}

void openFile(char* filename, ifstream& file) {
  // Input:  A filename, and an ifstream (passed by reference).
  // Output: Opens the file, if it exists, and attaches it to the ifstream,
  //         exiting with an error if the file does not exist, or cannot
  //         be read.

  file.open(filename, ios::in | ios::nocreate,
                filebuf::sh_read);            // Open file if and only if it
                                              //   already exists (nocreate)
                                              //   but allow others to read
                                              //   from the file (sh_read).
  if (!file.good()) {                         // File doesn't exist
    cerr << "Error: Unable to read from file '" << filename;
    cerr << "'." << endl;
    exit(1);
  }
}

void readData(char* dataFilename) {
  // Input:  A data file name.
  // Outupt: Checks the validity of the file. If it is valid, the file is read
  //         into the data struct. Otherwise, an error is output and the
  //         program exits.

  ifstream dataFile;

  openFile(dataFilename, dataFile);           // Open the data file

  while(dataFile.peek() == '/') {             // Ignore comment lines by reading
    dataFile.getline((char*) NULL, 1000);     //   into NULL.
  }

  data.polarization   = -1;                   // Set values to something out
  data.wavelength     = -1;                   //   of range, so after reading
  data.relativeIndex  = -1;                   //   file, we know they were ok.
  data.ambientIndex   = -1;
  data.scattererIndex = -1;

  char* item = new char[100];                 // Next item to be read and set

  dataFile >> item;                           // Read first item
  while(isdigit(item[0]) == 0) {              // As long as we haven't gotten
                                              //   to the angle-intensity data
    switch(item[0]) {                         //   yet, process items.
    case 'p':                                 // Assume it's "polarization"
      dataFile >> item;                       // Ignore = sign
      dataFile >> item;
      switch(item[1]) {                       // Check second character.
      case 'a':
        data.polarization = PARALLEL;
        break;
      case 'e':
        data.polarization = PERPENDICULAR;
        break;
      case 'n':
        data.polarization = UNPOLARIZED;
        break;
      default:
        cerr << "Error: Format of '" << dataFilename << "' is invalid." << endl;
        cerr << "Unknown polarization '" << item << "'." << endl;
        exit(1);
      }
      break;
    case 'w':
      dataFile >> item;
      dataFile >> data.wavelength;
      break;
    case 's':
      dataFile >> item;
      dataFile >> data.scattererIndex;
      break;
    case 'a':
      dataFile >> item;
      dataFile >> data.ambientIndex;
      break;
    case 'i':
    case 'r':
      dataFile >> item;
      dataFile >> data.relativeIndex;
      break;
    default:
      cerr << "Error: Format of '" << dataFilename << "' is invalid." << endl;
      cerr << "Unknown item '" << item << "'." << endl;
      exit(1);
    }
    dataFile >> item;                         // Read next item (or possibly the
  }                                           //   first angle in the data).

  if((data.relativeIndex == -1) && (data.scattererIndex != -1)) {
    if(data.ambientIndex == -1) {             // No relative index, so use the
      data.ambientIndex = N_AIR;              //   index of scatterers over the
    }                                         //   ambient index (or air).
    data.relativeIndex = data.scattererIndex/data.ambientIndex;
  } else if((data.relativeIndex != -1) && (data.scattererIndex == -1)) {
    if(data.ambientIndex == -1) {             // Compute scatterer index from
      data.ambientIndex = N_AIR;              //   relative and ambient indices.
      data.scattererIndex = data.relativeIndex;
    } else {
      data.scattererIndex = data.relativeIndex * data.ambientIndex;
    }
  } else if((data.relativeIndex != -1) && (data.scattererIndex != -1)) {
    if(data.ambientIndex == -1) {             // Infer ambient index
      data.ambientIndex = data.scattererIndex/data.relativeIndex;
    } else {                                  // Redo the division to be sure
      data.relativeIndex = data.scattererIndex/data.ambientIndex;
    }
  }

  if((data.wavelength    == -1) ||            // Make sure all items were given
     (data.relativeIndex == -1) ||
     (data.polarization  == -1)) {
    cerr << "Error: Format of '" << dataFilename << "' is invalid." << endl;
    cerr << "File must include polarization, wavelength, and "
         << "relativeIndex or" << endl;
    cerr << "scattererIndex and ambientIndex. (See readme.txt)" << endl;
    exit(1);
  }

  dataFile >> item >> item;                   // Already read first angle. Now
  data.numPoints = 0;                         //   read rest of lines and count.
  while((isdigit(item[0]) != 0) && (!dataFile.eof())) {
    data.numPoints++;
    dataFile >> item >> item >> item;
  }

  dataFile.close();                           // At end of file, so close it

  data.angles = new double[data.numPoints];   // Allocate space for data
  data.intensities = new double[data.numPoints];
  data.uncertainties = new double[data.numPoints];

  openFile(dataFilename, dataFile);           // Open the data file again

  while(isdigit(dataFile.peek()) == 0) {      // Ignore lines that don't
    dataFile.getline((char*) NULL, 1000);     //   start with digits.
  }

  int i = 0;
  for(i = 0; i < data.numPoints; i++) {       // Read values from file
    dataFile >> data.angles[i] >> data.intensities[i] >> data.uncertainties[i];
    if((i > 0) && (data.angles[i-1] >= data.angles[i])) {
      cerr << "Error: Format of '" << dataFilename << "' is invalid." << endl;
      cerr << "Data points must be unique and in sorted order. Check ";
      cerr << " the point at angle " << data.angles[i] << "." << endl;
      exit(1);
    }
  }

  dataFile.close();
  delete[] item;                              // Free memory
}

void readNecessaryMieTables(char* tablesFilename) {
  // Input:  A filename for the Mie tables, the struct in which to store
  //         them, and the experimental data already read in.
  // Output: Reads only those angles present in the theory tables, storing
  //         data in the theory struct.

  ifstream tablesFile;

  openFile(tablesFilename, tablesFile);       // Try to open the file

  char ch;
  tablesFile.get(ch);
  if(ch != '*') {
    cerr << "Error: Mie tables in file '" << tablesFilename << "' do not begin";
    cerr << " with * to signal first size parameter." << endl;
    exit(1);
  }
  tablesFile.getline((char*) NULL, 1000);     // Throw out rest of the * line
  
  double firstAngle;
  double secondAngle;
  tablesFile >> firstAngle;                   // Get first two angles and
  tablesFile.getline((char*) NULL, 1000);     //   compute the angle step so
  tablesFile >> secondAngle;                  //   angle range can be error
  tablesFile.getline((char*) NULL, 1000);     //   checked.
  theory.mieTablesAngleStep = secondAngle-firstAngle;

  int skipCount = 2;                          // Already read two angles
  while(tablesFile.peek() != '*') {
    skipCount++;
    tablesFile.getline((char*) NULL, 1000);
  }

  theory.mieTablesMaxAngle = firstAngle + 
                             (skipCount-1)*theory.mieTablesAngleStep;

  if(data.angles[data.numPoints-1] > theory.mieTablesMaxAngle) {
                                              // skipCount is one more than
                                              //   the greatest angle in the
                                              //   tables file.
    cerr << "Error: Mie tables in '" << tablesFilename << "' do not";
    cerr << " contain angles past " << theory.mieTablesMaxAngle << "." << endl;
    cerr << "Data file contains angles up to " << data.angles[data.numPoints-1];
    cerr << ". Please use larger tables." << endl;
    exit(1);
  }

  int i = 0;
  theory.numSizeParameters = 1;               // Initialize number of sizes. One
                                              //   was already read above.
  while(!tablesFile.eof()) {                  // Count *'s to determine number
    tablesFile.get(ch);                       //   of size params present.
    if(ch == '*') {
      theory.numSizeParameters++;
      tablesFile.get((char*) NULL, 1000);
    }
    for(i = 0; i <= skipCount; i++) {     // Ignore rest of this * line,
      tablesFile.get((char*) NULL, 1000);     //   and all actual data.
      tablesFile.ignore();                    // Ignore \n
    }
  }

  tablesFile.close();                         // Number of sizes is now known

  openFile(tablesFilename, tablesFile);       // Reopen the tables for reading

  theory.numAngles = data.numPoints;          // Will read only those angles
  theory.angles    = data.angles;             //   present in the data.
  theory.sizeParameters = new double[theory.numSizeParameters];

  double percentComplete = 0;                 // Percent of load complete
  double currentAngle;                        // Current angle in tables file
  char* garbage = new char[1000];             // For undesired strings & lines

  theory.intensities = new double*[theory.numSizeParameters];
  int j = 0;
  for(i = 0; i < theory.numSizeParameters; i++) {
    theory.intensities[i] = new double[theory.numAngles];

    tablesFile.get(ch);                       // Get * and ignore space
    tablesFile.ignore();
    tablesFile >> theory.sizeParameters[i];

    cerr << "Loading: ";
    percentComplete = (int)(100*((double) i/(double) theory.numSizeParameters));
    if(percentComplete < 10) {
      cerr << "0";
    }
    cerr << percentComplete << "%\b\b\b\b\b\b\b\b\b\b\b\b" << flush;

    for(j = 0; j < theory.numAngles; j++) {
      tablesFile >> currentAngle;
      while(currentAngle < theory.angles[j]) {
        tablesFile >> garbage >> garbage >> garbage >> currentAngle;
      }
      if(currentAngle != theory.angles[j]) {
        cerr << "Error: Angle " << theory.angles[j] << " in data not present ";
        cerr << "in Mie tables." << endl;
        exit(1);
      }
      switch(data.polarization){              // Read appropriate column for the
      case PERPENDICULAR:                     //   chosen polarization.
        tablesFile >> theory.intensities[i][j] >> garbage >> garbage;
        break;
      case PARALLEL:
        tablesFile >> garbage >> theory.intensities[i][j] >> garbage;
        break;
      case UNPOLARIZED:
        tablesFile >> garbage >> garbage >> theory.intensities[i][j];
        break;
      }
    }

    while((tablesFile.peek() != '*') && (!tablesFile.eof())) {
      tablesFile.getline(garbage, 1000);      // Skip the rest of the angles
    }
  }

  cerr << "Loading: 100%  Range: x = " << theory.sizeParameters[0] << " to "
       << theory.sizeParameters[theory.numSizeParameters-1] << "." << endl;
  tablesFile.close();                         // Completed load, so output table
  delete[] garbage;                           //   range and free memory.
}

void readFit(char* fitFilename) {
  // Input:  A filename and the fitting setup struct, along with the data
  //         so that wavelength and index of refraction can be used in
  //         the conversion from microns to size parameters.
  // Output: Reads information from the fit file into the fitting setup,
  //         including types of functions and ranges for parameters.


  delete[] setup.fittingFunctions;            // Delete old fitting functions
  setup.fittingFunctions = NULL;              //   and set the array to NULL.

  ifstream fitFile;

  openFile(fitFilename, fitFile);             // Open fit file.

  char ch;                                    // First, count how many fitting
  setup.fittingFunctionCount = 0;             //   functions are present.
  while(!fitFile.eof()) {                     
    fitFile.get(ch);
    if(ch == '[') {
      setup.fittingFunctionCount++;
    }
  }
  fitFile.close();
  setup.fittingFunctions = new fittingFunction*[setup.fittingFunctionCount];


  openFile(fitFilename, fitFile);
  
  char* item = new char[100];                 // Buffer for reading from file

  while((!fitFile.eof()) && (fitFile.peek() != '[')) {
    fitFile.getline((char*) NULL, 1000);      // Ignore lines until reaching a
  }                                           //   [FittingFunc] item.
  

  fittingFunction* f;

  int i = 0;
  for(i = 0; i < setup.fittingFunctionCount; i++) {
    fitFile >> item;

    switch(tolower(item[1])) {
    case 'g':                                 // Gaussian
      f = new gaussianFunction();
      f->readFromFile(fitFile);
      break;
    case 'd':                                 // Delta
      if(setup.fitMethod != SIMPLE_FIT) {
        cerr << "Warning: Delta functions incompatible with steepest "
             << "descent fitting." << endl;
        cerr << "         Switching to simple fit." << endl;
        setup.fitMethod = SIMPLE_FIT;         // Can't take derivatives very
      }                                       //   effectively, so don't try.
      f = new deltaFunction();
      f->readFromFile(fitFile);
      break;
    case 'l':                                 // Logonormal
      f = new logonormalFunction();
      f->readFromFile(fitFile);
      break;
    default:                                  // Unknown
      cerr << "Error: Fitting file '" << fitFilename << "' contains unknown";
      cerr << "function type " << item << "." << endl;
      exit(1);
    }

    f->toSizeParams(data.ambientIndex, data.wavelength);
                                              // Convert to size parameters in
                                              //   order to perform range check.
    double minX = theory.sizeParameters[0];
    double maxX = theory.sizeParameters[theory.numSizeParameters-1];

    if(!(f->inRange(minX, maxX))) {           // Check range. If out of range,
                                              //   minX and maxX are given the
                                              //   needed bounds for the tables.
      cerr << "Warning: Mie tables may not cover range of sizes needed." << endl;
      cerr << "         95% of the " << i+1;
      switch(i+1) {
      case 1:
        cerr << "st";
        break;
      case 2:
        cerr << "nd";
        break;
      case 3:
        cerr << "rd";
        break;
      default:
        cerr << "th";
      }
      cerr << " function lies in the range " << floor(100*minX)/100.0 << " to ";
      cerr << ceil(100.0*maxX)/100.0 << "." << endl;
    }

    setup.fittingFunctions[i] = f;            // Store the new function

    while((!fitFile.eof()) && (fitFile.peek() != '[')) {
      fitFile.getline((char*) NULL, 1000);    // Move to next [FittingFunc]
    }
  }
  fitFile.close();                            // Close and clean up
  delete[] item;
}

void outputSummary() {
  // Input:  None.
  // Output: Results of the fit, in abbreviated form.

  cerr << endl << "Results: " << endl;
  int i = 0;
  for(i = 0; i < setup.fittingFunctionCount; i++) {
    setup.fittingFunctions[i]->toMicronsR(data.ambientIndex, data.wavelength);
    cerr << "  " << *setup.fittingFunctions[i] << endl;
  }
}


bool fullOutputRequested() {
  // Input:  None.
  // Ouptut: Asks the user if another fit is desired.

  char ch;
  cerr << "Output full results to file? (y/n): " << flush;
  cin >> ch;
  while(cin.peek() != '\n') {
    cin.ignore();
  }
  cin.ignore();

  if(tolower(ch) == 'y') {
    return true;
  } else {
    return false;
  }
}


char* constructDefaultOutputName() {
  // Input:  None.
  // Output: A default output name based on the data input filename.

  ifstream file;

  char path[256];                             // Arrays to store the path,
  char filename[256];                         //   filename, and extension of
  char extension[256];                        //   the data file.
  strncpy(path, setup.dataFilename, 256);
  char* lastSlash = strrchr(path, '\\');      // Find end of path part

  if(lastSlash == NULL) {                     // No path was specified
    strncpy(filename, path, 256);
    path[0] = '\0';
  } else {                                    // copy everything after the last
    strncpy(filename, &(lastSlash[1]), 256);  //   slash into the filename, then
    lastSlash[1] = '\0';                      //   terminate the path after the
  }                                           //   last slash.

  char* lastPeriod = strrchr(filename, '.');  // Find beginning of extension

  if(lastPeriod == NULL) {                    // No extension
    extension[0] = '\0';
  } else {                                    // Copy the extension, then end
    strncpy(extension, lastPeriod, 256);      //   the filename before it.
    lastPeriod[0] = '\0';
  }

  int i = 0;                                  // Now check all files of the
  char digit[2];                              //   format filename-fit##.ext
  char* fullFilename = new char[270];         //   until an unused filename is
  bool existsOpenFile = false;                //   located. Make sure to close
  do {                                        //   files properly as they're
    if(existsOpenFile) {                      //   checked.
      file.close();
    }

    if(lastSlash != NULL) {                   // Construct candidate filename
      strncpy(fullFilename, path, 256);
    } else {
      strncpy(fullFilename, "\0", 1);
    }
    strncat(fullFilename, filename, 256);
    strncat(fullFilename, "-fit", 5);         // Add on the -fit
    if(i < 10) {
      strncat(fullFilename, "0", 2);
    } else {
      digit[0] = (char) ((i % 100)/10 + '0');
      digit[1] = '\0';
      strncat(fullFilename, digit, 2);
    }
    char digit[2];                            // Compute next number using the
    digit[0] = (char) ((i % 10) + '0');       //   mod function and some clever
    digit[1] = '\0';                          //   character arithmetic.
    strncat(fullFilename, digit, 2);

    strncat(fullFilename, extension, 256);    // Tack on extension

    file.open(fullFilename, ios::in | ios::nocreate,
                  filebuf::sh_read);          // Open file if and only if it
                                              //   already exists (nocreate)
                                              //   but allow others to read
                                              //   from the file (sh_read).
    existsOpenFile = true;                    // So it can be closed later
    i++;
  } while(file.good() && (i < 100));          // Continue to loop as long as
                                              //   the file exists and we
  file.close();                               //   haven't gotten 100 yet.
  return fullFilename;                        // Return first unused filename
}


char* getFilename() {
  // Input:  None.
  // Output: Prompts the user for a filename, and returns it.

  char* filename = new char[200];

  cerr << "Filename: " << flush;
  cin >> filename;

  return filename;
}

char* getFilename(char* defaultFilename) {
  // Input:  A default filename to be used if the user doesn't input a
  //         different one.
  // Output: Prompts the user for a filename, and returns it.

  char* filename = new char[250];

  cerr << "Filename (" << defaultFilename << "): " << flush;
  cin.getline(filename, 250);

  if(strncmp(filename, "", 1)==0) {           // If the user didn't input a
    return defaultFilename;                   //   filename, use the default.
  } else {
    return filename;
  }
}

void outputFullResults() {
  // Input:  The fitting setup should contain the best fit found.
  // Output: Outputs full fitting results to the file in setup.outputFilename.

  ofstream outputFile(setup.outputFilename);  // Initialize output file stream
  if(outputFile.bad()) {
    cerr << "Error: Unable to open output file " << setup.outputFilename;
    cerr << "." << endl;
    exit(1);
  }

  outputFile << "Distributed fit using fitdist " << VERSION << "." << endl;
  outputFile << "  Using ";
  switch(setup.fitMethod) {
  case SIMPLE_FIT:
    outputFile << "simple chi^2 exploration fit with "
               << setup.pointsPerParameter
               << " points per parameter.";
    break;
  case LEVMAR_FIT:
    outputFile << "Levenberg-Marquardt fit.";
    break;
  default:
    outputFile << "unknown fit method.";
  }
  outputFile << endl << endl;

  outputFile << "Data:" << endl;              // First copy data file into the
                                              //   output, so that this fit
  char* buffer = new char[1000];              //   can be reproduced.
  ifstream dataFile;
  openFile(setup.dataFilename, dataFile);
  while(dataFile.good()) {
    dataFile.getline(buffer, 1000);
    outputFile << buffer << endl;
  }
  dataFile.close();

  outputFile << endl << endl;

  outputFile << "Parameters:" << endl;        // Also output the fit file

  ifstream fitFile;
  openFile(setup.fitFilename, fitFile);
  while(fitFile.good()) {
    fitFile.getline(buffer, 1000);
    outputFile << buffer << endl;
  }
  fitFile.close();
  delete[] buffer;

  outputFile << endl << endl;

  outputFile << "Mie Tables:" << endl;        // Report on Mie tables used
  outputFile << "  Size Range: " << theory.sizeParameters[0];
  outputFile << " to " << theory.sizeParameters[theory.numSizeParameters-1];
  if(theory.numSizeParameters > 1) {
    outputFile << " in steps of "
               << theory.sizeParameters[1]-theory.sizeParameters[0];
  }
  outputFile << "." << endl;
  outputFile << "  Angles:     0 to " << theory.mieTablesMaxAngle << " degrees "
             << "in steps of " << theory.mieTablesAngleStep << "." << endl;
  
  outputFile << endl << endl;

  outputFile << "Summary of Results:    " << endl;
  outputFile << "  Chi^2:      " << setup.bestChiSquared << endl;

  int i = 0;
  for(i = 0; i < setup.fittingFunctionCount; i++) {
    setup.fittingFunctions[i]->toMicronsR(data.ambientIndex, data.wavelength);
    outputFile << "  " << *setup.fittingFunctions[i] << endl;
  }                                           // Convert to microns and output
                                              //   results summary.
  outputFile << endl << endl;

  outputFile << "Size Distribution:" << endl; // Sizes and abundances
  outputFile << "  radius (um)\trelative abundance" << endl;

  double* theoryDistribution = new double[theory.numAngles];
                                             // Properly scaled Mie theory
  for(i = 0; i < theory.numAngles; i++) {
    theoryDistribution[i] = 0;
  }

  for(i = 0; i < setup.fittingFunctionCount; i++) {
    setup.fittingFunctions[i]->toSizeParams(data.ambientIndex, data.wavelength);
  }

  int j = 0;
  for(i = 0; i < theory.numSizeParameters; i++) {
    double currentCoefficient = 0;
    for(j = 0; j < setup.fittingFunctionCount; j++) {
      currentCoefficient += setup.fittingFunctions[j]->
        evaluate(theory.sizeParameters[i]);
    }

    if(currentCoefficient != 0) {
      outputFile.setf(ios::fixed);
      outputFile << setw(5) << setprecision(2) << theory.sizeParameters[i]*
                    (data.wavelength*0.001/(2*PI*data.ambientIndex));
      outputFile.unsetf(ios::fixed);
      outputFile.setf(ios::floatfield);
      outputFile << setprecision(3) << setw(10) << "\t" << currentCoefficient << endl;
      outputFile.unsetf(ios::floatfield);
      for(j = 0; j < theory.numAngles; j++) {
        theoryDistribution[j] += currentCoefficient*theory.intensities[i][j];
      }
    }
  }

  outputFile << endl << endl;                 // Now, output scaled theory

  outputFile << "Theoretical Intensity Profile:" << endl;
  outputFile << "  Angle\t\tIntensity" << endl;

  for(j = 0; j < theory.numAngles; j++) {
    outputFile << theory.angles[j];
    outputFile.setf(ios::fixed);
    outputFile << setw(9) << setprecision(6) << "\t"
               << theoryDistribution[j] << endl;
    outputFile.unsetf(ios::fixed);
  }

  delete[] theoryDistribution;                // Free memory

  outputFile << endl;

  outputFile.close();
}

bool additionalFitRequested() {
  // Input:  None.
  // Output: Asks the user whether or not to perform another fit.

  char ch;
  cerr << "Perform another fit? (y/n): " << flush;
  cin >> ch;
  while(cin.peek() != '\n') {
    cin.ignore();
  }
  cin.ignore();

  if(tolower(ch) == 'y') {
    return true;
  } else {
    return false;
  }
}
// END Support Functions *******************************************************