// CalcMainB5.cc
// This programs reads in the harmonic amplitude data
// for a particular position and computes the spatial velocity waveforms.
// It then writes this waveform data to files with the same
// root filename as the harmonic amplitude data.  The pathnames for the 
// input and output data are specified within the program but the user
// may input the filenames while the code is running.
// 
// Original author:  Yurii Il'inskii
// Current author :  Ron Kumon
//
// Modifications:
// 08 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 program slightly to read files produced by RKint2
// 03 Feb 1997: Changed program so that the user needs only input the
//              root filename.  The program then appends the suffixes listed
//              in the "suffix" array.
// 27 Feb 1997: Updated to run with RKint3.  Now the positions for which
//              the waveforms are to be computed are input from a file
//              (with the filename of that file input by the user).
//              Also, I undid the previous change so that the user may 
//              input a separate filename root for the waveform output 
//              file set.
// 01 Mar 1997: Added code to write harmonic data out to a file for the 
//              harmonic propagation curves in a convenient form for
//              gnuplot to read in.
// 10 Apr 1997: Added code to generalize the way the user can write 
//              the harmonic propagation data out a file.  The user 
//              can now specify the starting harmonic, the skip interval
//              and the final harmonic.
// 16 Oct 1997: Added code to allow the user to set the maximum harmonic
//              number used to reconstruct the waveforms.  Note that the
//              any harmonics greater than the maximum number specified
//              will be set to zero and the total number of harmonics will
//              still be the number entered at the begin of the program.
//              Hence the bandwidth and sampling rate remain the same
//              as the full data set.
// 21 Oct 1998: Modified to produce a single waveform by executing a
//              simple inverse Fourier transform.

#include <iostream.h>
#include <iomanip.h>
#include <fstream.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include "Nrutil.h"
#include "Calc.h"

main()
{	
  int nm;                     // Number of harmonics
  int nm2;                    // Total number of harmonic comp. (real and imag)
  double *V;                  // Wavevector domain amplitude factors
  double *Vx,*Vxh,*Vy,*Vz;    // Spatial domain amplitude factors
  double v0;                  // Linear wave speed
  double Bxr,Bxi;             // Sum of surface components in x-dir
  double Byr,Byi;             // Sum of surface components in y-dir
  double Bzr,Bzi;             // Sum of surface components in z-dir
  double Zr[4],Zi[4];         // Real and imag parts of eigenvalues
  double FrMatrix[4][4][4];   // F Matrix values
  double FiMatrix[4][4][4];
  int nminput;                // Number of harmonics input
  int num;                    // Dummy integer for file input
  double *norm, *normdB;      // Dummy variables for file input
  double *dnormdBdx;          // Derivative of amplitudes in dB
  int yes=1;                  // Continuation flag
  double x0,x;                // Position coordinate
  double BxrN,BxiN,BxN;       // Surface x-dir comp (real, imag, norm)
  double ByrN,ByiN,ByN;       // Surface y-dir comp (real, imag, norm)
  double BzrN,BziN,BzN;       // Surface z-dir comp (real, imag, norm)
  const int STRLEN=81;        // String length for filenames
  char infilename[STRLEN];    // "Root" input filenames
  char infilename2[STRLEN];
  char infilename3[STRLEN];
  char inputfile[STRLEN];     // Pathname + "Root" input filename 
  char filename[STRLEN];      // Case name read from input file 
  char outfilename[STRLEN];   // "Root" input filename
  char outputfile[STRLEN];    // Pathname + "Root" output filename 
  char outputfile2[STRLEN];   
  char temp[STRLEN];          // Temporary string
  char positionstr[STRLEN];   // Position for data output
  char suffix[STRLEN];        // Suffix for position data files  
  double position;            // Position as number
  int nmout;                  // Number of harmonics to output for prop curves
  int nminit;                 // Initial harmonic number to output
  int nmoutskip;              // Skip interval for harmonic numbers to be output

  // Read in filename for input data
  // Note that the length of the pathname+filename+suffix must be
  // less than STRLEN-1 otherwise unpredictable results may occur.

  //   cout << "Input filename for parameter data (max 20 char): ";
  //   cin  >> setw(STRLEN) >> infilename; 
  //   const char inpathname[]="/home/kumon/crystal/nonlin/NLCrystal/data/";
  //   strcpy(inputfile,inpathname);    // Copy pathname to string
  //   strcat(inputfile,infilename);    // Concatenate filename to pathname
  //   strcat(inputfile,".out");        // Concatenate suffix to pathname+filename
  //   cout << "Will open file: " << inputfile << endl;
	
  //   ifstream InputFileParam(inputfile, ios::in);   // Open input file
  //   if (!InputFileParam) {
  //     cerr << "Input Data File could not be opened" << endl;
  //     exit(1);
  //   }

  //   InputFileParam >> filename;
  //   InputFileParam >> nm;
  //   InputFileParam >> v0;
  //   InputFileParam >> Bxr >> Bxi;
  //   InputFileParam >> Byr >> Byi;
  //   InputFileParam >> Bzr >> Bzi;
  //   InputFileParam >> Zr[1] >> Zi[1];
  //   InputFileParam >> Zr[2] >> Zi[2];
  //   InputFileParam >> Zr[3] >> Zi[3];
  //   for (int s1=1;s1<=3;s1++) {
  //     for (int s2=1;s2<=3;s2++) {
  //       for (int s3=1;s3<=3;s3++) {
  // 	InputFileParam >> FrMatrix[s1][s2][s3] >> FiMatrix[s1][s2][s3];
  //       }
  //     }
  //   }
  //   InputFileParam.close();          // Close input file
	
  //   BxN=sqrt(Bxr*Bxr+Bxi*Bxi);  // Compute surface component in x-dir
  //   if (BxN<1e-4)  BxN=1e-4;    
  //   BxrN=Bxr/BxN;
  //   BxiN=Bxi/BxN;
	
  //   ByN=sqrt(Byr*Byr+Byi*Byi);  // Compute surface component in y-dir
  //   if (ByN<1e-4)  ByN=1e-4;
  //   ByrN=Byr/ByN;
  //   ByiN=Byi/ByN;
	
  //   BzN=sqrt(Bzr*Bzr+Bzi*Bzi);  // Compute surface component in z-dir
  //   if (BzN<1e-4)  BzN=1e-4;
  //   BzrN=Bzr/BzN;
  //   BziN=Bzi/BzN;
		
//   BxrN=0.0;
//   BxiN=-1.0;
//   ByrN=0.0;
//   ByiN=0.0;
//   BzrN=1.0;
//   BziN=0.0;

  //  cout << "BxN=" << BxrN << " + " << BxiN << "i" << endl;
  //  cout << "ByN=" << ByrN << " + " << ByiN << "i" << endl;
  //  cout << "BzN=" << BzrN << " + " << BziN << "i" << endl;  

  cout << "Input the number of harmonics: ";
  cin >> nm;
  cout << "Input the number of harmonics to be used for waveform reconstruction: ";
  cin >> nminput;
  cout << endl << "Number of harmonics= " << nm << endl;
  cout <<         "Recon. harmonics   = " << nminput << endl;
  //  cout <<         "Linear wave speed  = " << v0 << endl;
	
  const int PRECISION=12;     // Define output formatting parameters
  const int COLW=20;         
  cout.precision(PRECISION);  // Set precision
  cout.setf(ios::scientific); // Set scientific notation output
   
  nm2=nm*2;                   // Total number of harmonic comp. (real and imag)
  V=dvector(1,nm2);           // Allocate memory for harmonic data
  norm=dvector(1,nm);
  normdB=dvector(1,nm);
  dnormdBdx=dvector(1,nm);
  Vx=dvector(0,nm2);          // Allocate memory for waveform data
  Vxh=dvector(0,nm2);
  Vy=dvector(0,nm2);
  Vz=dvector(0,nm2);
   
  for (int n=0;n<=nm2;n++)  { // Initialize vectors
    Vx[n]=Vxh[n]=Vy[n]=Vz[n]=0.0;
  }

  for (int n=0;n<=nm2;n++)  { // Initialize vectors
    V[n]=0.0;
  }
	
  // Input positions for data file output
  cout << endl;
  cout << "Input filename for positions for data output (max 80 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 << endl;
	
  ifstream InputFilePositions(inputfile, ios::in);   // Open input file
  if (!InputFilePositions) {
    cerr << "Input Data File could not be opened" << endl;
    exit(1);
  }

  // Read in Fourier component files and process waveforms
  cout << "Input \"root\" filename for harmonic amplitude data: ";
  cin  >> setw(STRLEN) >> infilename3; 
  char inpathname3[STRLEN]="/home/kumon/crystal/nonlin/RK4/main/testBurg/";  
  char outpathname[STRLEN]="/home/kumon/crystal/nonlin/RK4/main/testBurg/";

  strcpy(outputfile2,outpathname);      // Copy pathname to string
  strcat(outputfile2,infilename3);      // Concatenate filename to pathname
  strcat(outputfile2,".prop");          // Concatenate suffix to filename
  ofstream OutProp;                     // Declaration of output file object 
  OutProp.open(outputfile2, ios::out);  // Open data file
  if (!OutProp) {                       // Error checking
    cerr << "File could not be opened" << endl;
    exit(1);
  }
  OutProp.precision(PRECISION);         // Set output precision
  OutProp.setf(ios::scientific);        // Set scientific notation
  cout << endl << "Propagation data file: " << outputfile2 << endl;

  // Input parameters for output of propagation data
  cout << "Propagation data file specifications:" << endl;
  cout << endl << "Input the first value of the harmonic number to be written     : ";
  cin  >> setw(STRLEN) >> nminit;       // Set first harm. number to be written
  cout << "Input the skip interval for the harmonic numbers to be written : ";
  cin  >> setw(STRLEN) >> nmoutskip;    // Set harmonic number skip increment
  cout << "Input the maximum number of the harmonic number to be written  : ";
  cin  >> setw(STRLEN) >> nmout;        // Set number of harmonics to output

  while (InputFilePositions >> positionstr) {
    // Create filenames that contain Fourier components
    createSuffix(positionstr,suffix);   // Create appropriate suffix
    strcpy(inputfile,inpathname3);      // Copy pathname to string
    strcat(inputfile,infilename3);      // Concatenate filename to pathname
    strcat(inputfile,suffix);           // Concatenate filename to suffix1
    strcat(inputfile,".dat");           // Concatenate filename to suffix2  
    cout << endl << "Input File  :  " << inputfile << endl;

    ifstream InputFileComponents(inputfile, ios::in);   // Open input file
    if (!InputFileComponents) {
      cerr << "Input Data File could not be opened" << endl;
      exit(1);
    }		

    for (int n=1;n<=nminput;n++) {  // Read in harmonic amplitude data 
      InputFileComponents >> num >> V[2*n-1] >> V[2*n] 
                          >> norm[n] >> normdB[n]; //>> dnormdBdx[n];
    }
    InputFileComponents.close();          // Close input file
    		
    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;
    }    

    CalcVx(nm,V,BxrN,BxiN,Vx);        // Calculate waveform data for each 
    //    CalcVxh(nm,V,BxrN,BxiN,Vxh);// direction (x,y,z).  Note that Vxh
    //    CalcVy(nm,V,ByrN,ByiN,Vy);  // is the Hilbert transformed x-dir
    //    CalcVz(nm,V,BzrN,BziN,Vz);  // waveform.

    // Output waveform data to files
    strcpy(outputfile,outpathname);      // Copy pathname to string
    strcat(outputfile,infilename3);      // Concatenate filename to pathname
    strcat(outputfile,suffix);           // Concatenate suffix
    cout << "Output Files:  " << outputfile << ".*" << endl;
  
    //    ofstream OutVx,OutVxh,OutVy,OutVz;  // Declaration of output file objects
    ofstream OutVx;                      // Declaration of output file objects

    OutVx.open(strcat(strcpy(temp,outputfile), ".vx"), ios::out);
    //    OutVxh.open(strcat(strcpy(temp,outputfile), ".vxh"), ios::out);
    //    OutVy.open(strcat(strcpy(temp,outputfile), ".vy"), ios::out);
    //    OutVz.open(strcat(strcpy(temp,outputfile), ".vz"), ios::out);
		
    OutVx.precision(PRECISION);  // Set precision
    //    OutVxh.precision(PRECISION); 
    //    OutVy.precision(PRECISION);  
    //    OutVz.precision(PRECISION);  
    OutVx.setf(ios::scientific); // Set scientific notation output
    //    OutVxh.setf(ios::scientific);
    //    OutVy.setf(ios::scientific); 
    //    OutVz.setf(ios::scientific); 
  
    x0=1.0/nm2;                // Spatial step size
		
    for (int n=0;n<=nm2;n++) { // Output waveform data to files
      x=x0*n;
      OutVx  << setw(COLW) << x << setw(COLW) << Vx[n]  << endl;
      //      OutVxh << setw(COLW) << x << setw(COLW) << Vxh[n] << endl;
      //      OutVy  << setw(COLW) << x << setw(COLW) << Vy[n]  << endl;
      //      OutVz  << setw(COLW) << x << setw(COLW) << Vz[n]  << endl;
    }

    OutVx.close();            // Close files
    //    OutVxh.close();
    //    OutVy.close();
    //    OutVz.close();

    // Write out harmonic propagation curve data to file
    position=atof(positionstr);           // Convert string to double
    OutProp << setw(COLW) << position;    // Write out position
    for (int n=nminit; n<=nmout; n+=nmoutskip) {  // Write out |V_n|
      OutProp << setw(COLW) << norm[n];
    }
    for (int n=nminit; n<=nmout; n+=nmoutskip) {  // Write out 20*\log_10|V_n|
      OutProp << setw(COLW) << normdB[n];
    }
    OutProp << endl;
  }
  
  OutProp.close();           // Close data file

  free_dvector(V,1,nm2);     // Deallocate memory for vectors
  free_dvector(dnormdBdx,1,nm2);
  free_dvector(norm,1,nm);
  free_dvector(normdB,1,nm);
  free_dvector(Vx,0,nm2);
  free_dvector(Vxh,0,nm2);
  free_dvector(Vy,0,nm2);
  free_dvector(Vz,0,nm2);
	
  return 0;
}