---

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"

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];

    if(first){
      float xlo = 0.;
      float xhi = 50.;
      funlxp_(Cf252NeutronSpectrum,fspace,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;
    }

    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 energy
      int len = 1;
      funlux_(fspace,neutronE[n],len);
      //cout << "   Cf Neutron Energy = " << neutronE[n] << endl;

    } // done looping over sources

    // zero out remaining energy and position arrays
    for(int n=Nsources;n<k.Ndim;n++){
      neutronE[n] = 0.;
      R[n] = 0.;
      zR[n] = 0.;
      phiR[n] = 0.;
    }

  } // done with Cf 252

  if(Isotope == 255){

    Nsources = 0;
    for(int n=Nsources;n<k.Ndim;n++){
      
      neutronE[n] = 0.;
      R[n] = 0.;
      zR[n] = 0.;
      phiR[n] = 0.;
    }
  } // done with Cf 255

}

MyCfSource::~MyCfSource(){
  ;
}

MyCfSource::MyCfSource(const MyCfSource& CfSource){

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

  for(int n=0;n<Sdim;n++){
    neutronEnergy[n] = CfSource.neutronEnergy[n];
    neutronProb[n] = CfSource.neutronProb[n];
  }

  Nsources = CfSource.Nsources;

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

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

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

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

  for(int n=0;n<Sdim;n++){
    neutronEnergy[n] = rhs.neutronEnergy[n];
    neutronProb[n] = rhs.neutronProb[n];
  }

  Nsources = rhs.Nsources;

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

  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;
}

---