---

Photin.h

Go to the documentation of this file.

// -*- C++ -*-  $Id: Photin.h,v 1.1 2000/09/01 20:20:17 burnett Exp $
//
// This file is part of Gismo 2

#ifndef __PHOTIN_H
#define __PHOTIN_H

#include "facilities/error.h"

#include <iostream>
#include <strstream>

#include <cmath>
class Photin  // photon cross sections in this material
{
   friend class PEGSData;
public:

    ~Photin();
    // dtor
        
   float rsct(float xxx)const;
   // Rayleigh scattering angle function.

private:
   float ebinda;   // K edge energy (used in routine photo)

    float volatile eold;   // last energy (in MeV)
    float volatile loge;   // log of E
    int   volatile index;  // index

    float ge0, ge1;  // interpolating constants

    int   mge;  // length of following arrays
    float *gmfp0, *gmfp1,
         *gbr10, *gbr11,
         *gbr20, *gbr21,
         *cohe0, *cohe1;

    Photin():mge(0){eold = -1.;} // constructor

    // calculate new index if energy changes
    void updateIndex(float enew) const
    {  if (enew == eold) return;
#if 1
        Photin * self = const_cast<Photin*>(this);
        self->loge = log(enew);  // convert to MeV
        self->index= (int)(ge1*loge + ge0 -1); // note -1 for C arrays
        if (index>=mge) {
            std::ostrstream msg;
            msg << "Photin: photon energy, " << enew
                << "MeV, index "<< index << ", exceeds ap, mge= " << mge << '\0';
            FATAL(msg.str());
        }
        self->eold = enew;
#else // use volatile
      loge = log(enew);  // convert to MeV
      index= (int)(ge1*loge + ge0 -1); // note -1 for C arrays
      if (index>=mge)
      {
          FATAL("Photin: photon energy exceeds parameter ap");
      }
      eold = enew;
#endif
   }
   // the interpolating functions (should have designed a macro)

   float  gmfp()const{return  gmfp1[index]*loge + gmfp0[index];}
   float  gbr1()const{return  gbr11[index]*loge + gbr10[index];}
   float  gbr2()const{return  gbr21[index]*loge + gbr20[index];}
   float  cohe()const{return  cohe1[index]*loge + cohe0[index];}

   float rco0, rco1;

   int   ngr;   // length of following table
   float *rsct0, *rsct1;


   void readPEGS(std::istream&, int, int);
   void rescale(float);

};

#endif

---