---

RMPflux.cpp

#include <math.h>
#include <iostream.h>
#include "ReactorFortran.hh"

double RMPflux(double enu){
  //
  // from Vogel, Gratta RMP via Jon Link
  //
  int dim = 22;
  float etab[22];
  for(int n=0;n<dim;n++)etab[n]=1.85+n*0.3;
  float ftab[22]= {49.80,39.66,31.43,24.79,19.45,
                   15.19,11.80, 9.12, 7.01, 5.35,
                   4.06,  3.06, 2.28, 1.67, 1.21,
                   0.85,  0.58, 0.38, 0.23, 0.13,
                   0.06,0.03};
                 
  bool static first=1;
  float static ptab[22];
  
  if(first){
    double sum;
    for(int n=0;n<dim;n++)sum+=ftab[n];
    for(int n=0;n<dim;n++)ptab[n]=ftab[n]/sum;
    first = 0;
    cout << "RMPflux initialization complete" << endl;
  }
  
  int kpoly = 2;
  float e=enu;
  //
  // a CERNLIB routine to interpolate the table
  //
  return  divdif_(ptab,etab,dim,e,kpoly);

}
  

---