/*
   ___________           ______          _________
  |           |         |      |        |         |
  | PWR_UOX*4 |  =====> | Pool | =====> | Storage | 
  |___________|         |______|        |_________|
                                             ||
                                             ||
                                             \/
        __________       _________        __________________
       |         |      |         |      |                  |
       | Storage |<=====| PWR_URE |<=====| FabricationPlant |
       |_________|      |_________|      |__________________|
 
*/

// \rm CLASS* ; g++ -o CLASS_exec SimpleExample.cxx EQM_Simple.o -I $CLASS_include -L $CLASS_lib -lCLASSpkg `root-config --cflags` `root-config --libs` -fopenmp -lgomp -Wunused-result -lTMVA

#include "CLASSHeaders.hxx"
#include <sstream>
#include <iomanip>
#include <math.h>
#include <string>
#include "EQM_Simple.hxx"
#include "XS/XSM_CLOSEST.hxx"
#include "Irradiation/IM_RK4.hxx"
using namespace std;

int main(int argc, char** argv)
{  
  //All parameters;
  cSecond year=3600*24.*365.25;
  int Std_output_level=0;  // Only error are shown in terminal
  int File_output_level=2; // Error + Warning + Info are shown in the file CLASS_OUTPUT.log
  cSecond start=0.0;//Scenario start time (year);
  cSecond time=20.0;//Whole time of the scenario (year);
  cSecond coolingtime=5.0;//(year);
  cSecond fabricationtime=1.0;//(year);
  cSecond startreactor[]={1.,1.,1.,1.,15.};//(year);
  cSecond life[]={40.,40.,40.,40.,40.};//(year);
  double HMM[]={63.833491,63.833491,63.833491,63.833491,63.833491};//(ton);
  double power[]={40.17,40.17,40.17,40.17,40.17};//Specific power (MW/tHM);
  double burnup[]={33.5,33.5,33.5,33.5,28.790592};//(GWd/tHM);
  double loadfactor[]={0.85,0.85,0.85,0.85,0.85};
  int n=5;//Number of reactors;
  int i;
  
  //Change all parameters into the good units;
  start=start*year;
  time=time*year;
  coolingtime=coolingtime*year;
  fabricationtime=fabricationtime*year;
  for(i=0;i<=n-1;i++)
  {
    startreactor[i]=startreactor[i]*year;
    life[i]=life[i]*year;
    power[i]=power[i]*HMM[i]*1.0e6;
  }
  
  //Log management;
  CLASSLogger* Logger=new CLASSLogger("CLASS_OUTPUT.log",Std_output_level,File_output_level);

  //Scenario;
  Scenario* gCLASS=new Scenario(start,Logger);
  gCLASS->SetStockManagement(true);
  gCLASS->SetTimeStep(year/12.0);
  gCLASS->SetOutputFileName("CLASS.root");

  //Decay database;
  DecayDataBank* DecayDB=new DecayDataBank(gCLASS->GetLog(),"/home/yang/CLASS_4.1.0/DATA_BASES/DECAY/ALL/Decay.idx");
  gCLASS->SetDecayDataBase(DecayDB);
  
  //Fixed fuel reactors database;
  EvolutionData *database=new EvolutionData(Logger,"/home/yang/Results/RefFuel/UOX_3.2_40.17_33.5.dat");
  
  //Physics models;
  EQM_Simple *eqm=new EQM_Simple();
  XSM_CLOSEST *xsm=new XSM_CLOSEST(gCLASS->GetLog(),"/home/yang/Results/SymbioticScenarios/SimpleExample/SimpleXSM.idx");
  IM_RK4 *im=new IM_RK4(Logger);
  PhysicsModels *phym=new PhysicsModels(xsm,eqm,im);
  
  //Storages;
  Storage *stock=new Storage(Logger);
  stock->SetName("Stock_UOX");
  Storage *stockru=new Storage(Logger);
  stockru->SetName("Stock_RU");
  gCLASS->Add(stock);
  gCLASS->Add(stockru);
/*  Storage *Fissile=new Storage(Logger);
  Fissile->SetName("FissileStorage");
  gCLASS->Add(Fissile);
  Storage *Fertile=new Storage(Logger);
  Fertile->SetName("FertileStorage");
  gCLASS->Add(Fertile);
*/    
/*  //SeparationPlant;
  SeparationPlant *sp=new SeparationPlant(Logger);
  IsotopicVector fissile;
  fissile.Add(92,234,0,1.);
  fissile.Add(92,235,0,1.);
  fissile.Add(92,236,0,1.);
  IsotopicVector fertile;
  fertile.Add(92,238,0,1.);
  sp->SetBackEndDestination(Fissile,fissile,5*year);
  sp->SetBackEndDestination(Fertile,fertile,5*year);
*/  
  //Spent fuel pool;  
  Pool *pool=new Pool(Logger,stock,coolingtime);
  gCLASS->Add(pool);
  
  //FabricationPlant;
  FabricationPlant *fp=new FabricationPlant(gCLASS->GetLog(),fabricationtime);
  fp->SetFiFo(false);
  fp->AddFissileStorage(stock);
//  fp->AddFertileStorage(Fertile);
  fp->SetName("FP_RU");
  gCLASS->Add(fp);
   
  //Fixed fuel reactors;
  Reactor *uox1=new Reactor(Logger,
                            database,
                            pool,
                            startreactor[0],
                            life[0],
                            power[0],
                            HMM[0],
                            burnup[0],
                            loadfactor[0]);
  uox1->SetName("PWR_UOX_1");                            
  gCLASS->AddReactor(uox1);
  
  Reactor *uox2=new Reactor(Logger,
                            database,
                            pool,
                            startreactor[1],
                            life[1],
                            power[1],
                            HMM[1],
                            burnup[1],
                            loadfactor[1]);
  uox2->SetName("PWR_UOX_2");
  gCLASS->AddReactor(uox2);

  Reactor *uox3=new Reactor(Logger,
                            database,
                            pool,
                            startreactor[2],
                            life[2],
                            power[2],
                            HMM[2],
                            burnup[2],
                            loadfactor[2]);
  uox3->SetName("PWR_UOX_3");
  gCLASS->AddReactor(uox3);

  Reactor *uox4=new Reactor(Logger,
                            database,
                            pool,
                            startreactor[3],
                            life[3],
                            power[3],
                            HMM[3],
                            burnup[3],
                            loadfactor[3]);
  uox4->SetName("PWR_UOX_4");
  gCLASS->AddReactor(uox4);
  
  //Reprocessed fuel;
  Reactor *ru=new Reactor(Logger,
                          phym,
                          fp,
                          stockru,
                          startreactor[4],
                          life[4],
                          power[4],
                          HMM[4],
                          burnup[4],
                          loadfactor[4]);
  ru->SetName("PWR_RU");
  gCLASS->AddReactor(ru);                          
                  
  gCLASS->Evolution((double)(start+time));
  delete gCLASS;          

  return 1;
}