/*
Analysis tool for SPO-DVR data.
Last modification: 22-01-18

author: E.W.F. Smeets

Data is specified in a datalist file that also holds the dimension of all the splines needed as well as the state numbers for internal consistency checks
Rovib dat is specified in seperate file togeter with the quantun numbers for internal consistency checks

Input file 'input' has the neccesary parameters
--> Beam parameters are specified as commandline arguments
--> If no, or weird, input is given on the commandline, help is printed to the terminal

output files hold all relevant data: based on the Energy resolution specified.
--> degenracy averaged probabilities
--> rotational quadrupole alignment parameter
--> partifion function with state occupencies
--> fully state resolved reaction probabilities 

--code structure:
Each fully resolved state (v,j,mj) is a class
 -- initialization handels the interpolation of the (multiple) wavepacket(s) data 
 -- class can be evaluated for a particular energy by member function of class
 -- class can report state quantum numbers for internal checking

All other quantities, degeneracy averaged probabilities, rotational quadrupole alignming parameter,
are calculated using functions in the main program.
--> based on the findSTATE function that for a given v,j returns the integer
    of that state, with mj=0, in the array off classes that holds all states.

UPDATE: 13-7-2018
  --fixed double counting in RotQuad. 
  According to the degeneracy averaging, mj=0 should be counted once, all the others twice
*/

#include<iostream>
#include<fstream>
#include<string.h>
#include<sstream>
#include<math.h>
#include<iomanip>
#include<stdlib.h>
#include<stdio.h>
#include<time.h>
#include<unistd.h>
#include "ops-spline.cc"

extern  void spline_init(double*, double*, int, double*);
extern  double spline_interpolation(double*, double*, double*, int, double );


using namespace std;

//extern  void spline_init(double*, double*, int, double*);
//extern  void spline_interpolation(double*, double*, double*, int, double, double*, double* );

//Read file into file buffer -- ignore comments
  int rem_com(char* filename, char* streamstring, int string_length){ 
    const char com_B = '#';
    const char com_E = '\n';
    int pos = 0;
    char cc;
    ifstream inf(filename);
    while(inf.get(cc)&& pos < string_length-1){
      if(cc != com_B) 
	streamstring[pos++] = cc;
      else{
	while(cc != com_E && inf.get(cc));
	streamstring[pos++] = com_E;
      }
    }
    streamstring[pos] = 0;
    if(pos == string_length-1){
      cerr << "Buffer size exceeded !\n"; exit(0);
    }
    return(strlen(streamstring));
  }
int TotStates(int v0startj, int v0endj, int v1startj, int v1endj){
  int v, j, mj, startj, endj, TotStates=0;  
  for ( v = 0; v <= 1; v++){
    if  ( v == 0 ) {startj = v0startj; endj = v0endj;}
    if  ( v == 1 ) {startj = v1startj; endj = v1endj;}
    j = startj;
    while ( j <= endj ) {
      mj = 0;
      while ( mj <= j ){
      //cout << TotStates << "  v" << v << "j" << j << "mj" << mj <<"\n";  
      mj++; TotStates++;
      } //while mj
      j++;
      //cout << "\n";
    } //while j
    //cout << "\n";
  } //for v
  return TotStates;
} //TotStates

//PLAN: 
//make class wich encompasses a state
//two public functions: 1 initCLASS input: WPper and array with file names --> setup the neccesary splines
//			2 evalCLASS input: E in eV  out: ReacProb for that state in that energy
class STATE {
  private:
    int WPper, state;
    int qnu, qj, qmj;
    double rovib;
    string DataList;

    struct splines
    {
      string filename;
      double* values_X;
      double* values_FUNC_X;
      double* spl_cf;
      int DIM;
    };//struct splines


    splines* ALLsplines;
  public:
    void initSTATE(int instate, int inWPper, char* DataList);
    double evalSTATE(double valX );
    void print( );
    int getQNU( );
    int getQJ( );
    int getQMJ( );
    void findErange(double Erange[] );
}; //class STATE

int STATE::getQNU( )
{
  return ( qnu );
};//int STATE::getQNU
int STATE::getQJ( )
{
  return ( qj );
};//int STATE:getQJ
int STATE::getQMJ( )
{
  return ( qmj );
};//int STATE:getQMJ
void STATE::initSTATE(int instate, int inWPper, char* inDataList)
{

  WPper = inWPper;
  state = instate;
  DataList = inDataList;
  int buff_length  =  65536;
  char*  file_buff =  new char[buff_length];
  string dummy;
  double ddummy;
  int idummy;

  //make structs for the splines
  ALLsplines = new splines[WPper];

  //read file with filenames
  rem_com(inDataList, file_buff,  buff_length);
  istringstream ist(file_buff);
  //extract the WPper filenames we need
  //datalist structure                               qnu       qj       qmj       DIM     filename
  for (int i = 0; i < state * WPper; i++) { ist >> idummy >> idummy >> idummy >> idummy >> dummy;}
  for (int i = 0; i < WPper; i++) { ist >> qnu >> qj >> qmj >> ALLsplines[i].DIM >> ALLsplines[i].filename;}

  //cout << "I'm state "<<state<<"  WPper = " << WPper << "  Datalist= "<< DataList <<"\n";
  //cout << "     with splines:\n";
  //for (int i = 0; i < WPper; i++)
  //{
  //  cout << "spline " << i << "  DIM= " << ALLsplines[i].DIM << "   file= " << ALLsplines[i].filename << "\n";
  //}//for i

  //read file and initialize splines
  ifstream fin;
  for (int i = 0; i < WPper; i++)
  {
    //open data file 
    fin.open(ALLsplines[i].filename.c_str()); 
    if (fin.fail())
    {
      cout << "file operation failed --> " << ALLsplines[i].filename << " \n";
      exit(1);
    }//if
    
    //allocate spline arrays
    ALLsplines[i].values_X = 		new double[ALLsplines[i].DIM];
    ALLsplines[i].values_FUNC_X = 	new double[ALLsplines[i].DIM];
    ALLsplines[i].spl_cf = 		new double[ALLsplines[i].DIM];

    for (int j = 0; j < ALLsplines[i].DIM; j++)
    {
      //Read in values_X and values_FUNC_X from ALLsplines[i].filename
      fin >> ddummy >> ALLsplines[i].values_X[j] >> ALLsplines[i].values_FUNC_X[j] >> ddummy;
      //If smaller then 0 --> set to zero
      if ( ALLsplines[i].values_FUNC_X[j] < 0 ) { ALLsplines[i].values_FUNC_X[j] = 0.; }
      //cout << "j= " << j << "  values_x= " << ALLsplines[i].values_X[j] << " values_FUNC_X= " << ALLsplines[i].values_FUNC_X[j] << "  file= " << ALLsplines[i].filename << "\n";
    }//for j
    fin.close(); 
    
    //initialize splines
    spline_init( ALLsplines[i].values_X, ALLsplines[i].values_FUNC_X, ALLsplines[i].DIM, ALLsplines[i].spl_cf );
  }//for i
  //File IO done -- splines initiated 
  
}; //STATE::initSTATE

double STATE::evalSTATE( double valX )
{
  int check = 0;
  double splineX = 0., fr;
  
  //check for out of bounds
  if ( valX < ALLsplines[0].values_X[0] )
  {
    //cout << "CAUTION valX ("<<valX<<") is lower then lowest known energy ("<<ALLsplines[0].values_X[0]<<") .... \n";
    valX = ALLsplines[0].values_X[0];
  }
  if ( valX > ALLsplines[WPper-1].values_X[ALLsplines[WPper-1].DIM-1] )
  {
    cout << "CAUTION valX ("<<valX<<") is greater then highest known energy ("<<ALLsplines[WPper-1].values_X[ALLsplines[WPper-1].DIM-1]<<") .... \n";
    valX = ALLsplines[WPper-1].values_X[ALLsplines[WPper-1].DIM-1];
  }

  //eval splines
  for (int i = 0; i < WPper; i++)
  {
    if ( valX >= ALLsplines[i].values_X[0]  &&  valX <= ALLsplines[i].values_X[ALLsplines[i].DIM-1] )
    {
      splineX = spline_interpolation( ALLsplines[i].values_X, ALLsplines[i].values_FUNC_X, ALLsplines[i].spl_cf, ALLsplines[i].DIM, valX );
      check++;
    }//if
  }//for i
  
  if ( check != 1 )
  {
    splineX = 0.;
    for (int i = 0; i < WPper - 1; i++)
    {
      if ( valX > ALLsplines[i+1].values_X[0]  &&  valX < ALLsplines[i].values_X[ALLsplines[i].DIM-1] )
      {
        fr = ( valX - ALLsplines[i+1].values_X[0] ) / ( ALLsplines[i].values_X[ALLsplines[i].DIM-1] - ALLsplines[i+1].values_X[0] );
        //cout << "fr= "<< fr << "   cos(fr)= "<< cos( fr * 0.5 * 3.14159265) << "\n";
        splineX = (      cos( fr * 0.5 * 3.14159265) ) * spline_interpolation( ALLsplines[i].values_X, ALLsplines[i].values_FUNC_X, ALLsplines[i].spl_cf, ALLsplines[i].DIM, valX );
        splineX +=( 1.-  cos( fr * 0.5 * 3.14159265) ) * spline_interpolation( ALLsplines[i+1].values_X, ALLsplines[i+1].values_FUNC_X, ALLsplines[i+1].spl_cf, ALLsplines[i+1].DIM, valX );
      }//if
    }//for i
  }//if

  return( splineX );
}//STATE::evalSTATE

void STATE::findErange(double Erange[])
{
  Erange[0] = ALLsplines[0].values_X[0];
  Erange[1] = ALLsplines[WPper - 1].values_X[ALLsplines[WPper -1].DIM -1];
};//void STATE::findErange()

void STATE::print()
{
  cout << "print function\n";
  cout << "I'm state "<<state<<"  WPper = " << WPper << "  Datalist= "<< DataList <<"\n";
  cout << "     with splines:\n";
  for (int i = 0; i < WPper; i++)
  {
    cout << "spline " << i << "  DIM= " << ALLsplines[i].DIM << "   file= " << ALLsplines[i].filename << "\n";
  }//for i

};
//function to find the correct number belonging to a state nu = x j=0
//this in order to generate degeneracy averaged reaction probabilities and others
int findSTATE(int qnu, int qj, int v0endj){
  int state = 0, endj, v, j, mj;
  //cout << "\n\n find state   input qnu="<<qnu<<" qj="<<qj<<" v0endj="<<v0endj<<" \n";
  v = 0;
  while ( v <= qnu )
  {
    if ( qnu == 1 && v == 0 ) { endj = v0endj; } 
    if ( qnu == 1 && v == 1 ) { endj = qj;     }
    if ( qnu == 0 && v == 0 ) { endj = qj;     }
    j = 0;
    while ( j <= endj )
    {
      mj = 0;
      while ( mj <= j)
      {
        //cout << "nu="<<v<<" j="<<j<<" mj="<<mj<<" state="<<state<<"\n";
        mj++;
        state++; 
      }//while mj
      j++;
    }//while j
    v++;
  }//while v
  //returns the state number of ALLSTATES of state nu, j , mj=0
  return ( state - qj - 1 );
}// int findSTATE

//Degeneracy averaged reaction probabilities
double DegAveP( int qnu, int qj, int v0endj, double En, STATE ALLSTATES[])
{
  //Thesis M. Wijzenbroek p41 eq. 2.32
  //P_deg(qnu, qj, En) = /Sum^qj_mj=0 (2 -/d_mj) P_r(qnu, qj, qmj, En) / (2J+1)
  int numst;
  double P;

  //find state number 
  numst = findSTATE( qnu, qj, v0endj );

  //internal consistency check  
  if ( ALLSTATES[numst].getQNU() != qnu || ALLSTATES[numst].getQJ() != qj || ALLSTATES[numst].getQMJ() != 0 )
  {
    cout << "Error DegAveP  -- internal bookkeeping states incorrect\n"
         << "querry      qnu=" << qnu << " qj= "<< qj << " state number=" << numst << "\n"
         << "reply state qnu=" << ALLSTATES[numst].getQNU() << " qj=" << ALLSTATES[numst].getQJ() << " qmj=" << ALLSTATES[numst].getQMJ() << "\n";
  }//if
  P = 0.;
  //cout << "DegAveP ------------- qnu="<< qnu << " qj=" << qj << "\n";
  for (int mj = 0; mj <= qj; mj++)
  {
    if ( mj == 0 )
    {
      P += ALLSTATES[numst + mj].evalSTATE( En ) / ( 2. * qj + 1. );
    }//if
    else
    {
      P += 2. * ALLSTATES[numst + mj].evalSTATE( En ) / ( 2. * qj + 1. );
    }//else
    //cout << "mj=" << mj << " statemj="<< ALLSTATES[numst+mj].getQMJ() << " P=" << P << "\n";
  }//for mj
  return ( P ); 
}//double DegAveP()

double RotQuad( int qnu, int qj, int v0endj, double En, STATE ALLSTATES[])
{
  //Thesis M. Wijzenbroek p41 eq. 2.34
  //A_0^(2)(qnu, qj, En)= /Sum_mj^J P_r(qnu, qj, qmj, En) * ( (3*mj*mj / (J*(J+1))) - 1)
  //                       -------------------------------------------------------------
  //                                    /Sum_mj^J P_r(qnu, qj, qmj, En)

  int numst;
  double A0, sum;
 
  //find state number 
  numst = findSTATE( qnu, qj, v0endj );

  //internal consistency check  
  if ( ALLSTATES[numst].getQNU() != qnu || ALLSTATES[numst].getQJ() != qj || ALLSTATES[numst].getQMJ() != 0 )
  {
    cout << "Error RotQuad  -- internal bookkeeping states incorrect\n"
         << "querry      qnu=" << qnu << " qj= "<< qj << " state number=" << numst << "\n"
         << "reply state qnu=" << ALLSTATES[numst].getQNU() << " qj=" << ALLSTATES[numst].getQJ() << " qmj=" << ALLSTATES[numst].getQMJ() << "\n";
  }//if

  A0 = 0.; sum = 0.;
  //cout << "RotQuad ------------- qnu="<< qnu << " qj=" << qj << "\n";
  sum += ALLSTATES[numst + 0].evalSTATE( En );
  A0  += ALLSTATES[numst + 0].evalSTATE( En ) * (  ((3. * 0 * 0)/(qj * (qj +1.))) -1. ) ;
  for (int mj = 1; mj <= qj; mj++)
  {
    sum += ALLSTATES[numst + mj].evalSTATE( En ) * 2.;
    A0  += ALLSTATES[numst + mj].evalSTATE( En ) * (  ((3. * mj * mj)/(qj * (qj +1.))) -1. ) * 2.;
    //cout << "rotquad qnu=" << qnu << "  qj=" << qj << "  mj=" << mj << " statemj="<< ALLSTATES[numst+mj].getQMJ() << " A0=" << A0  << " Sum=" << sum << "\n";
  }//for mj
  A0 = A0 / sum;
  
  return ( A0 );
}//double RotQuad()

double PartFunc(int v0endj, int v1endj, int RVvmax, int RVjmax, double TN, double RotCool, double RoVib[], bool flag)
{
  //Calculate partition function based on nozzle temperature
  //flag = true --> for actual states defined by v0endj and v1endj
  //flag = false -> for all states in RoVib --> provide check to see if enough states are calculated to describe molbeam.
  int vmax, jmax;
  double kB = 3.1668152e-06; //boltzmann factor Hartree/K
  double partition_func = 0.;
  if ( flag == true  ) { vmax = 1; 	}
  if ( flag == false ) { vmax = RVvmax; }
  for (int v = 0; v <= vmax; v++)
  {
    if ( flag == true && v == 0 ) { jmax = v0endj; }
    if ( flag == true && v == 1 ) { jmax = v1endj; }
    if ( flag == false ) 	  { jmax = RVjmax; }
    for (int j = 0; j <= jmax; j++)
    {
      double FacVib = (RoVib[  v * (RVjmax + 1)     ] - RoVib[0                  ]) / (kB * TN);
      double FacRot = (RoVib[ (v * (RVjmax + 1)) + j] - RoVib[ (v * (RVjmax + 1))]) / (RotCool * kB * TN);
      double val    = (2. * j + 1.) * exp( -FacVib ) * exp( -FacRot );
      if( j%2 == 0   ) partition_func += val;
      if((j+1)%2 == 0) partition_func += (3*val);
    }//for j 
  }//for v
  return ( partition_func );
}//double PartFunc()

double BoltzW(int qnu, int qj, int RVjmax, double TN, double RotCool, double RoVib[], double PF)
{
  double kB = 3.1668152e-06;
  double FacVib = (RoVib[  qnu * (RVjmax + 1)     ] - RoVib[0                    ]) / (kB * TN);
  double FacRot = (RoVib[ (qnu * (RVjmax + 1)) +qj] - RoVib[ (qnu * (RVjmax + 1))]) / (RotCool * kB * TN);
  double val    = (2. * qj + 1.) * exp( -FacVib ) * exp( -FacRot );
  if((qj+1)%2 == 0) val = val*3.0;
  val = val / PF;
  return (val);
}//double BoltzW()

//mono energetic reaction probability
double Rmono(double En, int v0endj, int v1endj, int RVjmax, double TN, double RotCool, double RoVib[], double PF, STATE ALLSTATES[]){
  int endj;
  double sum = 0.;
  for (int v = 0; v <= 1; v++)
  {
    if ( v == 0 ) { endj = v0endj; }
    if ( v == 1 ) { endj = v1endj; }
    for (int j = 0; j <= endj; j++)
    {
      //loop over all v,j states
      sum += BoltzW(  v, j, RVjmax, TN, RotCool, RoVib, PF ) * 
	     DegAveP( v, j, v0endj, En, ALLSTATES);
      //cout << "v= " << v << " j= " << j << "  sum= "<< sum <<"\n"; 
    }//for j   
  }//for v
  return ( sum );
}//double Rmono()

double eVtoms(double eV, double MassH)
{
  double ms;
  ms = sqrt( 2. * eV / ( 2 * MassH) ) * 419382.9;
  return ( ms );
}//double eVtoms()

bool does_file_exist(char fileName[] )
{
    std::ifstream infile(fileName);
    return infile.good();
}
int main(int argc, char* argv[]){
  if ( argc < 5 || argc > 6 || atof( argv[2] ) > 3 || atof( argv[1] )> 3000 )
  {
    cout << "Program to compute molecular beam reaction probabilities from \n"
         << "fully state-resolved data. Designed for the output of SPO-DVR.\n"
         << "\n"
         << "command line arguments:\n"
         << "   1. Nozzle temperature [K]\n"
         << "   2. Stream energy [eV] <-- average energy of the molecular beam\n"
         << "   3. spread/alpha [m/s]\n"
         << "   4. energy resolution [eV]\n" 
         << "   5. Angle of incidence [degrees] (optional) default = normal incidence\n"
         << "\n"
         << "The input file: 'input'\n"
         << "2                       #number of wavepackets/datasets per v,j,mj state \n"
	 << "0 2                     #v0-jstart v0-jend \n"
         << "0 1                     #v1-jstart v1-jend \n"
         << "0.8                     #rotational cooling factor \n"
         << "RoVibLevels-Cu211.dat   #file with RoVib eigenvalues for H2\n"
	 << "DataList                #list with all QD reaction probabilities\n"
	 << "\n"
         << "The Rovibrational data file:\n"
         << "#Rovib Levels in eV for system from SPO-DVR -- each v should have j_max entries\n"
	 << "3 #v_max \n"
	 << "35 #j_max\n"
         << "144 #total values\n"
	 << "[v j Erovib[v,j](eV)] <-- in this example 144 lines \n"
	 << "\n"
	 << "The DataList file:\n"
	 << "#The datalist file holds the filenames for all wavepackets.\n"
	 << "#The files should be ordered according to ascending quantum numbers\n"
	 << "v j mj [number of points] filename\n";

    return 0;
  }//if help 
  char inputfile[] = "input";		//inputfile 
  int states;				//number of states
  int WPper;				//WP's per state
  int v0startj;		int v0endj;	//for the data
  int v1startj;		int v1endj;
  double RotCool;			//Rotational cooling fraction of nozzle temp
  double Erange[2];			//array for the minimum and maximum energy for which there is data. 

  char RoVibname[128];			//RoVib energy levels filename
  int RVvmax, RVjmax, RVnum;		//Rovib max nu, Rovib max J, Rovib tot. states
  double* RoVib;			//array of rovib level energies
  int dummy;

  double PF, PFfull;			//Partition function for available states, PF, and based on all RoVib states, PFfull

  char DataList[128];			//list with the the file names of all the (sub)states
 
  int buff_length  =  65536;			
  char*  file_buff =  new char[buff_length];

  //commandline arguments for molecular beam parameters
  double TN = atof( argv[1] );			//Nozzle Temp [k]
  double ES = atof( argv[2] );			//Estream [eV]
  double AL = atof( argv[3] );			//alpha [m/s]
  double Eres=atof( argv[4] );			//deltaE [eV]
  double ANG;
  if ( argc != 6 )
  {
    ANG = 0.; 
  }//if
  else 
  {
    ANG = atof( argv[5] );
    if ( ANG < 0 || ANG >= 90 )
    {
      cout << "Incedence angle needs to be between [0,90) interval....\n";
      return 0;
    }//if check angle
  }//else
  double MassH = 1837.3622; 			//ATOMICMASS OF H^1
  double mu = MassH*MassH/(MassH+MassH);	//reduced mass
  double VS = eVtoms( ES, MassH );		//Vstream [m/s]

  //check if input file 
  if ( does_file_exist( inputfile ) != true )
  {
    cout << "inputfile: 'input' was not found....\n"
         << "exiting...\n";
    return 0;
  }//if
  cout << "inputfile \n";

  //Read input file with parameters
  rem_com(inputfile, file_buff,  buff_length);
  istringstream ist(file_buff);
  ist >> WPper >> v0startj >> v0endj >> v1startj >> v1endj >> RotCool >> RoVibname >> DataList;
  delete file_buff; 
  ist.str(""); 
  //-----------------------------------------  

  //read RoVib lebels
  file_buff = new char[buff_length];
  rem_com( RoVibname, file_buff,  buff_length);
  istringstream ist2(file_buff);
  ist2 >> RVvmax >> RVjmax >> RVnum;
  if ( (RVvmax + 1)*(RVjmax + 1) != RVnum )
  {
    cout << "ERROR in reading RoVib data --> RVvmax * RVjmax != RVnum \n"
         << "RVvmax=" << RVvmax << " RVjmax=" << RVjmax << " RVnum=" << RVnum << "\n";
  }//if
  RoVib = new double[RVnum];
  for (int i = 0; i < RVnum; i++)
  {
    ist2 >> dummy >> dummy >> RoVib[i]; 
    RoVib[i] = RoVib[i] / 27.211383; //convert to hartree
  }//for i  
  delete file_buff;
  ist2.str("");
  //-----------------------------------------
  cout << "RoVib \n";

  ofstream out;		out.open("state.occupencies");
  //compute and check partition function
  //and compute state occupencies based on nozzle temp
  PF     = PartFunc( v0endj, v1endj, RVvmax, RVjmax, TN, RotCool, RoVib, true);
  PFfull = PartFunc( v0endj, v1endj, RVvmax, RVjmax, TN, RotCool, RoVib, false);
  out << "#Partition functions for available data:\n"
       << "#    Rotational cooling factor: "<<RotCool<<"\n"
       << "#    PF for nu in [0,1] with j in [0,"<<v0endj<<"/"<<v1endj<<"]                   --> PF    ="<<PF<<"\n"
       << "#    PF on provided RoVib data: nu in [0,"<<RVvmax<<"] with j in [0,"<<RVjmax<<"] --> PFfull="<<PFfull<<"\n"
       << "#    Difference="<<PFfull-PF<<"    --> "<< 100-((PF/PFfull)*100) << "%\n#\n"
       << "#Beam Parameters:\n"
       << "#   Tnozzle=        "<<TN<<" [K]\n"
       << "#   EStream=        "<<ES<<" [eV]            VStream= "<<VS<<" [m/s]\n"
       << "#   Alpha=          "<<AL<<" [m/s]\n"
       << "#   Eresolution=    "<<Eres<<" [eV]\n"
       << "#   Angle=          "<<ANG<<" [degrees]\n#\n"
       << "#Fractional occupencies per state:\n"
       << "#nu | J | Boltzman weight \n";
  for (int i = 0; i <=RVvmax; i++)
  {
    for (int j = 0; j <=RVjmax; j++)
    {
      out << i << "   " << j << "   " << BoltzW(i, j, RVjmax, TN, RotCool, RoVib, PF) << "\n";
    }
    out << "\n";
  }
  out.close();
  //-----------------------------------------

  //compute total number of states
  states = TotStates( v0startj, v0endj, v1startj, v1endj);

  //input method works -- can read strings and values 
  //cout << "states=   " << states << "\n"
  //     << "WPper=    " << WPper << "\n"
  //     << "v0startj= " << v0startj << "		v0endj= " << v0endj << "\n"
  //     << "v1startj= " << v1startj << "         	v1endj= " << v1endj << "\n"
  //     << "RoVibname= "<< RoVibname << "\n"
  //     << "DataList= " << DataList << "\n";
  
  //initialize classes
  STATE ALLSTATES[states];
  for ( int i = 0; i < states; i++) ALLSTATES[i].initSTATE( i, WPper, DataList );
  ALLSTATES[0].findErange( Erange );

  //now for finding the normalization factor for the velocity distribution
  //and checking if the defined energy interval is sufficient. 
  double sum1 = 0., sum2 = 0., v, sumV = 0., dvde;
  for (double Energy = 0. + Eres; Energy < Erange[1]; Energy += Eres)
  {
    Energy = Energy*cos(ANG*3.14159265/180); //normal energy scaling
    v = eVtoms( Energy, MassH);
    //sum1 = sum from 0 to Erange[1]
    //sum2 = sum from Erange[0] to Erange[1]
    
    //sum1 += v*v*v * exp( (-(v-VS)*(v-VS)) / ( AL*AL ));
    dvde = 1. / sqrt( (2* Energy * 1.6021766e-19) / 3.3474472e-27);
    //dvde = 1.;
    sum1 += v*v*v * exp( (-(v-VS)*(v-VS)) / ( AL*AL )) *dvde ;
    sumV += v*v*v * exp( (-(v-VS)*(v-VS)) / ( AL*AL )) ;
    //cout << "sum1= "<<sum1<< " v="<<v<<" En= "<<Energy<<"   "<< v*v*v * exp( (-(v-VS)*(v-VS)) / ( AL*AL ))<<"  " << v << "  "<< VS << "  " << AL << "\n";
    if ( Energy >= Erange[0] ) 
    { 
      sum2 += v*v*v * exp( (-(v-VS)*(v-VS)) / ( AL*AL )) *dvde; 
      //cout << "sum2=" <<sum2<<" v="<<v<<" En="<<Energy<<"   "<< v*v*v * exp( (-(v-VS)*(v-VS)) / ( AL*AL ))<<"\n";
    }
  }//for energy

  out.open("out.info");
  out  << "Analysis script --> Corrected for dv/dE jacobian when calculating mol beam S0\n"
       << "#BeamParameters:\n"
       << "#   Tnozzle=            "<<TN<<" [K]\n"
       << "#   EStream=            "<<ES<<" [eV]		VStream= "<<VS<<" [m/s]\n"
       << "#   Alpha=              "<<AL<<" [m/s]\n"
       << "#   Eresolution=        "<<Eres<<" [eV]\n"
       << "#   Rotational cooling= "<<RotCool<<"\n"
       << "#   Angle=          	   "<<ANG<<" [degrees]\n#\n"
       << "#Input:\n"
       << "#   Total number of states=          " << states << "\n"
       << "#   Number of wavepackets per state= " << WPper << "\n"
       << "#   v0startj=                        " << v0startj << "	v0endj= " << v0endj << "\n"
       << "#   v1startj=                        " << v1startj << "	v1endj= " << v1endj << "\n"
       << "#   RoVib data filename=             " << RoVibname << "\n"
       << "#   Datalist filename=               " << DataList << "\n#\n"
       << "#Partition function:\n"
       << "#    PF for nu in [0,1] with j in [0,"<<v0endj<<"/"<<v1endj<<"]                   --> PF    ="<<PF<<"\n"
       << "#    PF on provided RoVib data: nu in [0,"<<RVvmax<<"] with j in [0,"<<RVjmax<<"] --> PFfull="<<PFfull<<"\n"
       << "#    Difference="<<PFfull-PF<<"    --> "<< 100-((PF/PFfull)*100) << "%\n#\n"
       << "#Integration of velocity distribution:\n"
       << "#   sum1="<<sum1<<" on range [0:"<<Erange[1]<<"] [eV]\n"
       << "#   sum2="<<sum2<<" on range ["<<Erange[0]<<":"<<Erange[1]<<"] [eV]\n"
       << "#fraction that can be represented by current data: "<<(sum2/sum1)*100<<"%\n#\n"
       << "#Energy [eV] | E distro normalized | v [m/s] | v distro normalized | Rmono(Energy,Tnozzle)\n";
  //now calculate the molecular beam
  double Rbeam = 0.;
  for (double Energy = 0. + Eres; Energy < Erange[1]; Energy += Eres)
  {
    Energy = Energy*cos(ANG*3.14159265/180); //normal energy scaling
    v = eVtoms( Energy, MassH);
    dvde = 1. / sqrt( (2* Energy * 1.6021766e-19) / 3.3474472e-27);
    //dvde = 1.;
    Rbeam +=((v*v*v * exp( (-(v-VS)*(v-VS)) / ( AL*AL )) *dvde)/sum2) *  Rmono( Energy, v0endj, v1endj, RVjmax, TN, RotCool, RoVib, PF, ALLSTATES);

    out  << Energy                                                  << "          " 
         << (v*v*v * exp( (-(v-VS)*(v-VS)) / ( AL*AL )) *dvde)/sum2 << "          "
         << v                                                       << "          " 
         << (v*v*v * exp( (-(v-VS)*(v-VS)) / ( AL*AL )) ) / sumV    << "          "
         << Rmono( Energy, v0endj, v1endj, RVjmax, TN, RotCool, RoVib, PF, ALLSTATES) << "\n";
  }//for Energy
  out << "Rbeam = " << Rbeam << "\n";
  out.close();

  //print to terminal
  cout << "Analysis script --> Corrected for dv/dE jacobian when calculating mol beam S0\n"
       << "BeamParameters:\n"
       << "   Tnozzle            = "<<TN<<" [K]\n"
       << "   EStream            = "<<ES<<" [eV]            VStream = "<<VS<<" [m/s]\n"
       << "   Alpha              = "<<AL<<" [m/s]\n"
       << "   Eresolution        = "<<Eres<<" [eV]\n"
       << "   Rotational cooling = "<<RotCool<<"\n"
       << "   Angle              = "<<ANG<<" [degrees]\n\n"
       << "Input:\n"
       << "   Total number of states          = " << states << "\n"
       << "   Number of wavepackets per state = " << WPper << "\n"
       << "   v0startj                        = " << v0startj << "      v0endj= " << v0endj << "\n"
       << "   v1startj                        = " << v1startj << "      v1endj= " << v1endj << "\n"
       << "   RoVib data filename             = " << RoVibname << "\n"
       << "   Datalist filename               = " << DataList << "\n\n"
       << "Partition function:\n"
       << "    PF for nu in [0,1] with j in [0,"<<v0endj<<"/"<<v1endj<<"]                   --> PF    ="<<PF<<"\n"
       << "    PF on provided RoVib data: nu in [0,"<<RVvmax<<"] with j in [0,"<<RVjmax<<"] --> PFfull="<<PFfull<<"\n"
       << "    Difference="<<PFfull-PF<<"    --> "<< 100-((PF/PFfull)*100) << "%\n"
       << "    Fractional occupencies per state  --> file: state.occupencies \n\n"
       << "Integration of velocity distribution:\n"
       << "   sum1="<<sum1<<" on range [0:"<<Erange[1]<<"] [eV]\n"
       << "   sum2="<<sum2<<" on range ["<<Erange[0]<<":"<<Erange[1]<<"] [eV]\n"
       << "   Distribtuion in file: out.info \n"
       << "fraction that can be represented by current data: "<<(sum2/sum1)*100<<"%\n\n"
       << "   --------------------------------------------------------\n"
       << "   |   __                                                 |\n"
       << "   |   | \\                                                |\n"
       << "   |   |_/                                                |\n"
       << "   |   | \\                                                |\n"
       << "   |   |  \\ BEAM = "<<Rbeam<<"                               |\n"
       << "   |                                                      |\n"
       << "   --------------------------------------------------------\n";

  //output degenerace averaged reaction probabilities
  out.open("ALL.degeneracy.averaged");
  int endj;
  //generate comment line
  out << "#Degeneracy averaged reaction probabilities -- Energy resolution: "<<Eres<<"\n"
      << "# energy [eV] | ";
  for (int i = 0; i <=1; i++)
  {
    if ( i == 0 ) { endj = v0endj; }
    if ( i == 1 ) { endj = v1endj; }
    for (int j = 0; j <= endj; j++)
    {
      out << "v"<<i<<"j"<<j<<" | ";
    }//for j
  }//for i
  out << "\n";
  for (double Energy = Erange[0]; Energy < Erange[1]; Energy += Eres)
  {
    out << Energy << "  ";
    for (int i = 0; i <= 1; i++)
    {
      if ( i == 0 ) { endj = v0endj; }
      if ( i == 1 ) { endj = v1endj; }
      for (int j = 0; j <= endj; j++)
      { 
        out << DegAveP( i, j, v0endj, Energy, ALLSTATES ) << "  ";      
      }//for j
    }//for i
    out << "\n";
  }//for Energy
  out.close();

  //Quadrupole Alignment parameters
  out.open("ALL.quadrupole");
  //generate comment line
  out << "#quadrupole alignment parameters -- Energy resolution: "<<Eres<<"\n"
      << "# energy [eV] | ";
  for (int i = 0; i <= 1; i++)
  {
    if ( i == 0 ) { endj = v0endj; }
    if ( i == 1 ) { endj = v1endj; }
    for (int j = 1; j <= endj; j++)
    {
      out << "A_0^(2) (v"<<i<<"j"<<j<<") | ";
    }//for j
  }//for i
  out << "\n";
  for (double Energy = Erange[0]; Energy < Erange[1]; Energy += Eres)
  {
    out << Energy << "  ";
    for (int i = 0; i <= 1; i++)
    {
      if ( i == 0 ) { endj = v0endj; }
      if ( i == 1 ) { endj = v1endj; }
      for (int j = 1; j <= endj; j++)
      {
        out << RotQuad( i, j, v0endj, Energy, ALLSTATES ) << "  ";
      }//for j  
    }//for i
    out <<"\n";
  }//for Energy
  out.close(); 
 
  out.open("fully.state.resolved");
  out << "#fully state resolved data: according to ascending quantum numbers\n";
  for (double Energy = Erange[0]; Energy < Erange[1]; Energy += Eres)
  {
    out << Energy << "  ";
    for (int i = 0; i < states; i++)
    {
       out << ALLSTATES[i].evalSTATE( Energy ) << "  ";
    }//for i
    out << "\n";
  }//for Energy

  return 0;
}//main