/////////////////////// stable_current.cpp ///////////////////////
// Purpose: Stable Current Finder
// Author: Hassan Ibrahim, 2005
// Compile: g++ -o stable_current stable_current.cpp
// Use: ./stable_current 
//  or  ./stable_current [kinematics input_file output_file]

# include <cstdlib>
# include <cmath>
# include <string>
# include <iostream>
# include <fstream>
# include <iomanip>

using namespace std;

double bcmcur(string k, double u, double d, double t); // function prototype

const double curerr = 5.0;
const double edgetime = 20.0;
const double goodcur = 5.0;
const double goodtime = 60;

const double uconst[2] = {4092.4, 4101.6};
const double dconst[2] = {4188.5, 4165.65};
const double uoff[2] = {163.8, 163.8};
const double doff[2] = {110.0, 110.0};

const int nscal = 13; 

main(int argc, char* argv[])
{
  double ref[nscal];
  double tmp[nscal];
  double low[nscal];
  double pre[nscal];
  double high[nscal];

  double refcur;
  double tmpcur;
  double lowcur;
  double ctime;
  double cbcmu3;
  double cbcmd3;
  double avecur;
  double avechar;

  double ct1;
  double ct3;
  double ct5;
  double ce1;
  double ce3;
  double ce5;

  int n = 0;
  int s = 0;

  bool lowset = false;
  bool highset = false;

  string dummy;

  char* kin = new char [20];        // kinematics
  char* infile  = new char [20];    // input filename string variable
  char* outfile = new char [20];    // output filename string variable

  fstream fin;     // define fin as a file input object
  fstream fout;    // define fout as a file output object

  cout << '\n';
  cout << " -----------------------\n";
  cout << "| Stable Current Finder |\n";
  cout << " -----------------------\n";
  cout << '\n';

  // ------------- Inputs -------------

  if (argc == 4)
    {
      kin = argv[1];
      infile = argv[2];
      outfile = argv[3];
    }
  else
    {
      cout << "Kinematics (q1,q2,q3) : "; cin >> kin;
      cout << "Input Filename        : "; cin >> infile;
      cout << "Output Filename       : "; cin >> outfile;
      cout << '\n';
    }

  //------------- Calculations ------------------

  fin.open(infile, ios_base::in); // open the input file

  if ( !fin )
    {
      cerr << ">>> Error opening input file!\n\n";
      exit(1);
    }

  getline(fin, dummy);   // event, ctime, time, cbcmu3, bcmu3, cbcmd3, bcmd3

  for (int i=0; i<nscal; i++) fin >> ref[i];

  refcur = bcmcur(kin, ref[4], ref[6], ref[2]);

  for (int i=0; i<nscal; i++) fin >> tmp[i];

  fout.open(outfile, ios_base::out);       // open the output file

  fout.precision(10);
  fout.setf(ios_base::left);

  fout << setw(24) << "From (Event) To"  << setw(24) << "From (Time/s) To" << setw(12) << "Q (uC)" << setw(12) << "I (uA)" << setw(12) << "T1" << setw(12) << "T3" << setw(12) << "T5" << setw(12) << "E1" << setw(12) << "E3" << setw(12) << "E5" << '\n';

  while (!fin.eof())
    {

      n++;

      tmpcur = bcmcur(kin, tmp[4], tmp[6], tmp[2]);

      if (!lowset)
        {
          if (tmpcur > goodcur)
            {
              if (abs(refcur - tmpcur) < curerr)
                {
                  if ((tmp[1] - ref[1]) / 1024 > edgetime)
                    {
                      for (int i=0; i<nscal; i++) low[i] = tmp[i];
                      for (int i=0; i<nscal; i++) pre[i] = tmp[i];
                      lowcur = tmpcur;
                      lowset = true;
                    }
                }
              else
                {
                  for (int i=0; i<nscal; i++) ref[i] = tmp[i];
                  refcur = bcmcur(kin, ref[4], ref[6], ref[2]);
                }
            }
        }
      else
        {
          if (abs(lowcur - tmpcur) > curerr)
            {
              for (int i=0; i<nscal; i++) high[i] = pre[i];
              highset = true;
            }
          else
            {
              for (int i=0; i<nscal; i++) pre[i] = tmp[i];
            }
        }

      for (int i=0; i<nscal; i++) fin >> tmp[i];

      if (lowset && (highset || fin.eof()))
        {
          if (fin.eof())
            {
              for (int i=0; i<nscal; i++) high[i] = pre[i];    
            }
          ctime = high[1] - low[1];
          if (ctime / 1024 > goodtime)
            {
              s++;

              cbcmu3 = high[3] - low[3];
              cbcmd3 = high[5] - low[5];
              avecur = bcmcur(kin, cbcmu3, cbcmd3, ctime);
              avechar = avecur * ctime / 1024;

              ct1 = high[7] - low[7];
              ct3 = high[8] - low[8];
              ct5 = high[9] - low[9];
              ce1 = high[10] - low[10];
              ce3 = high[11] - low[11];
              ce5 = high[12] - low[12];

              fout << setw(12) << low[0] << setw(12) << high[0] << setw(12) << low[1] / 1024 << setw(12) << high[1] / 1024 << setw(12) << avechar << setw(12) << avecur << setw(12) << ct1 << setw(12) << ct3 << setw(12) << ct5 << setw(12) << ce1 << setw(12) << ce3 << setw(12) << ce5 << '\n';
            }
          lowset = false;
          highset = false;
        }
    }

  fin.close();     // close the input file
  fout.close();    // close the output file

  cout << "The input file " << infile << " includes " << n << " lines.\n";
  cout << "Number of stable current periods is " << s << ".\n";
  cout << "The output file " << outfile << " was saved.\n\n";     
}

// Function bcmcur

double bcmcur(string k, double u, double d, double t)
{

  double c;
  double uc;
  double dc;
  int m;

  m = (k == "q3") ? 1 : 0;

  uc = (u / t * 1024 - uoff[m]) / uconst[m];
  dc = (d / t * 1024 - doff[m]) / dconst[m];

  c = ( uc + dc ) / 2.0;

  return c;

}