// RKint3test.cc
// This program reads in parameters for the integration of a set of 
// differential equations, and integrates the amplitude equations starting
// with conditions input from a file.  It then writes the harmonics to 
// to a set of files.  Each file has a name of the form 
//              filename-#.#.dat
// where the filename is input by the user 
// 
// Original author:  Yurii Il'inskii
// Current author :  Ronald Kumon
// Original filename:  RKMmain.cc
//
// Global variables:
// nm, g, hi, eps, tiny, **Sr, **Si, *V, statflag, stepfilename, statfilename
//
// Modifications:
// 05 Dec 1996: Modified by Ronald Kumon; realigned text and added comments;
//              included introductory comments.  Changed C I/O functions to 
//              C++ I/O functions and added I/O error checking.
// 08 Dec 1996: Changed the way V was allocated with "dvector" so indices ran
//              from 1 to 2*nm; previously they ran from -1 to 2*nm+2
// 02 Feb 1997: Changed the way is data is written to files.  The actual
//              process is contained in a new function "freqproc" which 
//              writes the harmonic amplitudes, norm, and norm in dB.
//              Also changed the way the output "root" filename is set
//              so that it is the same as the input "root" filename.
// 08 Feb 1997: Added section to read in initial Fourier transformed velocity
//              components from a file.  Also changed the way the output
//              "root" file is set so that it back the way it used to be
//              (i.e., input by the user).
// 22 Feb 1997: Changed program so that the dissipation factor and
//              step size can be input by the user.
// 26 Feb 1997: Changed so that the user has the option to output 
//              integration statistics to files.  General statistics are
//              are output to "statfilename" while specific step size 
//              values for each step are written to "stepfilename" from
//              within the integration routine code (odeint.c).  The
//              step size output is optional.
// 23 Mar 1997: Changed program to allow the minimum scale factor
//              "tiny" to be input by user and set as a global variable.
// 11 Apr 1997: Added DiagnosticsOn as a boolean variable to turn on
//              integration diagnostics.  
// 24 Apr 1997: Modified code to integrate an arbitrary differential
//              equation using a fixed step size routine.
// 28 Sep 1997: Modified code to plot spectrum via the C binding library
//              to the PGPLOT FORTRAN library.

#include <iostream.h> 
#include <iomanip.h>
#include <fstream.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include "Nrutil.h"  
#include "RKtest.h"  
#include "RKinttest.h"
  
// Global variables
int nm;                       // Number of harmonics
int nm2;                      // Total number of harmonic comp. (real and imag)
int DiagnosticsOn;            // Boolean variable for integration diagnostics
double A;                     // Dissipation factor
double D;                     // Dispersion factor
double kxbar;                 // Wavenumber*shock formation distance
double *V;                    // Array of harmonic amplitudes
int statflag;                 // Flag for integration statistics
char stepfilename[81];        // Step file name for var. step size routine
char statfilename[81];        // Full filename for integration statistics
float *nvec;                  // Harmonic number vector
float *normdB;                // Vector for plotting amplitude norms in dB
float *normdBold;             // Vector for erasing old normdB values
double hi;                    // Initial step size
double eps;                   // Accuracy required for each step
double tiny;                  // Minimum scale factor for integration

void main()  
{  
  const int STRLEN=81;        // String length for filenames
  char infilename[STRLEN];    // "Root" input filename
  char infilename2[STRLEN];
  char inputfile[STRLEN];     // Pathname + "Root" input filename 
  char outfilename[STRLEN];   // "Root" input filename
  char outputfile[STRLEN];    // Pathname + "Root" output filename 
  char datafile[STRLEN];      // Output file + suffix for position data
  double x1, x2;              // Starting and ending times for integration
  double dx;                  // Step size for integration (fixed)
  double totaltime=0.0;       // Total time for all integrations
  int nminput;                // Number of harmonics initially input
  int num;                    // Initial cond. harmonic number (temp.)
  double Vamp,VampdB;         // Initial cond. vel. amplitude data (temp.)
  double dVampdBdx;           // Initial cond. ampl. deriv. data (temp.)
  char positionstr[STRLEN];   // Position for data output
  double position;            // Position as number
  char suffix[STRLEN];        // Suffix for position data files  	  
                              // Print data from input file    
  const int PRECISION=12;     // Define output formatting parameters
  cout.precision(PRECISION);  // Set precision
  cout.setf(ios::scientific); // Set scientific notation output

  // User-input parameters
  cout << "Dimensionless dissipation factor (A)   : ";
  cin  >> A;
  cout << "Dimensionless dispersion factor (D)    : ";
  cin  >> D;
  cout << "Wavenumber*shock formation distance    : ";
  cin  >> kxbar;
  cout << "Step size (fixed)                      : ";
  cin  >> dx; 
  cout << "Total number of harmonics in simulation: ";
  cin  >> nm;
  cout << "Integration diagnostics [1=yes 0=no]   : ";
  cin  >> DiagnosticsOn;  
  cout << endl;
  cout << "Dissipation factor     = " << A     << endl;
  cout << "Dispersion factor      = " << D     << endl;
  cout << "k_0*bar{x}             = " << kxbar << endl;
  cout << "Step size              = " << dx    << endl;
  cout << "Integration diagnostics= " << DiagnosticsOn << endl;
  cout << endl;

  // Allocate and initialize variables
  nm2=nm*2;                // Number of harmonics (positive & negative)
                           // Note that positive harmonics have odd indices
                           // and negative, even indices
  nvec=vector(1,nm);       // Allocate memory dynamically using NR routine
  normdB=vector(1,nm);
  normdBold=vector(1,nm);
  V=dvector(1,nm2);        
  for(int n=1; n<=nm;n++){  
    nvec[n]=n;             // Initialize harmonic number vector
  }
  for(int n=1; n<=nm2;n++){  
    V[n]=0.0;              // Initialize harmonic amplitudes to zero
  }

  // Input initial conditions
  cout << "Input filename for initial conditions data (max 20 char): ";
  cin  >> setw(STRLEN) >> infilename; 
  const char inpathname[]="/home/kumon/crystal/nonlin/RK4/main/testBurg/";
  strcpy(inputfile,inpathname);   // Copy pathname to string
  strcat(inputfile,infilename);   // Concatenate filename to pathname
  strcat(inputfile,".dat");       // Concatenate suffix to pathname+filename
  cout << "Will open file: " << inputfile << endl << endl;
	
  ifstream InputFileInitCond(inputfile, ios::in);   // Open input file
  if (!InputFileInitCond) {
    cerr << "Input Data File could not be opened" << endl;
    exit(1);
  }
  cout  << "Initial conditions:  " << endl;
  InputFileInitCond >> nminput;
  if (nminput > nm) {
    cerr << "Number of harmonics of data exceeds number in simulation!" << endl;
    cerr << "Please increase the number and try again." << endl;
    exit(1);
  }
  for (int n=1;n<=nminput;n++) {
    InputFileInitCond >> num >> V[2*n-1] >> V[2*n] >> Vamp >> VampdB >> dVampdBdx;
  }
  InputFileInitCond.close();          // Close input file

  cout  << "First ten harmonics:  " << endl;
  for (int n=1;n<=10;n++) {
    cout << "V[" << setw(2) << 2*n-1  << "]= " << V[2*n-1]   << "  "
         << "V[" << setw(2) << 2*n    << "]= " << V[2*n] << endl;
  }
 
  // Input "root" filename for the set of output data files 
  // Note that the length of the pathname+filename+suffix must be
  // less than STRLEN-1 otherwise unpredictable results may occur.
  cout << endl << "Input \"root\" filename for output data: ";
  cin  >> setw(STRLEN) >> outfilename; 
  const char outpathname[STRLEN]="/home/kumon/crystal/nonlin/RK4/main/testBurg/";
  strcpy(outputfile,outpathname);          // Copy pathname to string
  strcat(outputfile,outfilename);          // Concatenate filename to pathname
  cout << "Output File Root= " << outputfile << endl << endl;
  
  // First, configure the general statistics file
  strcpy(statfilename,outputfile);         // Copy root filename to string
  strcat(statfilename,".stat");            // Concatenate suffix to filename
  cout << endl << "General Stat. File= " << statfilename << endl;
  ofstream StatFile;                       // Decl. of output file object
  StatFile.open(statfilename, ios::out);   // Open data file
  if (!StatFile) {                         // Error checking
    cerr << "File could not be opened" << endl;
    exit(1);      
   }
  StatFile << statfilename << endl
           << "Initial conditions  : " << infilename2 << endl
           << "Dissipation         : " << A           << endl
           << "Dispersion          : " << D           << endl
           << "k_0*bar{x}          : " << kxbar       << endl
           << "Step size           : " << dx          << endl
           << "Total harmonics     : " << nm          << endl;
  StatFile.close();
   
  // Input positions for data file output
  cout << endl;
  cout << endl;
  cout << "Input filename for positions for data output (max 20 char): ";
  cin  >> setw(STRLEN) >> infilename2;
  const char inpathname2[]="/home/kumon/crystal/nonlin/RK4/main/testBurg/";
  strcpy(inputfile,inpathname2);  // Copy pathname to string
  strcat(inputfile,infilename2);  // Concatenate filename to pathname
  strcat(inputfile,".dat");       // Concatenate suffix to pathname+filename
  cout << "Will open file: " << inputfile << endl;
  ifstream InputFilePositions(inputfile, ios::in);   // Open input file
  if (!InputFilePositions) {
    cerr << "Input Data File could not be opened" << endl;
    exit(1);
  }

  // Perform integrations and write resulting data to files
  // Data at initial position
  InputFilePositions >> positionstr;      // Read in first position
  createSuffix(positionstr,suffix);       // Create appropriate suffix
  position=atof(positionstr);             // Convert from string to double
  x1=position; x2=position;   
  strcpy(datafile,outputfile);            // Copy root name to string   
  strcat(datafile,suffix);                // Add datafile suffix
  freqproc(datafile, nm, V, normdB);      // Write data to file
  openGraph(nm);                          // Open graphing window

  // Data at subsequent positions
  while (InputFilePositions >> positionstr) {
    createSuffix(positionstr,suffix);      // Create appropriate suffix
    position=atof(positionstr);            // Convert from string to double
    x1=x2; x2=position;                    // Integrate to next position
    integrateFixed(outputfile, suffix, nm, nm2, x1, x2, dx, 
                   V, &totaltime, nvec, normdB, normdBold);
  }

  // Print closing comments and clean up
  cout << "\a" << endl;                    // Sound system alarm
  cout << "Total time for all integrations was approximately "
       << totaltime << " s. " << endl;
  do { cout << "\b"; } while ( (cin.get()) != '\n' );  // Empty input stream
  closeGraph();                            // Close graphing window
  InputFilePositions.close();              // Close position input file
  free_dvector(V,1,nm2);                   // Deallocate memory for V 
  free_vector(normdBold,1,nm);
  free_vector(normdB,1,nm);
  free_vector(nvec,1,nm);

}