From a26b34674f2d74b710a9db59a7fb510bbd06bbe2 Mon Sep 17 00:00:00 2001
From: Baptiste Mouginot <mouginot.baptiste@gmail.com>
Date: Thu, 19 Jun 2014 08:13:09 +0000
Subject: [PATCH] Start Adding Irradiation Model, Equivalence Model and RK4
add HeavyMetalMass in EvolutionData
git-svn-id: svn+ssh://svn.in2p3.fr/class@261 0e7d625b-0364-4367-a6be-d5be4a48d228
---
.../CLASSV3/include/EquivalenceModel.hxx | 56 ++
.../CLASSV3/include/EvolutionData.hxx | 3 +
.../branches/CLASSV3/include/FuelDataBank.hxx | 3 +-
source/branches/CLASSV3/include/IM_RK4.hxx | 94 ++++
.../CLASSV3/include/IrradiationModel.hxx | 209 +++++++
.../branches/CLASSV3/src/EquivalenceModel.cxx | 0
source/branches/CLASSV3/src/EvolutionData.cxx | 17 +
source/branches/CLASSV3/src/IM_RK4.cxx | 441 +++++++++++++++
.../branches/CLASSV3/src/IrradiationModel.cxx | 512 ++++++++++++++++++
source/branches/CLASSV3/src/Makefile | 4 +-
10 files changed, 1336 insertions(+), 3 deletions(-)
create mode 100644 source/branches/CLASSV3/include/EquivalenceModel.hxx
create mode 100644 source/branches/CLASSV3/include/IM_RK4.hxx
create mode 100644 source/branches/CLASSV3/include/IrradiationModel.hxx
create mode 100644 source/branches/CLASSV3/src/EquivalenceModel.cxx
create mode 100644 source/branches/CLASSV3/src/IM_RK4.cxx
create mode 100644 source/branches/CLASSV3/src/IrradiationModel.cxx
diff --git a/source/branches/CLASSV3/include/EquivalenceModel.hxx b/source/branches/CLASSV3/include/EquivalenceModel.hxx
new file mode 100644
index 000000000..d44e12b87
--- /dev/null
+++ b/source/branches/CLASSV3/include/EquivalenceModel.hxx
@@ -0,0 +1,56 @@
+#ifndef _EQUIVALENCEMODEL_HXX
+#define _EQUIVALENCEMODEL_HXX
+
+
+/*!
+ \file
+ \brief Header file for EquivalenceModel class.
+
+
+ @author BaM
+ @version 2.0
+ */
+
+
+using namespace std;
+
+//-----------------------------------------------------------------------------//
+/*!
+ Define a EquivalenceModel.
+ The aim of these class is synthetyse all the commum properties to all Equivalence Model.
+
+
+ @author BaM
+ @version 3.0
+ */
+//________________________________________________________________________
+
+class IsotopicVector;
+
+
+class EquivalenceModel : public TObject
+{
+ public :
+
+ /// virtueal method called to build a reprocessed fuel as a function of the burnup requierement the stock, mass....
+ /*!
+ Build the fuel following the equivalance model with the proper requierment in term of mass burnup....
+ \param double BurnUp desireted burnup reached by the fuel at the end of irradiation
+ \param double HMMass, needed Heavy metal mass needed
+ \param vector<double> &lambda, fraction of the stock to take (initialy should be 0)
+ \param vector<IsotopicVector> FissilArray, isotopicvectors to use to get the fissil part of the fuel
+ \param vector<IsotopicVector> FertilArray, isotopicvectors to use to get the fertil part of the fuel (if empty take it from the god)
+ */
+
+ virtual IsotopicVector BuildFuel(double BurnUp, double HMMass, vector<double> &lambda, vector<IsotopicVector> FissilArray, vector<IsotopicVector> FertilArray = 0);
+ //}
+
+
+
+
+ private :
+
+};
+
+#endif
+
diff --git a/source/branches/CLASSV3/include/EvolutionData.hxx b/source/branches/CLASSV3/include/EvolutionData.hxx
index dd3d15f12..6ecab70e6 100755
--- a/source/branches/CLASSV3/include/EvolutionData.hxx
+++ b/source/branches/CLASSV3/include/EvolutionData.hxx
@@ -115,6 +115,7 @@ public :
\name Set Method
*/
//@{
+ void SetHeavyMetalMass(double Mass) {fHeavyMetalMass = Mass;} //!< Set the HeavyMetal Mass
void SetReactorType(string reactortype) { fReactorType = reactortype; } ///< Set the reactor Type (string)
void SetFuelType(string fueltype) { fFuelType = fueltype; } ///< Set the fuel Type (string)
@@ -158,6 +159,7 @@ public :
IsotopicVector GetIsotopicVectorAt(double t); ///< Return the Product IsotopicVector at t time
+ double GetHeavyMetalMass() const { return fHeavyMetalMass; } //!< Return the HeavyMetal Mass in the Core at the begining of the cycle
//{
@@ -222,6 +224,7 @@ protected :
string fFuelType; ///< Type of fuel
double fPower; ///< Power in W
double fCycleTime; ///< Cycle time of the DataBase
+ double fHeavyMetalMass; ///< Cycle time of the DataBase
double fNormFactor; ///< Normalisation factor needed to represent to full core (unsless)
diff --git a/source/branches/CLASSV3/include/FuelDataBank.hxx b/source/branches/CLASSV3/include/FuelDataBank.hxx
index e3e3134be..741501bcd 100644
--- a/source/branches/CLASSV3/include/FuelDataBank.hxx
+++ b/source/branches/CLASSV3/include/FuelDataBank.hxx
@@ -198,8 +198,7 @@ class FuelDataBank : public CLASSObject, DynamicalSystem
*/
//@{
-
- IsotopicVector Evolution(IsotopicVector IV, double dt); ///< Return the Product IsotopicVector evolution from zai during a dt time
+
void CalculateDistanceParameter(); ///< Calcul of the weight for each ZAI in the distance calculation from the mean XS of the FuelDataBank
void BuildDecayMatrix(); ///w Build the Decay Matrix for the futur evolution...
diff --git a/source/branches/CLASSV3/include/IM_RK4.hxx b/source/branches/CLASSV3/include/IM_RK4.hxx
new file mode 100644
index 000000000..875847bcd
--- /dev/null
+++ b/source/branches/CLASSV3/include/IM_RK4.hxx
@@ -0,0 +1,94 @@
+#ifndef _IMRK4_HXX
+#define _IMRK4_HXX
+
+
+/*!
+ \file
+ \brief Header file for IrradiationModel class.
+
+
+ @author BaM
+ @version 2.0
+ */
+
+
+using namespace std;
+
+//-----------------------------------------------------------------------------//
+/*!
+ Define a IM_RK4.
+ The aim of these class is perform the calculation of the evolution of a fuel trough irradiation solving numericaly the Bateman using RK4.
+
+
+ @author BaM
+ @version 3.0
+ */
+//________________________________________________________________________
+
+class EvolutionData;
+
+class IM_RK4 : public TObject, DynamicalSystem
+{
+
+ public :
+
+ /// virtueal method called to perform the irradiation calculation using a set of cross section.
+ /*!
+ Perform the Irradiation Calcultion using the XSSet data
+ \param IsotopicVector IV isotopic vector to irradiate
+ \param EvolutionData XSSet set of corss section to use to perform the evolution calculation
+ */
+ EvolutionData GenerateEvolutionData(IsotopicVector IV, EvolutionData XSSet, double Power);
+ //}
+
+ //********* RK4 Method *********//
+
+ //@}
+ /*!
+ \name RK4 Method
+ */
+ //@{
+
+ void UseRK4EvolutionMethod(bool usemethod = true) {fUseRK4EvolutionMethod = usemethod;}
+
+
+ using DynamicalSystem::RungeKutta;
+ //! Pre-treatment Runge-Kutta method.
+ /*!
+ // This method does initialisation and then call DynamicalSystem::RungeKutta
+ // \param t1: initial time
+ // \param t2: final time
+ */
+
+
+ void BuildEqns(double t, double *N, double *dNdt);
+ void SetTheMatrixToZero(); //!< Initialize the evolution Matrix
+ void ResetTheMatrix();
+ void SetTheMatrix(TMatrixT<double> BatemanMatrix); //!< Set the Evolution Matrix (Bateman equations)
+ TMatrixT<double> GetTheMatrix(); //!< return the Evolution Matrix (Bateman equations)
+
+ void SetTheNucleiVectorToZero(); //!< Initialize the evolution Matrix
+ void ResetTheNucleiVector();
+ void SetTheNucleiVector(TMatrixT<double> NEvolutionMatrix); //!< Set the Evolution Matrix (Bateman equations)
+ TMatrixT<double> GetTheNucleiVector(); //!< return the Evolution Matrix (Bateman equations)
+ //@}
+
+
+
+ private :
+
+ double *fTheNucleiVector; //!< The evolving atoms copied from Material proportions.
+ double **fTheMatrix; //!< The evolution Matrix
+
+ double fPrecision; //!< Precision of the RungeKutta
+ double fHestimate; //!< RK Step estimation.
+ double fHmin; //!< RK minimum Step.
+ double fMaxHdid; //!< store the effective RK max step
+ double fMinHdid; //!< store the effective RK min step
+ bool fIsNegativeValueAllowed; //!< whether or not negative value are physical.
+
+
+};
+
+#endif
+
diff --git a/source/branches/CLASSV3/include/IrradiationModel.hxx b/source/branches/CLASSV3/include/IrradiationModel.hxx
new file mode 100644
index 000000000..138b83950
--- /dev/null
+++ b/source/branches/CLASSV3/include/IrradiationModel.hxx
@@ -0,0 +1,209 @@
+#ifndef _IRRADIATIONMODEL_HXX
+#define _IRRADIATIONMODEL_HXX
+
+
+/*!
+ \file
+ \brief Header file for IrradiationModel class.
+
+
+ @author BaM
+ @version 2.0
+ */
+
+
+using namespace std;
+
+//-----------------------------------------------------------------------------//
+/*!
+ Define a IrradiationModel.
+ The aim of these class is synthetyse all the commum properties to all Irradiation Model.
+
+
+ @author BaM
+ @version 3.0
+ */
+//________________________________________________________________________
+
+class EvolutionData;
+
+class IrradiationModel : public TObject
+{
+
+ public :
+
+ /// virtueal method called to perform the irradiation calculation using a set of cross section.
+ /*!
+ Perform the Irradiation Calcultion using the XSSet data
+ \param IsotopicVector IV isotopic vector to irradiate
+ \param EvolutionData XSSet set of corss section to use to perform the evolution calculation
+ \param double Power, constant power to use for irradation
+ \param double irradiationtime, time of the irradiation
+ */
+ virtual EvolutionData GenerateEvolutionData(IsotopicVector IV, EvolutionData XSSet, double Power);
+ //}
+
+
+
+
+ //********* Get Method *********//
+ /*!
+ \name Get Method
+ */
+ //@{
+
+ double GetShorstestHalflife() const { return fShorstestHalflife; }
+
+ //@}
+
+
+
+
+ //********* Set Method *********//
+
+ /*!
+ \name Set Method
+ */
+ //@{
+
+ void SetFuelDataBank(map< IsotopicVector ,EvolutionData > mymap) { fFuelDataBank = mymap; } //!< Set the FuelDataBank map
+
+ void SetDataBaseIndex(string database) { fDataBaseIndex = database;; ReadDataBase(); } //!< Set the Name of the database index
+
+
+
+ //{
+ /// set Fission Energy using a file
+ /*!
+ // This method fill the Fission Energy map using a file
+ // \param FissionEnergyFile: filename containing the Fission Energy of some nuclei (form : Z A I Energy)
+ */
+ void SetFissionEnergy(string FissionEnergyFile);
+ //}
+
+ //{
+ /// set Fission Energy for a ZAI using ZAI(Z,A,I)
+ /*!
+ // This method fill the Fission Energy map of a set ZAI
+ // \param zai ZAI
+ // \param E Fission energy of the ZAI
+ */
+ void SetFissionEnergy(ZAI zai, double E);
+ //}
+
+ //{
+ /// set Fission Energy for a ZAI using the Z, A, I
+ /*!
+ // This method fill the Fission Energy map of a set ZAI
+ // \param Z Z of the ZAI
+ // \param A A of the ZAI
+ // \param I I of the ZAI
+ // \param E Fission energy of the ZAI
+ */
+ void SetFissionEnergy(int Z, int A, int I, double E ) { SetFissionEnergy(ZAI(Z,A,I), E);}
+ //}
+
+ void SetShortestHalfLife(double halflife) { fShorstestHalflife = halflife;} ///< Set the Half Life cut
+ void LoadFPYield(string SponfaneusYield, string ReactionYield); ///< Build Fision Yields maps;
+
+
+
+
+
+
+ //********* Evolution Method *********//
+
+ //@}
+ /*!
+ \name Evolution Method
+ */
+ //@{
+
+
+ void BuildDecayMatrix(); ///w Build the Decay Matrix for the futur evolution...
+
+ //@}
+
+
+ //********* Other Method *********//
+ /*!
+ \name Other Method
+ */
+ //@{
+ void Print() const;
+
+ //@}
+
+
+
+
+
+ protected :
+
+ double fShorstestHalflife;
+ int fZAIThreshold; //!< Lowest Mass deal by the evolution (default 90)
+
+
+ TMatrixT<double> fDecayMatrix; ///< Matrix with half life of each nuclei
+ map<ZAI, double > fFissionEnergy; ///< Store the Energy per fission use for the flux normalisation.
+ map<ZAI, map<ZAI, double> > fFastDecay; ///< Store the cut decay
+ map<ZAI, IsotopicVector> fSpontaneusYield; ///< Store the Spontaneus fission yield
+ map<ZAI, IsotopicVector> fReactionYield; ///< Store the reaction fission yield
+
+
+ int fNVar; //!< The size of the composition vector and /or number of ZAIs involved.
+
+
+ map<ZAI, int> findex_inver; ///< correspondance matrix from ZAI to the column (or line) of the different Reaction/Decay matrix
+ map<int, ZAI> findex; ///< correspondance matrix from the column (or line) of the different Reaction/Decay matrix to the ZAI
+
+ //{
+ /// Return the Fission XS Matrix at the time TStep
+ /*!
+ // This Method extract the Fission Cross section of an EvolutionData at the set time
+ // \param EvolutionDataStep: EvolutionData
+ // \param TStep: time
+ */
+ TMatrixT<double> GetFissionXsMatrix(EvolutionData EvolutionDataStep,double TStep);
+ //}
+
+ //{
+ /// Return the Capture XS Matrix at the time TStep
+ /*!
+ // This Method extract the capture Cross section of an EvolutionData at the set time
+ // \param EvolutionDataStep: EvolutionData
+ // \param TStep: time
+ */
+ TMatrixT<double> GetCaptureXsMatrix(EvolutionData EvolutionDataStep,double TStep);
+ //}
+
+ //{
+ /// Return the n2n XS Matrix at the time TStep
+ /*!
+ // This Method extract the (n,2n) Cross section of an EvolutionData at the set time
+ // \param EvolutionDataStep: EvolutionData
+ // \param TStep: time
+ */
+ TMatrixT<double> Getn2nXsMatrix(EvolutionData EvolutionDataStep,double TStep);
+ //}
+
+
+ //{
+ //! Returns a particular decay mode.
+ /*!
+ \param DecayModes : a list of decay modes with their branching ratios and isomeric state of the Daughters.
+ \param BR : branching ratio of the current decay mode
+ \param Iso : isomeric state of the Daughter of the current decay mode.
+ \param StartPos : the current decay mode to extract.
+ */
+ string GetDecay(string DecayModes, double &BR,int &Iso, int &StartPos);
+ //}
+
+ map< ZAI,IsotopicVector > ReadFPYield(string Yield); ///< Read a CLASSYield file and return the correpsponding map
+
+ private :
+
+};
+
+#endif
+
diff --git a/source/branches/CLASSV3/src/EquivalenceModel.cxx b/source/branches/CLASSV3/src/EquivalenceModel.cxx
new file mode 100644
index 000000000..e69de29bb
diff --git a/source/branches/CLASSV3/src/EvolutionData.cxx b/source/branches/CLASSV3/src/EvolutionData.cxx
index fd44aca89..dfce0182b 100755
--- a/source/branches/CLASSV3/src/EvolutionData.cxx
+++ b/source/branches/CLASSV3/src/EvolutionData.cxx
@@ -660,6 +660,23 @@ void EvolutionData::ReadDB(string DBfile, bool oldread)
}while ( !DecayDB.eof() );
+ double M = 0;
+ {
+ map<ZAI, double >::iterator it ;
+
+
+ IsotopicVector IVtmp = GetIsotopicVectorAt(0.).GetActinidesComposition();
+ map<ZAI, double >isotopicquantity = IVtmp.GetIsotopicQuantity();
+
+ for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
+ {
+ M += isotopicvector.GetActinidesComposition().GetZAIIsotopicQuantity( (*it).first )*cZAIMass.fZAIMass.find( (*it).first )->second/Na*1e-6;
+ }
+ isotopicquantity.clear();
+
+ }
+ fHeavyMetalMass = M;
+
DecayDB.close();
diff --git a/source/branches/CLASSV3/src/IM_RK4.cxx b/source/branches/CLASSV3/src/IM_RK4.cxx
new file mode 100644
index 000000000..f40faa5a9
--- /dev/null
+++ b/source/branches/CLASSV3/src/IM_RK4.cxx
@@ -0,0 +1,441 @@
+//
+// IM_RK4.cxx
+// CLASSSource
+//
+// Created by BaM on 04/05/2014.
+// Copyright (c) 2014 BaM. All rights reserved.
+//
+
+#include "IM_RK4.hxx"
+
+#include "IsotopicVector.hxx"
+#include "CLASSHeaders.hxx"
+#include "LogFile.hxx"
+#include "StringLine.hxx"
+
+#include <TGraph.h>
+#include <TString.h>
+
+
+#include <sstream>
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <algorithm>
+#include <cmath>
+
+
+
+
+
+//________________________________________________________________________
+IM_RK4::IM_RK4(): CLASSObject(), DynamicalSystem()
+{
+ fTheNucleiVector = 0;
+ fTheMatrix = 0;
+
+ fWeightedDistance = false;
+ fEvolutionDataInterpolation = false;
+
+
+
+ fOldReadMethod = true;
+ fUseRK4EvolutionMethod = true;
+ fDistanceType = 0;
+ fShorstestHalflife = 3600.*24*2.;
+ fZAIThreshold = 90;
+
+
+ fDataDirectoryName = getenv("CLASS_PATH");
+ fDataDirectoryName += "/source/data/";
+ fDataFileName = "chart.JEF3T";
+
+ SetForbidNegativeValue();
+
+}
+
+
+IM_RK4::IM_RK4(LogFile* Log):DynamicalSystem()
+{
+ SetLog(Log);
+ fWeightedDistance = false;
+ fEvolutionDataInterpolation = false;
+
+
+ fTheNucleiVector = 0;
+ fTheMatrix = 0;
+
+ fDataDirectoryName = getenv("CLASS_PATH");
+ fDataDirectoryName += "/source/data/";
+ fDataFileName = "chart.JEF3T";
+
+
+ fShorstestHalflife = 3600.*24*2.;
+ fZAIThreshold = 90;
+
+
+ SetForbidNegativeValue();
+
+}
+
+
+
+
+
+
+
+
+
+//________________________________________________________________________
+/* Evolution Calculation */
+//________________________________________________________________________
+EvolutionData IM_RK4::GenerateEvolutionData(IsotopicVector isotopicvector, EvolutionData XSSet, double Power);
+{
+
+ if(fFastDecay.size() == 0)
+ {
+ BuildDecayMatrix();
+ fNVar = findex_inver.size();
+ }
+
+ SetTheMatrixToZero();
+ SetTheNucleiVectorToZero();
+
+ string ReactorType;
+
+
+ vector< TMatrixT<double> > NMatrix ;// TMatrixT<double>(decayindex.size(),1))
+ { // Filling the t=0 State;
+ map<ZAI, double > isotopicquantity = isotopicvector.GetIsotopicQuantity();
+ TMatrixT<double> N_0Matrix = TMatrixT<double>( findex.size(),1) ;
+ for(int i = 0; i < (int)findex.size(); i++)
+ N_0Matrix[i] = 0;
+
+ map<ZAI, double >::iterator it ;
+ for(int i = 0; i < (int)findex.size(); i++)
+ N_0Matrix[i] = 0;
+
+ for(it = isotopicquantity.begin(); it != isotopicquantity.end(); it++)
+ {
+ /// Need TO change with FP managment
+ map<ZAI, int >::iterator it2;
+
+ if( (*it).first.Z() < fZAIThreshold )
+ it2 = findex_inver.find( ZAI(-2,-2,-2) );
+ else it2 = findex_inver.find( (*it).first );
+
+ if(it2 == findex_inver.end() ) //If not in index should be TMP, can't be fast decay for new Fuel !!!
+ it2 = findex_inver.find( ZAI(-3,-3,-3) );
+ N_0Matrix[ (*it2).second ][0] = (*it).second ;
+
+
+ }
+
+ isotopicquantity.clear();
+
+ NMatrix.push_back(N_0Matrix);
+ N_0Matrix.Clear();
+
+ }
+
+
+ //-------------------------//
+ //--- Perform Evolution ---//
+ //-------------------------//
+ ReactorType = XSSet.GetReactorType();
+
+ double Na = 6.02214129e23; //N Avogadro
+ double M_ref = XSSet.GetHeavyMetalMass();
+ double M = 0;
+ double Power_ref = XSSet.GetPower();
+
+ // Get the mass of the fuel to irradiate in order to perform the evolution at fixed burnup (the burnup step of the calculation will match the burnup step of the XSSet
+ {
+ map<ZAI, double >::iterator it ;
+
+
+ IsotopicVector IVtmp = isotopicvector.GetActinidesComposition();
+ map<ZAI, double >isotopicquantity = IVtmp.GetIsotopicQuantity();
+
+ for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
+ M += isotopicvector.GetActinidesComposition().GetZAIIsotopicQuantity( (*it).first )*cZAIMass.fZAIMass.find( (*it).first )->second/Na*1e-6;
+ isotopicquantity.clear();
+
+ }
+
+ int NStep = = XSSet.GetFissionXS().begin()->second->GetN();
+ double* DBTimeStep = XSSet.GetFissionXS().begin()->second->GetX();
+
+ int InsideStep = 10;
+
+ double timevector[NStep];
+ timevector[0] = 0;
+
+ double Flux[NStep];
+
+ TMatrixT<double> FissionEnergy = TMatrixT<double>(findex.size(),1);
+ for(int i = 0; i < (int)findex.size(); i++)
+ FissionEnergy[i] = 0;
+
+ {
+ map< ZAI, int >::iterator it;
+ for(it = findex_inver.begin(); it != findex_inver.end(); it++)
+ {
+ map< ZAI, double >::iterator it2 = fFissionEnergy.find(it->first);
+ if(it2 == fFissionEnergy.end())
+ {
+ if(it->first.Z() > fZAIThreshold)
+ FissionEnergy[it->second][0] = 1.9679e6*it->first.A()-2.601e8; // //simple linear fit to known values ;extrapolation to unknown isotopes
+ else FissionEnergy[it->second][0] = 0;
+ }
+ else
+ FissionEnergy[it->second][0] = it2->second;
+
+ }
+ }
+
+ vector< TMatrixT<double> > FissionXSMatrix; // Store The Fisison XS Matrix
+ vector< TMatrixT<double> > CaptureXSMatrix; // Store The Capture XS Matrix
+ vector< TMatrixT<double> > n2nXSMatrix; // Store The n2N XS Matrix
+
+ for(int i = 0; i < NStep-1; i++)
+ {
+ double TStepMax = ( (DBTimeStep[i+1]-DBTimeStep[i] ) ) * Power_ref/M_ref / Power*M ; // Get the next Time step
+
+
+ TMatrixT<double> BatemanMatrix = TMatrixT<double>(findex.size(),findex.size());
+ TMatrixT<double> BatemanReactionMatrix = TMatrixT<double>(findex.size(),findex.size());
+
+ TMatrixT<double> NEvolutionMatrix = TMatrixT<double>(findex.size(),1);
+ NEvolutionMatrix = NMatrix.back();
+
+
+
+ FissionXSMatrix.push_back(GetFissionXsMatrix(XSSet, DBTimeStep[i])); //Feel the fission reaction Matrix
+ CaptureXSMatrix.push_back(GetCaptureXsMatrix(XSSet, DBTimeStep[i])); //Feel the capture reaction Matrix
+ n2nXSMatrix.push_back(Getn2nXsMatrix(XSSet, DBTimeStep[i])); //Feel the (n,2n) reaction Matrix
+
+ // ---------------- Evolution
+
+ BatemanReactionMatrix = FissionXSMatrix[i];
+ BatemanReactionMatrix += CaptureXSMatrix[i];
+ BatemanReactionMatrix += n2nXSMatrix[i];
+
+ if(fUseRK4EvolutionMethod)
+ {
+ for(int k=0; k < InsideStep; k++)
+ {
+ double ESigmaN = 0;
+ for (int j = 0; j < (int)findex.size() ; j++)
+ ESigmaN -= FissionXSMatrix[i][j][j]*NEvolutionMatrix[j][0]*1.6e-19*FissionEnergy[j][0];
+ // Update Flux
+ double Flux_k = Power/ESigmaN;
+
+ if(k==0)
+ Flux[i]=Flux_k;
+
+ BatemanMatrix = BatemanReactionMatrix;
+ BatemanMatrix *= Flux_k;
+ BatemanMatrix += fDecayMatrix ;
+
+ SetTheMatrixToZero();
+ SetTheNucleiVectorToZero();
+
+ SetTheMatrix(BatemanMatrix);
+ SetTheNucleiVector(NEvolutionMatrix);
+
+
+ RungeKutta(fTheNucleiVector, timevector[i]+TStepMax/InsideStep*k, timevector[i]+TStepMax/InsideStep*(k+1), fNVar);
+ NEvolutionMatrix = GetTheNucleiVector();
+
+ }
+ NEvolutionMatrix = GetTheNucleiVector();
+ NMatrix.push_back(NEvolutionMatrix);
+ }
+
+ timevector[i+1] = timevector[i] + TStepMax;
+
+ BatemanMatrix.Clear();
+ BatemanReactionMatrix.Clear();
+ NEvolutionMatrix.Clear();
+
+
+ }
+ FissionXSMatrix.push_back(GetFissionXsMatrix(XSSet, DBTimeStep[NStep-1])); //Feel the reaction Matrix
+ CaptureXSMatrix.push_back(GetCaptureXsMatrix(XSSet, DBTimeStep[NStep-1])); //Feel the reaction Matrix
+ n2nXSMatrix.push_back(Getn2nXsMatrix(XSSet, DBTimeStep[NStep-1])); //Feel the reaction Matrix
+
+
+ EvolutionData GeneratedDB = EvolutionData(GetLog());
+
+ double ESigmaN = 0;
+ for (int j = 0; j < (int)findex.size() ; j++)
+ ESigmaN -= FissionXSMatrix.back()[j][j]*NMatrix.back()[j][0]*1.6e-19*FissionEnergy[j][0];
+
+ Flux[NStep-1] = Power/ESigmaN;
+
+ GeneratedDB.SetFlux( new TGraph(NStep, timevector, Flux) );
+
+ for(int i = 0; i < (int)findex.size(); i++)
+ {
+ double ZAIQuantity[NMatrix.size()];
+ double FissionXS[NStep];
+ double CaptureXS[NStep];
+ double n2nXS[NStep];
+ for(int j = 0; j < (int)NMatrix.size(); j++)
+ ZAIQuantity[j] = (NMatrix[j])[i][0];
+
+ for(int j = 0; j < NStep; j++)
+ {
+ FissionXS[j] = FissionXSMatrix[j][i][i];
+ CaptureXS[j] = CaptureXSMatrix[j][i][i];
+ n2nXS[j] = n2nXSMatrix[j][i][i];
+ }
+
+ GeneratedDB.NucleiInsert(pair<ZAI, TGraph*> (findex.find(i)->second, new TGraph(NMatrix.size(), timevector, ZAIQuantity)));
+ GeneratedDB.FissionXSInsert(pair<ZAI, TGraph*> (findex.find(i)->second, new TGraph(NStep, timevector, FissionXS)));
+ GeneratedDB.CaptureXSInsert(pair<ZAI, TGraph*> (findex.find(i)->second, new TGraph(NStep, timevector, CaptureXS)));
+ GeneratedDB.n2nXSInsert(pair<ZAI, TGraph*> (findex.find(i)->second, new TGraph(NStep, timevector, n2nXS)));
+ }
+
+ GeneratedDB.SetPower(Power );
+ GeneratedDB.SetFuelType(fFuelType );
+ GeneratedDB.SetReactorType(ReactorType );
+ GeneratedDB.SetCycleTime(cycletime);
+
+ ResetTheMatrix();
+ ResetTheNucleiVector();
+
+ for (int i = 0; i < (int) FissionXSMatrix.size(); i++)
+ {
+ FissionXSMatrix[i].Clear();
+ CaptureXSMatrix[i].Clear();
+ n2nXSMatrix[i].Clear();
+ }
+ FissionXSMatrix.clear();
+ CaptureXSMatrix.clear();
+ n2nXSMatrix.clear();
+
+ return GeneratedDB;
+
+}
+
+
+void IM_RK4::ResetTheMatrix()
+{
+
+ if(fTheMatrix)
+ {
+ for(int i= 0; i<fNVar; i++)
+ delete [] fTheMatrix[i];
+ delete [] fTheMatrix;
+ }
+ fTheMatrix = 0;
+}
+
+void IM_RK4::SetTheMatrixToZero()
+{
+ ResetTheMatrix();
+
+ fNVar = findex.size();
+ fTheMatrix = new double*[fNVar];
+
+#pragma omp parallel for
+ for(int i= 0; i < fNVar; i++)
+ fTheMatrix[i] = new double[fNVar];
+
+ for(int i = 0; i < fNVar; i++)
+ for(int k = 0; k < fNVar; k++)
+ {
+ fTheMatrix[i][k]=0.0;
+ }
+
+}
+
+//________________________________________________________________________
+void IM_RK4::ResetTheNucleiVector()
+{
+ if(fTheNucleiVector)
+ delete [] fTheNucleiVector;
+ fTheNucleiVector = 0;
+}
+
+//________________________________________________________________________
+void IM_RK4::SetTheNucleiVectorToZero()
+{
+ ResetTheNucleiVector();
+ fTheNucleiVector = new double[fNVar];
+
+#pragma omp parallel for
+ for(int i = 0; i < fNVar; i++)
+ fTheNucleiVector[i]=0.0;
+
+}
+
+//________________________________________________________________________
+void IM_RK4::BuildEqns(double t, double *N, double *dNdt)
+{
+ double sum=0;
+ // pragma omp parallel for reduction(+:sum)
+ for(int i = 0; i < fNVar; i++)
+ {
+ sum=0;
+ for(int k = 0; k < fNVar; k++)
+ {
+ sum += fTheMatrix[i][k]*N[k];
+ }
+ dNdt[i] = sum;
+ }
+}
+
+//________________________________________________________________________
+void IM_RK4::SetTheMatrix(TMatrixT<double> BatemanMatrix)
+{
+ for (int k = 0; k < (int)fNVar; k++)
+ for (int l = 0; l < (int)findex_inver.size(); l++)
+ fTheMatrix[l][k] = BatemanMatrix[l][k];
+}
+
+//________________________________________________________________________
+TMatrixT<double> IM_RK4::GetTheMatrix()
+{
+ TMatrixT<double> BatemanMatrix = TMatrixT<double>(findex.size(),findex.size());
+ for (int k = 0; k < (int)fNVar; k++)
+ for (int l = 0; l < (int)findex_inver.size(); l++)
+ BatemanMatrix[l][k] = fTheMatrix[l][k];
+
+ return BatemanMatrix;
+}
+
+//________________________________________________________________________
+void IM_RK4::SetTheNucleiVector(TMatrixT<double> NEvolutionMatrix)
+{
+ for (int k = 0; k < (int)fNVar; k++)
+ fTheNucleiVector[k] = NEvolutionMatrix[k][0];
+}
+
+//________________________________________________________________________
+TMatrixT<double> IM_RK4::GetTheNucleiVector()
+{
+ TMatrixT<double> NEvolutionMatrix = TMatrixT<double>(findex.size(),1);
+ for (int k = 0; k < (int)fNVar; k++)
+ NEvolutionMatrix[k][0] = fTheNucleiVector[k];
+
+ return NEvolutionMatrix;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/source/branches/CLASSV3/src/IrradiationModel.cxx b/source/branches/CLASSV3/src/IrradiationModel.cxx
new file mode 100644
index 000000000..a9acb3e9b
--- /dev/null
+++ b/source/branches/CLASSV3/src/IrradiationModel.cxx
@@ -0,0 +1,512 @@
+//
+// IrradiationModel.cxx
+// CLASSSource
+//
+// Created by BaM on 04/05/2014.
+// Copyright (c) 2014 BaM. All rights reserved.
+//
+
+#include "IrradiationModel.hxx"
+
+#include "IsotopicVector.hxx"
+#include "CLASSHeaders.hxx"
+#include "LogFile.hxx"
+#include "StringLine.hxx"
+
+#include <TGraph.h>
+#include <TString.h>
+
+
+#include <sstream>
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <algorithm>
+#include <cmath>
+
+
+
+
+
+
+
+
+
+
+
+
+
+//________________________________________________________________________
+//________________________________________________________________________
+/* Physics */
+//________________________________________________________________________
+//________________________________________________________________________
+
+
+
+//________________________________________________________________________
+/* Decay Stuff */
+//________________________________________________________________________
+string IrradiationModel::GetDecay(string DecayModes, double &BR,int &Iso, int &StartPos)
+{
+ string header;
+
+ BR=0;
+ //extraction of the decay mode and the BR
+ string DecayBR=StringLine::NextWord(DecayModes,StartPos,',');
+ //extraction of the decay
+ int ss=0;
+ string Decay=StringLine::NextWord(DecayBR,ss,':');
+ //extraction of the BR if exist (i.e. for non stable isotop)
+ if(ss<int(DecayBR.size()))
+ BR=atof(DecayBR.substr(ss+1).c_str());
+ //BR in % -> BR
+ BR/=100.;
+ //find the Isomeric state of Daughter
+ Iso=0;
+ if(Decay.find("/",0)<string::npos)
+ {
+ Iso=atoi(Decay.substr(Decay.find("/")+1).c_str());
+ Decay=Decay.substr(0,Decay.find("/"));
+ }
+ return Decay;
+}
+
+//________________________________________________________________________
+void IrradiationModel::BuildDecayMatrix()
+{
+ fDecayMatrix.Clear();
+
+ // List of Decay Time and Properties
+ map<ZAI, pair<double, map< ZAI, double > > > ZAIDecay;
+
+ { // TMP
+ map< ZAI, double > toAdd;
+ toAdd.insert(pair<ZAI, double> ( ZAI(-3,-3,-3) , 1) );
+ ZAIDecay.insert( pair< ZAI, pair<double, map< ZAI, double > > >( ZAI(-3,-3,-3), pair<double, map< ZAI, double > > ( 1e28 ,toAdd )) ) ;
+ }
+ { // PF
+ map< ZAI, double > toAdd;
+ toAdd.insert(pair<ZAI, double> ( ZAI(-2,-2,-2), 1) );
+ ZAIDecay.insert( pair< ZAI, pair<double, map< ZAI, double > > >( ZAI(-2,-2,-2), pair<double, map< ZAI, double > > ( 1e28 ,toAdd )) ) ;
+ }
+
+
+ string DataFullPathName = GetDataDirectoryName()+ GetDataFileName();
+ ifstream infile(DataFullPathName.c_str());
+
+ if(!infile)
+ {
+ cout << "!!Warning!! !!!IrradiationModel!!! \n Can't open \"" << DataFullPathName << "\"\n" << endl;
+ GetLog()->fLog << "!!Warning!! !!!IrradiationModel!!! \n Can't open \"" << DataFullPathName<< "\"\n" << endl;
+ }
+
+
+
+
+
+ do
+ {
+ int A = -1;
+ int Z = -1;
+ int I = -1;
+ string zainame;
+ string Iname;
+ int unkown;
+ double HalfLife;
+ string DecayModes;
+
+ infile >> A >> Z >> zainame >> Iname >> unkown >> HalfLife >> DecayModes;
+ if(Z >= fZAIThreshold )
+ {
+ // Get Isomeric State;
+
+ if(Iname=="gs")
+ I = 0;
+ else
+ if(Iname[0]=='m')
+ {
+ if( atoi( Iname.substr(1).c_str() )==0 )
+ I = 1;
+ else
+ I = atoi( Iname.substr(1).c_str() );
+
+ }
+
+
+ int start=0;
+ double branch_test=0;
+ double branch_test_f=0;
+
+ ZAI ParentZAI = ZAI(Z,A,I);
+ IsotopicVector DaughtersMap;
+ bool stable = true;
+
+ while(start<int(DecayModes.size()))
+ {
+ ZAI DaughterZAI;
+ double BR;
+ int daughter_A=0;
+ int daughter_N=0;
+ int daughter_Z=0;
+ int Iso=0;
+ int DM=-1;
+ // FPDistribution *FP=0;
+ string decay_name = GetDecay(DecayModes, BR, Iso, start);
+
+ if (decay_name == "s") {DM=0;daughter_N=(A-Z); daughter_Z=Z;BR=1;}
+ if (decay_name == "b-") {DM=1;stable=false; daughter_N=(A-Z)-1; daughter_Z=Z+1;}
+ if (decay_name == "n") {DM=2;stable=false; daughter_N=(A-Z)-1; daughter_Z=Z;}
+ if (decay_name == "nn") {DM=3;stable=false; daughter_N=(A-Z)-2; daughter_Z=Z;}
+ if (decay_name == "b-n"){DM=4;stable=false; daughter_N=(A-Z)-2; daughter_Z=Z+1;}
+ if (decay_name == "p") {DM=5;stable=false; daughter_N=(A-Z); daughter_Z=Z-1;}
+ if (decay_name == "b-a"){DM=6;stable=false; daughter_N=(A-Z)-3; daughter_Z=Z-1;}
+ if (decay_name == "pp") {DM=7;stable=false; daughter_N=(A-Z); daughter_Z=Z-2;}
+ if (decay_name == "ce") {DM=8;stable=false; daughter_N=(A-Z)+1; daughter_Z=Z-1;}
+ if (decay_name == "a") {DM=9;stable=false; daughter_N=(A-Z)-2; daughter_Z=Z-2;}
+ if (decay_name == "cen"){DM=10;stable=false; daughter_N=(A-Z); daughter_Z=Z-1;}
+ if (decay_name == "cep"){DM=11;stable=false; daughter_N=(A-Z)+1; daughter_Z=Z-2;}
+ if (decay_name == "it") {DM=12;stable=false; daughter_N=(A-Z); daughter_Z=Z;Iso = I-1;}
+ if (decay_name == "db-"){DM=13;stable=false; daughter_N=(A-Z)-2; daughter_Z=Z+2;}
+ if (decay_name == "db+"){DM=14;stable=false; daughter_N=(A-Z)+2; daughter_Z=Z-2;}
+ if (decay_name == "ita"){DM=15;stable=false; daughter_N=(A-Z)-2; daughter_Z=Z-2;}
+ if (decay_name == "sf") {DM=16;stable=false; daughter_N=0; daughter_Z=-2; Iso = -2;}
+ if (decay_name == "cesf"){DM=17;stable=false; daughter_N=0; daughter_Z=-2; Iso = -2;}
+ if (decay_name == "b-sf"){DM=18;stable=false; daughter_N=0; daughter_Z=-2; Iso = -2;}
+
+ daughter_A = daughter_Z + daughter_N;
+ {
+ if( daughter_Z < fZAIThreshold && daughter_Z!=-2 )
+ daughter_A = daughter_Z = Iso = -3;
+ // not spontaneous fission
+ ZAI DaughterZAI = ZAI(daughter_Z,daughter_A,Iso);
+
+ if((BR>1e-10) && (!stable))
+ {
+ if(DM <= 15)
+ {
+ DaughtersMap += BR * DaughterZAI;
+ branch_test+=BR;
+ }
+ else if( DM <= 18)
+ {
+ if(fSpontaneusYield.size() == 0 || fReactionYield.size() == 0 || DM == 17 || DM == 18)
+ {
+ DaughtersMap += 2*BR * ZAI(-2,-2,-2);
+
+ branch_test_f += BR;
+ }
+ else
+ {
+
+ map<ZAI, IsotopicVector>::iterator it_yield = fSpontaneusYield.find(ParentZAI);
+ if(it_yield != fSpontaneusYield.end())
+ {
+ DaughtersMap += (BR* (*it_yield).second );
+ branch_test_f += BR* (*it_yield).second.GetSumOfAll() / 2.;
+
+ }
+ else
+ {
+ DaughtersMap += 2*BR * ZAI(-2,-2,-2);
+ branch_test_f += BR;
+ }
+ }
+ }
+
+ }
+
+ }
+ if (DM !=0)
+ stable = false;
+ // End of While loop
+ }
+
+ double btest = fabs(branch_test + branch_test_f-1.0);
+ if ( btest > 1e-8 && !stable )
+ if(branch_test+branch_test_f > 0)
+ DaughtersMap = DaughtersMap /(branch_test+branch_test_f);
+
+
+
+ if (HalfLife < fShorstestHalflife && !stable)
+ fFastDecay.insert( pair< ZAI, map<ZAI, double> > ( ParentZAI, DaughtersMap.GetIsotopicQuantity() ) );
+ else if (stable)
+ {
+ IsotopicVector StableIV = ParentZAI *1;
+ ZAIDecay.insert( pair< ZAI, pair<double, map< ZAI, double > > >
+ (ParentZAI, pair<double, map< ZAI, double > >
+ ( 1e36, StableIV.GetIsotopicQuantity()) ) );
+ }
+ else
+ ZAIDecay.insert( pair< ZAI, pair<double, map< ZAI, double > > >
+ (ParentZAI, pair<double, map< ZAI, double > >
+ ( HalfLife, DaughtersMap.GetIsotopicQuantity()) ) );
+
+
+ }
+
+ } while (!infile.eof());
+
+ {
+ int i = 0;
+ map<ZAI, pair<double, map< ZAI, double > > >::iterator it;
+ for(it = ZAIDecay.begin() ; it != ZAIDecay.end(); it++)
+ {
+ findex.insert( pair<int, ZAI > ( i, (*it).first ) );
+ findex_inver.insert( pair<ZAI, int > ( (*it).first , i ));
+ i++;
+ }
+ }
+
+ // Fill the Decay Part of the Bateman Matrix Always the same !
+ bool FastDecayValidation = false;
+ while (!FastDecayValidation)
+ {
+ map<ZAI, map<ZAI, double> > FastDecayCopy = fFastDecay;
+
+ map<ZAI, map<ZAI, double> >::iterator FD_it;
+ map<ZAI, map<ZAI, double> >::iterator FD_it_Origin;
+
+
+ for(FD_it = FastDecayCopy.begin(); FD_it != FastDecayCopy.end(); FD_it++)
+ {
+
+ FD_it_Origin = fFastDecay.find(FD_it->first);
+
+ map<ZAI, double> BR = (*FD_it).second;
+ map<ZAI, double>::iterator BR_it;
+
+ for(BR_it = BR.begin(); BR_it != BR.end(); BR_it++)
+ {
+
+ map<ZAI, int>::iterator it_index = findex_inver.find( (*BR_it).first );
+
+ if( it_index == findex_inver.end() )
+ {
+ map<ZAI, map<ZAI, double> >::iterator FD2_it = FastDecayCopy.find((*BR_it).first);
+ if( FD2_it != FastDecayCopy.end() )
+ {
+ map<ZAI, double>::iterator BR2_it;
+ (*FD_it_Origin).second.erase((*BR_it).first);
+
+ for (BR2_it = (*FD2_it).second.begin(); BR2_it != (*FD2_it).second.end(); BR2_it++)
+ {
+
+ pair<map<ZAI, double>::iterator, bool> IResult;
+ IResult = (*FD_it_Origin).second.insert( pair<ZAI, double> ( (*BR2_it).first, (*BR_it).second * (*BR2_it).second ) );
+
+ if( !IResult.second)
+ (*IResult.first).second += (*BR_it).second * (*BR2_it).second ;
+
+
+ }
+
+ }
+ else
+ {
+ (*FD_it_Origin).second.erase( (*BR_it).first );
+ pair< map<ZAI, double>::iterator, bool> IResult;
+ IResult = (*FD_it_Origin).second.insert( pair<ZAI, double> ( ZAI(-3,-3,-3), (*BR_it).second) );
+ if( !IResult.second)
+ (*IResult.first).second += (*BR_it).second ;
+ }
+
+
+
+ }
+ }
+
+ }
+
+
+ FastDecayValidation = true;
+ for(FD_it = fFastDecay.begin(); FD_it != fFastDecay.end(); FD_it++)
+ {
+ map<ZAI, double>::iterator BR_it;
+ for (BR_it = (*FD_it).second.begin(); BR_it != (*FD_it).second.end(); BR_it++)
+ {
+ map<ZAI, int>::iterator Index_it = findex_inver.find( (*BR_it).first );
+ map<ZAI, map<ZAI, double> >::iterator FD2_it = fFastDecay.find( (*BR_it).first );
+ if(Index_it == findex_inver.end() && FD2_it == fFastDecay.end())
+ FastDecayValidation = false;
+ }
+ }
+ }
+
+
+ fDecayMatrix.ResizeTo(findex.size(),findex.size());
+ for(int i = 0; i < (int)findex.size(); i++)
+ for(int j = 0; j < (int)findex.size(); j++)
+ fDecayMatrix[i][j] = 0;
+
+
+
+
+ {
+ int i = 0;
+ map<ZAI, pair<double, map< ZAI, double > > >::iterator it;
+ for(it = ZAIDecay.begin() ; it != ZAIDecay.end(); it++)
+ {
+ map< ZAI, double >::iterator it2;
+ map< ZAI, double > decaylist = (*it).second.second;
+ for(it2 = decaylist.begin(); it2!= decaylist.end(); it2++)
+ {
+
+ map<ZAI, int >::iterator it3 = findex_inver.find( (*it2).first );
+ if( it3 != findex_inver.end() )
+ {
+ fDecayMatrix[(*it3).second][i] += log(2.)/(*it).second.first * (*it2).second;
+ }
+ else
+ {
+ map<ZAI, map<ZAI, double> >::iterator it4 = fFastDecay.find( (*it2).first );
+
+ if( it4 == fFastDecay.end() )
+ {
+
+
+ fDecayMatrix[0][i] += log(2.)/(*it).second.first * (*it2).second;
+ }
+ else
+ {
+ map< ZAI, double >::iterator it5;
+ map< ZAI, double > decaylist2 = (*it4).second;
+ for(it5 = decaylist2.begin(); it5!= decaylist2.end(); it5++)
+ {
+ it3 = findex_inver.find( (*it5).first );
+ if( it3 == findex_inver.end() )
+ fDecayMatrix[0][i] += log(2.)/(*it).second.first * (*it2).second * (*it5).second;
+
+ else
+ {
+ fDecayMatrix[(*it3).second][i] += log(2.)/(*it).second.first * (*it2).second * (*it5).second;
+
+ }
+ }
+ }
+
+ }
+ }
+ fDecayMatrix[i][i] += -log(2.)/(*it).second.first;
+
+ i++;
+
+
+ }
+ }
+
+}
+
+//________________________________________________________________________
+/* Fission Stuff */
+//________________________________________________________________________
+void IrradiationModel::SetFissionEnergy(ZAI zai, double E)
+{
+ pair<map<ZAI, double>::iterator, bool> IResult;
+ IResult = fFissionEnergy.insert( pair<ZAI ,double>(zai, E));
+ if(!IResult.second)
+ IResult.first->second = E;
+
+}
+
+//________________________________________________________________________
+void IrradiationModel::SetFissionEnergy(string FissionEnergyFile)
+{
+ ifstream FissionFile(FissionEnergyFile.c_str()); // Open the File
+ if(!FissionFile) //check if file is correctly open
+ {
+ cout << "!!Warning!! !!!IrradiationModel!!! \n Can't open \"" << FissionFile << "\"\n" << endl;
+ GetLog()->fLog << "!!Warning!! !!!IrradiationModel!!! \n Can't open \"" << FissionFile << "\"\n" << endl;
+ }
+
+ do {
+ int Z = 0;
+ int A = 0;
+ int I = 0;
+ double E = 0;
+ FissionFile >> Z >> A >> I >> E;
+ SetFissionEnergy(Z, A, I, E);
+ } while (!FissionFile.eof());
+}
+
+//________________________________________________________________________
+map< ZAI,IsotopicVector > IrradiationModel::ReadFPYield(string Yield)
+{
+ IsotopicVector EmptyIV;
+ map< ZAI,IsotopicVector > Yield_map;
+
+ ifstream infile(Yield.c_str());
+ if(!infile)
+ {
+ cout << "!!Warning!! !!!IrradiationModel!!! \n Can't open \"" << Yield << "\"\n" << endl;
+ GetLog()->fLog << "!!Warning!! !!!IrradiationModel!!! \n Can't open \"" << Yield<< "\"\n" << endl;
+ }
+
+
+ string line;
+ int start = 0;
+
+ getline(infile, line);
+
+ do
+ {
+ int Z = atof(StringLine::NextWord(line, start, ' ').c_str());
+ int A = atof(StringLine::NextWord(line, start, ' ').c_str());
+ int I = 0;
+
+ // if(Z!=0 && A!=0)
+ {
+ pair<map<ZAI, IsotopicVector>::iterator, bool> IResult;
+ IResult = Yield_map.insert(pair<ZAI,IsotopicVector>(ZAI(Z,A,I),EmptyIV) );
+ if(!IResult.second)
+ {
+ cout << "!!Error!! !!!IrradiationModel!!! Many accurance of ZAI " << Z << " " << A;
+ cout << " in " << Yield << " file!! Please Check it !!!" << endl;
+ exit(1);
+
+ }
+ }
+ }while(start < (int)line.size()-1);
+
+ do
+ {
+ start = 0;
+
+ getline(infile, line);
+ int Z = atof(StringLine::NextWord(line, start, ' ').c_str());
+ int A = atof(StringLine::NextWord(line, start, ' ').c_str());
+ int I = atof(StringLine::NextWord(line, start, ' ').c_str());
+ map<ZAI, IsotopicVector>::iterator it = Yield_map.begin();
+ do
+ {
+ if (it == Yield_map.end())
+ {
+ cout << "!!Error!! !!!IrradiationModel!!! Many accurance of the PF " << Z << " " << A;
+ cout << " in " << Yield << " file!! Please Check it";
+ cout << "(Number of yield does not match the number of ZAI that fission !!!" << endl;
+ exit(1);
+
+ }
+
+ double Yield_values = atof(StringLine::NextWord(line, start, ' ').c_str());
+ (*it).second += Yield_values * ZAI(Z,A,I);
+
+ it++;
+ }while(start < (int)line.size()-1);
+
+
+
+
+ } while (!infile.eof());
+ return Yield_map;
+}
+
+//________________________________________________________________________
+void IrradiationModel::LoadFPYield(string SponfaneusYield, string ReactionYield)
+{
+
+ fSpontaneusYield = ReadFPYield(SponfaneusYield);
+ fReactionYield = ReadFPYield(ReactionYield);
+ fZAIThreshold = 0;
+}
diff --git a/source/branches/CLASSV3/src/Makefile b/source/branches/CLASSV3/src/Makefile
index b13a3c235..651df1500 100755
--- a/source/branches/CLASSV3/src/Makefile
+++ b/source/branches/CLASSV3/src/Makefile
@@ -29,7 +29,9 @@ OBJS = CLASS.o \
FuelDataBank.o DecayDataBank.o \
DynamicalSystem.o\
EvolutionData.o EvolutionDataDict.o \
- LogFile.o
+ IrradiationModel.o\
+ EquivalenceModel.o\
+ LogFile.o
ROOTOBJS = CLASSObject.o CLASSObjectDict.o\
CLASSFacility.o CLASSFacilityDict.o\
--
GitLab