00001
00008 #include <math.h>
00009 #include <iostream.h>
00010 #include <fstream.h>
00011 #include <iomanip.h>
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
00024
00025
00026
00027 if(Isotope == 8){
00028
00029
00030
00031 Nneutron = 1;
00032 Ngamma = 0;
00033
00034
00035
00036 float static br[4] = {0.080,0.160,0.169,1.0000};
00037
00038
00039 float static betaQ[4] = {7.442,5.252,1.852,0.000};
00040
00041
00042 float static totQ = 8.619;
00043
00044
00045 float static neuE[4] = {0.697,2.887,6.727,0.000};
00046
00047
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
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 rv[0] = 0.167;
00070 for(int i=0;i<4; i++){
00071 if(rv[0] < br[i] && !passed){
00072
00073
00074 neutronE = neuE[i];
00075
00076 if(neutronE){
00077
00078
00079 len = 1;
00080 funlux_(fspace[i],betaE,len);
00081 passed = true;
00082 }
00083 }
00084 }
00085 ntrial++;
00086 }
00087 Weight = 1./float(ntrial);
00088
00089
00090
00091
00092
00093
00094 double HalfLife = 1.2e8;
00095
00096 int len=1;
00097 float rv[len];
00098 ranlux_(rv,len);
00099 DecayTime = HalfLife*log(1/rv[0]);
00100
00101 }
00102 }
00103
00104 MyHeSource::~MyHeSource(){
00105 ;
00106 }
00107
00108 MyHeSource::MyHeSource(const MyHeSource& HeSource){
00109
00110 Isotope = HeSource.Isotope;
00111 Rmaxgen = HeSource.Rmaxgen;
00112
00113 Nneutron = HeSource.Nneutron;
00114 Nbeta = HeSource.Nbeta;
00115 Ngamma = HeSource.Ngamma;
00116 neutronE = HeSource.neutronE;
00117 betaE = HeSource.betaE;
00118 gammaE = HeSource.gammaE;
00119
00120 Weight = HeSource.Weight;
00121 DecayTime = HeSource.DecayTime;
00122 }
00123
00124 MyHeSource& MyHeSource::operator=(const MyHeSource& rhs){
00125
00126 if(this == &rhs)return *this;
00127
00128 Isotope = rhs.Isotope;
00129 Rmaxgen = rhs.Rmaxgen;
00130
00131 Nneutron = rhs.Nneutron;
00132 Nbeta = rhs.Nbeta;
00133 Ngamma = rhs.Ngamma;
00134 neutronE = rhs.neutronE;
00135 betaE = rhs.betaE;
00136 gammaE = rhs.gammaE;
00137
00138 Weight = rhs.Weight;
00139 DecayTime = rhs.DecayTime;
00140
00141 return *this;
00142 }
00143
00144 float MyHeSource::BetaHighQ(float const& x){
00145
00146
00147 static double me = k.Melectron;
00148 float beta = 0.;
00149 int A = 8;
00150 int Z = 2;
00151 double Q = 7.442;
00152
00153
00154 if(x == 0.)
00155 return 4.222;
00156 else{
00157 beta = sqrt(x*x+2*x*me)*(Q-x)*(Q-x)*(x+me)*Fermi(Z,A,x);
00158 return beta;
00159 }
00160 }
00161
00162 float MyHeSource::BetaMiddleQ(float const& x){
00163
00164
00165 static double me = k.Melectron;
00166 float beta = 0.;
00167 int A = 8;
00168 int Z = 2;
00169 double Q = 5.252;
00170
00171
00172 if(x == 0.)
00173 return 1.1565;
00174 else{
00175 beta = sqrt(x*x+2*x*me)*(Q-x)*(Q-x)*(x+me)*Fermi(Z,A,x);
00176 return beta;
00177 }
00178 }
00179
00180 float MyHeSource::BetaLowQ(float const& x){
00181
00182
00183 static double me = k.Melectron;
00184 float beta = 0.;
00185 int A = 8;
00186 int Z = 2;
00187 double Q = 1.852;
00188
00189
00190 if(x == 0.)
00191 return 0.1944;
00192 else{
00193 beta = sqrt(x*x+2*x*me)*(Q-x)*(Q-x)*(x+me)*Fermi(Z,A,x);
00194 return beta;
00195 }
00196 }
00197
00198 double MyHeSource::Fermi(int const& Z, int const& A, float const& x){
00199
00200
00201 static double me = k.Melectron;
00202 static double pi = k.pi;
00203 static double alpha = 1./k.alphainv;
00204 double fermi = 0.;
00205
00206 double xx = x/me+1.;
00207 double pe = sqrt(xx*xx-1.);
00208 double g0 = 1.-(alpha*Z)*(alpha*Z);
00209 double rnu = alpha*Z*xx/pe;
00210 double r = pow(A,1./3.);
00211
00212 complex<double> arg(g0,rnu);
00213 complex<double> cg = wgamma_(arg);
00214 double abscgsq = real(conj(cg)*cg);
00215
00216 fermi = exp(pi*rnu)*abscgsq*pow((pe*r),(2.*(g0-1.)));
00217 return fermi;
00218 }