Main Page   Compound List   File List   Compound Members   File Members  

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

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