---

NeutronPropagator.cpp

#include <math.h>
#include <iostream.h>
#include <fstream.h>
#include <string.h>
#include "ReactorFortran.hh"
#include "ReactorConstants.hh"
#include "ReactorNeutrons.hh"
void NeutronPropagator$(double,double,double&,double&,double*,double*,double&);
void NeutronPropagator(double R0,double R2,
                        double& ke,double& A,int& passes,
                        double* pn,double* rn,double& tn){

  double r=sqrt(rn[0]*rn[0]+rn[1]*rn[1]+rn[2]*rn[2]);
  while(ke>0 && r<R2){
    passes++;
    NeutronPropagator$(R0,R2,ke,A,pn,rn,tn);
    r=sqrt(rn[0]*rn[0]+rn[1]*rn[1]+rn[2]*rn[2]);
  }
  return;
}
void NeutronPropagator$(double R0,double R2,
                        double& ke,double& A,
                        double* pn,double* rn,double& tn){
  //
  if(ke==0)return;
  //
  static ReactorConstants RC;
  //
  double r=sqrt(rn[0]*rn[0]+rn[1]*rn[1]+rn[2]*rn[2]);
  if(r>R2)return;
  double rstart=r;
  //
  // Set up cross section and materials data
  //
  static ReactorNeutrons XCdata;
  static bool first=1;
  static double a[4];
  static double z[4];
  static double nn[4];
  //
  // Initialize on first call
  // 
  if(first){
    //
    a[0] = RC.A0;
    a[1] = RC.A1;
    a[2] = RC.A2;
    a[3] = RC.A3;
    //
    z[0] = RC.Z0;
    z[1] = RC.Z1;
    z[2] = RC.Z2;
    z[3] = RC.Z3;
    //
    nn[0]=RC.Navogadro*RC.f0*RC.density/RC.A0;
    nn[1]=RC.Navogadro*RC.f1*RC.density/RC.A1;
    nn[2]=RC.Navogadro*RC.GdConcentration/100*RC.f2*RC.density/RC.A2;
    nn[3]=RC.Navogadro*RC.GdConcentration/100*RC.f3*RC.density/RC.A3;
    //
    ReactorNeutrons Init(0);
    XCdata=Init;
    first=0;
  }
  //
  // fetch cross section info
  //
  int nproc=8;
  double probproc[8];
  for(int i=0;i<nproc/2;i++){
    probproc[2*i]=nn[i]*XCdata.GetElastXC(z[i],a[i],1e6*ke);
    probproc[2*i+1]=nn[i]*XCdata.GetGammaXC(z[i],a[i],1e6*ke);
  }
  //
  // No Gd if outside inner volume
  //
  if(r>R0)for(int i=4;i<nproc;i++)probproc[i]=0;
  //
  // Normalize
  //
  double probtot=0;
  for(int i=0;i<nproc;i++)probtot+=probproc[i];
  for(int i=0;i<nproc;i++)probproc[i]/=probtot;
  for(int i=1;i<nproc;i++)probproc[i]+=probproc[i-1];
  //
  // Take steps until there is an interaction
  //
  double step=RC.SlowNeutronStep;
  double mn = RC.Mproton+RC.Delta;
  //
  double mfp=1/probtot/1e-24;
  double probint = 1-exp(-step/mfp);
  int len=5;
  float rv[5];
  double t[3];
  double p=sqrt(pn[0]*pn[0]+pn[1]*pn[1]+pn[2]*pn[2]);
  double v=p/mn*RC.c;
  for(int i=0;i<3;i++)t[i]=pn[i]/p;
  bool interacted=0;
  while(!interacted){
    ranlux_(rv,len);
    interacted=(rv[0]<probint);
    if(interacted){
      //
      // throw the exponential within the step 
      //
      double step2=-mfp*log(1-probint*rv[1]);
      for(int i=0;i<3;i++)rn[i]+=t[i]*step2;
      tn+=step2/v;
    }
    else{
      //
      // just step 
      //
      for(int i=0;i<3;i++)rn[i]+=t[i]*step;
      tn+=step/v;
    }
    //
    // Have to return if volume changes
    //
    r=sqrt(rn[0]*rn[0]+rn[1]*rn[1]+rn[2]*rn[2]);
    if(rstart<R0&&r>R0)return;
    if(rstart>R0&&r<R0)return;
  }
  //
  // figure out what kind of interaction
  // only consider elastic(odd) or absorb.(even)
  //
  bool stopped=0;
  if(rv[2]<probproc[0]){A=a[0];stopped=0;}
  else if(rv[2]<probproc[1]){A=a[0];stopped=1;}
  else if(rv[2]<probproc[2]){A=a[1];stopped=0;}
  else if(rv[2]<probproc[3]){A=a[1];stopped=1;}
  else if(rv[2]<probproc[4]){A=a[2];stopped=0;}
  else if(rv[2]<probproc[5]){A=a[2];stopped=1;}
  else if(rv[2]<probproc[6]){A=a[3];stopped=0;}
  else if(rv[2]<probproc[7]){A=a[3];stopped=1;}
  //
  // if stopped finished
  //
  if(stopped){
    ke=0;
    return;
  }
  //
  // otherwise generate an elastic scatter
  // throw target momentum via Max-Boltz. 
  //
  double Sqrt2mkT = sqrt(2*A*RC.Kboltzman*RC.Temperature*RC.Mproton);
  float rw[3];
  len=3;
  rnormx_(rw,len,ranlux_);
  double pt[3];
  for(int i=0;i<3;i++)pt[i]=Sqrt2mkT*rw[i];
  //
  double vr[3];
  for(int i=0;i<3;i++)vr[i]=pn[i]-pt[i]/A;
  double vvr=sqrt(vr[0]*vr[0]+vr[1]*vr[1]+vr[2]*vr[2]);
  double u[3];
  double cosu = -1+2*rv[3];
  double phiu = 2*RC.pi*rv[4];
  u[0] = sqrt(1-cosu*cosu)*cos(phiu);
  u[1] = sqrt(1-cosu*cosu)*sin(phiu);
  u[2] = cosu;
  //
  for(int i=0;i<3;i++)pn[i] = pn[i]/(A+1)+pt[i]/(A+1)+vvr*A/(A+1)*u[i];
  ke=0;
  for(int i=0;i<3;i++)ke+=pn[i]*pn[i]/2/mn;
  //
  return;
}

---