Main Page   Compound List   File List   Compound Members   File Members  

NeutronPropagator.cpp

00001 #include <math.h>
00002 #include <iostream.h>
00003 #include <fstream.h>
00004 #include <string.h>
00005 #include "ReactorFortran.hh"
00006 #include "ReactorConstants.hh"
00007 #include "ReactorNeutrons.hh"
00008 //void NeutronPropagator$(double,double,double&,double&,double*,
00009 //                      double*,double&,double,double*,double*,double*);
00010 void NeutronPropagator(double,double,double&,double&,double*,
00011                         double*,double&,double,double*,double*,double*);
00012 void NeutronPropagator(double R0,double R2,
00013                        double& ke,double& A,double& passes,
00014                        double* pn,double* rn,double& tn,
00015                        double Temperature,double* a,double* z,double* nn){
00016 
00017   double r=sqrt(rn[0]*rn[0]+rn[1]*rn[1]+rn[2]*rn[2]);
00018   while(ke>0 && r<R2){
00019     passes++;
00020     NeutronPropagator(R0,R2,ke,A,pn,rn,tn,Temperature,a,z,nn);
00021     r=sqrt(rn[0]*rn[0]+rn[1]*rn[1]+rn[2]*rn[2]);
00022   }
00023   return;
00024 }
00025 //void NeutronPropagator$(double R0,double R2,
00026 //                      double& ke,double& A,
00027 //                      double* pn,double* rn,double& tn,
00028 //                      double Temperature,double* a,double* z,double* nn){
00029 void NeutronPropagator(double R0,double R2,
00030                         double& ke,double& A,
00031                         double* pn,double* rn,double& tn,
00032                         double Temperature,double* a,double* z,double* nn){
00033   //
00034   if(ke==0)return;
00035   //
00036   //
00037   double r=sqrt(rn[0]*rn[0]+rn[1]*rn[1]+rn[2]*rn[2]);
00038   if(r>R2)return;
00039   double rstart=r;
00040   // removed following section so that when you reset gd concentration,
00041   // these numbers get updated as well, passed through from ReactorDetector
00042   //
00043   // Set up cross section and materials data
00044   //
00045   static ReactorNeutrons XCdata;
00046   static bool first=1;
00047   //static double a[4];
00048   //static double z[4];
00049   //static double nn[4];
00050   //
00051   // Initialize on first call
00052   // 
00053   if(first){
00054     /*
00055       a[0] = RC.A0;
00056       a[1] = RC.A1;
00057       a[2] = RC.A2;
00058       a[3] = RC.A3;
00059       //
00060       z[0] = RC.Z0;
00061       z[1] = RC.Z1;
00062       z[2] = RC.Z2;
00063       z[3] = RC.Z3;
00064       //
00065       nn[0]=RC.Navogadro*RC.f0*RC.density/RC.A0;
00066       nn[1]=RC.Navogadro*RC.f1*RC.density/RC.A1;
00067       nn[2]=RC.Navogadro*RC.GdConcentration/100*RC.f2*RC.density/RC.A2;
00068       nn[3]=RC.Navogadro*RC.GdConcentration/100*RC.f3*RC.density/RC.A3;
00069     */
00070     ReactorNeutrons Init(0);
00071     XCdata=Init;
00072     first=0;
00073   }
00074   //
00075   // fetch cross section info
00076   //
00077   int nproc=8;
00078   double probproc[8];
00079   for(int i=0;i<nproc/2;i++){
00080     probproc[2*i]=nn[i]*XCdata.GetElastXC(int(z[i]),int(a[i]),1e6*ke);
00081     probproc[2*i+1]=nn[i]*XCdata.GetGammaXC(int(z[i]),int(a[i]),1e6*ke);
00082   }
00083   //
00084   // No Gd if outside inner volume
00085   //
00086   if(r>R0)for(int i=4;i<nproc;i++)probproc[i]=0;
00087   //
00088   // Normalize
00089   //
00090   double probtot=0;
00091   for(int i=0;i<nproc;i++)probtot+=probproc[i];
00092   for(int i=0;i<nproc;i++)probproc[i]/=probtot;
00093   for(int i=1;i<nproc;i++)probproc[i]+=probproc[i-1];
00094   //
00095   // Take steps until there is an interaction
00096   //
00097   double step=SlowNeutronStep;
00098   double mn = MPROTON+DELTA;
00099   //
00100   double mfp=1/probtot/1e-24;
00101   double probint = 1-exp(-step/mfp);
00102   int len=5;
00103   float rv[5];
00104   double t[3];
00105   double p=sqrt(pn[0]*pn[0]+pn[1]*pn[1]+pn[2]*pn[2]);
00106   double v=p/mn*c;
00107   for(int i=0;i<3;i++)t[i]=pn[i]/p;
00108   bool interacted=0;
00109   while(!interacted){
00110     ranlux_(rv,len);
00111     interacted=(rv[0]<probint);
00112     if(interacted){
00113       //
00114       // throw the exponential within the step 
00115       //
00116       double step2=-mfp*log(1-probint*rv[1]);
00117       for(int i=0;i<3;i++)rn[i]+=t[i]*step2;
00118       tn+=step2/v;
00119     }
00120     else{
00121       //
00122       // just step 
00123       //
00124       for(int i=0;i<3;i++)rn[i]+=t[i]*step;
00125       tn+=step/v;
00126     }
00127     //
00128     // Have to return if volume changes
00129     //
00130     r=sqrt(rn[0]*rn[0]+rn[1]*rn[1]+rn[2]*rn[2]);
00131     if(rstart<R0&&r>R0)return;
00132     if(rstart>R0&&r<R0)return;
00133   }
00134   //
00135   // figure out what kind of interaction
00136   // only consider elastic(odd) or absorb.(even)
00137   //
00138   bool stopped=0;
00139   if(rv[2]<probproc[0]){A=a[0];stopped=0;}
00140   else if(rv[2]<probproc[1]){A=a[0];stopped=1;}
00141   else if(rv[2]<probproc[2]){A=a[1];stopped=0;}
00142   else if(rv[2]<probproc[3]){A=a[1];stopped=1;}
00143   else if(rv[2]<probproc[4]){A=a[2];stopped=0;}
00144   else if(rv[2]<probproc[5]){A=a[2];stopped=1;}
00145   else if(rv[2]<probproc[6]){A=a[3];stopped=0;}
00146   else if(rv[2]<probproc[7]){A=a[3];stopped=1;}
00147   //
00148   // if stopped finished
00149   //
00150   if(stopped){
00151     ke=0;
00152     return;
00153   }
00154   //
00155   // otherwise generate an elastic scatter
00156   // throw target momentum via Max-Boltz. 
00157   //
00158   double Sqrt2mkT = sqrt(2*A*KBOLTZMAN*Temperature*MPROTON);
00159   float rw[3];
00160   len=3;
00161   rnormx_(rw,len,ranlux_);
00162   double pt[3];
00163   for(int i=0;i<3;i++)pt[i]=Sqrt2mkT*rw[i];
00164   //
00165   double vr[3];
00166   for(int i=0;i<3;i++)vr[i]=pn[i]-pt[i]/A;
00167   double vvr=sqrt(vr[0]*vr[0]+vr[1]*vr[1]+vr[2]*vr[2]);
00168   double u[3];
00169   double cosu = -1+2*rv[3];
00170   double phiu = 2*PI*rv[4];
00171   u[0] = sqrt(1-cosu*cosu)*cos(phiu);
00172   u[1] = sqrt(1-cosu*cosu)*sin(phiu);
00173   u[2] = cosu;
00174   //
00175   for(int i=0;i<3;i++)pn[i] = pn[i]/(A+1)+pt[i]/(A+1)+vvr*A/(A+1)*u[i];
00176   ke=0;
00177   for(int i=0;i<3;i++)ke+=pn[i]*pn[i]/2/mn;
00178   //
00179   return;
00180 }
00181 

Generated on Mon Feb 21 16:11:19 2005 for ReactorFsim by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002