---

MyCfSource.cpp

Go to the documentation of this file.

#include <math.h>
#include <iostream.h>
#include <fstream.h>   // file I/O
#include <iomanip.h>   // format manipulation
#include <string>
#include "ReactorConstants.hh"
#include "ReactorFortran.hh"
#include "MyCfSource.hh"

long long factorial(int n);

MyCfSource::MyCfSource(int newIsotope,
                       int newRmaxgen){

  Isotope = newIsotope;
  Rmaxgen = newRmaxgen;

  //  cout << "MyCfSource::MyCfSource(" << Isotope << "," << Rmaxgen 
  //       << ")" << endl;

  // Cf252
  if(Isotope == 252){

    // setup the probability density as a function of energy
    bool static first=1;
    float static fspace[200];
    float static mspace[200];
    float static gspace[200];

    if(first){
      float xlo = 0.;
      float xhi = 50.;
      funlxp_(Cf252NeutronSpectrum,fspace,xlo,xhi);
      xhi = 25.;
      funlxp_(Cf252GammaMultiplicityFit,mspace,xlo,xhi);
      funlxp_(Cf252GammaSpectrum,gspace,xlo,xhi);
      first = 0;
    } 
 
    // input the source number and position
    fstream cf252_position("data/cf252_position.txt", ios::in);
    //cout << "Cf input file: data/cf252_position.txt" << endl;

    Nsources = 0;
    char paramLine[256];
    string dummy;
    bool done = false;
    bool numSet = false;
    int positionNum = 0;

    while(!done){

      // '#' specifies a comment line -- ignore it
      if(cf252_position.peek() == '#'){

        cf252_position.getline(paramLine, 256);
        continue;
      }

      // 'E'ND or 'e'nd specifies end of file -- set done flag to true
      if(cf252_position.peek() == 'e' || cf252_position.peek() == 'E'){

        if(Nsources == 0) cout << "   Cf252 ERROR: " 
                               << "number of cf252 sources = 0" << endl;
        done = true;
        continue;
      }

      // 'N'SOURCES = or 'n'sources = specifies number of sources 
      if(cf252_position.peek() == 'n' || cf252_position.peek() == 'N'){

        if(numSet){

          cf252_position.getline( paramLine, 256 );
          continue;
        }
        cf252_position >> dummy;
        while(cf252_position.peek() == ' ' || cf252_position.peek() == '=') 
          cf252_position.ignore(1);

        cf252_position >> Nsources;
        if(Nsources == 0){
          cout << "   Cf252 ERROR: number of cf252 sources = 0" << endl;
          done = true;
          continue;
        }

        numSet = true;
        continue;
      }

      // 'P'OSITION = or 'p'osition = specifies position of the sourse(s)
      if(cf252_position.peek() == 'p' || cf252_position.peek() == 'P'){

        if(Nsources == 0 || positionNum >= Nsources){

          cout << "   Cf252 ERROR: too many positions " << positionNum+1 
               << " given for the number of sources " << Nsources << endl;
          cf252_position.getline( paramLine, 256 );
          positionNum++;
          continue;
        }
        cf252_position >> dummy;
        while(cf252_position.peek() == ' ' || cf252_position.peek() == '=') 
          cf252_position.ignore(1);

        cf252_position >> R[positionNum] >> zR[positionNum] 
                       >> phiR[positionNum];

        // check for unphysical source positions
        if(R[positionNum] > Rmaxgen){
          cout << "   Cf252 ERROR: R[" << positionNum << "] of " 
               << R[positionNum] << " > Rmaxgen" << endl;
          done = true;
          continue;
        }
        if(zR[positionNum] < -1. || zR[positionNum] > 1.){
          cout << "   Cf252 ERROR: CosTheta[" 
               << positionNum << "] of " << zR[positionNum] << endl;
          done = true;
          continue;
        }
        if(phiR[positionNum] > 2*k.pi){
          cout << "   Cf252 ERROR: Phi[" << positionNum << "] of " 
               << phiR[positionNum] << " > 2*pi" << endl;
          done = true;
          continue;
        }

        positionNum++;
        continue;
      }

      // ignore extra whitespace
      if(cf252_position.peek() == ' ' || cf252_position.peek() == '\n'){

        cf252_position.ignore(1);
        continue;
      }
    }

    if(positionNum < Nsources){
      cout << "   Cf252 ERROR: not enough positions " << positionNum 
           << " given for the number of sources " << Nsources << endl;
    }
    //
    // neutron multiplicity distribution (via Los Alamo manual)
    //
    double p[maxNeutron+1] = {0.002,0.028,0.155,0.428,
                              0.732,0.917,0.983,0.998,1.000};
    //
    for(int n=0;n<Nsources;n++){
      //cout << "Source " << n+1 << endl << "   Cf Source Position = " 
      //           << R[n] << ", " << zR[n] << ", " << phiR[n] << endl;

      // pick a neutron multiplicity
      bool passed = false;
      int len;
      while(!passed){
        len = 1;
        float r[len];
        ranlux_(r,len);

        for(int i=0;i<maxNeutron+1; i++){

          if(p[i] > r[0] && !passed){
            Nneutron[n] = i;
            if(Nneutron[n])passed = true;
          }
        }
      }
      //cout << "   " << Nneutron[n] << " neutrons" << endl;
      //
      // pick an energy for each neutron
      //
      for(int nn=0;nn<Nneutron[n]; nn++){
        len = 1;
        funlux_(fspace,neutronE[n][nn],len);
        //cout << "   Cf neutron energy = " << neutronE[n][nn] << endl;
        Tneutron[n][nn] = 0.;
      }
      /*
        The total energy in the prompt gammas is a function of the
        material and the number of prompt neutrons produced; these 
        formulae come from T. Valentine's summary
      */
      float Z = 98.;
      float A = 252.;
      //
      double avge = -1.33 + 119.6*pow(Z,0.3333)/A;
      double phi = 2.51 - 1.13e-5*Z*Z*sqrt(A);
      double etot = phi*Nneutron[n]+4.0;
      //
      // use the Brunson multiplicity (also via Valentine)
      //
      len = 1;
      float ng = 0.;
      funlux_(mspace,ng,len);
      Ngamma[n] = (int)ng;
      //cout << "   " << Ngamma[n] << " photons" << endl;
      //
      // pick an energy for each gamma
      //
      double tote = 0.;
      for(int nn=0;nn<Ngamma[n]; nn++){
        len = 1;
        funlux_(gspace,gammaE[n][nn],len);
        //cout << "   252Cf photon energy = " << gammaE[n][nn] << endl;
        tote += gammaE[n][nn];

        len = 2;
        float rv[len];
        ranlux_(rv,len);
        //
        // 80% of gammas have T_1/2 = 0.01 ns and 20% have T_1/2 = 1 ns
        //
        double halflife;
        if(rv[0] < 0.8)
          halflife = 0.01;
        else
          halflife = 1.0;
        Tgamma[n][nn] = halflife*log(1/rv[1]);
      }
      //cout << "          total energy = " << tote << endl;
    } // done looping over sources

    // zero out remaining energy and position arrays
    for(int n=0;n<Nsources;n++){
      for(int nn=Nneutron[n];nn<maxNeutron;nn++){
        neutronE[n][nn] = 0.;
        Tneutron[n][nn] = 0.;
      }
      for(int nn=Ngamma[n];nn<maxGamma;nn++){
        gammaE[n][nn] = 0.;
        Tgamma[n][nn] = 0.;
      }
    }
    for(int n=Nsources;n<k.Ndim;n++){
      Nneutron[n] = 0;
      Ngamma[n] = 0;
      R[n] = 0.;
      zR[n] = 0.;
      phiR[n] = 0.;
      for(int nn=0;nn<maxNeutron;nn++){
        neutronE[n][nn] = 0.;
        Tneutron[n][nn] = 0.;
      }
      for(int nn=0;nn<maxGamma;nn++){
        gammaE[n][nn] = 0.;
        Tgamma[n][nn] = 0.;
      }
    }

  } // done with Cf 252

  if(Isotope == 255){

    Nsources = 0;
    for(int n=Nsources;n<k.Ndim;n++){
      Nneutron[n] = 0;
      Ngamma[n] = 0;
      R[n] = 0.;
      zR[n] = 0.;
      phiR[n] = 0.;
      for(int nn=0;nn<maxNeutron;nn++){
        neutronE[n][nn] = 0.;
        Tneutron[n][nn] = 0.;
      }
      for(int nn=0;nn<maxGamma;nn++){
        gammaE[n][nn] = 0.;
        Tgamma[n][nn] = 0.;
      }
    }
  } // done with Cf 255
}

MyCfSource::~MyCfSource(){
  ;
}

MyCfSource::MyCfSource(const MyCfSource& CfSource){

  Isotope = CfSource.Isotope;
  Rmaxgen = CfSource.Rmaxgen;

  Nsources = CfSource.Nsources;

  for(int n=0;n<k.Ndim;n++){
    Nneutron[n] = CfSource.Nneutron[n];
    Ngamma[n] = CfSource.Ngamma[n];
    R[n] = CfSource.R[n];
    zR[n] = CfSource.zR[n];
    phiR[n] = CfSource.phiR[n];

    for(int nn=0;nn<maxNeutron;nn++){ 
      neutronE[n][nn] = CfSource.neutronE[n][nn];
      Tneutron[n][nn] = CfSource.Tneutron[n][nn];
    }
    for(int nn=0;nn<maxGamma;nn++){
      gammaE[n][nn] = CfSource.gammaE[n][nn];
      Tgamma[n][nn] = CfSource.Tgamma[n][nn];
    }
  }
}    

MyCfSource& MyCfSource::operator=(const MyCfSource& rhs){

  if(this == &rhs)return *this;

  Isotope = rhs.Isotope;  
  Rmaxgen = rhs.Rmaxgen;

  Nsources = rhs.Nsources;

  for(int n=0;n<k.Ndim;n++){
    Nneutron[n] = rhs.Nneutron[n];
    Ngamma[n] = rhs.Ngamma[n];
    R[n] = rhs.R[n];
    zR[n] = rhs.zR[n];
    phiR[n] = rhs.phiR[n];

    for(int nn=0;nn<maxNeutron;nn++){
      neutronE[n][nn] = rhs.neutronE[n][nn];
      Tneutron[n][nn] = rhs.Tneutron[n][nn];
    }
    for(int nn=0;nn<maxGamma;nn++){
      gammaE[n][nn] = rhs.gammaE[n][nn];
      Tgamma[n][nn] = rhs.Tgamma[n][nn];
    }
  }
  return *this;
}    

float MyCfSource::Cf252NeutronSpectrum(const float& x){

  // return the neutron spectrum N(x)
  float N = 0.;

  float scale = 1/(2*k.pi*0.359);
  //cout << "scale " << scale << endl;

  float fminus = exp(-(sqrt(x)-sqrt(0.359))*(sqrt(x)-sqrt(0.359))/1.175);
  float fplus  = exp(-(sqrt(x)+sqrt(0.359))*(sqrt(x)+sqrt(0.359))/1.175);

  N = scale*(fminus-fplus);

  //cout << "N " << N << endl;
  return N;
}

float MyCfSource::Cf252GammaMultiplicity(const int& x){

  // return the gamma multiplicity M(x)
  double M = 0.;

  double C1 = 0.675;
  double C2 = 6.78;
  double C3 = 9.92;

  M = C1*pow(C2,x)*exp(-C2)/factorial(x) + (1-C1)*pow(C3,x)*exp(-C3)/factorial(x);

  //cout << "M(" << x << ") " << M << endl;

  return M;
}

long long factorial(int n){

  // cannot return above n = 20 because of memory
  if(n < 1)
    return 1;
  else
    return n * factorial(n-1);
}

float MyCfSource::Cf252GammaMultiplicityFit(const float& x){

  // return the gamma multiplicity M(x)
  double M = 0.;

  float gaussian = 0.13345*exp(-0.5*((x-7.0322)/2.6301)*((x-7.0322)/2.6301));
  float exponential = 0.11255*(exp(-(sqrt(x)-sqrt(7.5987))
                                   *(sqrt(x)-sqrt(7.5987))/0.56213));

  if(x <= 9.0){
    M = gaussian;
  }
  else{
    M = exponential;
  }

  //cout << "M(" << x << ") " << M << endl;

  return M;
}

float MyCfSource::Cf252GammaSpectrum(const float& x){

  // return the gamma spectrum N(x)
  float N = 0.;

  float gaussian = 1.8367*exp(-0.5*((x-0.45934)/0.31290)*((x-0.45934)/0.31290));
  float exponential = exp(0.84774-0.89396*x);

  if(x <= 0.744){
    N = gaussian;
  }
  else{
    N = exponential;
  }
  
  //cout << "N " << N << endl;
  return N;
}

---