//
#ifdef WIN32
#pragma warning(disable:4786)
/*Warning messages that are identified as Warning C4786: are created
  when MSVC generates extremely long names that it uses for debugging
  purposes. The long names are generated by the template expansion
  process and the warning messages normally can be ignored. Since these
  warnings tend to hide more interesting warning/error messages, you may
  wish to suppress the warning*/
#endif

#include <sstream>
#include <iostream>
#include <cantera/Cantera.h>
#include <cantera/onedim.h>
#include <cantera/IdealGasMix.h>
#include <cantera/equilibrium.h>
#include <cantera/transport.h>
#include "example_utils.h"

using namespace Cantera;

int flamespeed(ofstream& myresult , double& p) {
  try {
    int i;
    IdealGasMix gas("1STEP.cti","gas");
    double transl = 0.0;
    doublereal temp = 300.0; // K
    doublereal pressure = 1.0*OneAtm; //Atm
    doublereal uin=0.37; //m/sec
    double phi=p;
    int n_reac = gas.nReactions();

    int nsp = gas.nSpecies();
    double n_press = 1.0;

/*  cout<<"Temperature of Fresh Gas [K]:"<<endl;
    cin>>temp;*/
    cout<<"Pressure [Atm]:"<<endl;
    cin>>n_press;
    cout<<"Equivalence Ratio:"<<endl;
    cin>>phi;
    pressure = n_press*OneAtm;
 /*   cout<<"Translation"<<endl;
    cin>>transl;
    //cout<<report(gas);
  */  
    vector_fp x;
    x.resize(nsp);
//++++++++COMPUTION OF INITIAL COMPOSITION AS A FUNCTION OF PHI++++++++++++++++++++++//
    doublereal nu_F=1.0;
    doublereal nu_O=5.0;
    doublereal fa_stoic=nu_O/nu_F;
    doublereal X_F=1/(1+4.76*fa_stoic/phi);
    doublereal X_O2=1/((phi/fa_stoic)+4.76);
    doublereal X_N2=3.76*X_O2;
    for(int k=0;k<nsp;k++){
      if(k==gas.speciesIndex("C3H8")){ x[k]=X_F; }
      else if(k==gas.speciesIndex("O2")){ x[k]=X_O2; }
      else if(k==gas.speciesIndex("N2")){ x[k]=X_N2; }
      else{ x[k]=0.0; }
    }
    gas.setState_TPX(temp,pressure,DATA_PTR(x));
    doublereal rho_in=gas.density();
    double *yin=new double[nsp];
    gas.getMassFractions(yin);

//++++++++EQUILIBRIUM STATE CALCULATION++++++++++++++++++++++++++++
    try {
      equilibrate(gas,"HP");
    }
    catch (CanteraError) {
      showErrors(cout);
    }
    double *yout=new double[nsp];
    gas.getMassFractions(yout);
    doublereal rho_out = gas.density();
    doublereal Tad=gas.temperature();
    //=============  build each domain ========================


    //-------- step 1: create the flow -------------

    FreeFlame flow(&gas);

    // create an initial grid
    int nz=10; 
    doublereal lz=0.01;
    doublereal *z=new double[nz+1];
    doublereal dz=lz/((doublereal)(nz-1));
    for(int iz=0;iz<nz;iz++){
      z[iz]=((doublereal)iz)*dz;
    }
    //add one node onto end of domain to help with zero gradient at outlet
    z[nz]=lz*1.05;
    nz++;

    flow.setupGrid(nz, z);

    // specify the objects to use to compute kinetic rates and 
    // transport properties

    Transport* tr = newTransportMgr("AVBP", &gas);
    flow.setTransport(*tr);
    flow.setKinetics(gas);
    flow.setPressure(pressure);

    //------- step 2: create the inlet  -----------------------

    Inlet1D inlet;

    inlet.setMoleFractions(DATA_PTR(x));
    doublereal mdot=uin*rho_in;
    inlet.setMdot(mdot);
    inlet.setTemperature(temp);


    //------- step 3: create the outlet  ---------------------

    Outlet1D outlet;

    //=================== create the container and insert the domains =====

    vector<Domain1D*> domains;
    domains.push_back(&inlet);
    domains.push_back(&flow);
    domains.push_back(&outlet);

    //    OneDim flamesim(domains);

    Sim1D flame(domains);

    //----------- Supply initial guess----------------------

    vector_fp locs;
    vector_fp value;

    locs.resize(3);
    value.resize(3);

    //ramp values from inlet to adiabatic flame conditions 
    //  over 70% of domain and then level off at equilibrium
    double z1=0.7;

    double uout;
    uout=inlet.mdot()/rho_out;
    uin=inlet.mdot()/rho_in;
    locs[0]=0.0; locs[1]=z1; locs[2]=1.0;
    value[0]=uin; value[1]=uout; value[2]=uout;
    flame.setInitialGuess("u",locs,value);

    value[0]=temp; value[1]=Tad; value[2]=Tad;
    flame.setInitialGuess("T",locs,value);

    for(i=0;i<nsp;i++){
      value[0]=yin[i]; value[1]=yout[i]; value[2]=yout[i];
      flame.setInitialGuess(gas.speciesName(i),locs,value);
    }
//lecture solution initiale:
    flame.restore("flame_0.9_2.2atm.xml","2");
    flow.setPressure(pressure);

    inlet.setMoleFractions(DATA_PTR(x));
    inlet.setMdot(mdot);
    inlet.setTemperature(temp);

    flame.showSolution();

    int flowdomain=1;
    double ratio=1.1;
    double slope=0.05;
    double curve=0.05;
    double prune=-0.05;

    int npmax=600;
//    doublereal timestep=1.0e-7;
    int num_steps=4;
    integer* tsteps=new int[num_steps];
    tsteps[0]=1;
    tsteps[1]=2;
    tsteps[2]=5;
    tsteps[3]=10;
//    doublereal tmin=1.0e-8;
//    doublereal tmax=1.0e-2;

 
    flame.setMaxGridPoints(flowdomain,npmax);
    flame.setRefineCriteria(flowdomain,ratio,slope,curve,prune);



    int loglevel=1; 
    bool refine_grid = true;


    /* Solve freely propagating flame*/

    /* Linearally interpolate to find location where this
       temperature would exist. The temperature at this
       location will then be fixed for remainder of
       calculation.*/

    flow.fixTemperature();
    refine_grid=false;
    flame.setFixedTemperature(900.0);
    flame.solve(loglevel,refine_grid);
    refine_grid = true;
    flow.solveEnergyEqn();
    flame.solve(loglevel,refine_grid);

    flow.setTolerancesTS(1.0e-2,1.0e-3);
    flow.fixTemperature();
    refine_grid=false;
    flame.setFixedTemperature(900.0);
    //    flame.setAdiabaticFlame();
    flame.solve(loglevel,refine_grid);
    refine_grid = true;
    flow.solveEnergyEqn();
    flame.solve(loglevel,refine_grid);
    flame.setMaxGridPoints(flowdomain, npmax);
    flow.setTolerancesTS(1.0e-4,1.0e-5);
    flame.solve(loglevel,refine_grid);
    npmax=npmax+10;
    flame.setMaxGridPoints(flowdomain, npmax);
    flow.setTolerancesTS(5.0e-5,5.0e-6);
    flame.solve(loglevel,refine_grid);
    npmax=npmax+10;
    flame.setMaxGridPoints(flowdomain, npmax);
    flow.setTolerancesTS(1.0e-5,1.0e-6);
    flame.solve(loglevel,refine_grid);
    npmax=npmax+10;
    flame.setMaxGridPoints(flowdomain, npmax);
    flow.setTolerancesTS(5.0e-6,5.0e-7);
    flame.solve(loglevel,refine_grid);
    npmax=npmax+10;
    flame.setMaxGridPoints(flowdomain, npmax);
    flow.setTolerancesTS(1.0e-6,1.0e-7);
    flame.solve(loglevel,refine_grid);
    npmax=npmax+10;
    flame.setMaxGridPoints(flowdomain, npmax);
    flow.setTolerancesTS(5.0e-7,5.0e-8);




    int np=flow.nPoints();
// CALCULATE NEW SOLUTION WITH DIEFFERENT PHI, TEMP or PRESSION
    doublereal temp_new=temp;
    doublereal press_new=pressure;
    double phi_new=phi;
    myresult<<phi_new<<" "<<Tad<<" ";
    myresult<<flame.value(flowdomain,flow.componentIndex("u"),0)<<" ";
    myresult<<flame.value(flowdomain,flow.componentIndex("T"),np-1)<<" ";
    for(int n=0; n<nsp;n++)
      myresult<<flame.value(flowdomain,n+4,np-1)<<" ";
    myresult<<"\n";
	
    for(int k=0;k<0;k++){
      cout<<"++++++++++CALCULATE SOLUTION FOR DIFFERENT PHI+++++++++"<<endl;
      temp_new=temp_new;//*1.01;
      press_new=press_new;//+1.0*OneAtm;//*1.01;
      phi_new=phi_new-0.001;//
    
      nu_F=1.0;
      nu_O=5.0;
      fa_stoic=nu_O/nu_F;
      X_F=1/(1+4.76*fa_stoic/phi_new);
      X_O2=1/((phi_new/fa_stoic)+4.76);
      X_N2=3.76*X_O2;

      for(int k=0;k<nsp;k++){
        if(k==gas.speciesIndex("C3H8")){ x[k]=X_F; }
        else if(k==gas.speciesIndex("O2")){ x[k]=X_O2; }
        else if(k==gas.speciesIndex("N2")){ x[k]=X_N2; }
        else{ x[k]=0.0; }
      }



      inlet.setTemperature(temp_new);
      inlet.setMoleFractions(DATA_PTR(x));
      flow.setPressure(press_new);

      npmax=npmax+50;
      flame.setMaxGridPoints(flowdomain,npmax);
	
      flame.solve(loglevel,refine_grid);
      np=flow.nPoints();
      myresult<<phi_new<<" "<<Tad<<" ";
      myresult<<flame.value(flowdomain,flow.componentIndex("u"),0)<<" ";
      myresult<<flame.value(flowdomain,flow.componentIndex("T"),np-1)<<" ";
      for(int n=0; n<nsp;n++)
        myresult<<flame.value(flowdomain,n+4,np-1)<<" ";
      myresult<<"\n";

      Array2D soln(2*nsp+5+3*n_reac+2,1);
      for(int n=0;n<np;n++){
        flow.setGas(flame.solution()+flame.start(flowdomain),n);
        (n==0) ? saveSoln(0, flow.grid(n)+transl, flame.value(flowdomain, flow.componentIndex("u"),n), gas, soln, flow.AVBPthick(n)) : saveSoln(flow.grid(n)+transl, flame.value(flowdomain, flow.componentIndex("u"),n), gas, soln, flow.AVBPthick(n));
      }
      cout << endl<<"Adiabatic flame temperature from equilibrium is: "<<Tad<<endl;
      cout << "Flame speed for phi="<<phi_new<<" is "<< flame.value(flowdomain,flow.componentIndex("u"),0)<<" m/s."<<endl;

      char str_phi[4];
      sprintf(str_phi, "%1.2f", phi_new);
      string plotTitle = "1D laminar flame phi = ";
      plotTitle += str_phi;
      string filename, extension;

      filename = "soln_phi_new";
      extension = ".IGOR";                              // ".csv" ; ".IGOR"
      filename = filename + "-" + str_phi + extension;
      plotSoln(filename, "IGO", plotTitle, gas, soln); //  "XL"  ; "IGO"

    }
    np=flow.nPoints();
    vector<doublereal> zvec,Tvec,H2Ovec,CO2vec,Uvec;

    Array2D soln(2*nsp+5+3*n_reac+2,1);
    for(int n=0;n<np;n++){
      flow.setGas(flame.solution()+flame.start(flowdomain),n);
      (n==0) ? saveSoln(0, flow.grid(n)+transl, flame.value(flowdomain, flow.componentIndex("u"),n), gas, soln, flow.AVBPthick(n)) : saveSoln(flow.grid(n)+transl, flame.value(flowdomain, flow.componentIndex("u"),n), gas, soln, flow.AVBPthick(n));
    }


    //char str_phi[4];
    stringstream pn (stringstream::in | stringstream::out);
    pn << press_new;
    string str_press = pn.str();
    stringstream phn (stringstream::in | stringstream::out);
    phn << phi_new;
    string str_phi = phn.str();
    //sprintf(str_phi, "%1.2f", phi_new);
    string plotTitle = "1D laminar flame phi = ";
    plotTitle += str_phi;
    string filename, extension;
    string filename2, extension2;
    string filename3, extension3;

// sauvegarde solution finale pour restart
    flame.save("flame_" + str_phi + "_" + str_press + "atm.xml","2","Italiachampionnedumonde");

    filename = "soln_phi";
    filename3 = filename;
    extension3 = ".csv";
    extension = ".IGOR";                              // ".csv" ; ".IGOR"
    filename = filename + "-" + str_phi + extension;
    filename3 = filename3 + "-" + str_phi + extension3;
    plotSoln(filename, "IGO", plotTitle, gas, soln); //  "XL"  ; "IGO"
    plotSoln(filename3, "XL", plotTitle, gas, soln); //  "XL"  ; "IGO"

    filename2 = "Initial_solution_AVBP";
    extension2=".h5";
    filename2 = filename2 + "-" + str_phi + extension2;
    
     cout<<endl<<"+++++++++++++++Creating can2av.choices+++++++++++++++++++"<<endl;
     ofstream choices;
     choices.open("can2av.choices");
     choices<<"'./"<<filename3<<"'                  ! CANTERA last or restart file (1)"<<endl;
     choices<<"'./"<<filename2<<"'      ! The AVBP solution"<<endl;
     choices<<"'../MESH/flame.mesh.h5'                ! The AVBP mesh coordinate"<<endl;
     choices<<"'./input_species.dat'                  ! The AVBP species file"<<endl;
     choices<<"'./input_premix.dat'                   ! The AVBP premix file(the correct order of species is preserved automatically) "<<endl;
     choices<<"'./input_ifp.dat'                      ! The AVBP ifp file "<<endl;
     choices<<gas.pressure()<<"                                 ! Pressure of CANTERA flame"<<endl;
     choices<<nsp<<"                                      ! Number of species in CANTERA solution"<<endl;
     choices<<np<<"                                    ! Number of the grid points of CANTERA solution"<<endl;
     choices.close();
     cout<<"File can2av.choices created. Please check it out before run moulinette CAN2AV"<<endl;

     cout<<endl<<"Results for fresh gas temperature= "<<temp_new<<", pressure= "<<press_new<<" and phi= "<<phi_new<<":"<<endl;
     cout << "Adiabatic flame temperature from equilibrium is: "<<flame.value(flowdomain,flow.componentIndex("T"),np-1)<<endl;
     cout << "Flame speed is: "<<flame.value(flowdomain,flow.componentIndex("u"),0)<<" m/s."<<endl;
   
     return 0;    
   }
   catch (CanteraError) {
     showErrors(cerr);
     cerr << "program terminating." << endl;
     return -1;
   }
}

#ifndef CXX_DEMO
int main() {
    double p=1;    
    ofstream myresult;
    myresult.open("my_result.dat");
    int err;
    err= flamespeed(myresult,p);
   myresult.close();
   return err;

}
#endif