// Filename: fitmie.cpp
// Author:   Mike Rust
// Modified: Mike Schubmehl
// Date(s):  June 1998, 07/13/2000
// Purpose:  Fits theoretical Mie scattering curves to experimental data using
//           one of a variety of fitting methods, including:
//
//           1. A standard alpha regression. In this method, the theory curve
//              is first scaled so that it has the same amplitude as the data,
//              taking into account the uncertainty at each point in the data
//              set. This new theory curve is then compared to the data set by
//              computing a chi^2 (i.e., the sum of the squares of the
//              residuals at each point in the data set, divided by the number
//              of points in the set.) A lower chi^2 indicates a better fit.
//           2. An alpha regression disregarding uncertainties. Same as above,
//              but without taking into account uncertainty. This could be
//              useful for data without uncertainties, or suspect uncertainty.
//           3. A divide-by-point normalization. This method, suggested by
//              Bohren and Huffman, uses a scaling factor for both the data
//              and theory. Both scaling factors are computed by dividing all
//              of the intensities present by the intensity at some chosen
//              angle. A chi^2 is then computed as above, taking into account
//              uncertainties.
//           4. A front-back fit. This fitting method takes the average of a
//              given number of points from the front and back of the fitting
//              range in the data and in the theory. It then compares the ratio
//              of these to averages in the data and theory, and uses the
//              squared difference between the theory ratio and data ratio to
//              measure fit quality.
//
//             The program takes as input a data table in the format:
//
//                  4    24782.229965156796   1526.0627957332767
//                  6    4436.597110754414    575.72747321778
//                  8    1076.4872521246461   158.36175478043836
//               <angle>     <intensity>        <uncertainty>
//                  .             .                    .
//                  .             .                    .
//                  .             .                    .
//
//             It also takes a theory table in MIEV0 format:
//
//                 *	0.500000
//                      0.000000	0.000269	0.000269	0.000269
//                      1.000000	0.000269	0.000269	0.000269
//                      2.000000	0.000269	0.000269	0.000269
//                      3.000000	0.000269	0.000268	0.000269
//                       <angle>   <perp>     <par>     <unp>
//                          .         .         .         .
//                          .         .         .         .
//                          .         .         .         .
//                 * <size param>...
//
//             The program reads in the theory tables one at a time, computes
//           the appropriate scaling factor according to the fit metohd, and
//           then computes a chi^2 (or other measure of fit quality). It keeps
//           track of the minimum chi^2 encountered, the size parameter that
//           produced it, and the corresponding alpha (scaling) factor. It then
//           prints the results of the fit, dumps the data, dumps the chi^2 at
//           each point, or dumps the scaled theory at a particular intensity,
//           depending on the command line switches specified.
//             Command-line syntax is:
//
//  fitmie datafile theoryfile [-pol {1,2,3}] [-nouncert]
//         [-dividebypoint double] [-frontback int] [-dump {data, chi, double}]
//         [-range double double] [-help | -?]
//
// datafile                  - Experimental data (angle-intensity-sigma table).
// theoryfile                - Theoretical Mie tables in MIEV0 format.
// [-pol {1,2,3}]            - Polarization. 1 = perp, 2 = par, 3 = unp.
// [-nouncert]               - Fit by alpha regression without uncertainties.
// [-dividebypoint double]   - Fit by Bohren and Huffman's divide by point
//                             method, using the specified angle.
// [-frontback int]          - Fit by back-to-front ratio using averages of
//                             a number of points from front and back of data.
// [-dump {data,chi,double}] - Dump input data, chi-square, or the fit curve
//                             at the size parameter nearest to the double.
// [-range double double]    - Fit only data in the range from the first angle
//                             specified to the second angle specified.
// [-help | -?]              - Print help message.
//
// Notes:      Requires stdlib.h, iostream.h, fstream.h, string.h. Uses
//           dataset.hh for the dataSet class in which the theory data is
//           stored for use during the correction.

#include <stdlib.h>
#include <iostream.h>
#include <fstream.h>
#include <string.h>
#include "dataset.hh"                         // dataSet class

                                              //       -- FIT MODES --
#define ALPHA         1                       // Alpha regression with uncert
#define ALPHANOUNCERT 2                       // Disregard uncertainties
#define DIVIDEBYPOINT 3                       // Bohren & Huffman div by point
#define FRONTBACK     4                       // Front-back fit

#define PERPENDICULAR 1                       // Polarization constants. These
#define PARALLEL      2                       //   should be as they appear in
#define UNPOLARIZED   3                       //   dataset.hh.

#define EXTENDED      true                    // For use with the printUsage
#define SIMPLE        false                   //   routine.

struct configuration {                        // Parameters for the correction
  char*  dataFilename;                        // Filenames for data and theory
  char*  theoryFilename;                      //   files to be opened.

  int    polarization;                        // Polarization of incident light
  int    fitMode;                             // Fitting mode to use

  double pointToDivide;                       // Used for the DIVIDEBYPOINT
                                              //   fitting method.
  int    frontBackPoints;                     // Used for FRONTBACK fitting

  double startFit;                            // The angle range over which to
  double endFit;                              //   perform the fit.

  bool   dumpData;                            // True if the corresponding
  bool   dumpChi;                             //   object is to be output.
  bool   dumpTheory;                          //   Only one can be true.
  double dumpTheoryX;                         // Size parameter for which to
                                              //   output scaled theory curve.
};


// Table of Contents **********************************************************
// -- Mathematical Support Functions --
bool   inFitRange(int i, dataSet& data,       // Returns true if the ith point
                  configuration& config);     //   of the dataSet is in the
                                              //   angle range config.startFit
                                              //   to config.endFit.
double scaleFactor(dataSet& data,             // Returns the scale factor by
                   dataSet& theory,           //   which to multiply the theory
                   configuration& config);    //   so it has the same amplitude
                                              //   as the data in given mode.
double chiSquared(dataSet& data,              // Computes the chi^2 between 
                  dataSet& theory,            //   data and scaled theory,
                  configuration& config,      //   based on user-selected mode.
                  double alpha);

// -- Support Functions --
void setDefaults(configuration& config);      // Setup default configuration
void processArguments(int argc, char* argv[], // Process command line arguments
                      configuration& config); //   and modify configuration.
void checkConfiguration(configuration config);// Verify that all values are
                                              //   acceptable or print errors.
void openFiles(char* theoryFilename,          // Opens the files specified by
               ifstream& theoryFile,          //   the given filenames into
               char* dataFilename,            //   the ifstreams. Prints errors
               ifstream& dataFile);           //   if the files do not exist.
void outputScaledTheory(dataSet& theory,      // Write the theory dataSet to
                        double alpha);        //   cout with values multiplied
                                              //   by alpha.
void printResults(double sizeParameter,       // Write the results of the fit
                  double chiSquare,           //   to cout.
                  double alpha);
void printUsage(bool extended);               // Prints useage message to cerr,
                                              //   with detailed help if
                                              //   extended is set to true.
// END Table of Contents ******************************************************


// Main Program ***************************************************************
int main(int argc, char* argv[])  {
  configuration config;                       // Parameters for this fit

  setDefaults(config);                        // Setup default configuration
  processArguments(argc, argv, config);       // Process command line arguments
  checkConfiguration(config);                 // Make sure all values are ok

  ifstream theoryFile, dataFile;              // Create ifstreams and attach
  openFiles(config.theoryFilename, theoryFile,//   the theory and data files.
            config.dataFilename, dataFile);

  dataSet data;
  if (!data.loadData(dataFile)) {
    cerr << "Error: Unable to read data from file " << config.dataFilename;
    cerr << "." << endl;
    exit(1);
  }

  if (config.dumpData) {                      // Output the data file to cout
    cout << data;                             //   if the user wants data dump.
  }

  double bestFitChiSquared = -1;              // Set flag for first pass. These
  double bestFitAlpha;                        //   variables hold information
  double bestFitSizeParameter;                //   on the best fit.

  while(!theoryFile.eof()) {                  // As long as there are more size
    dataSet theory;                           //   parameters remaining, loop
    double sizeParameter;                     //   through the theory file and
                                              //   compute chi^2 at each param.
    while(!theoryFile.eof() && theoryFile.peek() != '*') {
      theoryFile.ignore();                    // Read file until end or a *
    }                                         //   is found, indicating the
    theoryFile.ignore();                      //   size parameter is next.
                                              //   Then skip the * and read
    theoryFile >> sizeParameter;              //   the size parameter.

    if(!theory.loadTheory(theoryFile, config.polarization)) { // Load theory
      cerr << "Error: Unable to read Mie theory from file '";
      cerr << config.theoryFilename << "'." << endl;
      exit(1);
    }

    double currentAlpha = scaleFactor(data, theory, config);
    double currentChiSquared = chiSquared(data, theory, config, currentAlpha);

    if ((bestFitChiSquared == -1) ||          // Compare best values to these
        (currentChiSquared < bestFitChiSquared)) {
      bestFitChiSquared = currentChiSquared;  // If the current values are
      bestFitAlpha = currentAlpha;            //   better, or if there are no
      bestFitSizeParameter = sizeParameter;   //   best values (best chi^2 is
    }                                         //   still -1), record current
                                              //   values as the best so far.

    if (config.dumpChi) {                     // If the user requested that
      cout << sizeParameter << "\t";          //   chi^2 be dumped, output
      cout << currentChiSquared << endl;      //   the current chi^2 and the
    }                                         //   size parameter.

    if (config.dumpTheory &&
        ((sizeParameter >= config.dumpTheoryX) || (theoryFile.eof()))) {
      cout << "* " << sizeParameter << endl;
      outputScaledTheory(theory, currentAlpha);// If user wanted theory dumped
      return 0;                                //   at this size parameter, do
    }                                          //   it and exit. If eof is true
  }                                            //   dump last theory.

  if (config.dumpChi) {                       // Now that all chi^2 have been
    return 0;                                 //   printed, exit without error.
  }

  if (!config.dumpChi && !config.dumpData && !config.dumpTheory) {
    printResults(bestFitSizeParameter,        // If no special output options
                 bestFitChiSquared,           //   were specified, just print
                 bestFitAlpha);               //   the results of the fit.
  }

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


// Mathematical Support Functions *********************************************
bool inFitRange(int i, dataSet& data, configuration& config) {
  // Input:  An index i, the dataSet and configuration to be used in
  //         determining whether the ith point in the dataSet is in the
  //         angle range config.startFit .. config.endFit.
  // Output: True if the ith point of the dataSet has an angle that is in
  //         the range config.startFit .. config.endFit.

  return ((i < data.getSize()) &&
          (data.getPoint(i)->getAngle() <= config.endFit) &&
          (data.getPoint(i)->getAngle() >= config.startFit));
}

double scaleFactor(dataSet& data, dataSet& theory,
                   configuration& config) {
  // Input:  Experimental and theory data, and a configuration.
  // Output: Returns the scale factor by which to multiply the theory in
  //         order to give it the same amplitude as the data. Takes into
  //         account the fit mode being used, and the fit range.

  double dataValue;                           // Values of data and theory at
  double theoryValue;                         //   current angle.

  double numerator   = 0;                     // The scale factor alpha is
  double denominator = 0;                     //   numerator/denominator.

  int i = 0;                                  // First move to start of fit
  while(!inFitRange(i, data, config)) {       //   range.
    i++;
  }

  switch(config.fitMode) {
    case ALPHA:
      while(inFitRange(i, data, config)) {    // Sum over entire fit range
        dataValue = data.getPoint(i)->getValue();
        double theta = data.getPoint(i)->getAngle();
        double sigma = data.getPoint(i)->getSigma();
        theoryValue = theory.interpolateIntensity(theta);

        numerator   += (dataValue*theoryValue)/(sigma*sigma);
        denominator += (theoryValue*theoryValue)/(sigma*sigma);
        i++;
      }

      return numerator/denominator;
      break;
    case ALPHANOUNCERT:
      while(inFitRange(i, data, config)) {
        dataValue = data.getPoint(i)->getValue();
        double theta = data.getPoint(i)->getAngle();
        theoryValue = theory.interpolateIntensity(theta);

        numerator   += (dataValue*theoryValue);
        denominator += (theoryValue*theoryValue);
        i++;
      }

      return numerator/denominator;

      break;
    case DIVIDEBYPOINT:
      return 1/theory.interpolateIntensity(config.pointToDivide);
      break;
    case FRONTBACK:
      return 1;                               // No scaling for front-back
      break;
  }

  return 1;
}

double chiSquared(dataSet& data, dataSet& theory,
                  configuration& config, double alpha) {
  // Input:  Experimental data and theory, a configuration containing
  //         information on fitting mode, ranges, etc., and the scaling
  //         factor to make the theory fit the data.
  // Output: Returns chi-squared (or some other measure of fit quality,
  //         depending on the fit mode).

  double result = 0;                          // Holds partial chi^2 values
  int    count  = 0;                          // Number of points in fit range

  int i = 0;
  while(!inFitRange(i, data, config)) {       // Keep incrementing i until it's
    i++;                                      //   at the start of the fit
  }                                           //   range.

  switch(config.fitMode) {
    case ALPHA:
      while (inFitRange(i, data, config)) {    // Sum residuals over fit range
        double difference = alpha *
          theory.interpolateIntensity(data.getPoint(i)->getAngle()) -
          data.getPoint(i)->getValue();
        double sigma = data.getPoint(i)->getSigma();
        result += (difference*difference)/(sigma*sigma);
        count++;
        i++;
      }

      result /= count;                        // Want chi-squared per degree
                                              //   of freedom.
      break;
    case ALPHANOUNCERT:
      while(inFitRange(i, data, config)) {    // Sum over fit range
        double difference = alpha *
          theory.interpolateIntensity(data.getPoint(i)->getAngle()) -
          data.getPoint(i)->getValue();
        result += (difference*difference);
        count++;
        i++;
      }

      result /= count;                        // Per degree of freedom
      break;
    case DIVIDEBYPOINT:
      double experimentalAlpha;               // Get experimental scale factor
      bool foundDividePoint;
      int j;

      foundDividePoint = false;
      for(j = 0; j < data.getSize(); j++) {
        if (data.getPoint(j)->getAngle() == config.pointToDivide) {
          experimentalAlpha = 1/data.getPoint(j)->getValue();
          foundDividePoint = true;
        }
      }

      if (!foundDividePoint) {
        cerr << "Error: Unable to divide by point at " << config.pointToDivide;
        cerr << " degrees, because that point is not" << endl;
        cerr << "in the data set." << endl;
        exit(1);
      }

      while(inFitRange(i, data, config)) {    // Sum over fit range
        double difference = alpha *
          theory.interpolateIntensity(data.getPoint(i)->getAngle()) -
          data.getPoint(i)->getValue() * experimentalAlpha;
        double sigma = data.getPoint(i)->getSigma() * experimentalAlpha;
        result += (difference*difference)/(sigma*sigma);
        count++;
        i++;
      }

      result /= count;                        // Per degree of freedom
      break;
    case FRONTBACK:
      double frontDataAvg, backDataAvg;       // Front and back data averages
      double frontTheoryAvg, backTheoryAvg;   // Front and back theory averages

      count = 0;
      frontDataAvg = frontTheoryAvg = 0;      // Initialize front averages
      while(count < config.frontBackPoints) { // Sum values
        frontDataAvg += data.getPoint(i)->getValue();
        double theta = data.getPoint(i)->getAngle();
        frontTheoryAvg += theory.interpolateIntensity(theta);
        i++;
        count++;
      }
      frontDataAvg /= config.frontBackPoints; // Turn sums into averages
      frontTheoryAvg /= config.frontBackPoints;

      i--;                                    // Step back in case we've passed
                                              //   end of fit range (averaging
                                              //   all points in range).
      while(inFitRange(i+1, data, config)) {  // Advance to end of fit range to
        i++;                                  //   compute back average.
      }
      count = 0;
      backDataAvg = backTheoryAvg = 0;
      while(count < config.frontBackPoints) { // Sum values
        backDataAvg += data.getPoint(i)->getValue();
        double theta = data.getPoint(i)->getAngle();
        backTheoryAvg += theory.interpolateIntensity(theta);
        i--;
        count++;
      }
      backDataAvg /= config.frontBackPoints;  // Convert to averages
      backTheoryAvg /= config.frontBackPoints;

      double difference = (frontDataAvg/backDataAvg 
                           - frontTheoryAvg/backTheoryAvg);
      result = difference*difference;
      break;
  }

  return result;
}
// END Mathematical Support Functions *****************************************


// Support Functions ********************************************************** 
void setDefaults(configuration& config) {
  // Input:  A configuration.
  // Output: Sets the values in that configuration to the defaults.
  
  config.polarization    = UNPOLARIZED;
  config.fitMode         = ALPHA;

  config.startFit        = 0;
  config.endFit          = 180;
  config.frontBackPoints = 1;
  config.pointToDivide   = 45;

  config.dumpChi         = false;
  config.dumpData        = false;
  config.dumpTheory      = false;
  config.dumpTheoryX     = 1;
}

void processArguments(int argc, char* argv[], configuration& config) {
  // Input:  The argument count and vector from the main program, and a
  //         configuration to parse that information into.
  // Output: Looks at the command line arguments one by one and sets up
  //         configuration accordingly.

  if (argc < 3) {                             // If the user didn't at least
    printUsage(EXTENDED);                     //   give filenames, print
    exit(1);                                  //   extended help.
  }

  config.dataFilename   = argv[1];            // Read filenames first.
  config.theoryFilename = argv[2];

  int i = 3;
  while(i < argc) {
    char* argument;

    argument = argv[i];      
    if (argument[0] != '-') {                 // Switch doesn't begin with '-'
      cerr << "Error: Invalid switch '" << argument << "'." << endl;
      printUsage(SIMPLE);                     // Print command line syntax
      exit(1);
    }

    switch(argument[1]) {                     // Check first character of
      case 'p':                               // Switch is actually -pol
        if ((argc <= i + 1) || (argv[i+1][0] == '-')) {
          cerr << "Error: Switch -pol requires a value." << endl;
          printUsage(SIMPLE);                 // No value included with switch
          exit(1);                            //   so print error message.
        }
        config.polarization = atoi(argv[++i]);
        break;
      case 'n':                               // Switch is -nouncert
        config.fitMode = ALPHANOUNCERT;
        break;
      case 'd':                               // Switch is -dividebypoint or
                                              //   -dump, so check which one.

        if ((argc <= i + 1) || (argv[i+1][0] == '-')) {
          cerr << "Error: Switch " << argument << " requires a value." << endl;
          printUsage(SIMPLE);                 // No value included with switch
          exit(1);                            //   so print error message.
        }

        if (argument[2] == 'i') {             // This is -dividebypoint
          config.fitMode = DIVIDEBYPOINT;
          config.pointToDivide = atof(argv[++i]);
        }
        else {                                // This is -dump
          if (argv[i+1][0] == 'd') {          // Value is 'data', so dump data
            config.dumpData = true;
          }
          else if (argv[i+1][0] == 'c') {     // Value is 'chi', so dump chi
            config.dumpChi = true;
          }
          else {                              // Dumping theory at x=value
            config.dumpTheory = true;
            config.dumpTheoryX = atof(argv[i+1]);
          }
          i++;
        }
        break;
      case 'f':                               // Switch is -frontback
        if ((argc <= i + 1) || (argv[i+1][0] == '-')) {
          cerr << "Error: Switch -frontback requires a value." << endl;
          printUsage(SIMPLE);                 // No value included with switch
          exit(1);                            //   so print error message.
        }
        config.fitMode = FRONTBACK;
        config.frontBackPoints = atoi(argv[++i]);
        break;
      case 'r':                               // Switch is -range
        if ((argc <= i + 2) || (argv[i+1][0] == '-') ||
            (argv[i+2][0] == '-')) {
          cerr << "Error: Switch -range requires two values." << endl;
          printUsage(SIMPLE);                 // No value included with switch
          exit(1);                            //   so print error message.
        }
        config.startFit = atof(argv[++i]);
        config.endFit   = atof(argv[++i]);
        break;
      case 'h':                               // In either case, help is being
      case '?':                               //   requested, so print extended
        printUsage(EXTENDED);                 //   help message and exit.
        exit(1);
        break;
      default:                                // Unknown switch
        cerr << "Error: Unknown switch '" << argument << "'." << endl;
        printUsage(SIMPLE);
        exit(1);
        break;
    }
    i += 1;                                   // Increment i to the next switch
  }
}

void checkConfiguration(configuration config) {
  // Input:  A configuration.
  // Output: Checks to see if that configuration is valid, and prints error
  //         messages if necessary.

  if ((config.fitMode < ALPHA) || (config.fitMode > FRONTBACK)) {
    cerr << "Error: Invalid fit mode specified." << endl;
    exit(1);
  }
  if (config.frontBackPoints < 1) {
    cerr << "Error: Must specify at least 1 point for a front-back fit." << endl;
    exit(1);
  }

  if ((config.startFit > config.endFit) || (config.endFit < 0)) {
    cerr << "Error: Invalid or reversed range, [" << config.startFit;
    cerr << " .. " << config.endFit << "]." << endl;
    exit(1);
  } else {                                    // Open files and check ranges
    ifstream dataFile(config.dataFilename, ios::in | ios::nocreate,
                      filebuf::sh_read);      // Open data file only if it
                                              //   already exists (nocreate)
                                              //   but allow others to read
                                              //   from the file (sh_read).
    dataSet data;
    if (dataFile.good() && data.loadData(dataFile)) {
      int i = 0;
      while(data.getPoint(i)->getAngle() < config.startFit) {
        i++;                                  // Find start of fit range
        if(i >= data.getSize()) {
          cerr << "Error: Encountered end of data before reaching beginning";
          cerr << " of specified fitting range." << endl;
          exit(1);
        }
      }

      if (data.getPoint(i)->getAngle() > config.endFit) { // See if there is data in range
        cerr << "Error: No data points fall within the specified fitting";
        cerr << " range." << endl;
        exit(1);
      }

      int count = 0;                          // Count how many points are in
      while(inFitRange(i, data, config)) {    //   the fitting range.
        i++;
        count++;
      }

      if (count < config.frontBackPoints) {   // Are there enough data points?
        cerr << "Error: Attempted front-back fit with ";
        cerr << config.frontBackPoints << " points in a fitting range with ";
        cerr << " only " << count << " points." << endl;
        exit(1);
      }

      i = 0;                                  // Return to start of fitting
      while(!inFitRange(i, data, config)) {   //   range and look for the
        i++;                                  //   point to divide by in a
      }                                       //   DIVIDEBYPOINT fit. If it
                                              //   isn't there (or isn't in
      bool foundDivideByPoint = false;        //   range), output an error.
      while(inFitRange(i, data, config)) {
        if (data.getPoint(i)->getAngle() == config.pointToDivide) {
          foundDivideByPoint = true;
        }
        i++;
      }

      if (!foundDivideByPoint && (config.fitMode == DIVIDEBYPOINT)) {
        cerr << "Error: Cannot divide by point at angle ";
        cerr << config.pointToDivide << " because that angle is not present";
        cerr << endl;
        cerr << "in the fitting range. (Use -dump data to examine data).";
        cerr << endl;
        exit(1);
      }
    }
    dataFile.close();
  }

  if ((config.polarization < 1) || (config.polarization > 3)) {
    cerr << "Error: Invalid polarization " << config.polarization << ".";
    cerr << " Use 1 = perp, 2 = par, 3 = unp." << endl;
    printUsage(SIMPLE);                       // Print short error message
    exit(1);
  }
  if ((config.dumpChi && config.dumpTheory) ||
      (config.dumpChi && config.dumpData)   ||
      (config.dumpData && config.dumpTheory)) {
    cerr << "Error: Multiple -dump switches specified." << endl;
    exit(1);
  }
}

void openFiles(char* theoryFilename, ifstream& theoryFile,
               char* dataFilename, ifstream& dataFile) {
  // Input:  Data and theory filenames, and two ifstreams (passed by reference)
  //         into which those files are to be opened.
  // Output: Opens theoryFilename (if it exists) and attaches it to theoryFile.
  //         Similarly attaches dataFilename to dataFile.

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

  theoryFile.open(theoryFilename, ios::in | ios::nocreate,
                  filebuf::sh_read);          // Open theory file
  if (!theoryFile.good()) {                   // File doesn't exist
    cerr << "Error: Unable to read theory file '" << theoryFilename;
    cerr << "'." << endl;
    exit(1);
  }
}

void outputScaledTheory(dataSet& theory, double alpha) {
  // Input:  A theory dataSet to be output with scaling given by alpha.
  // Output: Writes the theory set to cout with scaling alpha.

  for(int i = 0; i < theory.getSize(); i++) {
    cout << theory.getPoint(i)->getAngle() << "\t";
    cout << theory.getPoint(i)->getValue() * alpha << endl;
  }
}

void printResults(double sizeParameter, double chiSquare, double alpha) {
  // Input:  A size parameter, chiSquare and alpha corresponding to the best
  //         fit according to the method given.
  // Output: Writes the results of the fit to cout.

  cout << "Best fit is x = " << sizeParameter << " with chi^2 = ";
  cout << chiSquare << "." << endl;
  cout << "Theory scaling factor is alpha = " << alpha << "." << endl;
}
  
void printUsage(bool extended) {
  // Input:  A boolean, indicating whether or not extended help is to be
  //         printed. True indicates that extended help is to be printed.
  // Output: Prints the appropriate error message.
  
  cerr << "Usage: fitmie datafile theoryfile [-pol {1,2,3}] [-nouncert]" << endl;
  cerr << " [-dividebypoint double] [-frontback int] [-dump {data, chi, double}]" << endl;
  cerr << " [-range double double] [-help | -?]" << endl;

  if (extended) {
    cerr << endl;
    cerr << "datafile                  - Experimental data in angle-intensity-sigma format." << endl;
    cerr << "theoryfile                - Theoretical Mie tables in MIEV0 format." << endl;
    cerr << "[-pol {1,2,3}]            - Polarization. 1 = perp, 2 = par, 3 = unp (default)" << endl;
    cerr << "[-nouncert]               - Fit by alpha regression without uncertainties." << endl;
    cerr << "[-dividebypoint double]   - Fit by Bohren and Huffman's divide by point method," << endl;
    cerr << "                            using the specified angle." << endl;
    cerr << "[-frontback int]          - Fit by back-to-front ratio using averages of given" << endl;
    cerr << "                            number of points from front and back of data." << endl;
    cerr << "[-dump {data,chi,double}] - Dump input data, chi-square, or the fit curve" << endl;
    cerr << "                            at the size parameter nearest to the double given." << endl;
    cerr << "[-range double double]    - Fit only data in the range from the first angle" << endl;
    cerr << "                            specified to the second angle specified." << endl;
    cerr << "[-help | -?]              - Print this message." << endl;
  }
}