---

ReactorNeutrons.cpp

#include <math.h>
#include <iostream.h>
#include <fstream.h>
#include <string.h>
#include "ReactorNeutrons.hh"
#include "ReactorConstants.hh"

ReactorNeutrons::ReactorNeutrons(){
  Data = new XCdata[0];
}

ReactorNeutrons::ReactorNeutrons(int option){
  Data = new XCdata[DataSize];
  DataPtr = Data;
  AddHdata();
  AddCdata();
  AddGd155data();
  AddGd157data();
}
ReactorNeutrons::~ReactorNeutrons(){
  delete[] Data;
}
ReactorNeutrons::ReactorNeutrons(const ReactorNeutrons& N){

  Data = new XCdata[DataSize];
  XCdata* pold = N.Data+DataSize;
  XCdata* pnew = Data+DataSize;
  while (pnew>Data)*--pnew=*--pold;
  
  DataPtr=Data;
}
ReactorNeutrons& ReactorNeutrons::operator=(const ReactorNeutrons& rhs){
  if(this==&rhs)return *this;

  delete[] Data;
  Data = new XCdata[DataSize];
  XCdata* pold = rhs.Data+DataSize;
  XCdata* pnew = Data+DataSize;
  while (pnew>Data)*--pnew=*--pold;

  DataPtr=Data;

  return *this;

}
//
//
//
ReactorNeutrons::XCdata::XCdata(){
  Size = 0;
  KE = new double[0];
  XC = new double[0];
  LastKE=KE;
  LastXC=XC;
}
ReactorNeutrons::XCdata::XCdata(char* filename){

  int n;
  double ke[MaxSize];
  double xc[MaxSize];

  strcpy(FileName,filename);

  GetData(filename,n,ke,xc);

  Size = n;
  KE = new double[n];
  XC = new double[n];

  for(int i=0;i<n;i++){
    *(KE+i)=ke[i];
    *(XC+i)=xc[i];
  }
  LastKE = KE+n-1;
  LastXC = XC+n-1;

  strncpy(Process,FileName+11,5);
  Process[5] = '\0';

  Z = 100*(FileName[5]-'0')+10*(FileName[6]-'0')+(FileName[7]-'0');
  A = 100*(FileName[8]-'0')+10*(FileName[9]-'0')+(FileName[10]-'0');

  /*
  cout << " Cross Section data for " << Process 
       <<" from file " << FileName 
       << " for A= " << A
       << " and Z= " << Z << endl;
  */
}
ReactorNeutrons::XCdata::~XCdata(){
  delete[] KE;
  delete[] XC;
}
ReactorNeutrons::XCdata::XCdata(const ReactorNeutrons::XCdata& X){
  
  A=X.A;
  Z=X.Z;
  //
  strcpy(FileName,X.FileName);
  strcpy(Process,X.Process);
  //
  Size = X.Size;
  //
  KE = new double[Size];
  XC = new double[Size];
  //
  LastKE = KE+Size-1;
  LastXC = XC+Size-1;
  //
  double* pold = X.KE+Size;
  double* pnew = KE+Size;
  while(pnew>KE)*--pnew=*--pold;
  //
  pold = X.XC+Size;
  pnew = XC+Size;
  while(pnew>XC)*--pnew=*--pold;
}
//
//
//
ReactorNeutrons::XCdata& 
ReactorNeutrons::XCdata::operator=(const ReactorNeutrons::XCdata& rhs){
  
  if(this==&rhs)return *this;
  //
  A=rhs.A;
  Z=rhs.Z;
  //
  strcpy(FileName,rhs.FileName);
  strcpy(Process,rhs.Process);
  //
  Size = rhs.Size;
  delete[] KE;
  delete[] XC;
  //
  KE = new double[Size];
  XC = new double[Size];
  //
  LastKE = KE+Size-1;
  LastXC = XC+Size-1;
  //
  double* pold = rhs.KE+Size;
  double* pnew = KE+Size;
  while(pnew>KE)*--pnew=*--pold;
  //
  pold = rhs.XC+Size;
  pnew = XC+Size;
  while(pnew>XC)*--pnew=*--pold;
  //
  return *this;
}
//
//
//
void ReactorNeutrons::AddData(char* file1,char* file2,char* file3){
  ReactorNeutrons::XCdata d1(file1);
  *DataPtr++=d1;
  ReactorNeutrons::XCdata d2(file2);
  *DataPtr++=d2;
  ReactorNeutrons::XCdata d3(file3);
  *DataPtr++=d3;
}

void ReactorNeutrons::XCdata::GetData(const char* file,int& size,double* e,double* xc){
  size=0;
  double* eptr=e;
  double* xcptr=xc;
  ifstream infile(file);
  while(infile.good()){
    infile>> *eptr++ >> *xcptr++;size++;
  }  
  size-=1;
  infile.close();
}
//
//
//
int ReactorNeutrons::GetTotalSize(int Z,int A){
  for(ReactorNeutrons::XCdata* data = Data;data<Data+DataSize;data++){
    if(data->GetA()==A&&data->GetZ()==Z&&data->Total())return data->GetSize();
  }
}    
int ReactorNeutrons::GetElastSize(int Z,int A){
  for(ReactorNeutrons::XCdata* data = Data;data<Data+DataSize;data++){
    if(data->GetA()==A&&data->GetZ()==Z&&data->Elastic())return data->GetSize();
  }
}    
int ReactorNeutrons::GetGammaSize(int Z,int A){
  for(ReactorNeutrons::XCdata* data = Data;data<Data+DataSize;data++){
    if(data->GetA()==A&&data->GetZ()==Z&&data->Gamma())return data->GetSize();
  }
}    
//
//
//
double* ReactorNeutrons::GetTotalKE(int Z,int A){
  for(ReactorNeutrons::XCdata* data = Data;data<Data+DataSize;data++){
    if(data->GetA()==A&&data->GetZ()==Z&&data->Total())return data->GetKE();
  }
}    
double* ReactorNeutrons::GetElastKE(int Z,int A){
  for(ReactorNeutrons::XCdata* data = Data;data<Data+DataSize;data++){
    if(data->GetA()==A&&data->GetZ()==Z&&data->Elastic())return data->GetKE();
  }
}    
double* ReactorNeutrons::GetGammaKE(int Z,int A){
  for(ReactorNeutrons::XCdata* data = Data;data<Data+DataSize;data++){
    if(data->GetA()==A&&data->GetZ()==Z&&data->Gamma())return data->GetKE();
  }
}    
//
double* ReactorNeutrons::GetTotalXC(int Z,int A){
  for(ReactorNeutrons::XCdata* data = Data;data<Data+DataSize;data++){
    if(data->GetA()==A&&data->GetZ()==Z&&data->Total())return data->GetXC();
  }
}    
double* ReactorNeutrons::GetElastXC(int Z,int A){
  for(ReactorNeutrons::XCdata* data = Data;data<Data+DataSize;data++){
    if(data->GetA()==A&&data->GetZ()==Z&&data->Elastic())return data->GetXC();
  }
}    
double* ReactorNeutrons::GetGammaXC(int Z,int A){
  for(ReactorNeutrons::XCdata* data = Data;data<Data+DataSize;data++){
    if(data->GetA()==A&&data->GetZ()==Z&&data->Gamma())return data->GetXC();
  }
}    
//
//
//
double* ReactorNeutrons::GetLastKE(int Z,int A,char option){
  for(ReactorNeutrons::XCdata* data = Data;data<Data+DataSize;data++){
    if(data->GetA()==A&&data->GetZ()==Z&&option=='T'&&data->Total()){
      return data->GetLastKE();
    }
    if(data->GetA()==A&&data->GetZ()==Z&&option=='E'&&data->Elastic()){
      return data->GetLastKE();
    }
    if(data->GetA()==A&&data->GetZ()==Z&&option=='G'&&data->Gamma()){
      return data->GetLastKE();
    }
  } 
}
double* ReactorNeutrons::GetLastXC(int Z,int A,char option){
  for(ReactorNeutrons::XCdata* data = Data;data<Data+DataSize;data++){
    if(data->GetA()==A&&data->GetZ()==Z&&option=='T'&&data->Total()){
      return data->GetLastXC();
    }
    if(data->GetA()==A&&data->GetZ()==Z&&option=='E'&&data->Elastic()){
      return data->GetLastXC();
    }
    if(data->GetA()==A&&data->GetZ()==Z&&option=='G'&&data->Gamma()){
      return data->GetLastXC();
    }
  } 
}
//
void ReactorNeutrons::SetLastData(int Z,int A,char option,double* lastKE,double* lastXC){
  for(ReactorNeutrons::XCdata* data = Data;data<Data+DataSize;data++){
    if(data->GetA()==A&&data->GetZ()==Z&&option=='T'&&data->Total()){
      data->SetLastKE(lastKE);
      data->SetLastXC(lastXC);
      return;
    }
    if(data->GetA()==A&&data->GetZ()==Z&&option=='E'&&data->Elastic()){
      data->SetLastKE(lastKE);
      data->SetLastXC(lastXC);
      return;
    }
    if(data->GetA()==A&&data->GetZ()==Z&&option=='G'&&data->Gamma()){
      data->SetLastKE(lastKE);
      data->SetLastXC(lastXC);
      return;
    }
  } 
}
double ReactorNeutrons::GetXC(int Z,int A,double ke,char option){
  double xc=1e-6;
  //
  double* KE;
  double* XC;
  if(option=='T'){
    KE=GetTotalKE(Z,A);
    XC=GetTotalXC(Z,A);
  }
  if(option=='E'){
    KE=GetElastKE(Z,A);
    XC=GetElastXC(Z,A);
  }
  if(option=='G'){
    KE=GetGammaKE(Z,A);
    XC=GetGammaXC(Z,A);
  }
  //
  double* LastKE=GetLastKE(Z,A,option);
  double* LastXC=GetLastXC(Z,A,option);
  //
  if(ke<*KE){
    return *XC;
  }
  else if(ke<*LastKE){
    double* xcptr=LastXC;
    for(double* keptr=LastKE;keptr>KE;keptr--){
      if(ke<*keptr && ke>*(keptr-1)){
        xc=*(xcptr-1)+(*xcptr-*(xcptr-1))/(*keptr-*(keptr-1))*(ke-*(keptr-1));
        SetLastData(Z,A,option,keptr,xcptr);
        return xc;
      }
      xcptr--;
    }
  }
  else{
    int size=0;
    if(option=='T')size=GetTotalSize(Z,A);
    if(option=='E')size=GetElastSize(Z,A);
    if(option=='G')size=GetGammaSize(Z,A);
    if(size == 0) cout << " ERROR: ReactorNeutron size is zero!" << endl;
    //
    double *xcptr=LastXC;
    for(double* keptr=LastKE;keptr<KE+size-1;keptr++){
      if(ke<*(keptr+1) && ke>*keptr){
        xc=*xcptr+(*(xcptr+1)-*xcptr)/(*(keptr+1)-*keptr)*(ke-*keptr);
        SetLastData(Z,A,option,keptr,xcptr);
        return xc;
      }
      ++xcptr;
    }
  }
  return xc;
}

---