diff --git a/source/branches/CLASSV3/include/CLASSHeaders.hxx b/source/branches/CLASSV3/include/CLASSHeaders.hxx
index 4766ce48378be43c913d6417ebceefee02245ad4..46739c8be5355e832a68721bc2d9ddc3ba62fbd3 100755
--- a/source/branches/CLASSV3/include/CLASSHeaders.hxx
+++ b/source/branches/CLASSV3/include/CLASSHeaders.hxx
@@ -11,8 +11,6 @@
#include "Reactor.hxx"
#include "Pool.hxx"
#include "FabricationPlant.hxx"
-#include "PWR_THU_FabricationPlant.hxx"
-#include "PWR_THPU_FabricationPlant.hxx"
#include "Storage.hxx"
#include "IsotopicVector.hxx"
#include "ZAI.hxx"
@@ -20,7 +18,6 @@
#include "EvolutionData.hxx"
#include "DecayDataBank.hxx"
-#include "FuelDataBank.hxx"
#include "StringLine.hxx"
#include "ZAIMass.hxx"
diff --git a/source/branches/CLASSV3/src/EQM_LIN_PWR_MOX.cxx b/source/branches/CLASSV3/src/EQM_LIN_PWR_MOX.cxx
index c4443f35e2b72713b1c2fd1d5dcf551eeb282ebc..5fc357c4ff3b0fdf5d4776c6b8027710e44dac27 100644
--- a/source/branches/CLASSV3/src/EQM_LIN_PWR_MOX.cxx
+++ b/source/branches/CLASSV3/src/EQM_LIN_PWR_MOX.cxx
@@ -4,8 +4,8 @@
#include "StringLine.hxx"
#include "LogFile.hxx"
-
-
+#include "IsotopicVector.hxx"
+#include "CLASSHeaders.hxx"
@@ -74,7 +74,7 @@ vector<double> EQM_LIN_PWR_MOX::BuildFuel(double BurnUp, double HMMass,vector<Is
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
- // ADD ENrichment of the U check !!!!!!!!!!!!!!!!!!!!!!!!!!!! //
+ // ADD ENrichment of the U check ++ Check Un seul fertile !!!! //
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
@@ -82,6 +82,119 @@ vector<double> EQM_LIN_PWR_MOX::BuildFuel(double BurnUp, double HMMass,vector<Is
vector<double> lambda;
+ for(int i = 0; i < (int) (FissilArray.size() + FertilArray.size()); i++)
+ lambda.push_back(0);
+
+
+ double Na = 6.02214129e23; //N Avogadro
+
+ IsotopicVector FullUsedStock;
+ IsotopicVector stock;
+
+ bool FuelBuild = false;
+
+ int N_FissilStock_OnCheck = 0;
+
+ while(!FuelBuild)
+ {
+ double nPu_0 = 0;
+ double MPu_0 = 0;
+ {
+ map<ZAI ,double>::iterator it;
+
+ map<ZAI ,double> isotopicquantity = FullUsedStock.GetSpeciesComposition(94).GetIsotopicQuantity();
+ for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
+ nPu_0 += (*it).second;
+
+ isotopicquantity = (FullUsedStock.GetSpeciesComposition(94) + ZAI(94,241,0)*FullUsedStock.GetZAIIsotopicQuantity(95,241,0)).GetIsotopicQuantity(); //Add the 241Am as 241Pu... the Pu is not old in the Eq Model but is in the FissileArray....;
+ for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
+ MPu_0 += (*it).second*cZAIMass.fZAIMass.find( (*it).first )->second/Na*1e-6;
+ }
+
+ stock = FissilArray[N_FissilStock_OnCheck];
+ double nPu_1 = 0;
+ double MPu_1 = 0;
+ double Sum_AlphaI_nPuI = 0;
+ double Sum_AlphaI_nPuI0 = 0;
+ {
+ map<ZAI ,double>::iterator it;
+ map<ZAI ,double> isotopicquantity = stock.GetSpeciesComposition(94).GetIsotopicQuantity();
+
+ for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
+ {
+ if ((*it).first.A() >= 238 && (*it).first.A() <= 242)
+ {
+ nPu_1 += (*it).second;
+ Sum_AlphaI_nPuI += fFuelParameter[(*it).first.A() -237]*(*it).second;
+ }
+ }
+
+ isotopicquantity = (stock.GetSpeciesComposition(94) + ZAI(94,241,0)*stock.GetZAIIsotopicQuantity(95,241,0)).GetIsotopicQuantity(); //Add the 241Am as 241Pu... the Pu is not old in the Eq Model but is in the FissileArray....
+ for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
+ if ((*it).first.A() >= 238 && (*it).first.A() <= 242)
+ {
+ MPu_1 += (*it).second * (cZAIMass.fZAIMass.find( (*it).first )->second)/Na*1e-6;
+ }
+
+ isotopicquantity = FullUsedStock.GetSpeciesComposition(94).GetIsotopicQuantity();
+ for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
+ if ((*it).first.A() >= 238 && (*it).first.A() <= 242)
+ {
+ Sum_AlphaI_nPuI0 += fFuelParameter[(*it).first.A() -237]*(*it).second;
+ }
+ }
+
+ double StockFactionToUse = 0;
+
+ double NT = HMMass*1e6 * Na / (cZAIMass.fZAIMass.find( ZAI(92,238,0) )->second*0.997 + cZAIMass.fZAIMass.find( ZAI(92,235,0) )->second*0.003 );
+
+ double N1 = (BurnUp - fFuelParameter[6]) * NT;
+ double N2 = -Sum_AlphaI_nPuI0;
+ double N3 = -fFuelParameter[0] * Na / (cZAIMass.fZAIMass.find( ZAI(92,238,0) )->second*0.997 + cZAIMass.fZAIMass.find( ZAI(92,235,0) )->second*0.003 ) * (HMMass*1e6 - MPu_0*1e6);
+
+ double D1 = Sum_AlphaI_nPuI;
+ double D2 = -fFuelParameter[0] * MPu_1*1e6 * Na / (cZAIMass.fZAIMass.find( ZAI(92,238,0) )->second*0.997 + cZAIMass.fZAIMass.find( ZAI(92,235,0) )->second*0.003 ) ;
+
+ StockFactionToUse = (N1 + N2 + N3) / (D1 + D2);
+
+ if(StockFactionToUse < 0)
+ {
+ cout << "!!Bad Trouble!! !!!FabricationPlant!!! Oups Bug in calculating stock fraction to use "<< endl;
+ GetLog()->fLog << "!!Bad Trouble!! !!!FabricationPlant!!! Oups Bug in calculating stock fraction to use" << endl;
+ lambda[N_FissilStock_OnCheck] = 0.;
+ FuelBuild = false;
+ }
+ else if( StockFactionToUse > 1 )
+ {
+
+ FullUsedStock += stock;
+ lambda[N_FissilStock_OnCheck] = 1;
+ N_FissilStock_OnCheck++;
+ FuelBuild = false;
+
+
+ if( N_FissilStock_OnCheck > (int) lambda.size() ) // Check if the last Fissil stock has been tested... quit if so...
+ {
+ for(int i = 0; i < (int)lambda.size(); i++)
+ lambda[i] = -1;
+
+ FuelBuild = true;
+ }
+ }
+ else
+ {
+ lambda[N_FissilStock_OnCheck] = StockFactionToUse;
+
+ FuelBuild = true;
+
+ double U8_Quantity = (HMMass - (MPu_0+StockFactionToUse*MPu_1 ))/(cZAIMass.fZAIMass.find( ZAI(92,238,0) )->second*0.997 + cZAIMass.fZAIMass.find( ZAI(92,235,0) )->second*0.003 )*Na/1e-6;
+
+ lambda.back() = U8_Quantity / FertilArray[0].GetSumOfAll();
+ }
+
+ }
+
+
return lambda;
}