diff --git a/source/trunk/include/CLSSFacility.hxx b/source/trunk/include/CLSSFacility.hxx
index be28389958ce62f9cb436246d301afb6c9b62f52..77780ee81e41b50be9fa773b7a00ea9934151a90 100644
--- a/source/trunk/include/CLSSFacility.hxx
+++ b/source/trunk/include/CLSSFacility.hxx
@@ -33,7 +33,7 @@ public :
//********* Get Method *********//
int GetId() const { return fId; } //!< Return the Facility Parc'Is
- IsotopicVector GetInsideIV() const { return fInsideIV; } //!< Return the IV contain in the Facility
+ virtual IsotopicVector GetInsideIV() const { return fInsideIV; } //!< Return the IV contain in the Facility
cSecond GetInternalTime() const { return fInternalTime; } //!< Return Creation Time
diff --git a/source/trunk/include/DataBank.hxx b/source/trunk/include/DataBank.hxx
index 10e21a865a396d4fb7ec58a0a732378e9130bfa7..300f78195abcdd8b851642221ad184d6135ce807 100755
--- a/source/trunk/include/DataBank.hxx
+++ b/source/trunk/include/DataBank.hxx
@@ -13,6 +13,7 @@
#include "CLSSObject.hxx"
#include "TMatrix.h"
#include "IsotopicVector.hxx"
+#include "DynamicalSystem.hxx"
#include <map>
#include <vector>
@@ -30,7 +31,7 @@ double ReactionRateWeightedDistance(EvolutionData DB, IsotopicVector IV1 );
template <class T>
-class DataBank : public CLSSObject
+class DataBank : public CLSSObject, DynamicalSystem
{
public :
@@ -54,20 +55,33 @@ public :
map<double, EvolutionData> GetDistancesTo(IsotopicVector isotopicvector, double t = 0) const; //! Return a map containing the distance of each EvolutionData in the DataBase to the set IV at the t time
EvolutionData GetClosest(IsotopicVector isotopicvector, double t = 0) const; //! Return the closest EvolutionData from the DataBank.
+
+ string GetDataFileName() const { return fDataFileName; }
+ string GetDataDirectoryName() const { return fDataDirectoryName; }
+
+ double GetShorstestHalflife() const { return fShorstestHalflife; }
+
+
+
+
//********* Set Method *********//
- void SetDataBank(map<T ,EvolutionData > mymap) { fDataBank = mymap; }
+ void SetDataBank(map<T ,EvolutionData > mymap) { fDataBank = mymap; }
void SetDataBaseIndex(string database) { fDataBaseIndex = database; }
EvolutionData GenerateEvolutionData(IsotopicVector isotopicvector, double cycletime, double Power); //!< Genration of a New EvolutionData From the one already present
EvolutionData OldGenerateEvolutionData(IsotopicVector isotopicvector, double cycletime, double Power); //!< Genration of a New EvolutionData From the one already present
- void SetUpdateReferenceDBatEachStep(bool val) {fUpdateReferenceDBatEachStep = val;}
void SetOldReadMethod(bool val) { fOldReadMethod = val;}
void SetFissionEnergy(string FissionEnergyFile);
void SetFissionEnergy(ZAI zai, double E);
void SetFissionEnergy(int Z, int A, int I, double E ) { SetFissionEnergy(ZAI(Z,A,I), E);}
+ void SetDataFileName(string name) { fDataFileName = name;}
+ void SetDataDirectoryName(string name) { fDataDirectoryName = name;}
+ void SetShartestHalfLife(double halflife) { fShorstestHalflife = halflife; BuildDecayMatrix(); ReadDataBase();}
+
+
//********* Modification Method *********//
IsotopicVector Evolution(const T &key, double dt); ///< Return the Product IsotopicVector evolution from zai during a dt time
@@ -79,24 +93,51 @@ public :
///< 1 for each ZAI weighted with its XS,
///< 2 for each ZAI weighted with coefficient given by the user.
- void UseDBTimeStep(bool oldmethod = true) {fUseDBTimeStep = oldmethod;}
-
+ void BuildDecayMatrix();
+
void UseOldGeneration(bool oldmethod = true) {fUseOldGeneration = oldmethod;}
+ 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)
+
+
+
//********* Printing Method *********//
void Print() const;
protected :
+ double fShorstestHalflife;
+
+ string fDataFileName;
+ string fDataDirectoryName;
+
map<T, EvolutionData> fDataBank;
map<T, EvolutionData> fDataBankCalculated;
string fDataBaseIndex;
- bool fUpdateReferenceDBatEachStep;
+ bool fUseRK4EvolutionMethod;
bool fOldReadMethod;
bool fUseOldGeneration;
- bool fUseDBTimeStep;
-
+
string fFuelType;
pair<double,double> fBurnUpRange;
vector<double> fFuelParameter;
@@ -106,44 +147,33 @@ protected :
IsotopicVector fDistanceParameter; ///< weight for each ZAI in the distance calculation
- TMatrixT<double> fDecayMatrix; ///< Matrix with half life of each nuclei
- void BuildDecayMatrix();
+
TMatrixT<double> GetFissionXsMatrix(EvolutionData EvolutionDataStep,double TStep);
TMatrixT<double> GetCaptureXsMatrix(EvolutionData EvolutionDataStep,double TStep);
TMatrixT<double> Getn2nXsMatrix(EvolutionData EvolutionDataStep,double TStep);
-
-
TMatrixT<double> ExtractXS(EvolutionData EvolutionDataStep,double TStep);
+ void OldBuildDecayMatrix();
+
+ string GetDecay(string DecayModes, double &BR,int &Iso, int &StartPos);
+ 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;
+
+ double *fTheNucleiVector; //!< The evolving atoms copied from Material proportions.
+ double **fTheMatrix; //!< The evolution Matrix
+
+ int fNVar; //!< The size of the composition vector and /or number of ZAIs involved.
+ 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.
+
map<ZAI, int> findex_inver;
map<int, ZAI> findex;
- //0 TMP
- //1 PF
- //2 232Th
- //3 233U
- //4 234U
- //5 235U
- //6 236U
- //7 238U
- //8 237Np
- //9 238Pu
- //10 239Pu
- //11 240Pu
- //12 241Pu
- //13 242Pu
- //14 241Am
- //15 242Am*
- //16 243Am
- //17 242Cm
- //18 243Cm
- //19 244Cm
- //20 245Cm
- //21 246Cm
- //22 247Cm
- //23 248Cm
};
diff --git a/source/trunk/include/DynamicalSystem.hxx b/source/trunk/include/DynamicalSystem.hxx
new file mode 100755
index 0000000000000000000000000000000000000000..66880571e476dabbda3442999e8f7c3a18d307e1
--- /dev/null
+++ b/source/trunk/include/DynamicalSystem.hxx
@@ -0,0 +1,127 @@
+#ifndef _DynamicalSystem_
+#define _DynamicalSystem_
+/*!
+ \file
+ \brief Header file for DynamicalSystem class.
+*/
+
+#include <math.h>
+#include <vector>
+
+
+using namespace std;
+
+//-----------------------------------------------------------------------------//
+
+
+
+//! DynamicalSystem class solves system of differential equations.
+/*!
+// A DynamicalSystem is a base class thatsolves system of 1st order of differential equations.
+// \f[ \frac{d\vec{Y}}{dt}=A\vec{Y}\f]
+// The differential equations are built in DynamicalSystem::BuildEqns ; the method MUST be
+// defined in the derived classes.
+// In this first version only Runge-Kutta method is implemented, but the aim of this class
+// is to provide also other methods such as CRAM (Chebyshev rational approximation method).
+//
+//
+// @author PTO.
+// @version 1.0
+*/
+//________________________________________________________________________
+
+class DynamicalSystem
+{
+ public :
+
+ DynamicalSystem(); //!< Normal Constructor
+ DynamicalSystem(const DynamicalSystem & DS); //!< Copy Constructor
+ virtual ~DynamicalSystem(); //!< Destructor
+
+ /*!
+ \name Mains attributes of the DynamicalSystem
+ */
+ //@{
+ int GetNumberOfEquationSize(){return fNVar;} //!< return the number of equations.
+ void SetNumberOfEquationSize(int n){fNVar=n;} //!< set the number of equations.
+ //@}
+
+ /*!
+ \name Runge-Kutta related methods
+ Algorithms are taken from Numerical Receipes.
+ */
+ //@{
+
+ void SetPrecision(double eps=1e-5){fPrecision=eps;} //!< set RK precision to change the integration step
+
+ //! Forbid negative value during integration.
+ /*!
+ For some quantities (such as nuclei composition), negative values are forbidden.
+ But, due to integration step and very fast variation of the integrated variables
+ Runge-Kutta wil produce very small negative value. This method is used to force
+ negative value to be zero.
+ */
+ void SetForbidNegativeValue(){fIsNegativeValueAllowed=false;}
+
+ //! Runge Kutta calling method.
+ /*!
+ // \param YStart: input : the initial condition Y(t1) ; output the final value Y(t2)
+ // \param t1: initial time
+ // \param t2: final time
+ */
+ void RungeKutta(double *YStart, double t1, double t2, int EquationNumber);
+
+ //! Builds the equations for integration.
+ /*!
+ This method is an abstract method ; it MUST be overwritten by derived classes.
+ // \param t: time at which the equations are built
+ // \param Y: array of variable at time t
+ // \param dYdt: ode's variable.
+ */
+ virtual void BuildEqns(double t, double *Y, double *dYdt){}
+
+ //@}
+
+ /*!
+ \name Miscellaneous methods
+ */
+ //@{
+
+ //@}
+
+ protected :
+ //! Runge Kutta main method.
+ /*!
+ // Call by RungeKutta
+ // \param y: initial values to integrate
+ // \param dydx: ode's equations (variable is x)
+ // \param x: variable of integration
+ // \param h: step size for integration
+ // \param yout: result after integration
+ */
+ void RK4(double *y, double *dydx, double x, double h, double *yout);
+ //! Adaptative Step Size method for RK.
+ /*!
+ // Call by RK4
+ // \param y: initial values to integrate
+ // \param dydx: ode's equations (variable is x)
+ // \param x: new value of the variable after the adaptative step
+ // \param htry: try step size for integration
+ // \param eps: precision
+ // \param yscal: result after hdid step integration
+ // \param hdid: did step size for integration
+ // \param hnext: next step size for integration
+ */
+ void AdaptStepSize(double *y, double *dydx, double *x, double htry, double eps, double *yscal, double *hdid, double *hnext);
+
+ int fNVar; //!< The size of the composition vector and /or number of ZAIs involved.
+ 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/trunk/include/EvolutionData.hxx b/source/trunk/include/EvolutionData.hxx
index 241247eae59cf59edd569b67421cd66f1c1cf3c8..2fb3861c3c584b28332105985fadd725de303d2e 100755
--- a/source/trunk/include/EvolutionData.hxx
+++ b/source/trunk/include/EvolutionData.hxx
@@ -122,7 +122,8 @@ protected :
string fFuelType;
double fPower;
double fCycleTime;
-
+ double fNormFactor;
+
void OldReadDB(string DBfile);
void ReadDB(string DBfile, bool oldread = false);
diff --git a/source/trunk/include/IsotopicVector.hxx b/source/trunk/include/IsotopicVector.hxx
index 05b907b225383e00bda69034da53c323f2f2b144..d9513743ed521584a630b99875ef1abc0c50cb77 100755
--- a/source/trunk/include/IsotopicVector.hxx
+++ b/source/trunk/include/IsotopicVector.hxx
@@ -45,7 +45,9 @@ public :
vector<int> GetChemicalSpecies() const; //!< Return the Species Species contained
int GetZAIQuantity() {return fIsotopicQuantity.size(); }
-
+
+ double GetSumOfAll();
+
//********* Modification Method *********//
void Clear(); //!< Empty all the IV
void ClearNeed(); //!< Empty Need componant of the IV
diff --git a/source/trunk/include/ZAI.hxx b/source/trunk/include/ZAI.hxx
index 52a83356a777b0e201f30a97e70421d87770ea18..903bb134e51c9bb1f1592644c3185ce7ba2922d7 100755
--- a/source/trunk/include/ZAI.hxx
+++ b/source/trunk/include/ZAI.hxx
@@ -60,7 +60,8 @@ class ZAI : public TObject
(fZ < zai.Z()) : ( (fA != zai.A())?
(fA < zai.A()) : (fI < zai.I()) ); }
//!< ZAI Comparator
- bool operator !=(const ZAI& zai) const { return ( fZ != zai.Z() ) || ( fA != zai.A() ) || ( fI < zai.I() ); }
+ bool operator !=(const ZAI& zai) const { return ( fZ != zai.Z() ) || ( fA != zai.A() ) || ( fI != zai.I() ); }
+ bool operator ==(const ZAI& zai) const { return ( fZ == zai.Z() && fA == zai.A() && fI == zai.I()); }
//!< ZAI Comparator
void Print() const { cout << fZ << " " << fA << " " << fI << endl;}
protected :
diff --git a/source/trunk/src/DataBank.cxx b/source/trunk/src/DataBank.cxx
index 914bb25aa47ea2016b4d6fb063aefa68e68879f5..aa4e4a821209df940d79fc18876988cd15d9d14d 100755
--- a/source/trunk/src/DataBank.cxx
+++ b/source/trunk/src/DataBank.cxx
@@ -1,6 +1,7 @@
#include "DataBank.hxx"
#include "IsotopicVector.hxx"
+#include "CLASSHeaders.hxx"
#include "EvolutionData.hxx"
#include "LogFile.hxx"
#include "StringLine.hxx"
@@ -182,6 +183,9 @@ bool DataBank<ZAI>::IsDefine(const ZAI& zai) const
//________________________________________________________________________
//________________________________________________________________________
+template<>
+void DataBank<ZAI>:: BuildEqns(double t, double *N, double *dNdt)
+{}
@@ -194,48 +198,103 @@ bool DataBank<ZAI>::IsDefine(const ZAI& zai) const
//________________________________________________________________________
template<>
void DataBank<IsotopicVector>::ReadDataBase();
+
template<>
-TMatrixT<double> DataBank<IsotopicVector>::GetFissionXsMatrix(EvolutionData EvolutionDataStep,double TStep);
+TMatrixT<double> DataBank<IsotopicVector>::GetFissionXsMatrix(EvolutionData EvolutionDataStep,double BUStep);
+
template<>
-TMatrixT<double> DataBank<IsotopicVector>::GetCaptureXsMatrix(EvolutionData EvolutionDataStep,double TStep);
+TMatrixT<double> DataBank<IsotopicVector>::GetCaptureXsMatrix(EvolutionData EvolutionDataStep,double BUStep);
+
template<>
-TMatrixT<double> DataBank<IsotopicVector>::Getn2nXsMatrix(EvolutionData EvolutionDataStep,double TStep);
+TMatrixT<double> DataBank<IsotopicVector>::Getn2nXsMatrix(EvolutionData EvolutionDataStep,double BUStep);
+
template<>
-TMatrixT<double> DataBank<IsotopicVector>::ExtractXS(EvolutionData EvolutionDataStep,double TStep);
+TMatrixT<double> DataBank<IsotopicVector>::ExtractXS(EvolutionData EvolutionDataStep,double BUStep);
+
template<>
void DataBank<IsotopicVector>::BuildDecayMatrix();
+
template<>
EvolutionData DataBank<IsotopicVector>::OldGenerateEvolutionData(IsotopicVector isotopicvector, double cycletime, double Power);
+template<>
+string DataBank<IsotopicVector>::GetDecay(string DecayModes, double &BR,int &Iso, int &StartPos);
+
+template<>
+void DataBank<IsotopicVector>::SetTheMatrixToZero();
+
+template<>
+void DataBank<IsotopicVector>::ResetTheMatrix();
+
+template<>
+void DataBank<IsotopicVector>::SetTheNucleiVectorToZero();
+
+template<>
+void DataBank<IsotopicVector>::ResetTheNucleiVector();
+
+template<>
+void DataBank<IsotopicVector>::BuildEqns(double t, double *N, double *dNdt);
+
+template<>
+void DataBank<IsotopicVector>::SetTheMatrix(TMatrixT<double> BatemanMatrix);
+
+template<>
+void DataBank<IsotopicVector>::SetTheNucleiVector(TMatrixT<double> NEvolutionMatrix);
+
+template<>
+TMatrixT<double> DataBank<IsotopicVector>::GetTheNucleiVector();
+
//________________________________________________________________________
template<>
-DataBank<IsotopicVector>::DataBank()
+DataBank<IsotopicVector>::DataBank():DynamicalSystem()
{
- fUpdateReferenceDBatEachStep = false;
+ fTheNucleiVector = 0;
+ fTheMatrix = 0;
+
fOldReadMethod = true;
fUseOldGeneration = false;
- fUseDBTimeStep = false;
fDistanceType = 0;
+ fShorstestHalflife = 3600.;
+
+ fDataDirectoryName = getenv("CLASS_PATH");
+ fDataDirectoryName += "/source/data/";
+ fDataFileName = "chart.JEF3T";
+
BuildDecayMatrix();
+
+ fNVar = findex_inver.size();
+ SetForbidNegativeValue();
+
}
template<>
-DataBank<IsotopicVector>::DataBank(LogFile* Log, string DB_index_file, bool setlog, bool olfreadmethod)
+DataBank<IsotopicVector>::DataBank(LogFile* Log, string DB_index_file, bool setlog, bool olfreadmethod):DynamicalSystem()
{
SetLog(Log);
IsLog(setlog);
+ fTheNucleiVector = 0;
+ fTheMatrix = 0;
+
+ fDataDirectoryName = getenv("CLASS_PATH");
+ fDataDirectoryName += "/source/data/";
+ fDataFileName = "chart.JEF3T";
+
fDataBaseIndex = DB_index_file;
- fUpdateReferenceDBatEachStep = false;
fOldReadMethod = olfreadmethod;
fUseOldGeneration = false;
- fUseDBTimeStep = false;
fDistanceType = 0;
+
+ fShorstestHalflife = 3600.;
+
BuildDecayMatrix();
ReadDataBase();
+ fNVar = findex_inver.size();
+ SetForbidNegativeValue();
+
if(PrintLog())
{
// Warning
@@ -362,7 +421,7 @@ map<double, EvolutionData> DataBank<IsotopicVector>::GetDistancesTo(IsotopicVect
pair<map<double, EvolutionData>::iterator, bool> IResult;
double D = Distance(isotopicvector.GetActinidesComposition(), (*it).second.GetIsotopicVectorAt(t).GetActinidesComposition()/ Norme( (*it).second.GetIsotopicVectorAt(t).GetActinidesComposition() )*Norme(isotopicvector.GetActinidesComposition())
- ,fDistanceType, fDistanceParameter);
+ ,fDistanceType, fDistanceParameter);
IResult = distances.insert( pair<double, EvolutionData>( D , (*it).second ) );
}
@@ -378,10 +437,10 @@ EvolutionData DataBank<IsotopicVector>::GetClosest(IsotopicVector isotopicvector
map<IsotopicVector, EvolutionData > evolutiondb = fDataBank;
double distance = Distance(isotopicvector.GetActinidesComposition(),
- evolutiondb.begin()->second.GetIsotopicVectorAt(t).GetActinidesComposition()
- / Norme( evolutiondb.begin()->second.GetIsotopicVectorAt(t).GetActinidesComposition() )
- * Norme(isotopicvector.GetActinidesComposition()),
- fDistanceType, fDistanceParameter);
+ evolutiondb.begin()->second.GetIsotopicVectorAt(t).GetActinidesComposition()
+ / evolutiondb.begin()->second.GetIsotopicVectorAt(t).GetActinidesComposition().GetSumOfAll()
+ * isotopicvector.GetActinidesComposition().GetSumOfAll(),
+ fDistanceType, fDistanceParameter);
EvolutionData CloseEvolData = evolutiondb.begin()->second ;
@@ -390,17 +449,16 @@ EvolutionData DataBank<IsotopicVector>::GetClosest(IsotopicVector isotopicvector
{
pair<map<double, EvolutionData>::iterator, bool> IResult;
double D = Distance(isotopicvector.GetActinidesComposition(),
- (*it).second.GetIsotopicVectorAt(t).GetActinidesComposition()
- / Norme( (*it).second.GetIsotopicVectorAt(t).GetActinidesComposition() )
- * Norme( isotopicvector.GetActinidesComposition()),
- fDistanceType, fDistanceParameter);
+ (*it).second.GetIsotopicVectorAt(t).GetActinidesComposition()
+ / (*it).second.GetIsotopicVectorAt(t).GetActinidesComposition().GetSumOfAll()
+ * isotopicvector.GetActinidesComposition().GetSumOfAll(),
+ fDistanceType, fDistanceParameter);
if (D< distance)
{
distance = D;
CloseEvolData = (*it).second;
}
}
-
return CloseEvolData;
}
@@ -516,21 +574,11 @@ void DataBank<IsotopicVector>::SetDistanceType(int DistanceType)
-
template<>
EvolutionData DataBank<IsotopicVector>::GenerateEvolutionData(IsotopicVector isotopicvector, double cycletime, double Power)
{
- /*{
- map<IsotopicVector, EvolutionData>::iterator it;
- for( it = fDataBankCalculated.begin(); it != fDataBankCalculated.end(); it++)
- {
- if (Distance(isotopicvector, (*it).first) == 0
- && cycletime == (*it).second.GetCycleTime()
- && Power == (*it).second.GetPower() ) {
- return (*it).second;
- }
- }
- }*/
+ SetTheMatrixToZero();
+ SetTheNucleiVectorToZero();
if(fUseOldGeneration)
{
@@ -541,12 +589,7 @@ EvolutionData DataBank<IsotopicVector>::GenerateEvolutionData(IsotopicVector iso
}
string ReactorType;
- //-------------------------//
- //--- Perform Evolution ---//
- //-------------------------//
- int NStep = 80.;
- double timevector[NStep+1];
- timevector[0] = 0.;
+
vector< TMatrixT<double> > NMatrix ;// TMatrixT<double>(decayindex.size(),1))
{ // Filling the t=0 State;
map<ZAI, double > isotopicquantity = isotopicvector.GetIsotopicQuantity();
@@ -567,19 +610,52 @@ EvolutionData DataBank<IsotopicVector>::GenerateEvolutionData(IsotopicVector iso
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 ;
-
-
+
+
}
NMatrix.push_back(N_0Matrix);
+
}
- TMatrixT<double> SigmaPhi = TMatrixT<double>(findex.size()*3+1,NStep); // Store the XS and the flux trought the evolution calculation.
+
+ //-------------------------//
+ //--- Perform Evolution ---//
+ //-------------------------//
EvolutionData EvolutionDataStep = GetClosest(isotopicvector.GetActinidesComposition(), 0.); //GetCLosest at the begining of evolution
-
ReactorType = EvolutionDataStep.GetReactorType();
+
+ double Na = 6.02214129e23; //N Avogadro
+ double M_ref = 0;
+ double M = 0;
+ double Power_ref = EvolutionDataStep.GetPower();
+ {
+ map<ZAI, double >::iterator it ;
+
+
+ IsotopicVector IVtmp = isotopicvector.GetActinidesComposition() + EvolutionDataStep.GetIsotopicVectorAt(0.).GetActinidesComposition();
+ map<ZAI, double >isotopicquantity = IVtmp.GetIsotopicQuantity();
+
+ for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
+ {
+ M_ref += EvolutionDataStep.GetIsotopicVectorAt(0.).GetActinidesComposition().GetZAIIsotopicQuantity( (*it).first )*cZAIMass.fZAIMass.find( (*it).first )->second/Na*1e-6;
+ M += isotopicvector.GetActinidesComposition().GetZAIIsotopicQuantity( (*it).first )*cZAIMass.fZAIMass.find( (*it).first )->second/Na*1e-6;
+ }
+ }
+
+ int DBTimeStepN = EvolutionDataStep.GetFissionXS().begin()->second->GetN();
+ double* DBTimeStep = EvolutionDataStep.GetFissionXS().begin()->second->GetX();
+
+ int InsideStep = 10;
+
+ int NStep = (DBTimeStepN);
+ double timevector[NStep];
+ timevector[0] = 0;
+ double Flux[NStep];
+
+ TMatrixT<double> SigmaPhi = TMatrixT<double>(findex.size()*3+1,NStep); // Store the XS and the flux trought the evolution calculation.
+
TMatrixT<double> FissionEnergy = TMatrixT<double>(findex.size(),1);
{
map< ZAI, int >::iterator it;
@@ -588,137 +664,154 @@ EvolutionData DataBank<IsotopicVector>::GenerateEvolutionData(IsotopicVector iso
map< ZAI, double >::iterator it2 = fFissionEnergy.find(it->first);
if(it2 == fFissionEnergy.end())
{
- if(it->first.A() > 90)
+ if(it->first.Z() > 90)
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;
-
-
+ FissionEnergy[it->second][0] = it2->second;
+
}
}
-
-
-
+
vector< TMatrixT<double> > FissionXSMatrix; //The Fisison XS Matrix
vector< TMatrixT<double> > CaptureXSMatrix; //The Capture XS Matrix
- vector< TMatrixT<double> > n2nXSMatrix; //The n2N XS Matrix
- double Flux[NStep];
- for(int i = 0; i < NStep; i++)
+ vector< TMatrixT<double> > n2nXSMatrix; //The n2N XS Matrix
+
+ for(int i = 0; i < NStep-1; i++)
{
-
- double TStep = cycletime/NStep*i;
-
- if(fUpdateReferenceDBatEachStep) //Updated At each Step
- {
- IsotopicVector IVStep;
- for(int k = 0; k < (int)findex.size(); k++)
- IVStep += findex.find(k)->second * NMatrix.back()[k][0];
-
- if(i != 0) // Update closest
- EvolutionDataStep = GetClosest(IVStep, TStep);
-
- //Extraxt XS (if the refecence is not updated will get it from the Evolution DataBase
- TMatrixT<double> SigmaPhiTmp = ExtractXS(EvolutionDataStep, TStep);
-
- for(int j=0; j < (int)findex.size()*3; j++) // Store it
- SigmaPhi[j][i] = SigmaPhiTmp[j][0];
- }
+ double TStepMax = ( (DBTimeStep[i+1]-DBTimeStep[i] ) ) * Power_ref/M_ref / Power*M ;
+
+ 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();
- double TStepForXS = 0;
-
- if( fUseDBTimeStep)
- {
- double DBStepTime = EvolutionDataStep.GetFinalTime()/(EvolutionDataStep.GetFissionXS().begin()->second->GetN() -1);
- while( (TStepForXS + DBStepTime) < TStep)
- TStepForXS += DBStepTime;
-
- }
- else
- TStepForXS =TStep ;
-
-
-
- FissionXSMatrix.push_back(GetFissionXsMatrix(EvolutionDataStep, TStepForXS)); //Feel the reaction Matrix
- CaptureXSMatrix.push_back(Getn2nXsMatrix(EvolutionDataStep, TStepForXS)); //Feel the reaction Matrix
- n2nXSMatrix.push_back(GetCaptureXsMatrix(EvolutionDataStep, TStepForXS)); //Feel the reaction Matrix
-
- double ESigmaN = 0;
-
- for (int j = 0; j < (int)findex.size() ; j++)
- ESigmaN -= FissionXSMatrix.back()[j][j]*NEvolutionMatrix[j][0]*1.60217653e-19*FissionEnergy[j][0];
-
- // Update Flux
- Flux[i] = Power/ESigmaN ;///(1+ 4*(2.e-10*TStepForXS/3600./24.*TStepForXS/3600./24. + 2.e-7*TStepForXS/3600./24. + 4.e-6));
-// cout << Flux[i] << endl;
+ FissionXSMatrix.push_back(GetFissionXsMatrix(EvolutionDataStep, DBTimeStep[i])); //Feel the reaction Matrix
+ CaptureXSMatrix.push_back(GetCaptureXsMatrix(EvolutionDataStep, DBTimeStep[i])); //Feel the reaction Matrix
+ n2nXSMatrix.push_back(Getn2nXsMatrix(EvolutionDataStep, DBTimeStep[i])); //Feel the reaction Matrix
+
// ---------------- Evolution
- TMatrixT<double> BatemanMatrix = TMatrixT<double>(findex.size(),findex.size());
-
- TMatrixT<double> BatemanReactionMatrix = TMatrixT<double>(findex.size(),findex.size());
+
BatemanReactionMatrix = FissionXSMatrix[i];
BatemanReactionMatrix += CaptureXSMatrix[i];
BatemanReactionMatrix += n2nXSMatrix[i];
-
- BatemanMatrix = BatemanReactionMatrix;
- BatemanMatrix *= Flux[i];
-
- BatemanMatrix += fDecayMatrix ;
-
- double TStepMax = cycletime/(NStep);
- timevector[i+1] = timevector[i] + TStepMax;
-
- BatemanMatrix *= TStepMax;
- TMatrixT<double> IdMatrix = TMatrixT<double>(findex.size(),findex.size());
- for(int j = 0; j < (int)findex.size(); j++)
- for(int k = 0; k < (int)findex.size(); k++)
+
+ if(fUseRK4EvolutionMethod)
+ {
+ for(int k=0; k < InsideStep; k++)
{
- if(k == j) IdMatrix[j][k] = 1;
- else IdMatrix[j][k] = 0;
+ 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();
+
}
-
-
- TMatrixT<double> BatemanMatrixDL = TMatrixT<double>(findex.size(),findex.size()); // Order 0 Term from the DL : Id
- TMatrixT<double> BatemanMatrixDLTermN = TMatrixT<double>(findex.size(),findex.size()); // Addind it;
-
+ NEvolutionMatrix = GetTheNucleiVector();
+ NMatrix.push_back(NEvolutionMatrix);
+ }
+ else
{
- BatemanMatrixDLTermN = IdMatrix;
- BatemanMatrixDL = BatemanMatrixDLTermN;
- int j = 1;
- double NormN;
-
- do
+
+ for(int k=0; k < InsideStep; k++)
{
- TMatrixT<double> BatemanMatrixDLTermtmp = TMatrixT<double>(findex.size(),findex.size()); // Adding it;
- BatemanMatrixDLTermtmp = BatemanMatrixDLTermN;
+ 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;
- BatemanMatrixDLTermN.Mult(BatemanMatrixDLTermtmp, BatemanMatrix );
-
- BatemanMatrixDLTermN *= 1./j;
- BatemanMatrixDL += BatemanMatrixDLTermN;
-
- NormN = 0;
- for(int m = 0; m < (int)findex.size(); m++)
- for(int n = 0; n < (int)findex.size(); n++)
- NormN += BatemanMatrixDLTermN[m][n]*BatemanMatrixDLTermN[m][n];
- j++;
-
- } while ( NormN != 0 );
+ if(k==0)
+ Flux[i]=Flux_k;
+
+ BatemanMatrix = BatemanReactionMatrix;
+ BatemanMatrix *= Flux_k;
+ BatemanMatrix += fDecayMatrix ;
+ BatemanMatrix *= TStepMax/InsideStep ;
+
+
+ TMatrixT<double> IdMatrix = TMatrixT<double>(findex.size(),findex.size());
+ for(int j = 0; j < (int)findex.size(); j++)
+ for(int k = 0; k < (int)findex.size(); k++)
+ {
+ if(k == j) IdMatrix[j][k] = 1;
+ else IdMatrix[j][k] = 0;
+ }
+
+
+ TMatrixT<double> BatemanMatrixDL = TMatrixT<double>(findex.size(),findex.size()); // Order 0 Term from the DL : Id
+ TMatrixT<double> BatemanMatrixDLTermN = TMatrixT<double>(findex.size(),findex.size()); // Addind it;
+
+ {
+ BatemanMatrix *= TStepMax ;
+ BatemanMatrixDLTermN = IdMatrix;
+ BatemanMatrixDL = BatemanMatrixDLTermN;
+ int j = 1;
+ double NormN;
+
+ do
+ {
+ TMatrixT<double> BatemanMatrixDLTermtmp = TMatrixT<double>(findex.size(),findex.size()); // Adding it;
+ BatemanMatrixDLTermtmp = BatemanMatrixDLTermN;
+
+ BatemanMatrixDLTermN.Mult(BatemanMatrixDLTermtmp, BatemanMatrix );
+
+ BatemanMatrixDLTermN *= 1./j;
+ BatemanMatrixDL += BatemanMatrixDLTermN;
+
+ NormN = 0;
+ for(int m = 0; m < (int)findex.size(); m++)
+ for(int n = 0; n < (int)findex.size(); n++)
+ NormN += BatemanMatrixDLTermN[m][n]*BatemanMatrixDLTermN[m][n];
+ j++;
+
+ } while ( NormN != 0 );
+ }
+ NEvolutionMatrix = BatemanMatrixDL * NEvolutionMatrix ;
+ }
+ NMatrix.push_back(NEvolutionMatrix);
+
}
-
- NEvolutionMatrix = BatemanMatrixDL * NMatrix.back() ;
- NMatrix.push_back(NEvolutionMatrix);
+
+ timevector[i+1] = timevector[i] + TStepMax;
+
+
}
-
-
+ FissionXSMatrix.push_back(GetFissionXsMatrix(EvolutionDataStep, DBTimeStep[NStep])); //Feel the reaction Matrix
+ CaptureXSMatrix.push_back(GetCaptureXsMatrix(EvolutionDataStep, DBTimeStep[NStep])); //Feel the reaction Matrix
+ n2nXSMatrix.push_back(Getn2nXsMatrix(EvolutionDataStep, DBTimeStep[NStep])); //Feel the reaction Matrix
+
+
EvolutionData GeneratedDB = EvolutionData(GetLog());
-
+
+
+ double ESigmaN = 0;
+ for (int j = 0; j < (int)findex.size() ; j++)
+ ESigmaN -= FissionXSMatrix[NStep-1][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++)
@@ -729,7 +822,7 @@ EvolutionData DataBank<IsotopicVector>::GenerateEvolutionData(IsotopicVector iso
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];
@@ -738,20 +831,10 @@ EvolutionData DataBank<IsotopicVector>::GenerateEvolutionData(IsotopicVector iso
}
GeneratedDB.NucleiInsert(pair<ZAI, TGraph*> (findex.find(i)->second, new TGraph(NMatrix.size(), timevector, ZAIQuantity)));
- if(fUpdateReferenceDBatEachStep)
- {
- 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)));
- }
- else
- {
- GeneratedDB.SetFissionXS(EvolutionDataStep.GetFissionXS());
- GeneratedDB.SetCaptureXS(EvolutionDataStep.GetCaptureXS());
- GeneratedDB.Setn2nXS(EvolutionDataStep.Getn2nXS());
-
+ 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 );
@@ -759,266 +842,2065 @@ EvolutionData DataBank<IsotopicVector>::GenerateEvolutionData(IsotopicVector iso
GeneratedDB.SetCycleTime(cycletime);
fDataBankCalculated.insert( pair< IsotopicVector, EvolutionData > ( GeneratedDB.GetIsotopicVectorAt(0.), GeneratedDB) );
-
+
+ //exit(1);
+ ResetTheMatrix();
+ ResetTheNucleiVector();
return GeneratedDB;
}
-//________________________________________________________________________
+//________________________________________________________________________
+
+template<>
+TMatrixT<double> DataBank<IsotopicVector>::GetFissionXsMatrix(EvolutionData EvolutionDataStep,double TStep)
+{
+
+ map<ZAI ,TGraph* >::iterator it;
+ TMatrixT<double> BatemanMatrix = TMatrixT<double>(findex.size(),findex.size());
+
+ // ---------------- A(n,.) X+Y
+
+ map<ZAI ,TGraph* > FissionXS = EvolutionDataStep.GetFissionXS();
+
+ for(it = FissionXS.begin() ; it != FissionXS.end(); it++)
+ {
+ map<ZAI, int>::iterator findex_inver_it = findex_inver.find( (*it).first );
+ if( findex_inver_it != findex_inver.end() )
+ {
+ double y = (*it).second->Eval(TStep);
+ BatemanMatrix[ findex_inver_it->second ][ findex_inver_it->second ] += -y* 1e-24;
+ BatemanMatrix[1][ findex_inver_it->second ] += 2*y* 1e-24;
+ }
+
+ }
+
+ return BatemanMatrix;
+
+}
+//________________________________________________________________________
+template<>
+TMatrixT<double> DataBank<IsotopicVector>::GetCaptureXsMatrix(EvolutionData EvolutionDataStep,double TStep)
+{
+
+
+ map<ZAI ,TGraph* >::iterator it;
+ TMatrixT<double> BatemanMatrix = TMatrixT<double>(findex.size(),findex.size());
+
+ map<ZAI, map<ZAI, double> > Capture;
+ { // 241Am
+ map<ZAI, double> toAdd ;
+ toAdd.insert(pair<ZAI, double> ( ZAI(96,242,0) , 0.8733*0.827) ); //directly cut the Am242 as in MURE
+ toAdd.insert(pair<ZAI, double> ( ZAI(94,242,0) , 0.8733*0.173) ); //directly cut the Am242 as in MURE
+ toAdd.insert(pair<ZAI, double> ( ZAI(95,242,1) , 0.1267) );
+ Capture.insert( pair< ZAI, map<ZAI, double> > ( ZAI(95,241,0), toAdd ) );
+ }
+ { // 242Am*
+ map<ZAI, double> toAdd ;
+ toAdd.insert(pair<ZAI, double> ( ZAI(95,243,0) , 1) );
+ Capture.insert( pair< ZAI, map<ZAI, double> > ( ZAI(95,242,1), toAdd ) );
+ }
+
+
+ // ---------------- A(n,.)A+1
+ map<ZAI ,TGraph* > CaptureXS = EvolutionDataStep.GetCaptureXS();
+ for(it = CaptureXS.begin(); it != CaptureXS.end(); it++)
+ {
+ map<ZAI, int>::iterator Index_it = findex_inver.find( (*it).first );
+ if( Index_it != findex_inver.end() )
+ {
+ double y;
+ y = (*it).second->Eval(TStep);
+
+ BatemanMatrix[Index_it->second][ Index_it->second ] += -y* 1e-24 ;
+
+ map<ZAI, map<ZAI, double> >::iterator it3 = Capture.find( (*it).first );
+
+ if( it3 == Capture.end() )
+ {
+ map<ZAI, int >::iterator it6 = findex_inver.find( ZAI( (*it).first.Z(), (*it).first.A()+1, (*it).first.I()) );
+
+ if( it6 != findex_inver.end() )
+ {
+ BatemanMatrix[(*it6).second][Index_it->second] += y* 1e-24 ;
+ }
+ else
+ {
+ map<ZAI, map<ZAI, double> >::iterator it4 = fFastDecay.find( ZAI( (*it).first.Z(), (*it).first.A()+1, (*it).first.I()) );
+
+ if( it4 == fFastDecay.end() )
+ {
+ cout << "Capture Problem in FastDecay for nuclei " << (*it).first.Z() << " " << (*it).first.A()+1 << " " << (*it).first.I() << endl;
+
+ BatemanMatrix[0][Index_it->second] += y* 1e-24 ;
+ }
+ else
+ {
+
+ map< ZAI, double >::iterator it5;
+ map< ZAI, double > decaylist2 = (*it4).second;
+ for(it5 = decaylist2.begin(); it5!= decaylist2.end(); it5++)
+ {
+ it6 = findex_inver.find( (*it5).first );
+ if( it6 == findex_inver.end() )
+ BatemanMatrix[0][Index_it->second] += y* 1e-24 * (*it5).second;
+ else
+ BatemanMatrix[(*it6).second][Index_it->second] += y* 1e-24 * (*it5).second;
+ }
+ }
+ }
+ }
+ else
+ {
+ map<ZAI, double>::iterator it4;
+ map<ZAI, double> CaptureList = (*it3).second;
+ for(it4 = CaptureList.begin(); it4 != CaptureList.end() ; it4++)
+ {
+
+ map<ZAI, int >::iterator it6 = findex_inver.find( (*it4).first );
+ if( it6 != findex_inver.end() )
+ BatemanMatrix[(*it6).second][Index_it->second] += y* 1e-24 * (*it4).second ;
+ else
+ {
+ map<ZAI, map<ZAI, double> >::iterator it7 = fFastDecay.find( (*it4).first );
+
+ if( it7 == fFastDecay.end() )
+ {
+ cout << "CaptureList Problem in FastDecay for nuclei " << (*it7).first.Z() << " " << (*it7).first.A() << " " << (*it7).first.I() << endl;
+ exit(1);
+ }
+
+ map< ZAI, double >::iterator it5;
+ map< ZAI, double > decaylist2 = (*it7).second;
+ for(it5 = decaylist2.begin(); it5!= decaylist2.end(); it5++)
+ {
+
+ it6 = findex_inver.find( (*it5).first );
+ if( it6 == findex_inver.end() )
+ {
+ cout << "CaptureList Problem in FastDecay for nuclei " << (*it7).first.Z() << " " << (*it7).first.A() << " " << (*it7).first.I() << endl;
+ exit(1);
+ }
+
+ BatemanMatrix[(*it6).second][Index_it->second] += y * 1e-24 * (*it5).second * (*it4).second;
+ }
+ }
+
+ }
+ }
+
+
+ }
+ }
+ return BatemanMatrix;
+
+}
+
+
+//________________________________________________________________________
+
+template<>
+TMatrixT<double> DataBank<IsotopicVector>::Getn2nXsMatrix(EvolutionData EvolutionDataStep,double TStep)
+{
+
+
+ map<ZAI ,TGraph* >::iterator it;
+ TMatrixT<double> BatemanMatrix = TMatrixT<double>(findex.size(),findex.size());
+
+ map<ZAI, map<ZAI, double> > n2n;
+ { // 237Np
+ map<ZAI, double> toAdd ;
+ toAdd.insert(pair<ZAI, double> ( ZAI(93,236,0) , 0.2) );
+ toAdd.insert(pair<ZAI, double> ( ZAI(93,236,1) , 0.8) );
+ n2n.insert( pair< ZAI, map<ZAI, double> > ( ZAI(93,237,0), toAdd ) );
+ }
+ { // 242Am*
+ map<ZAI, double> toAdd ;
+ toAdd.insert(pair<ZAI, double> ( ZAI(95,241,0) , 1) );
+ n2n.insert( pair< ZAI, map<ZAI, double> > ( ZAI(95,242,1), toAdd ) );
+ }
+
+ // ---------------- A(n,2n)A-1
+ map<ZAI ,TGraph* > n2nXS = EvolutionDataStep.Getn2nXS();
+ for(it = n2nXS.begin(); it != n2nXS.end(); it++)
+ {
+ map<ZAI, int>::iterator Index_it = findex_inver.find( (*it).first );
+ if( Index_it != findex_inver.end() )
+ {
+ double y;
+ y = (*it).second->Eval(TStep);
+
+ BatemanMatrix[Index_it->second][ Index_it->second ] += -y* 1e-24 ;
+
+ map<ZAI, map<ZAI, double> >::iterator it3 = n2n.find( (*it).first );
+
+ if( it3 == n2n.end() )
+ {
+ map<ZAI, int >::iterator it6 = findex_inver.find( ZAI( (*it).first.Z(), (*it).first.A()-1, (*it).first.I()) );
+
+ if( it6 != findex_inver.end() )
+ {
+ BatemanMatrix[(*it6).second][Index_it->second] += y* 1e-24 ;
+ }
+ else
+ {
+ map<ZAI, map<ZAI, double> >::iterator it4 = fFastDecay.find( ZAI( (*it).first.Z(), (*it).first.A()-1, (*it).first.I()) );
+
+ if( it4 == fFastDecay.end() )
+ {
+ cout << "n2n Problem in FastDecay for nuclei " << (*it).first.Z() << " " << (*it).first.A()-1 << " " << (*it).first.I() << endl;
+
+ BatemanMatrix[0][Index_it->second] += y* 1e-24 ;
+ }
+ else
+ {
+
+ map< ZAI, double >::iterator it5;
+ map< ZAI, double > decaylist2 = (*it4).second;
+ for(it5 = decaylist2.begin(); it5!= decaylist2.end(); it5++)
+ {
+ it6 = findex_inver.find( (*it5).first );
+ if( it6 == findex_inver.end() )
+ BatemanMatrix[0][Index_it->second] += y* 1e-24 * (*it5).second;
+ else
+ BatemanMatrix[(*it6).second][Index_it->second] += y* 1e-24 * (*it5).second;
+ }
+ }
+ }
+ }
+ else
+ {
+ map<ZAI, double>::iterator it4;
+ map<ZAI, double> n2nList = (*it3).second;
+ for(it4 = n2nList.begin(); it4 != n2nList.end() ; it4++)
+ {
+
+ map<ZAI, int >::iterator it6 = findex_inver.find( (*it4).first );
+ if( it6 != findex_inver.end() )
+ BatemanMatrix[(*it6).second][Index_it->second] += y* 1e-24 * (*it4).second ;
+ else
+ {
+ map<ZAI, map<ZAI, double> >::iterator it7 = fFastDecay.find( (*it4).first );
+
+ if( it7 == fFastDecay.end() )
+ {
+ cout << "n2nList Problem in FastDecay for nuclei " << (*it7).first.Z() << " " << (*it7).first.A() << " " << (*it7).first.I() << endl;
+ exit(1);
+ }
+
+ map< ZAI, double >::iterator it5;
+ map< ZAI, double > decaylist2 = (*it7).second;
+ for(it5 = decaylist2.begin(); it5!= decaylist2.end(); it5++)
+ {
+
+ it6 = findex_inver.find( (*it5).first );
+ if( it6 == findex_inver.end() )
+ {
+ cout << "n2nList Problem in FastDecay for nuclei " << (*it7).first.Z() << " " << (*it7).first.A() << " " << (*it7).first.I() << endl;
+ exit(1);
+ }
+
+ BatemanMatrix[(*it6).second][Index_it->second] += y * 1e-24 * (*it5).second * (*it4).second;
+ }
+ }
+
+ }
+ }
+
+
+ }
+ }
+ return BatemanMatrix;
+}
+
+
+template<>
+TMatrixT<double> DataBank<IsotopicVector>::ExtractXS(EvolutionData EvolutionDataStep,double TStep)
+{
+
+
+ map<ZAI ,TGraph* >::iterator it;
+ // ---------------- A(n,.) X+Y
+
+ map<ZAI ,TGraph* > FissionXS = EvolutionDataStep.GetFissionXS();
+
+ TMatrixT<double> SigmaPhi = TMatrixT<double>(findex.size()*3+1,1);
+ for(it = FissionXS.begin() ; it != FissionXS.end(); it++)
+ {
+
+ if( findex_inver.find( (*it).first ) != findex_inver.end() )
+ {
+ double y;
+ y = (*it).second->Eval(TStep);
+ SigmaPhi[findex_inver.find( (*it).first )->second][0] = y ;
+ }
+
+ }
+
+ // ---------------- A(n,.)A+1
+ map<ZAI ,TGraph* > CaptureXS = EvolutionDataStep.GetCaptureXS();
+ for(it = CaptureXS.begin(); it != CaptureXS.end(); it++)
+ {
+ if( findex_inver.find( (*it).first ) != findex_inver.end() )
+ {
+ double y;
+ y = (*it).second->Eval(TStep);
+ SigmaPhi[findex_inver.find( (*it).first )->second + findex.size() ][0] = y ;
+
+ }
+ }
+
+ // ---------------- A(n,2n)A-1
+ map<ZAI ,TGraph* > n2nXS = EvolutionDataStep.Getn2nXS();
+ for(it = n2nXS.begin() ; it != n2nXS.end(); it++)
+ {
+ if( findex_inver.find( (*it).first ) != findex_inver.end() )
+ {
+ double y;
+ y = (*it).second->Eval(TStep);
+ SigmaPhi[findex_inver.find( (*it).first )->second + findex.size() + findex.size()][0] = y ;
+
+ }
+ }
+ return SigmaPhi;
+}
+//________________________________________________________________________
+//________________________________________________________________________
+template<>
+string DataBank<IsotopicVector>::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;
+}
+
+
+
+template<>
+void DataBank<IsotopicVector>::BuildDecayMatrix()
+{
+ // 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!! !!!DataBank!!! \n Can't open \"" << DataFullPathName << "\"\n" << endl;
+ GetLog()->fLog << "!!Warning!! !!!DataBank!!! \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 >= 90 )
+ {
+ // 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);
+ map< ZAI, double > 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;}
+ 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 < 90 && 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)
+ {
+ pair<map<ZAI, double>::iterator, bool> IResult;
+ IResult = DaughtersMap.insert(pair<ZAI, double> (DaughterZAI , BR) );
+ if( !IResult.second)
+ (*IResult.first).second += BR;
+
+ branch_test+=BR;
+ }
+ else if( DM <= 18)
+ {
+ pair<map<ZAI, double>::iterator, bool> IResult;
+ IResult = DaughtersMap.insert(pair<ZAI, double> (ZAI(-2,-2,-2) , 2*BR) );
+ if( !IResult.second)
+ (*IResult.first).second += 2*BR;
+ 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 )
+ {
+
+ map< ZAI, double >::iterator DM_it = DaughtersMap.begin();
+ if(branch_test+branch_test_f>0)
+ for( DM_it = DaughtersMap.begin(); DM_it != DaughtersMap.end(); DM_it++)
+ {
+ if ( (*DM_it).first != ZAI(-2,-2,-2) )
+ (*DM_it).second *= 1./(branch_test+branch_test_f);
+ else
+ (*DM_it).second *= 1./(branch_test+branch_test_f);
+ }
+ }
+
+
+
+ if (HalfLife < fShorstestHalflife)
+ {
+ fFastDecay.insert( pair< ZAI, map<ZAI, double> > ( ParentZAI, DaughtersMap ) );
+ }
+ else
+ {
+ ZAIDecay.insert( pair< ZAI, pair<double, map< ZAI, double > > >
+ ( ParentZAI, pair<double, map< ZAI, double > >
+ ( HalfLife, DaughtersMap) ) );
+ }
+
+
+ }
+
+ } 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 == fFastDecay.end() )
+ {
+ (*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 ;
+ }
+ else
+ {
+ 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 ;
+ }
+
+ }
+
+ }
+ }
+
+ }
+
+ 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++;
+
+
+ }
+ }
+ //exit(1);
+
+}
+
+
+template<>
+void DataBank<IsotopicVector>::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;
+ }
+ }
+
+}
+
+template<>
+void DataBank<IsotopicVector>::ResetTheMatrix()
+{
+
+ if(fTheMatrix)
+ {
+ for(int i= 0; i<fNVar; i++)
+ delete [] fTheMatrix[i];
+ delete [] fTheMatrix;
+ }
+ fTheMatrix = 0;
+}
+
+
+template<>
+void DataBank<IsotopicVector>::SetTheNucleiVectorToZero()
+{
+ ResetTheNucleiVector();
+ fTheNucleiVector = new double[fNVar];
+
+#pragma omp parallel for
+ for(int i = 0; i < fNVar; i++)
+ fTheNucleiVector[i]=0.0;
+
+}
+
+template<>
+void DataBank<IsotopicVector>::ResetTheNucleiVector()
+{
+ if(fTheNucleiVector)
+ delete [] fTheNucleiVector;
+ fTheNucleiVector = 0;
+}
+
+
+template<>
+void DataBank<IsotopicVector>::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;
+ }
+}
+
+template<>
+void DataBank<IsotopicVector>::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];
+}
+
+template<>
+TMatrixT<double> DataBank<IsotopicVector>::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;
+}
+
+
+template<>
+void DataBank<IsotopicVector>::SetTheNucleiVector(TMatrixT<double> NEvolutionMatrix)
+{
+ for (int k = 0; k < (int)fNVar; k++)
+ fTheNucleiVector[k] = NEvolutionMatrix[k][0];
+}
+
+template<>
+TMatrixT<double> DataBank<IsotopicVector>::GetTheNucleiVector()
+{
+ TMatrixT<double> NEvolutionMatrix = TMatrixT<double>(findex.size(),1);
+ for (int k = 0; k < (int)fNVar; k++)
+ NEvolutionMatrix[k][0] = fTheNucleiVector[k];
+
+ return NEvolutionMatrix;
+}
+
+
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+
+/*
+
+
+ OLD READ
+
+
+ */
+
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+//________________________________________________________________________
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-template<>
-TMatrixT<double> DataBank<IsotopicVector>::GetFissionXsMatrix(EvolutionData EvolutionDataStep,double TStep)
-{
-
- map<ZAI ,TGraph* >::iterator it;
- TMatrixT<double> BatemanMatrix = TMatrixT<double>(findex.size(),findex.size());
-
- // ---------------- A(n,.) X+Y
-
- map<ZAI ,TGraph* > FissionXS = EvolutionDataStep.GetFissionXS();
-
- for(it = FissionXS.begin() ; it != FissionXS.end(); it++)
- {
- map<ZAI, int>::iterator findex_inver_it = findex_inver.find( (*it).first );
- if( findex_inver_it != findex_inver.end() )
- {
- double y = (*it).second->Eval(TStep);
-
- BatemanMatrix[ findex_inver_it->second ][ findex_inver_it->second ] += -y* 1e-24;
- BatemanMatrix[1][ findex_inver_it->second ] += 2*y* 1e-24;
- }
-
- }
-
- return BatemanMatrix;
-
-}
-//________________________________________________________________________
-template<>
-TMatrixT<double> DataBank<IsotopicVector>::GetCaptureXsMatrix(EvolutionData EvolutionDataStep,double TStep)
-{
- map<ZAI ,TGraph* >::iterator it;
- TMatrixT<double> BatemanMatrix = TMatrixT<double>(findex.size(),findex.size());
-
- map<ZAI, map<ZAI, double> > Capture;
- { // 241Am
- map<ZAI, double> toAdd ;
- toAdd.insert(pair<ZAI, double> ( ZAI(95,242,0) , 0.086) );
- toAdd.insert(pair<ZAI, double> ( ZAI(95,242,1) , 0.914) );
- Capture.insert( pair< ZAI, map<ZAI, double> > ( ZAI(95,241,0), toAdd ) );
- }
- { // 242Am*
- map<ZAI, double> toAdd ;
- toAdd.insert(pair<ZAI, double> ( ZAI(95,243,0) , 1) );
- Capture.insert( pair< ZAI, map<ZAI, double> > ( ZAI(95,242,1), toAdd ) );
- }
-
-
- // ---------------- A(n,.)A+1
- map<ZAI ,TGraph* > CaptureXS = EvolutionDataStep.GetCaptureXS();
- for(it = CaptureXS.begin(); it != CaptureXS.end(); it++)
- {
- if( findex_inver.find( (*it).first ) != findex_inver.end() )
- {
- double y;
- y = (*it).second->Eval(TStep);
-
- BatemanMatrix[findex_inver.find( (*it).first )->second][ findex_inver.find( (*it).first )->second ] += -y* 1e-24 ;
-
- map<ZAI, map<ZAI, double> >::iterator it3 = Capture.find( (*it).first );
-
- if( it3 == Capture.end() )
- {
- map<ZAI, int >::iterator it6 = findex_inver.find( ZAI( (*it).first.Z(), (*it).first.A()+1, (*it).first.I()) );
-
- if( it6 != findex_inver.end() )
- {
- BatemanMatrix[(*it6).second][findex_inver.find( (*it).first )->second] += y* 1e-24 ;
- }
- else
- {
- map<ZAI, map<ZAI, double> >::iterator it4 = fFastDecay.find( ZAI( (*it).first.Z(), (*it).first.A()+1, (*it).first.I()) );
-
- if( it4 == fFastDecay.end() )
- {
- cout << "Problem in FastDecay for nuclei " << (*it).first.Z() << " " << (*it).first.A()+1 << " " << (*it).first.I() << endl;
- exit(1);
- }
-
- map< ZAI, double >::iterator it5;
- map< ZAI, double > decaylist2 = (*it4).second;
- for(it5 = decaylist2.begin(); it5!= decaylist2.end(); it5++)
- {
- it6 = findex_inver.find( (*it5).first );
- if( it6 == findex_inver.end() )
- {
- cout << "Problem in FastDecay for nuclei " << (*it).first.Z() << " " << (*it).first.A() << " " << (*it).first.I() << endl;
- exit(1);
- }
- BatemanMatrix[(*it6).second][findex_inver.find( (*it).first )->second] += y* 1e-24 * (*it5).second;
- }
- }
- }
- else
- {
- map<ZAI, double>::iterator it4;
- map<ZAI, double> CaptureList = (*it3).second;
- for(it4 = CaptureList.begin(); it4 != CaptureList.end() ; it4++)
- {
-
- map<ZAI, int >::iterator it6 = findex_inver.find( (*it4).first );
- if( it6 != findex_inver.end() )
- BatemanMatrix[(*it6).second][findex_inver.find( (*it).first )->second] += y* 1e-24 * (*it4).second ;
- else
- {
- map<ZAI, map<ZAI, double> >::iterator it7 = fFastDecay.find( (*it4).first );
-
- if( it7 == fFastDecay.end() )
- {
- cout << "Problem in FastDecay for nuclei " << (*it7).first.Z() << " " << (*it7).first.A() << " " << (*it7).first.I() << endl;
- exit(1);
- }
-
- map< ZAI, double >::iterator it5;
- map< ZAI, double > decaylist2 = (*it7).second;
- for(it5 = decaylist2.begin(); it5!= decaylist2.end(); it5++)
- {
-
- it6 = findex_inver.find( (*it5).first );
- if( it6 == findex_inver.end() )
- {
- cout << "Problem in FastDecay for nuclei " << (*it7).first.Z() << " " << (*it7).first.A() << " " << (*it7).first.I() << endl;
- exit(1);
- }
-
- BatemanMatrix[(*it6).second][findex_inver.find( (*it).first )->second] += y * 1e-24 * (*it5).second * (*it4).second;
- }
- }
-
- }
- }
-
-
- }
- }
- return BatemanMatrix;
-
-}
-//________________________________________________________________________
-template<>
-TMatrixT<double> DataBank<IsotopicVector>::Getn2nXsMatrix(EvolutionData EvolutionDataStep,double TStep)
-{
- map<ZAI ,TGraph* >::iterator it;
- TMatrixT<double> BatemanMatrix = TMatrixT<double>(findex.size(),findex.size());
-
- // ---------------- A(n,2n)A-1
- map<ZAI ,TGraph* > n2nXS = EvolutionDataStep.Getn2nXS();
- for(it = n2nXS.begin() ; it != n2nXS.end(); it++)
- {
- if( findex_inver.find( (*it).first ) != findex_inver.end() )
- {
- double y;
- y = (*it).second->Eval(TStep);
- BatemanMatrix[ findex_inver.find( (*it).first )->second ][findex_inver.find( (*it).first )->second] += -y* 1e-24 ;
-
-
- map<ZAI, int>::iterator it3 = findex_inver.find( ZAI( (*it).first.Z(), (*it).first.A()-1, 0) );
-
- if( it3 != findex_inver.end() )
- BatemanMatrix[(*it3).second][findex_inver.find( (*it).first )->second] += y* 1e-24 ;
- else
- {
-
- map<ZAI, map<ZAI, double> >::iterator it4 = fFastDecay.find( ZAI( (*it).first.Z(), (*it).first.A()-1, 0) );
-
- if( it4 == fFastDecay.end() )
- {
- it3 = findex_inver.find( ZAI( -3, -3, -3 ) );
- BatemanMatrix[(*it3).second][findex_inver.find( (*it).first )->second] += y* 1e-24;
- }
- 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() )
- {
- cout << "Problem in FastDecay for nuclei " << (*it4).first.Z() << " " << (*it4).first.A() << " " << (*it4).first.I() << endl;
- exit(1);
- }
- BatemanMatrix[(*it3).second][findex_inver.find( (*it).first )->second] += y* 1e-24 * (*it5).second ;
- }
- }
- }
- }
- }
- return BatemanMatrix;
-}
-template<>
-TMatrixT<double> DataBank<IsotopicVector>::ExtractXS(EvolutionData EvolutionDataStep,double TStep)
-{
-
- map<ZAI ,TGraph* >::iterator it;
- // ---------------- A(n,.) X+Y
-
- map<ZAI ,TGraph* > FissionXS = EvolutionDataStep.GetFissionXS();
-
- TMatrixT<double> SigmaPhi = TMatrixT<double>(findex.size()*3+1,1);
- for(it = FissionXS.begin() ; it != FissionXS.end(); it++)
- {
-
- if( findex_inver.find( (*it).first ) != findex_inver.end() )
- {
- double y;
- y = (*it).second->Eval(TStep);
- SigmaPhi[findex_inver.find( (*it).first )->second][0] = y ;
- }
-
- }
-
- // ---------------- A(n,.)A+1
- map<ZAI ,TGraph* > CaptureXS = EvolutionDataStep.GetCaptureXS();
- for(it = CaptureXS.begin(); it != CaptureXS.end(); it++)
- {
- if( findex_inver.find( (*it).first ) != findex_inver.end() )
- {
- double y;
- y = (*it).second->Eval(TStep);
- SigmaPhi[findex_inver.find( (*it).first )->second + findex.size() ][0] = y ;
-
- }
- }
-
- // ---------------- A(n,2n)A-1
- map<ZAI ,TGraph* > n2nXS = EvolutionDataStep.Getn2nXS();
- for(it = n2nXS.begin() ; it != n2nXS.end(); it++)
- {
- if( findex_inver.find( (*it).first ) != findex_inver.end() )
- {
- double y;
- y = (*it).second->Eval(TStep);
- SigmaPhi[findex_inver.find( (*it).first )->second + findex.size() + findex.size()][0] = y ;
-
- }
- }
- return SigmaPhi;
-}
-//________________________________________________________________________
-//________________________________________________________________________
template<>
-void DataBank<IsotopicVector>::BuildDecayMatrix()
+void DataBank<IsotopicVector>::OldBuildDecayMatrix()
{
// List of Decay Time and Properties
map<ZAI, pair<double, map< ZAI, double > > > ZAIDecay;
@@ -1299,9 +3181,14 @@ void DataBank<IsotopicVector>::BuildDecayMatrix()
}
}
+
+
+
+
}
+
template<>
EvolutionData DataBank<IsotopicVector>::OldGenerateEvolutionData(IsotopicVector isotopicvector, double cycletime, double Power)
{
@@ -1637,7 +3524,7 @@ EvolutionData DataBank<IsotopicVector>::OldGenerateEvolutionData(IsotopicVector
EvolutionData EvolutionDataStep = GetClosest(isotopicvector.GetActinidesComposition(), 0.); //GetCLosest at the begining of evolution
-
+
ReactorType = EvolutionDataStep.GetReactorType();
for(int i = 0; i < 16; i++)
@@ -1654,7 +3541,7 @@ EvolutionData DataBank<IsotopicVector>::OldGenerateEvolutionData(IsotopicVector
for(int k = 0; k < (int)index.size(); k++)
IVStep += index.find(k)->second * NMatrix.back()[k][0];
- if(fUpdateReferenceDBatEachStep && i != 0); //GetCLosest at the each of evolution step (begining already done...)
+ if(i != 0); //GetCLosest at the each of evolution step (begining already done...)
EvolutionDataStep = GetClosest(IVStep, TStep);
double NormFactor = 1;
@@ -1749,7 +3636,6 @@ EvolutionData DataBank<IsotopicVector>::OldGenerateEvolutionData(IsotopicVector
}
else
{
- //if( (*it3).first.Z() == 90 && (*it3).first.A() == 232) cout << y* 1e-24 *Flux << endl;
map<ZAI, double>::iterator it4;
map<ZAI, double> CaptureList = (*it3).second;
for(it4 = CaptureList.begin(); it4 != CaptureList.end() ; it4++)
@@ -1779,7 +3665,6 @@ EvolutionData DataBank<IsotopicVector>::OldGenerateEvolutionData(IsotopicVector
cout << "Problem in FastDecay for nuclei " << (*it7).first.Z() << " " << (*it7).first.A() << " " << (*it7).first.I() << endl;
exit(1);
}
- //if( (*it6).first.Z() == 92 && (*it6).first.A() == 233) cout << y* 1e-24 *Flux * (*it5).second << endl;
BatemanMatrix[(*it6).second][index_inver.find( (*it).first )->second] += y * 1e-24 * Flux * (*it5).second * (*it4).second;
}
}
diff --git a/source/trunk/src/DynamicalSystem.cxx b/source/trunk/src/DynamicalSystem.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..a2f0efad124a22b7d2254cce8571749665ef0ab3
--- /dev/null
+++ b/source/trunk/src/DynamicalSystem.cxx
@@ -0,0 +1,214 @@
+#include <cstdio>
+#include <iostream>
+#include <fstream>
+#include <cmath>
+#include <vector>
+#include <algorithm>
+
+#include "DynamicalSystem.hxx"
+
+using namespace std;
+
+//________________________________________________________________________
+DynamicalSystem::DynamicalSystem()
+{
+
+ SetPrecision();
+ fHestimate=1.; //this value is change in DynamicalSystem::RungeKutta
+ fHmin=0.;
+ fMaxHdid=-1e30;
+ fMinHdid=1e30;
+ fIsNegativeValueAllowed=true;
+}
+
+//________________________________________________________________________
+DynamicalSystem::DynamicalSystem(const DynamicalSystem & DS)
+{
+ fNVar=DS.fNVar;
+ fPrecision=DS.fPrecision;
+ fHestimate=DS.fHestimate;
+ fHmin=DS.fHmin;
+ fMaxHdid=DS.fMaxHdid;
+ fMinHdid=DS.fMinHdid;
+ fIsNegativeValueAllowed=DS.fIsNegativeValueAllowed;
+}
+
+//________________________________________________________________________
+DynamicalSystem::~DynamicalSystem()
+{
+
+}
+
+//________________________________________________________________________
+void DynamicalSystem::RungeKutta(double *YStart,double t1, double t2, int EquationNumber)
+{
+ //double shortestHalfLife=gMURE->GetShortestHalfLife();
+ fNVar = EquationNumber;
+ int nstp;
+ double t,hnext,hdid,h;
+ double *yscal=new double[fNVar];
+ double *y=new double[fNVar];
+ double *dydt=new double[fNVar];
+
+ const double MAXSTP = 10000;
+ const double TINY = 1.0e-30;
+
+
+
+ if(fabs(t1-t2)<=TINY)
+ cout << "Integration time is 0." << endl;
+
+ t=t1;
+ fHestimate=(t2-t1)/100;
+ //h=(t2 > t1) ? fabs(fHestimate) : -fabs(fHestimate);
+ h=fHestimate;
+
+ // pragma omp parallel for
+ for (int i = 0; i < fNVar; i++) y[i] = YStart[i];
+
+ for ( nstp = 1; nstp <= MAXSTP; nstp++)
+ {
+ BuildEqns(t,y,dydt);
+
+ // pragma omp parallel for
+ for ( int i = 0; i < fNVar; i++)
+ yscal[i] = fabs(y[i]) + fabs(dydt[i]*h)+TINY;
+
+ if ( (t+h-t2) * (t+h-t1) > 0.0 ) h=t2-t;
+
+ AdaptStepSize(y,dydt,&t,h,fPrecision,yscal,&hdid,&hnext);
+
+ if(fMaxHdid<hdid) fMaxHdid=hdid;
+ if(fMinHdid>hdid) fMinHdid=hdid;
+
+ if ((t-t2)*(t2-t1) >= 0.0)
+ {
+ // pragma omp parallel for
+ for (int i=0;i<fNVar;i++) YStart[i]=y[i];
+
+ delete [] dydt;
+ delete [] y;
+ delete [] yscal;
+ //cout << "The maximum step used in RK was "<<fMaxHdid<<" Step NUM in RK was "<<nstp << endl;
+ return;
+ }
+ if (fabs(hnext) <= fHmin)
+ cout << "Step size too small in RungeKutta" << endl;
+
+ h=hnext;
+ }
+ cout << "Too many steps in routine RungeKutta" << endl;
+}
+
+//________________________________________________________________________
+void DynamicalSystem::RK4(double *y, double *dydx, double x, double h, double *yout)
+{
+ //cout<<"Calling Function RK4"<<endl;
+ double xh,hh,h6;
+
+ double *dym=new double[fNVar];
+ double *dyt=new double[fNVar];
+ double *yt=new double[fNVar];
+ hh=h*0.5;
+ h6=h/6.0;
+ xh=x+hh;
+
+ // pragma omp parallel for
+ for (int i=0;i<fNVar;i++) yt[i]=y[i]+hh*dydx[i];
+
+ BuildEqns(xh,yt,dyt);
+
+ // pragma omp parallel for
+ for (int i=0;i<fNVar;i++) yt[i]=y[i]+hh*dyt[i];
+
+ BuildEqns(xh,yt,dym);
+
+ for (int i=0;i<fNVar;i++)
+ {
+ yt[i]=y[i]+h*dym[i];
+ dym[i] += dyt[i];
+ }
+
+ BuildEqns(x+h,yt,dyt);
+ // pragma omp parallel for
+ for (int i=0;i<fNVar;i++)
+ {
+ yout[i]=y[i]+h6*(dydx[i]+dyt[i]+2.0*dym[i]);
+ if(!fIsNegativeValueAllowed && yout[i]<0)
+ {
+ cout << "Material composition is negative "
+ <<"i="<<i<<" ("/*<<fEvolvingMaterial->GetComposition()[i]->GetZAI()->PrintName()
+ */<<") old proportion="<<y[i]<<" new="<<yout[i]
+ <<". Setting to 0." << endl;
+ yout[i]=0.;
+ }
+ }
+ delete [] yt;
+ delete [] dyt;
+ delete [] dym;
+}
+
+//________________________________________________________________________
+void DynamicalSystem::AdaptStepSize(double *y, double *dydx, double *x, double htry,
+ double eps, double *yscal, double *hdid, double *hnext)
+{
+ //cout<<"Calling Function AdaptStepSize()"<<endl;
+ double xsav,hh,h,temp,errmax;
+ double *dysav=new double[fNVar];
+ double *ysav=new double[fNVar];
+ double *ytemp=new double[fNVar];
+
+ const double PGROW =-0.20;
+ const double PSHRNK =-0.25;
+ const double FCOR =0.06666666 ;
+ const double SAFETY =0.9;
+ const double ERRCON =6.0e-4;
+
+ xsav=(*x);
+
+ // pragma omp parallel for
+ for (int i=0;i<fNVar;i++)
+ {
+ ysav[i]=y[i];
+ dysav[i]=dydx[i];
+ }
+ h=htry;
+ for (;;)
+ {
+ hh=0.5*h;
+ RK4(ysav,dysav,xsav,hh,ytemp);
+ *x=xsav+hh;
+
+ BuildEqns(*x,ytemp,dydx);
+ RK4(ytemp,dydx,*x,hh,y);
+ *x=xsav+h;
+ if (*x == xsav )
+ {
+ //cout << "Step size ("<<h<<") too small in routine AdaptStepSize" << endl;
+ }
+ RK4(ysav,dysav,xsav,h,ytemp);
+ errmax=0.0;
+
+ for (int i=0;i<fNVar;i++)
+ {
+ ytemp[i]=y[i]-ytemp[i];
+ temp=fabs(ytemp[i]/yscal[i]);
+ if (errmax < temp) errmax=temp;
+ }
+ errmax /= eps;
+ if (errmax <= 1.0)
+ {
+ *hdid=h;
+ *hnext=(errmax > ERRCON ? SAFETY*h*exp(PGROW*log(errmax)) : 4.0*h);
+ break;
+ }
+ h=SAFETY*h*exp(PSHRNK*log(errmax));
+ }
+
+ for (int i=0;i<fNVar;i++) y[i] += ytemp[i]*FCOR;
+
+ delete [] ytemp;
+ delete [] dysav;
+ delete [] ysav;
+}
+
diff --git a/source/trunk/src/EvolutionData.cxx b/source/trunk/src/EvolutionData.cxx
index 017f06c2b8a7a5a80a1de3529db3b87a16944544..cf240d812e3d38f26a4257e1270be4383e560ddb 100755
--- a/source/trunk/src/EvolutionData.cxx
+++ b/source/trunk/src/EvolutionData.cxx
@@ -98,6 +98,10 @@ EvolutionData::EvolutionData(LogFile* Log)
SetLog(Log);
fIsCrossSection = false;
+ fPower = 0;
+ fCycleTime = 0;
+ fNormFactor = 0;
+
}
@@ -108,6 +112,9 @@ EvolutionData::EvolutionData(LogFile* Log, string DB_file, bool oldread, ZAI zai
SetLog(Log);
fIsCrossSection = false;
fDB_file = DB_file;
+ fPower = 0;
+ fCycleTime = 0;
+ fNormFactor = 0;
if(zai != ZAI(0,0,0))
AddAsStable(zai);
@@ -265,6 +272,7 @@ void EvolutionData::ReadDB(string DBfile, bool oldread)
OldReadDB(DBfile);
return;
}
+ ReadInfo(); // Read the .info associeted
ifstream DecayDB(DBfile.c_str()); // Open the File
if(!DecayDB) //check if file is correctly open
@@ -348,7 +356,7 @@ void EvolutionData::ReadDB(string DBfile, bool oldread)
}while ( !DecayDB.eof() );
-
+
}
@@ -374,7 +382,7 @@ void EvolutionData::ReadKeff(string line, double* time, int NTimeStep)
i++;
}
- fFlux = new TGraph(NTimeStep, time, Keff); // Add the TGraph
+ fKeff = new TGraph(NTimeStep, time, Keff); // Add the TGraph
@@ -406,7 +414,6 @@ void EvolutionData::ReadFlux(string line, double* time, int NTimeStep)
fFlux = new TGraph(NTimeStep, time, Flux); // Add the TGraph
- ReadInfo(); // Read the .info associeted
}
@@ -435,7 +442,7 @@ void EvolutionData::ReadInv(string line, double* time, int NTimeStep)
int i = 0;
while(start < (int)line.size() && i < NTimeStep ) // Read the Data
{
- Inv[i] = atof(StringLine::NextWord(line, start, ' ').c_str()) ;
+ Inv[i] = atof(StringLine::NextWord(line, start, ' ').c_str())*fNormFactor ;
i++;
}
// Add the TGraph
@@ -587,8 +594,8 @@ void EvolutionData::ReadInfo()
getline(InfoDB, line);
if ( tlc(StringLine::NextWord(line, start, ' ')) == "normalizationfactor")
{
- double NormFactor = atof(StringLine::NextWord(line, start, ' ').c_str());
- fPower = fPower * NormFactor;
+ fNormFactor = atof(StringLine::NextWord(line, start, ' ').c_str());
+ fPower = fPower * fNormFactor;
}
}
diff --git a/source/trunk/src/FabricationPlant.cxx b/source/trunk/src/FabricationPlant.cxx
index ba8b7a96d419b2bd004a264a2c1fd5155371b5de..2c554ae228631ae9633d3f6184d41f52b57f0f0f 100644
--- a/source/trunk/src/FabricationPlant.cxx
+++ b/source/trunk/src/FabricationPlant.cxx
@@ -163,8 +163,7 @@ void FabricationPlant::FabricationPlantEvolution(cSecond t)
}
}
- fInsideIV = GetFullFabrication();
-
+
}
@@ -196,7 +195,6 @@ void FabricationPlant::BuildFuelForReactor(int ReactorId)
{
double nPu_0 = 0;
double MPu_0 = 0;
-
{
map<ZAI ,double>::iterator it;
@@ -330,14 +328,9 @@ void FabricationPlant::BuildFuelForReactor(int ReactorId)
for(int i = (int)fFractionToTake.size()-1; i >= 0; i--)
{
- IVBeginCycle += fStorage->GetStock()[fFractionToTake[i].first].GetSpeciesComposition(94)*( fFractionToTake[i].second );
- IsotopicVector UnusedIV = fStorage->GetStock()[fFractionToTake[i].first]*(fFractionToTake[i].second)
- - fStorage->GetStock()[fFractionToTake[i].first].GetSpeciesComposition(94)*(fFractionToTake[i].second) ;
-
- pair<IsotopicVector, IsotopicVector> SepatationIV = Separation(UnusedIV);
-
- fReUsable->AddToStock( SepatationIV.first );
- GetParc()->AddWaste( SepatationIV.second );
+ IVBeginCycle += fStorage->GetStock()[fFractionToTake[i].first].GetSpeciesComposition(94)*( fFractionToTake[i].second );
+ fReUsable->AddToStock(fStorage->GetStock()[fFractionToTake[i].first]*(fFractionToTake[i].second)
+ - fStorage->GetStock()[fFractionToTake[i].first].GetSpeciesComposition(94)*(fFractionToTake[i].second));
fStorage->TakeFractionFromStock(fFractionToTake[i].first,fFractionToTake[i].second);
@@ -418,9 +411,9 @@ EvolutionData FabricationPlant::BuildEvolutiveDB(int ReactorId,IsotopicVector is
EvolutionData EvolBuild;
- EvolBuild = evolutiondb->GenerateEvolutionData(isotopicvector,
- GetParc()->GetReactor()[ReactorId]->GetCycleTime(),
- GetParc()->GetReactor()[ReactorId]->GetPower());
+ EvolBuild = evolutiondb->GenerateEvolutionData(isotopicvector,
+ GetParc()->GetReactor()[ReactorId]->GetCycleTime(),
+ GetParc()->GetReactor()[ReactorId]->GetPower());
return EvolBuild;
}
@@ -494,6 +487,7 @@ void FabricationPlant::RecycleStock(double fraction)
}
+
//________________________________________________________________________
void FabricationPlant::DumpStock()
{
diff --git a/source/trunk/src/IsotopicVector.cxx b/source/trunk/src/IsotopicVector.cxx
index 60cac2675f3f7ae3749b74f504c6c2295b54d561..1a2eaff7a63277b5431e4cc452f51f58471d6a1a 100755
--- a/source/trunk/src/IsotopicVector.cxx
+++ b/source/trunk/src/IsotopicVector.cxx
@@ -284,6 +284,18 @@ void IsotopicVector::Multiply(double factor)
}
+//________________________________________________________________________
+
+double IsotopicVector::GetSumOfAll()
+{
+ double Sum = 0;
+ map<ZAI ,double >::iterator it;
+ for( it = fIsotopicQuantity.begin(); it != fIsotopicQuantity.end(); it++)
+ Sum += (*it).second;
+
+ return Sum;
+
+}
//________________________________________________________________________
void IsotopicVector::Add(const ZAI& zai, double quantity)
{
@@ -297,7 +309,7 @@ void IsotopicVector::Add(const ZAI& zai, double quantity)
{
pair<map<ZAI, double>::iterator, bool> IResult;
IResult = fIsotopicQuantity.insert( pair<ZAI ,double>(zai, quantity));
- if(IResult.second == false)
+ if(!IResult.second)
IResult.first->second += quantity;
}
@@ -337,10 +349,6 @@ void IsotopicVector::Remove(const ZAI& zai, double quantity)
map<ZAI ,double>::iterator it;
it = fIsotopicQuantity.find(zai);
- if( ceil(quantity*1e25) - quantity*1e25 > quantity*1e25 - floor(quantity*1e25) )
- quantity = floor(quantity*1e25)*1/1e25;
- else quantity = ceil(quantity*1e25)*1/1e25;
-
if(quantity > 0)
{
if ( it != fIsotopicQuantity.end() )
@@ -376,12 +384,6 @@ void IsotopicVector::Remove(const IsotopicVector& isotopicvector)
//________________________________________________________________________
void IsotopicVector::Need(const ZAI& zai, double quantity)
{
-
- if( ceil(quantity*1e25) - quantity*1e25 > quantity*1e25 - floor(quantity*1e25) )
- quantity = floor(quantity*1e25)*1/1e25;
- else quantity = ceil(quantity*1e25)*1/1e25;
-
-
pair<map<ZAI, double>::iterator, bool> IResult;
if(quantity > 0)
{
diff --git a/source/trunk/src/Makefile b/source/trunk/src/Makefile
index 9dde40f5a161538a9a965ddd80dcbe4ffc750a08..b672e8944308d92f4d758860432f3a13534194e5 100755
--- a/source/trunk/src/Makefile
+++ b/source/trunk/src/Makefile
@@ -23,9 +23,9 @@ OBJS = CLASS.o \
IsotopicVector.o IsotopicVectorDict.o \
ZAIMass.o \
ZAI.o ZAIDict.o \
- DataBank.o \
+ DataBank.o DynamicalSystem.o \
EvolutionData.o EvolutionDataDict.o \
- LogFile.o
+ LogFile.o
ROOTOBJS = CLSSObject.o CLSSObjectDict.o\
CLSSFacility.o CLSSFacilityDict.o\
@@ -36,7 +36,7 @@ ROOTOBJS = CLSSObject.o CLSSObjectDict.o\
IsotopicVector.o IsotopicVectorDict.o \
ZAIMass.o \
ZAI.o ZAIDict.o \
- DataBank.o \
+ DataBank.o DynamicalSystem.o\
EvolutionData.o EvolutionDataDict.o \
LogFile.o
@@ -83,7 +83,6 @@ EvolutionDataDict.cxx: $(LOCALINC)/EvolutionData.hxx
-
clean:
@rm -vf $(OBJS) *~ core *Dict.cxx *Dict.h
oclean :
diff --git a/source/trunk/src/Reactor.cxx b/source/trunk/src/Reactor.cxx
index 3c75dce2d9e327d0bb307d7c6818c5ec0684abd7..efa7bc86c2b822d599006ab01cb035f466f7df5b 100755
--- a/source/trunk/src/Reactor.cxx
+++ b/source/trunk/src/Reactor.cxx
@@ -213,7 +213,6 @@ Reactor::Reactor(LogFile* log, EvolutionData evolutivedb,
fFixedFuel = true;
fIsStorage = false;
- fFabricationPlant = 0;
fAssociedPool = Pool;
fInternalTime = 0;
@@ -350,7 +349,7 @@ void Reactor::Evolution(cSecond t)
if(fInternalTime == 0 && fIsStarted == false) // Start of the Reactor
{
fEndOfCycle = true;
- fInsideIV = fIVBeginCycle; ///useless ??
+ fInsideIV = fIVBeginCycle;
fInternalTime = t;
}