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 "MyCfSource.hh"
00016
00017 long long factorial(int n);
00018
00019 MyCfSource::MyCfSource(int newIsotope,
00020 double newRmaxgen){
00021
00022 Isotope = newIsotope;
00023 Rmaxgen = newRmaxgen;
00024
00025
00026
00027
00028
00029 if(Isotope == 252){
00030
00031
00032 bool static first=1;
00033 float static fspace[200];
00034 float static mspace[200];
00035 float static gspace[200];
00036
00037 if(first){
00038 float xlo = 0.;
00039 float xhi = 50.;
00040 funlxp_(Cf252NeutronSpectrum,fspace,xlo,xhi);
00041 xhi = 25.;
00042 funlxp_(Cf252GammaMultiplicityFit,mspace,xlo,xhi);
00043 funlxp_(Cf252GammaSpectrum,gspace,xlo,xhi);
00044 first = 0;
00045 }
00046
00047
00048 fstream cf252_position("data/cf252_position.txt", ios::in);
00049
00050
00051 Nsources = 0;
00052 char paramLine[256];
00053 string dummy;
00054 bool done = false;
00055 bool numSet = false;
00056 int positionNum = 0;
00057
00058 while(!done){
00059
00060
00061 if(cf252_position.peek() == '#'){
00062
00063 cf252_position.getline(paramLine, 256);
00064 continue;
00065 }
00066
00067
00068 if(cf252_position.peek() == 'e' || cf252_position.peek() == 'E'){
00069
00070 if(Nsources == 0) cout << " Cf252 ERROR: "
00071 << "number of cf252 sources = 0" << endl;
00072 done = true;
00073 continue;
00074 }
00075
00076
00077 if(cf252_position.peek() == 'n' || cf252_position.peek() == 'N'){
00078
00079 if(numSet){
00080
00081 cf252_position.getline( paramLine, 256 );
00082 continue;
00083 }
00084 cf252_position >> dummy;
00085 while(cf252_position.peek() == ' ' || cf252_position.peek() == '=')
00086 cf252_position.ignore(1);
00087
00088 cf252_position >> Nsources;
00089 if(Nsources == 0){
00090 cout << " Cf252 ERROR: number of cf252 sources = 0" << endl;
00091 done = true;
00092 continue;
00093 }
00094
00095 numSet = true;
00096 continue;
00097 }
00098
00099
00100 if(cf252_position.peek() == 'p' || cf252_position.peek() == 'P'){
00101
00102 if(Nsources == 0 || positionNum >= Nsources){
00103
00104 cout << " Cf252 ERROR: too many positions " << positionNum+1
00105 << " given for the number of sources " << Nsources << endl;
00106 cf252_position.getline( paramLine, 256 );
00107 positionNum++;
00108 continue;
00109 }
00110 cf252_position >> dummy;
00111 while(cf252_position.peek() == ' ' || cf252_position.peek() == '=')
00112 cf252_position.ignore(1);
00113
00114 cf252_position >> R[positionNum] >> zR[positionNum]
00115 >> phiR[positionNum];
00116
00117
00118 if(R[positionNum] > Rmaxgen){
00119 cout << " Cf252 ERROR: R[" << positionNum << "] of "
00120 << R[positionNum] << " > Rmaxgen" << endl;
00121 done = true;
00122 continue;
00123 }
00124 if(zR[positionNum] < -1. || zR[positionNum] > 1.){
00125 cout << " Cf252 ERROR: CosTheta["
00126 << positionNum << "] of " << zR[positionNum] << endl;
00127 done = true;
00128 continue;
00129 }
00130 if(phiR[positionNum] > 2*PI){
00131 cout << " Cf252 ERROR: Phi[" << positionNum << "] of "
00132 << phiR[positionNum] << " > 2*pi" << endl;
00133 done = true;
00134 continue;
00135 }
00136
00137 positionNum++;
00138 continue;
00139 }
00140
00141
00142 if(cf252_position.peek() == ' ' || cf252_position.peek() == '\n'){
00143
00144 cf252_position.ignore(1);
00145 continue;
00146 }
00147 }
00148
00149 if(positionNum < Nsources){
00150 cout << " Cf252 ERROR: not enough positions " << positionNum
00151 << " given for the number of sources " << Nsources << endl;
00152 }
00153
00154
00155
00156 double p[maxNeutron+1] = {0.002,0.028,0.155,0.428,
00157 0.732,0.917,0.983,0.998,1.000};
00158
00159 for(int n=0;n<Nsources;n++){
00160
00161
00162
00163
00164 bool passed = false;
00165 int len;
00166 while(!passed){
00167 const int len = 1;
00168 float r[len];
00169 ranlux_(r,len);
00170
00171 for(int i=0;i<maxNeutron+1; i++){
00172
00173 if(p[i] > r[0] && !passed){
00174 Nneutron[n] = i;
00175 if(Nneutron[n])passed = true;
00176 }
00177 }
00178 }
00179
00180
00181
00182
00183 for(int nn=0;nn<Nneutron[n]; nn++){
00184 len = 1;
00185 funlux_(fspace,neutronE[n][nn],len);
00186
00187 Tneutron[n][nn] = 0.;
00188 }
00189
00190
00191
00192
00193
00194 float Z = 98.;
00195 float A = 252.;
00196
00197 double avge;
00198 double phi;
00199 double etot;
00200
00201 avge = -1.33 + 119.6*pow(Z,0.3333)/A;
00202 phi = 2.51 - 1.13e-5*Z*Z*sqrt(A);
00203 etot = phi*Nneutron[n]+4.0;
00204
00205
00206
00207 len = 1;
00208 float ng = 0.;
00209 funlux_(mspace,ng,len);
00210 Ngamma[n] = (int)ng;
00211
00212
00213
00214
00215 double tote = 0.;
00216 for(int nn=0;nn<Ngamma[n]; nn++){
00217 int len = 1;
00218 funlux_(gspace,gammaE[n][nn],len);
00219
00220 tote += gammaE[n][nn];
00221
00222 const int len2 = 2;
00223 float rv[len2];
00224 ranlux_(rv,len2);
00225
00226
00227
00228 double halflife;
00229 if(rv[0] < 0.8)
00230 halflife = 0.01;
00231 else
00232 halflife = 1.0;
00233 Tgamma[n][nn] = halflife*log(1/rv[1]);
00234 }
00235
00236 }
00237
00238
00239 for(int n=0;n<Nsources;n++){
00240 for(int nn=Nneutron[n];nn<maxNeutron;nn++){
00241 neutronE[n][nn] = 0.;
00242 Tneutron[n][nn] = 0.;
00243 }
00244 for(int nn=Ngamma[n];nn<maxGamma;nn++){
00245 gammaE[n][nn] = 0.;
00246 Tgamma[n][nn] = 0.;
00247 }
00248 }
00249 for(int n=Nsources;n<Ndim;n++){
00250 Nneutron[n] = 0;
00251 Ngamma[n] = 0;
00252 R[n] = 0.;
00253 zR[n] = 0.;
00254 phiR[n] = 0.;
00255 for(int nn=0;nn<maxNeutron;nn++){
00256 neutronE[n][nn] = 0.;
00257 Tneutron[n][nn] = 0.;
00258 }
00259 for(int nn=0;nn<maxGamma;nn++){
00260 gammaE[n][nn] = 0.;
00261 Tgamma[n][nn] = 0.;
00262 }
00263 }
00264
00265 }
00266
00267 if(Isotope == 255){
00268
00269 Nsources = 0;
00270 for(int n=Nsources;n<Ndim;n++){
00271 Nneutron[n] = 0;
00272 Ngamma[n] = 0;
00273 R[n] = 0.;
00274 zR[n] = 0.;
00275 phiR[n] = 0.;
00276 for(int nn=0;nn<maxNeutron;nn++){
00277 neutronE[n][nn] = 0.;
00278 Tneutron[n][nn] = 0.;
00279 }
00280 for(int nn=0;nn<maxGamma;nn++){
00281 gammaE[n][nn] = 0.;
00282 Tgamma[n][nn] = 0.;
00283 }
00284 }
00285 }
00286 }
00287
00288 MyCfSource::~MyCfSource(){
00289 ;
00290 }
00291
00292 MyCfSource::MyCfSource(const MyCfSource& CfSource){
00293
00294 Isotope = CfSource.Isotope;
00295 Rmaxgen = CfSource.Rmaxgen;
00296
00297 Nsources = CfSource.Nsources;
00298
00299 for(int n=0;n<Ndim;n++){
00300 Nneutron[n] = CfSource.Nneutron[n];
00301 Ngamma[n] = CfSource.Ngamma[n];
00302 R[n] = CfSource.R[n];
00303 zR[n] = CfSource.zR[n];
00304 phiR[n] = CfSource.phiR[n];
00305
00306 for(int nn=0;nn<maxNeutron;nn++){
00307 neutronE[n][nn] = CfSource.neutronE[n][nn];
00308 Tneutron[n][nn] = CfSource.Tneutron[n][nn];
00309 }
00310 for(int nn=0;nn<maxGamma;nn++){
00311 gammaE[n][nn] = CfSource.gammaE[n][nn];
00312 Tgamma[n][nn] = CfSource.Tgamma[n][nn];
00313 }
00314 }
00315 }
00316
00317 MyCfSource& MyCfSource::operator=(const MyCfSource& rhs){
00318
00319 if(this == &rhs)return *this;
00320
00321 Isotope = rhs.Isotope;
00322 Rmaxgen = rhs.Rmaxgen;
00323
00324 Nsources = rhs.Nsources;
00325
00326 for(int n=0;n<Ndim;n++){
00327 Nneutron[n] = rhs.Nneutron[n];
00328 Ngamma[n] = rhs.Ngamma[n];
00329 R[n] = rhs.R[n];
00330 zR[n] = rhs.zR[n];
00331 phiR[n] = rhs.phiR[n];
00332
00333 for(int nn=0;nn<maxNeutron;nn++){
00334 neutronE[n][nn] = rhs.neutronE[n][nn];
00335 Tneutron[n][nn] = rhs.Tneutron[n][nn];
00336 }
00337 for(int nn=0;nn<maxGamma;nn++){
00338 gammaE[n][nn] = rhs.gammaE[n][nn];
00339 Tgamma[n][nn] = rhs.Tgamma[n][nn];
00340 }
00341 }
00342 return *this;
00343 }
00344
00345 float MyCfSource::Cf252NeutronSpectrum(const float& x){
00346
00347
00348 float N = 0.;
00349
00350 float scale = 1/(2*PI*0.359);
00351
00352
00353 float fminus = exp(-(sqrt(x)-sqrt(0.359))*(sqrt(x)-sqrt(0.359))/1.175);
00354 float fplus = exp(-(sqrt(x)+sqrt(0.359))*(sqrt(x)+sqrt(0.359))/1.175);
00355
00356 N = scale*(fminus-fplus);
00357
00358
00359 return N;
00360 }
00361
00362 float MyCfSource::Cf252GammaMultiplicity(const int& x){
00363
00364
00365 double M = 0.;
00366
00367 double C1 = 0.675;
00368 double C2 = 6.78;
00369 double C3 = 9.92;
00370
00371 M = C1*pow(C2,x)*exp(-C2)/factorial(x) + (1-C1)*pow(C3,x)*exp(-C3)/factorial(x);
00372
00373
00374
00375 return M;
00376 }
00377
00378 long long factorial(int n){
00379
00380
00381 if(n < 1)
00382 return 1;
00383 else
00384 return n * factorial(n-1);
00385 }
00386
00387 float MyCfSource::Cf252GammaMultiplicityFit(const float& x){
00388
00389
00390 double M = 0.;
00391
00392 float gaussian = 0.13345*exp(-0.5*((x-7.0322)/2.6301)*((x-7.0322)/2.6301));
00393 float exponential = 0.11255*(exp(-(sqrt(x)-sqrt(7.5987))
00394 *(sqrt(x)-sqrt(7.5987))/0.56213));
00395
00396 if(x <= 9.0){
00397 M = gaussian;
00398 }
00399 else{
00400 M = exponential;
00401 }
00402
00403
00404
00405 return M;
00406 }
00407
00408 float MyCfSource::Cf252GammaSpectrum(const float& x){
00409
00410
00411 float N = 0.;
00412
00413 float gaussian = 1.8367*exp(-0.5*((x-0.45934)/0.31290)*((x-0.45934)/0.31290));
00414 float exponential = exp(0.84774-0.89396*x);
00415
00416 if(x <= 0.744){
00417 N = gaussian;
00418 }
00419 else{
00420 N = exponential;
00421 }
00422
00423
00424 return N;
00425 }