Main Page   Compound List   File List   Compound Members   File Members  

MyLiSource.cpp

Go to the documentation of this file.
00001 
00008 #include <math.h>
00009 #include <iostream.h>
00010 #include <fstream.h>   // file I/O
00011 #include <iomanip.h>   // format manipulation
00012 #include <complex>
00013 #include <string>
00014 #include "ReactorConstants.hh"
00015 #include "ReactorFortran.hh"
00016 #include "MyLiSource.hh"
00017 
00018 MyLiSource::MyLiSource(int newIsotope,
00019                        int newRmaxgen){
00020 
00021   Isotope = newIsotope;
00022   Rmaxgen = newRmaxgen;
00023 
00024   //cout << "MyLiSource::MyLiSource(" << Isotope << "," << Rmaxgen 
00025   //     << ")" << endl;
00026 
00027   // 9Li
00028   if(Isotope == 9){
00029 
00030     // setup the decay parameters (from Jon Link)
00031     // Nneutron is always 1 and Ngamma always 0
00032     Nneutron = 1;
00033     Ngamma = 0;
00034 
00035     // branching fraction (the 1.0000 bin contains all decays 
00036     // which do nothing: neutron E = 0, beta Q = 0)
00037     float static br[8] = {0.0238,0.3400,0.4400,0.4475,
00038                           0.4550,0.4894,0.4950,1.0000};
00039 
00040     // Q-value for beta decay
00041     float static betaQ[8] = {11.176,11.176,10.826,5.666,
00042                              5.666,2.323,2.323,0.000};
00043 
00044     // Q-value for the total decay chain
00045     float static totQ = 12.033;
00046 
00047     // neutron energy
00048     float static neuE[8] = {0.764,0.856,1.115,3.235,
00049                             5.473,8.816,6.578,0.000};
00050 
00051     // prepare the beta decay functions
00052     bool static first=1;
00053     float static fspace[7][200];
00054 
00055     if(first){
00056       float xlo = 0.;
00057       funlxp_(BetaHighestQ,fspace[0],xlo,betaQ[0]);
00058       funlxp_(BetaHighestQ,fspace[1],xlo,betaQ[1]);
00059       funlxp_(BetaHighQ,fspace[2],xlo,betaQ[2]);
00060       funlxp_(BetaMiddleQ,fspace[3],xlo,betaQ[3]);
00061       funlxp_(BetaMiddleQ,fspace[4],xlo,betaQ[4]);
00062       funlxp_(BetaLowQ,fspace[5],xlo,betaQ[5]);
00063       funlxp_(BetaLowQ,fspace[6],xlo,betaQ[6]);
00064 
00065       first = 0;
00066     }
00067     //
00068     // let the 9Li decay
00069     //
00070     int ntrial = 0;
00071     bool passed = false;
00072     while(!passed){
00073       int len=1;
00074       float rv[len];
00075       ranlux_(rv,len);
00076       //
00077       for(int i=0;i<8; i++){
00078         if(rv[0] < br[i] && !passed){
00079 
00080           // pick an energy for each neutron
00081           neutronE = neuE[i];
00082 
00083           if(neutronE){
00084 
00085             // get the energy of the beta
00086             len = 1;
00087             funlux_(fspace[i],betaE,len);
00088             passed = true;
00089           }
00090         }
00091       }
00092       ntrial++;
00093     }
00094     Weight = 1./float(ntrial);
00095 
00096     //cout << "   9Li neutron energy = " << neutronE << " MeV" << endl;
00097     //cout << "   9Li beta energy = " << betaE << " MeV" << endl;
00098     //cout << "   9Li weight = " << Weight << endl;
00099 
00100     // pick a decay time in ns
00101     double HalfLife = 1.8e8;
00102 
00103     int len=1;
00104     float rv[len];
00105     ranlux_(rv,len);
00106     DecayTime = HalfLife*log(1/rv[0]);
00107 
00108   } // done with 9Li
00109 }
00110 
00111 MyLiSource::~MyLiSource(){
00112   ;
00113 }
00114 
00115 MyLiSource::MyLiSource(const MyLiSource& LiSource){
00116 
00117   Isotope = LiSource.Isotope;
00118   Rmaxgen = LiSource.Rmaxgen;
00119 
00120   Nneutron = LiSource.Nneutron;
00121   Nbeta    = LiSource.Nbeta;
00122   Ngamma   = LiSource.Ngamma;
00123   neutronE = LiSource.neutronE;
00124   betaE    = LiSource.betaE;
00125   gammaE   = LiSource.gammaE;
00126 
00127   Weight = LiSource.Weight;
00128   DecayTime = LiSource.DecayTime;
00129 } 
00130    
00131 
00132 MyLiSource& MyLiSource::operator=(const MyLiSource& rhs){
00133 
00134   if(this == &rhs)return *this;
00135 
00136   Isotope = rhs.Isotope;  
00137   Rmaxgen = rhs.Rmaxgen;
00138 
00139   Nneutron = rhs.Nneutron;
00140   Nbeta    = rhs.Nbeta;
00141   Ngamma   = rhs.Ngamma;
00142   neutronE = rhs.neutronE;
00143   betaE    = rhs.betaE;
00144   gammaE   = rhs.gammaE;
00145 
00146   Weight = rhs.Weight;
00147   DecayTime = rhs.DecayTime;
00148 
00149   return *this;
00150 }
00151 
00152 float MyLiSource::BetaHighestQ(float const& x){
00153 
00154   // beta decay function from Hannah and Jon
00155   static double me = k.Melectron;
00156   float beta = 0.;
00157   int A = 9;
00158   int Z = 3;
00159   double Q = 11.176;
00160 
00161   // x -> 0, beta -> 4.5, but Fermi -> inf so don't call it
00162   if(x == 0.)
00163     return 4.5;
00164   else{
00165     beta = sqrt(x*x+2*x*me)*(Q-x)*(Q-x)*(x+me)*Fermi(Z,A,x);
00166     return beta;
00167   }
00168 }
00169 
00170 float MyLiSource::BetaHighQ(float const& x){
00171 
00172   // beta decay function from Hannah and Jon
00173   static double me = k.Melectron;
00174   float beta = 0.;
00175   int A = 9;
00176   int Z = 3;
00177   double Q = 10.826;
00178 
00179   // x -> 0, beta -> 4.222, but Fermi -> inf so don't call it
00180   if(x == 0.)
00181     return 4.222;
00182   else{
00183     beta = sqrt(x*x+2*x*me)*(Q-x)*(Q-x)*(x+me)*Fermi(Z,A,x);
00184     return beta;
00185   }
00186 }
00187 
00188 float MyLiSource::BetaMiddleQ(float const& x){
00189 
00190   // beta decay function from Hannah and Jon
00191   static double me = k.Melectron;
00192   float beta = 0.;
00193   int A = 9;
00194   int Z = 3;
00195   double Q = 5.666;
00196 
00197   // x -> 0, beta -> 1.1565, but Fermi -> inf so don't call it
00198   if(x == 0.)
00199     return 1.1565;
00200   else{
00201     beta = sqrt(x*x+2*x*me)*(Q-x)*(Q-x)*(x+me)*Fermi(Z,A,x);
00202     return beta;
00203   }
00204 }
00205 
00206 float MyLiSource::BetaLowQ(float const& x){
00207 
00208   // beta decay function from Hannah and Jon
00209   static double me = k.Melectron;
00210   float beta = 0.;
00211   int A = 9;
00212   int Z = 3;
00213   double Q = 2.323;
00214 
00215   // x -> 0, beta -> 0.1944, but Fermi -> inf so don't call it
00216   if(x == 0.)
00217     return 0.1944;
00218   else{
00219     beta = sqrt(x*x+2*x*me)*(Q-x)*(Q-x)*(x+me)*Fermi(Z,A,x);
00220     return beta;
00221   }
00222 }
00223 
00224 double MyLiSource::Fermi(int const& Z, int const& A, float const& x){
00225 
00226   // Fermi function for beta decay from Hannah and Jon
00227   static double me    = k.Melectron;
00228   static double pi    = k.pi;
00229   static double alpha = 1./k.alphainv;
00230   double fermi = 0.;
00231 
00232   double xx = x/me+1.;
00233   double pe = sqrt(xx*xx-1.);
00234   double g0 = 1.-(alpha*Z)*(alpha*Z);
00235   double rnu = alpha*Z*xx/pe;
00236   double r = pow(A,1./3.);
00237 
00238   complex<double> arg(g0,rnu);
00239   complex<double> cg = wgamma_(arg);
00240   double abscgsq = real(conj(cg)*cg);
00241 
00242   fermi = exp(pi*rnu)*abscgsq*pow((pe*r),(2.*(g0-1.)));
00243   return fermi;
00244 }

Generated on Mon Dec 27 11:10:13 2004 for ReactorFsim by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002