// Filename: snell.cpp
// Author:   Mike Rust and Ian Weiner
// Modified: Mike Schubmehl
// Date(s):  July 1999, 07/07/2000
// Purpose:    Corrects theoretical Mie tables for the Snell's law and
//           foreshortening effects present in an experimental setup.
//
//            This program reads in a theoretical table of scattering
//           efficiency as a function of angle in the format produced by
//           the MIEV0 program. Each size parameter is flagged by an
//           asterisk. Angle is printed in the first column, followed by
//           perpendicular-polarized theory, parallel-polarized theory,
//           and finally unpolarized theory. All fields are seperated by tabs.
//           For example,
//
//                 *	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
//
//             After reading in one polarization from the file, specified on
//           the command line, the program performs a correction for Snell's
//           law and foreshortening effects present in the experimental setup.
//           It outputs these new theory values to the output file in MIEV0
//           format, with the two unused columns (polarizations) zeroed out.
//             The command line syntax for the program is
//
// snell inputfile outputfile [-R int] [-L int] [-n double] [-pol {1,2,3}]
//       [-w double] [-N int] [-max double] [-front {1,0}] [-side {1,0}] [-help | -?]
// 
// inputfile      - Uncorrected Mie tables in MIEV0 format.
// outputfile     - File to output corrected Mie tables to.
// [-R int]       - Radius from detector to cuvette in mm (default = 245).
// [-L int]       - Half the length of cuvette side in mm (default = 6).
// [-n double]    - Relative index of refraction (e.g., water to air = 1.33
//                  when using particles suspended in water) (default = 1.33).
// [-pol {1,2,3}] - Polarization. 1 = perp, 2 = par, 3 = unp (default).
// [-w double]    - Angular resolution of detector, in degrees (default = 0.3).
// [-N int]       - Number of scatterers to use in correction (default = 20).
// [-max double]  - Maximum angle to output corrected data at (default = 70.0).
// [-front {1,0}] - Use 0 to disable front correction, 1 to enable (default).
// [-side {1,0}]  - Use 0 to disable side correction, 1 to enable (default).
// [-help | -?]   - Print this message.
//
// Notes:      Requires stdlib.h, iostream.h, fstream.h, iomanip.h, string.h,
//           math.h, and ctype.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 <math.h>
#include <ctype.h>
#include "dataset.hh"                         // dataSets used to hold data

#define Pi            3.1415926535897932385   // Numerical approximation of Pi
#define DEG           0.0174532925199432958   // Multiply by 180/Pi to convert
                                              //   degrees to radians.
const double rootTolerance = 0.000001;        // Tolerance for finding roots of
                                              //   equations.
const double deltaTheta    = 0.000001;        // Small angle difference used to
                                              //   keep endpoints in bounds.
#define PERPENDICULAR 1                       // These polarization constants
#define PARALLEL      2                       //   should be exactly as they
#define UNPOLARIZED   3                       //   appear in dataset.cpp.

#define EXTENDED      true                    // Constants for use with
#define SIMPLE        false                   //   printUsage(bool extended).

struct configuration {                        // Parameters for the correction
  char*  inputFilename;                       // Filenames for input and output
  char*  outputFilename;                      //   files to be opened.

  int    polarization;                        // Polarization of incident light
  int    R;                                   // Photodetector arm radius
  int    L;                                   // Half side length of cuvette
  double nr;                                  // Relative index of refraction;
                                              //   water to air = 1.33 for an
                                              //   experiment using particles
                                              //   suspended in water.
  double omega;                               // Angular resolution of detector
  int    numScatterers;                       // Number of scatterers used to
                                              //   approximate infinitely many.
  double maxOutputAngle;                      // Maximum angle to output
  bool   frontFace;                           // True if front and side faces
  bool   sideFace;                            //   are to be included.
};


// Table of Contents **********************************************************
// -- Mathematical Support Functions --
double Sin(double x);                         // Returns the sin of an angle x,
                                              //   specified in degrees.
double Cos(double x);
double Tan(double x);
int    sign(double x);                        // Returns the sign of x

double scatterAngleFront(double h,              // Find angle in original table
                         double detectorAngle,  //   corresponding to detector
                         configuration config); //   angle and scatterer
                                                //   position, using the values
                                                //   in config. FRONT face.
double scatterAngleSide(double h,               // Same, for SIDE face.
                        double detectorAngle,
                        configuration config);

void outputCorrectedTheory(dataSet& theory,     // Computes corrections and
                           configuration config,//   outputs the result to
                           ofstream& file);     //   file in MIEVO format.

// -- 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 printUsage(bool extended);               // Prints useage message to cerr,
                                              //   with detailed help if
                                              //   extended is set to true.
int  initProgressIndicator(configuration& config); // Find size of file and
                                                   //   create a status bar.
// END Table of Contents ******************************************************


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

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

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

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

  int starCount = initProgressIndicator(config);  // Get total # of stars

  int currentCount = 0;                       // Count the stars as they pass
  while(!inputFile.eof()) {
    dataSet theory;                           // dataSet for uncorrected theory
    double  sizeParameter;                    // Size parameter being computed

    while(!inputFile.eof() && inputFile.peek() != '*') {
      inputFile.ignore();                     // Read file until end or a *
    }                                         //   is found, indicating the
    inputFile.ignore();                       //   size parameter is next.
                                              //   Then skip the * and read
    inputFile >> sizeParameter;               //   the size parameter.

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

    outputFile << "*\t" << sizeParameter << endl;
    outputCorrectedTheory(theory, config, outputFile);
    currentCount++;                           // Increment stars passed so far
    int percentCompleted = (int) (100*((float) currentCount/ (float) starCount));
    if(percentCompleted < 10) {
      cout << "0";
    }
    cout << percentCompleted << "%\b\b\b" << flush;  // Output % complete
  }
  cout << endl;

  return 0;
}
// END Main Program ***********************************************************


// Mathematical Support Functions *********************************************
double Sin(double x) {
  // Input:  An angle x, in degrees.
  // Output: The sin of that angle.

  return sin(x * DEG);
}

double Cos(double x) {
  // Input:  An angle x, in degrees.
  // Output: The cos of that angle.

  return cos(x * DEG);
}

double Tan(double x) {
  // Input:  An angle x, in degrees.
  // Output: The tan of that angle.

  return tan(x * DEG);
}

int sign(double x) {
  // Input:  A real number x.
  // Output: The sign of x (i.e., 1 if x > 0, 0 if x = 0, and -1 if x < 0).

  if (x < 0) {
    return -1;
  }
  else if (x > 0) {
    return 1;
  }
  else {
    return 0;
  }
}

double scatterAngleFront(double h, double detectorAngle,
                         configuration config) {
  // Input:  The distance of the scatterer from the front of the cuvette, h,
  //         and the angle at which the detector is sitting.
  // Output: Uses Snell's law to figure out what scattered angle this detector
  //         angle corresponds to. That is, figures out what angle to look up
  //         in the uncorrected Mie tables to determine what intensity we ought
  //         to be seeing at this detector angle through the FRONT face of the
  //         cuvette. The routine does this by numerically solving the equation
  //
  //                                          nr Sin[theta]
  // h Tan[theta] + (R Cos[phi] - L) ----------------------------- = R Sin[phi]
  //                                  Sqrt[1 - (nr Sin[theta])^2]
  //
  //         where theta is the angle the scattered ray makes with the normal,
  //         phi is the angle of the detector with respect to the beam, and nr
  //         is the relative index of refraction between the medium inside and
  //         outside the cuvette. This equation is discussed in greater detail
  //         in Weiner, Rust, and Donnelly's paper. The equation is solved here
  //         for theta using bisection.


  int    R  = config.R;                       // Copy values from config to
  int    L  = config.L;                       //   simplify formulas below.
  double nr = config.nr;

  double lower;                               // Endpoints and midpoint of the
  double upper;                               //   interval containing the root
  double middle;                              //   in the bisection algorithm.

  double lowerValue;                          // Values at the endpoints and
  double upperValue;                          //   midpoint of the interval.
  double middleValue;
  
  lower = 0;                                  // Initialize endpoints.
  upper = asin(1/nr)/DEG - deltaTheta;        // If we evaluate the function at
                                              //   higher angles, it becomes
                                              //   imaginary. Want real roots.

  lowerValue = h * Tan(lower) + (R * Cos(detectorAngle) - L) *
               (nr * Sin(lower))/sqrt(1 - pow(nr * Sin(lower), 2))
               - R * Sin(detectorAngle);      // Compute the difference between
                                              //   the right and left sides of
                                              //   the equation above, at the
                                              //   lower end of the interval.

  upperValue = h * Tan(upper) + (R * Cos(detectorAngle) - L) *
               (nr * Sin(upper))/sqrt(1 - pow(nr * Sin(upper), 2))
               - R * Sin(detectorAngle);      // Difference at upper endpoint

  while(fabs(lowerValue - upperValue) > rootTolerance) {
    if(sign(lowerValue) == sign(upperValue)) {
      cerr << "Error: Unable to find scatter angle corresponding to a";
      cerr << " detector angle of " << detectorAngle << "." << endl;
      exit(1);
    }

    middle = (upper + lower)/2;
    middleValue = h * Tan(middle) + (R * Cos(detectorAngle) - L) *
                 (nr * Sin(middle))/sqrt(1 - pow(nr * Sin(middle), 2))
                 - R * Sin(detectorAngle);    // Difference at midpoint

    if(fabs(middleValue) < rootTolerance) {   // Midpoint is an acceptably
      return middle;                          //   accurate solution.
    }
    else {                                    // Cut interval in half
      if(sign(middleValue) != sign(upperValue)) {
        lower = middle;                       // Use midpoint as the new lower
        lowerValue = middleValue;             //   endpoint, keeping the sign
      }                                       //   change in the new interval.
      else {                                  // Midpoint becomes the new
        upper = middle;                       //   upper endpoint.
        upperValue = middleValue;
      }
    }
  }

  return (upper + lower)/2;                    // This is where the function has
}                                              //   a zero.

double scatterAngleSide(double h, double detectorAngle,
                        configuration config) {
  // Input:  The distance of the scatterer from the front of the cuvette, h,
  //         and the angle at which the detector is sitting.
  // Output: Uses Snell's law to figure out what scattered angle this detector
  //         angle corresponds to. That is, figures out what angle to look up
  //         in the uncorrected Mie tables to determine what intensity we ought
  //         to be seeing at this detector angle through the SIDE face of the
  //         cuvette. The routine does this by numerically solving the equation
  //
  //                                          nr Cos[theta]
  // L/Tan[theta] + (R Sin[phi] - L) ----------------------------- = R Cos[phi] + L - h
  //                                  Sqrt[1 - (nr Cos[theta])^2]
  //
  //         where theta is the angle the scattered ray makes with the normal,
  //         phi is the angle of the detector with respect to the beam, and nr
  //         is the relative index of refraction between the medium inside and
  //         outside the cuvette. This equation is discussed in greater detail
  //         in Weiner, Rust, and Donnelly's paper. The equation is solved here
  //         for theta using bisection.


  int    R  = config.R;                       // Copy values from config to
  int    L  = config.L;                       //   simplify formulas below.
  double nr = config.nr;

  double lower;                               // Endpoints and midpoint of the
  double upper;                               //   interval containing the root
  double middle;                              //   in the bisection algorithm.

  double lowerValue;                          // Values at the endpoints and
  double upperValue;                          //   midpoint of the interval.
  double middleValue;
  
  lower = acos(1/nr)/DEG + deltaTheta;        // Initialize endpoints, being
  upper = 180 - lower - deltaTheta;           //   careful to stay above the
                                              //   angle below which it is
                                              //   imaginary. Want real roots.

  lowerValue = L/Tan(lower) + (R * Sin(detectorAngle) - L) *
               (nr * Cos(lower))/sqrt(1 - pow(nr * Cos(lower), 2))
               - R * Cos(detectorAngle) + h - L;

  upperValue = L/Tan(upper) + (R * Sin(detectorAngle) - L) *
               (nr * Cos(upper))/sqrt(1 - pow(nr * Cos(upper), 2))
               - R * Cos(detectorAngle) + h - L;

  while(fabs(lowerValue - upperValue) > rootTolerance) {
    if(sign(lowerValue) == sign(upperValue)) {
      cerr << "Error: Unable to find scatter angle corresponding to a";
      cerr << " detector angle of " << detectorAngle << "." << endl;
      exit(1);
    }

    middle = (upper + lower)/2;
    middleValue = L/Tan(middle) + (R * Sin(detectorAngle) - L) *
               (nr * Cos(middle))/sqrt(1 - pow(nr * Cos(middle), 2))
               - R * Cos(detectorAngle) + h - L;

    if(fabs(middleValue) < rootTolerance) {   // Midpoint is an acceptably
      return middle;                          //   accurate solution.
    }
    else {                                    // Cut interval in half
      if(sign(middleValue) != sign(upperValue)) {
        lower = middle;                       // Use midpoint as the new lower
        lowerValue = middleValue;             //   endpoint, keeping the sign
      }                                       //   change in the new interval.
      else {                                  // Midpoint becomes the new
        upper = middle;                       //   upper endpoint.
        upperValue = middleValue;
      }
    }
  }

  return (upper + lower)/2;                    // This is where the function has
}                                              //   a zero.

void outputCorrectedTheory(dataSet& theory, configuration config,
                           ofstream& file) {
  // Input:  An uncorrected dataSet called theory, a configuration, and
  //         an ofstream to output the results to.
  // Output: Computes the snell correction for the theory and outputs the
  //         result to the output stream in MIEVO format (see opening
  //         comments). The correction is made by summing the scattered
  //         intensity from a number (config.numScatterers) of individual,
  //         evenly spaced scatterers in the cuvette. This also corrects for
  //         foreshortening effects.
  //           At a given angle psi and scatterer position h, the routine
  //         uses some simple trigonometry to check whether the detector
  //         can see the front and side faces of the cuvette. If it can, and
  //         the user wishes to make the relevant correction, the contribution
  //         from each face is computed and added into the total theoretical
  //         signal at angle psi (result). Each contribution is computed by
  //         integrating the interpolated theoretical intensity over the
  //         angular resolution of the detector.

  double omega = config.omega;                // Copy values from config to
  double R     = config.R;                    //   simplify the formulas below.
  double L     = config.L;
  
  double h;                                   // Distance of scatterer from
                                              //   front of cuvette.
  double hStep = 2 * L/(double)(config.numScatterers + 1);

  double psi;                                 // Detector angle
  double psiStep = theory.getStepSize();      // Use steps the same size as
                                              //   those in the theory file,
                                              //   to avoid over-interpolating.
  
  double theta;                               // Scatter angle that produced
                                              //   the signal at angle psi.
  double theta_0;                             // Scatter angles that produced
  double theta_f;                             //   signal at psi +/- half the
                                              //   detector's resolution.
  double result;                              // Total corrected theory at psi
  double frontContribution;                   // Contribution from font face
  double sideContribution;                    // Contribution from side face

  for(psi = 1; psi <= config.maxOutputAngle; psi += psiStep) {
    result = 0;
    for(h = hStep; (int)(h/hStep) < config.numScatterers; h += hStep) {
                                              // Sum contributions from all
                                              //   scatterers in the cuvette.
      frontContribution = 0;                  // Is front face visible?
      if ((R * Cos(psi + (omega/2)) > L) && config.frontFace) {
        theta = scatterAngleFront(h, psi + (omega/2), config);                             /* changed from psi */
        double x = h * Tan(theta);                                                         /* changed from Tan(theta + omega/2) */
        if(fabs(x) < L) {                     // Can scatter even hit cuvette?
          theta_0 = scatterAngleFront(h, psi - (omega/2), config);
          theta_f = theta;
          frontContribution = theory.interpolateIntegral(theta_0, theta_f);
        }
      }
      result += frontContribution/config.numScatterers;
     
      sideContribution = 0;                   // Is side face visible?
      if ((R * Sin(psi - (omega/2)) > L) && config.sideFace) {
		    theta = scatterAngleSide(h, psi - (omega/2), config);
		    double y = h - L/Tan(theta);
        if((y > 0) && (y < L)) {              // Can scatter even hit cuvette?
 		      theta_0 = theta;
          theta_f = scatterAngleSide(h, psi + (omega/2), config);
          sideContribution = theory.interpolateIntegral(theta_0, theta_f);
        }
      }
      result += sideContribution/config.numScatterers;
    }

    switch(config.polarization) {             // Output in MIEV0 format, with
      case PERPENDICULAR:                     //   results in proper column.
        file << "\t" << (double)psi << "\t" << result << "\t" << "0" << "\t" << "0" << endl;
        break;
      case PARALLEL:
        file << "\t" << (double)psi << "\t" << "0" << "\t" << result << "\t" << "0" << endl;
        break;
      case UNPOLARIZED:
        file << "\t" << (double)psi << "\t" << "0" << "\t" << "0" << "\t" << result << endl;
        break;
    }
  }
}
// END Mathematical Support Functions *****************************************


// Support Functions **********************************************************
void setDefaults(configuration& config) {
  // Input:  A configuration.
  // Output: Sets that configuration to the default values given below, and in
  //         the printUsage message.

  config.R              = 245;
  config.L              = 6;
  config.nr             = 1.33;
  config.polarization   = 3;
  config.omega          = 0.3;
  config.numScatterers  = 20;
  config.maxOutputAngle = 70.0;
  config.frontFace      = true;
  config.sideFace       = true;
}

void processArguments(int argc, char* argv[], configuration& config) {
  // Input:  The argc (argument count) and argv (argument vector) from the
  //         main() program, and a configuration to parse them into.
  // Output: Checks the validity of the command line arguments and sets options
  //         in config accordingly.

  if (argc < 3) {                             // If the user didn't at least
    printUsage(EXTENDED);                     //   give filenames, print
    exit(1);                                  //   extended help.
  }
    
  config.inputFilename  = argv[1];            // Read filenames first.
  config.outputFilename = 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);
    }

    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.
    }

    char* value = argv[i + 1];                // Get value for argument
    switch(argument[1]) {                     // Check first character of
      case 'R':                               //   argument to determine what
        config.R = atoi(value);               //   switch was specified.
        break;
      case 'L':
        config.L = atoi(value);
        break;
      case 'n':
        config.nr = atof(value);
        break;
      case 'p':                               // Switch is actually -pol
       config.polarization = atoi(value);
        break;
      case 'w':
        config.omega = atof(value);
        break;
      case 'N':
        config.numScatterers = atoi(value);
        break;
      case 'm':                               // Switch is -max
        config.maxOutputAngle = atof(value);
        break;
      case 'f':                               // Switch is -front
        config.frontFace = (atoi(value) != 0);
        break;
      case 's':                               // Switch is -side
        config.sideFace = (atoi(value) != 0);
        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;
        break;
    }
    i += 2;                                   // Increment i to the next switch
  }                                           //   and argument pair.
}

void checkConfiguration(configuration config) {
  // Input:  A configuration.
  // Output: Checks all values of config to be sure they are valid, and outputs
  //         error messages if necessary.

  if (config.R <=0) {
    cerr << "Error: Value of R must be positive." << endl;
    exit(1);
  }
  if (config.L <= 0) {
    cerr << "Error: Value of L must be positive." << endl;
    exit(1);
  }
  if (config.nr <= 0) {
    cerr << "Error: Invalid relative index of refraction " << config.nr;
    cerr << "." << endl;
    exit(1);
  }
  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.omega <= 0) {
    cerr << "Error: Value of omega must be positive." << endl;
    exit(1);
  }
  if (config.numScatterers < 1) {
    cerr << "Error: Too few scatterers specified (" << config.numScatterers;
    cerr << ")." << endl;
    exit(1);
  }
  if (config.maxOutputAngle < 1) {
    cerr << "Warning: Maximum angle to be output is extremely";
    cerr << " low (" << config.maxOutputAngle << ")." << endl;
  }
  if (!config.frontFace && !config.sideFace) {
    cerr << "Error: Must specify front face and/or";
    cerr << " side face corrections." << endl;
    exit(1);
  }
}

void printUsage(bool extended) {
  // Input:  A boolean telling whether to print extended help.
  // Output: Prints useage to cerr.

  cerr << "Usage: snell inputfile outputfile [-R int] [-L int] [-n double] [-pol {1,2,3}]" << endl;;
  cerr << " [-w double] [-N int] [-max double] [-front {1,0}] [-side {1,0}] [-help | -?]" << endl;

  if (extended) {
    cerr << endl;
    cerr << "inputfile      - Uncorrected Mie tables in MIEV0 format." << endl;
    cerr << "outputfile     - File to output corrected Mie tables to." << endl;
    cerr << "[-R int]       - Radius from detector to cuvette in mm (default = 245)." << endl;
    cerr << "[-L int]       - Half the length of cuvette side in mm (default = 6)." << endl;
    cerr << "[-n double]    - Relative index of refraction (e.g., water to air = 1.33" << endl;
    cerr << "                 when using particles suspended in water) (default = 1.33)." << endl;
    cerr << "[-pol {1,2,3}] - Polarization of light. 1 = perp, 2 = par, 3 = unp (default)." << endl;
    cerr << "[-w double]    - Angular resolution of detector, in degrees (default = 0.3)." << endl;
    cerr << "[-N int]       - Number of scatterers to use in correction (default = 20)." << endl;
    cerr << "[-max double]  - Maximum angle to output corrected data for (default = 70.0)." << endl;
    cerr << "[-front {1,0}] - Use 0 to disable front face correction, 1 to enable (default)." << endl;
    cerr << "[-side {1,0}]  - Use 0 to disable side face correction, 1 to enable (default)." << endl;
    cerr << "[-help | -?]   - Print this message." << endl;
  }
}

int initProgressIndicator(configuration& config) {
  // Input:  None.
  // Output: Returns the number of *'s in config.inputFilename, outputs a 00%,
  //         then repositions the cursor at the beginning of the 00%.

  ifstream inputFile(config.inputFilename);   // Open file specified by
                                              //   config's inputFilename,
                                              //   which contains MIEVO theory.

  cout << "Performing correction on file '" << config.inputFilename;
  cout << "': " << flush;
  
  char ch;                                    // Character read from inputFile
  int count = 0;                              // Number of *'s in the file
  while(inputFile.good()) {                   // As long as there is input,
    inputFile.get(ch);                        //   get a character. If it is a
    if (ch == '*') {                          //   *, increment count.
      count++;
    }
  }
  cout << "00%" << flush;                     // Output 00% to screen and then
  cout << "\b\b\b" << flush;                  //   move back three spaces.

  inputFile.close();                          // Done with file, so close it
  return count;                               // Return number of *'s.
}
// END Support Functions ******************************************************