// Filename: dataset.cpp
// Author:   Mike Rust and Ian Weiner
// Modified: Michael P. Schubmehl
// Date(s):  07/03/2000
// Purpose:    Implements dataset.hh, the interface to the dataSet class used
//           in the Mie scattering programs snell and fitmie.
//
// Notes:      Implements dataset.hh, and requires iostream.h  and fstream.h
//           for stream-handling routines. Uses stdlib.h for the NULL constant,
//           and ctype.h for the isdigit() function.

#include <stdlib.h>
#include <iostream.h>
#include <fstream.h>
#include <ctype.h>                            // For isdigit() function
#include "polynomial.hh"                      // Polynomials for interpolation
#include "list.hh"                            // Templated List class
#include "dataset.hh"                         // Header file for this class

#define PERPENDICULAR 1                       // Define polarizations
#define PARALLEL      2
#define UNPOLARIZED   3

// Table of Contents **********************************************************
//
// dataPoint:
// -- public --
// dataPoint();                               // Constructors
// dataPoint(double angle_, double value_,
//           double sigma_);
//
// dataPoint& operator=(const dataPoint& source); // Assignment operator
//
// void set(double angle_, double value_,     // Sets the dataPoint to the
//          double sigma_);                   //   values provided.
// void set(double angle_, double value_);
//
// double getAngle();                         // Accessor functions
// double getValue();
// double getSigma();
//
// dataSet: 
// -- public --
// dataSet();                                 // Default constructor
// ~dataSet();                                // Destructor
//
// void setPolarization(int polarization_);   // Set dataSet's polarization
// void setMaxOutputAngle(double maxOutputAngle_); // Set max angle to output
//
// double stepSize();                         // Accessor function for stepSize
//
// bool loadData(ifstream &file);             // Load experimental data file
// bool loadTheory(ifstream &file,            // Load the theoretical file with
//                int polarization);          //   specified polarization.
//
//
// double getMinAngle();                      // Return minimum and maximum
// double getMaxAngle();                      //   angle in this set.
// int getSize();                             // Returns size of this set
// dataPoint* getPoint(int i);                // Returns a pointer to the ith
//                                            //   dataPoint in the set.
// double interpolateIntensity(double theta); // Interpolates intensity using
//	                                          //   Lagrange's formula (cubic).
// double interpolateIntegral(double theta0,  // Integrates the interpolating
//		                       double theta1);  //   polynomial between points.
//
// other:
// ostream& operator<<(ostream& stream,       // Used to output dataSets
//                     const dataSet& source)
// END Table of Contents ******************************************************


// class dataPoint ************************************************************
dataPoint::dataPoint() {
  // Input:  None.
  // Output: Creates a new, blank data point.

  angle = 0;
  value = 0;
  sigma = 0;
}

dataPoint::dataPoint(double angle_, double value_, double sigma_) {
  // Input:  Values for angle, value, and sigma.
  // Output: Creates a new dataPoint with those values.

  angle = angle_;
  value = value_;
  sigma = sigma_;
}

void dataPoint::set(double angle_, double value_, double sigma_) {
  // Input:  Values for angle, value, and sigma.
  // Output: Sets the values of this dataPoint to those provided.

  angle = angle_;
  value = value_;
  sigma = sigma_;
}

void dataPoint::set(double angle_, double value_) {
  // Input:  Values for angle and value.
  // Output: Sets the values of this dataPoint to those provided.

  set(angle_, value_, 0);
}

double dataPoint::getAngle() {
  // Input:  None.
  // Output: Returns the angle of this dataPoint.

  return angle;
}

double dataPoint::getValue() {
  // Input:  None.
  // Output: Returns the value of this dataPoint.

  return value;
}

double dataPoint::getSigma() {
  // Input:  None.
  // Output: Returns the sigma of this dataPoint.

  return sigma;
}

dataPoint& dataPoint::operator=(const dataPoint& source) {
  // Input:  An source dataPoint.
  // Output: Sets the values of this data point equal to those of the old one.

  angle = source.angle;
  value = source.value;
  sigma = source.sigma;

  return *this;
}
// END class dataPoint ********************************************************


// class dataSet **************************************************************
dataSet::dataSet() {
  // Input:  None.
  // Output: Creates a new, empty dataSet.

  size           = 0;                         // Initially contains no data
  stepSize       = 1;                         // Most common angle step size
  minAngle       = 0;                         // Minimum and maximum angle
  maxAngle       = 180;                       //   present in this dataSet.
  maxOutputAngle = 70;                        // Output commonly used angles
  polarization   = UNPOLARIZED;               // Polarization of theory

  data           = NULL;                      // Allocate array on load
  interpolatingPolynomials = NULL;            // Allocate for polynomials
}

dataSet::~dataSet() {
  // Input:  None.
  // Output: Deletes the dataSet.

  delete[] data;                              // Delete the whole data array
  delete[] interpolatingPolynomials;          // Delete polynomial array
}

void dataSet::setPolarization(int polarization_) {
  // Input:  A polarization. 1 = perp, 2 = par, 3 = unp.
  // Output: Sets this file's polarization to that value.

  switch(polarization_) {
    case PERPENDICULAR:
      polarization = PERPENDICULAR;
      break;
    case PARALLEL:
      polarization = PARALLEL;
      break;
    case UNPOLARIZED:
      polarization = UNPOLARIZED;
      break;
    default:
      cerr << "Error: Invalid polarization specified in setPolarization.";
      cerr << endl;
      exit(1);
      break;
  }
}

void dataSet::setMaxOutputAngle(double maxOutputAngle_) {
  // Input:  A value for the maximum angle to output.
  // Output: Sets the value of maxOutputAngle as specified.

  maxOutputAngle = maxOutputAngle_;
}

double dataSet::getStepSize() {
  // Input:  None.
  // Output: The value of stepSize for this dataSet.

  return stepSize;
}

bool dataSet::loadData(ifstream &file) {
  // Input:  An input stream containing the experimental data to be read.
  // Output: Reads data form file, setting the value of size and allocating
  //         an array to hold all of the dataPoints form the file. Extends
  //         the array by 1 for each point read.

  double angle;                               // Local variables to store
  double value;                               //   values read from file.
  double sigma;
  List<dataPoint*> dataList;                  // Temporary list to hold data
                                              //   until an appropriate array
                                              //   size can be determined.
  size = 0;
  file >> angle >> value >> sigma;            // Read first data point
  while(file.good())  {                       // As long as there is more data
    size++;                                   //   in the file, read in angle,
                                              //   value, and sigma from file.
    dataPoint* temp = new dataPoint(angle, value, sigma);
    dataList.pushAtTail(temp);                // Store in list
    file >> angle >> value >> sigma;          // Read next point
  }

  delete[] data;                              // Delete any old data
  data = new dataPoint[size];                 // Allocate space for data array

  int i = 0;
  while(!dataList.isEmpty() && i < size) {    // As long as there is data left
    dataPoint* temp = dataList.pop();         //   in the list and space in the
    data[i] = *temp;                          //   array, move data into the
    delete temp;                              //   array and delete from list.
    i++;
  }

  minAngle = data[0].angle;                   // Record values of minAngle and
  maxAngle = data[size-1].angle;              //   maxAngle.
  
  isExperimental = true;                      // This is experimental data
  interpolatingPolynomials = NULL;            // No interpolation allowed

  return (size != 0);                         // Operation succeeded only if
                                              //   some data was read.
}

bool dataSet::loadTheory(ifstream &file, int polarization_) {
  // Input:  An input stream and the polarization to read from that input
  //         stream. 1 = perp, 2 = par, 3 = unpolarized.
  // Output: Reads the specified column from the file, replacing the data
  //         array with a new array containing the information from the file.
  //         Also sets size appropriately. Note that this routine expects the
  //         * sentinel immediately preceeding the data to have been read, as
  //         it has been when the size parameter has already been recorded.

  double angle;                               // Local variables to store
  double value;                               //   values from file.
  double perp, par, unp;                      // Different polarizations read
                                              //   from file's three columns.
  List<dataPoint*> dataList;                  // Temporary list to hold the
                                              //   data while it is read from
                                              //   file and counted.
  bool done = false;                          // True if eof has been hit or
                                              //   if a * indicates end of set.

  polarization = polarization_;               // Set dataSet's polarization to
                                              //   that specified.
  size = 0;                                   // Initialize size
  while(!done) {
    if((!file.good()) || (file.peek() == '*')) {
      done = true;                            // See if data is finished
    }
    else {
      file >> angle >> perp;                  // Read angle and perpendicular,
      if(file.peek() == 35) {                 //   then check to see if there
        perp = 0;                             //   was an error (1.#INF000). If
        while(file.peek() != 10) {            //   so, skip to end of line.
          file.ignore();
        }
      }
      else {
        file >> par;                          // Read parallel.
        if(file.peek() == 35) {               // Check to see if there
          par = 0;                            //   was an error. If so, skip
          while(file.peek() != 10) {          //   to end of this line.
            file.ignore();
          }
        }
        else {  
          file >> unp;                        // Read unpolarized.
          if(file.peek() == 35) {                
            unp = 0;                      
            while(file.peek() != 10) {      
            file.ignore();
            }
          }
        }
      }

      size++;                                 // Increment number of points.

      switch(polarization)  {                 // Select appropriate value from
        case PERPENDICULAR:                   //   file, based on specified
          value = perp;                       //   polarization.
          break;
        case PARALLEL:
          value = par;
          break;
        case UNPOLARIZED:
          value = unp;
          break;
        default:
          cerr << "Error: Illegal polarization in call to loadTheory" << endl;
          exit(1);
          break;
      }
      
      dataPoint* temp = new dataPoint(angle, value, 0);
      dataList.pushAtTail(temp);              // Store in list

      while(file.good() && (!isdigit(file.peek())) && (file.peek()!='*'))
        file.ignore();                        // Skip newline before next read
    }
  }

  delete[] data;                              // Delete any old data
  data = new dataPoint[size];                 // Allocate space for data array

  int i = 0;
  while(!dataList.isEmpty() && i < size) {    // As long as there is data left
    dataPoint* temp = dataList.pop();         //   in the list and space in the
    data[i] = *temp;                          //   array, move data into the
    delete temp;                              //   array and delete from list.
    i++;
  }

  if (size > 1) {                             // Compute step size if possible
    stepSize = data[1].angle - data[0].angle;
    minAngle = data[0].angle;
    maxAngle = data[size-1].angle;
  }
  else {
    cerr << "Warning: Theory data set has size less than 2." << endl;
  }

  isExperimental = false;                     // This is theoretical data

  delete[] interpolatingPolynomials;          // Delete old polynomials
  interpolatingPolynomials = new polynomial*[size]; // Allocate poly array
  for(int j = 0; j < size; j++) {             // Initialize all polynomials
    interpolatingPolynomials[j] = NULL;       //   to NULL (not computed).
  }

  return (size != 0);                         // Operation succeeded only if
                                              //   some data was read.
}

double dataSet::getMinAngle() {
  // Input:  None.
  // Output: Returns the smallest angle in this dataSet.

  return minAngle;
}

double dataSet::getMaxAngle() {
  // Input:  None.
  // Output: Returns the largest angle in this dataSet.

  return maxAngle;
}

int dataSet::getSize() {
  // Input:  None.
  // Output: Returns the size of this dataSet.

  return size;
}

dataPoint* dataSet::getPoint(int i) {
  // Input:  An integer i, whose value is less than the size of this dataSet.
  // Output: Reutrns a pointer to the ith dataPoint in this set.

  if ((i < 0) || (i >= size)) {
    cerr << "Error: Cannot get dataPoint " << i << " in a set";
    cerr << " of size " << size << "." << endl;
    exit(1);
  }

  return &(data[i]);
}

double dataSet::interpolateIntensity(double theta) {
  // Input:  An angle at which to interpolate the intensity.
  // Output: The intensity at that angle, found by using a Lagrange
  //         interpolating cubic, stored in the interpolatingPolynomials
  //         array. If this polynomial has not already been computed,
  //         it is computed and then stored.

  int index;                                  // Index of next lower angle
                                              //   present in table.

  index = (int) ((theta - minAngle)/stepSize);    // Compute the index of the
  if (interpolatingPolynomials[index] == NULL) {  //   next lower angle in the
    constructInterpPoly(index);                   //   table and build an
  }                                               //   interpolating poly, if
                                                  //   needed.

  return interpolatingPolynomials[index]->evaluate(theta);
}

double dataSet::interpolateIntegral(double theta0, double theta1) {
  // Input:  An angle range over which to integrate the interpolating
  //         polynomial used at theta1.
  // Output: The intensity integrated over those angles. If the interpolating
  //         polynomial has not already been computed, it is computed and
  //         then stored.

  int index;                                  // Index of next lower angle
                                              //   present in table.

  index = (int) ((theta1 - minAngle)/stepSize);   // Compute the index of the
  if (interpolatingPolynomials[index] == NULL) {  //   next lower angle in the
    constructInterpPoly(index);                   //   table and build an
  }                                               //   interpolating poly, if
                                                  //   needed.

  return interpolatingPolynomials[index]->evaluateIntegral(theta0, theta1);
}

void dataSet::constructInterpPoly(int index) {
  // Input:  An index in the array for which to build a Lagrange interpolating
  //         cubic. Uses the values in the data array to build this polynomial.
  // Output: Uses Lagrange's formula for an interpolating cubic:
  //
  //                  (x-x2)(x-x3)(x-x4)              (x-x1)(x-x3)(x-x4)
  //  y[x] = y[x1] * --------------------- + y[x2] * --------------------- +
  //                 (x1-x2)(x1-x3)(x1-x4)           (x2-x1)(x2-x3)(x2-x4)
  //
  //                  (x-x1)(x-x2)(x-x4)              (x-x1)(x-x2)(x-x3)
  //         y[x3] * --------------------- + y[x4] * ---------------------
  //                 (x3-x1)(x3-x2)(x3-x4)           (x4-x1)(x4-x2)(x4-x3)
  //
  //  with x2 chosen as the index provided, and x1, x3, and x4 the surrounding
  //  points on either side. Stores the result (in a form with terms expanded
  //  and collected by powers of x) in the polynomial at index in the array
  //  interpolatingPolynomials.

  int index1 = index - 1;                     // The indices of the points
  int index2 = index;                         //   to use in generating the
  int index3 = index + 1;                     //   Lagrange polynomial.
  int index4 = index + 2;

  double x1 = index1 * stepSize + minAngle;   // The four points for the
  double x2 = index2 * stepSize + minAngle;   //   interpolating polynomial,
  double x3 = index3 * stepSize + minAngle;   //   as described above.
  double x4 = index4 * stepSize + minAngle;


  if (index1 < 0) {                           // Negative indices illegal
    if (-x1 < minAngle) {                     // See if we can use symmetry
      cerr << "Error: Can't interpolate at theta = " << x1 << endl;
      exit(1);
    }
    else {
      index1 = -index1;                       // Use the intensity at |x1|
    }
  }

  if (index3 >= size) {                       // At end of array => no data
    if ((size - 1) * stepSize + minAngle < 180) { // Can't use symmetry
      cerr << "Error: Not enough data to interpolate at theta = ";
      cerr << x2 << endl;
      exit(1);
    }
    else {                                    // Use data at 360 - x3,
                                              //   since data is symmetric.
      index3 = (int) ((360 - x3 - minAngle)/stepSize);
    }
  }

  if (index4 >= size) {                       // See if past end of array
    if ((size - 1) * stepSize + minAngle < 180) { // Can't use symmetry
      cerr << "Error: Not enough data to interpolate at theta = ";
      cerr << x2 << endl;
      exit(1);
    }
    else {                                    // Use data at 360 - x4
      index4 = (int) ((360 - x4 - minAngle)/stepSize);
    }
  }

  double y1 = data[index1].value;             // Get y[x1], as described above
  double y2 = data[index2].value;
  double y3 = data[index3].value;
  double y4 = data[index4].value;

  double a[4] = {0, 0, 0, 0};                 // Coefficients of polynomial,
                                              //   initialized to zero.
  double temp;

  temp = y1/((x1-x2) * (x1-x3) * (x1-x4));    // Compute the Lagrange
  a[0] += temp * -(x2 * x3 * x4);             //   coefficients, collected in
  a[1] += temp * (x2*x3 + x2*x4 + x3*x4);     //   terms with like power of x.
  a[2] += temp * (-x2 - x3 - x4);
  a[3] += temp;

  temp = y2/((x2-x1) * (x2-x3) * (x2-x4));
  a[0] += temp * -(x1 * x3 * x4);
  a[1] += temp * (x1*x3 + x1*x4 + x3*x4);
  a[2] += temp * (-x1 - x3 - x4);
  a[3] += temp;

  temp = y3/((x3-x1) * (x3-x2) * (x3-x4));
  a[0] += temp * -(x1 * x2 * x4);
  a[1] += temp * (x1*x2 + x1*x4 + x2*x4);
  a[2] += temp * (-x1 - x2 - x4);
  a[3] += temp;

  temp = y4/((x4-x1) * (x4-x2) * (x4-x3));
  a[0] += temp * -(x1 * x2 * x3);
  a[1] += temp * (x1*x2 + x1*x3 + x2*x3);
  a[2] += temp * (-x1 - x2 - x3);
  a[3] += temp;

  interpolatingPolynomials[index] = new polynomial(3, a); // Create polynomial
}
// END class dataSet **********************************************************

ostream& operator<<(ostream& stream, const dataSet& source) {
  // Input:  A dataSet to be written to the output stream.
  // Output: Writes the dataSet to that stream in the form
  //
  //         <tab> angle <tab> perp <tab> par <tab> unp
  //
  //         For all points in the dataSet less than maxOuputAngle.
  //         The two polarizations not relevant to this dataSet
  //         will be printed as zeroes.


  if(source.isExperimental) {                 // Output in experimental format
    for(int i = 0; i < source.size; i++) {
      stream << source.data[i].angle << "\t" << source.data[i].value;
      stream << "\t" << source.data[i].sigma << endl;
    }
  }
  else {                                      // Output in theoretical (MIEVO)
    int i = 0;
    while((i < source.size) && (source.data[i].angle <= source.maxOutputAngle)) {
      stream << "\t" << source.data[i].angle << "\t";
      switch(source.polarization) {
        case PERPENDICULAR:
          stream << source.data[i].value << "\t" << 0 << "\t" << 0 << endl;
          break;
        case PARALLEL:
          stream << 0 << "\t" << source.data[i].value << "\t" << 0 << endl;
          break;
        case UNPOLARIZED:
          stream << 0 << "\t" << 0 << "\t" << source.data[i].value << endl;
          break;
      }
      i++;
    }
  }

  return stream;
}