/RTFunctions.cc

Go to the documentation of this file.

#include "RTFunctions.h"

#include <utl/config.h>

#include <utl/Photon.h>
#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/RandomEngine.h>
#include <utl/MathConstants.h>
#include <utl/AugerException.h>
#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>

#include <CLHEP/Random/Randomize.h>

#include <iostream>
#include <vector>
#include <utility>


using namespace utl;
using namespace std;

using CLHEP::RandFlat;
using CLHEP::RandGauss;

namespace TelescopeSimulatorKG2 {

  namespace RTFunctions {

    void
    Reflection(const utl::Photon& photonIn,
               const Vector& normal,
               utl::Photon& photonOut) 
    {
      //This asumes that normal is normalized check?
      photonOut = photonIn;
      const Vector& nIn = photonIn.GetDirection();
      const double dot = nIn * normal;   // projection
      photonOut.SetDirection(nIn - 2. * dot * normal); // -> new direction
    }

    //a phenomenological description of the absorption by the dust layer on the outside of the filter
    //the weight of the photon is reduced by the absorption factor set in TelescopeSimulatorKG2.xml.in
    //for 10% absorption, set FilterDustAbsorptionOutside = 0.1
    //if not specified, the default value is 0.0
    void
    Absorption(const utl::Photon& photonIn,
               const double filterAbsorptionFactor,
               utl::Photon& photonOut) 
    {
      const double weightIn = photonIn.GetWeight();
      const double weightOut = weightIn * ( 1.0 - filterAbsorptionFactor);
      photonOut = photonIn;
      photonOut.SetWeight(weightOut);
    }

    //a phenomenological description of the absorption by the dust layer on the mirror, default: no absorption
    //the weight of the photon is reduced by an elevation-dependent absorption factor set in TelescopeSimulatorKG2.xml.in
    //for e.g. 3% absorption at the top of the mirror, set MirrorAbsorptionTop = 0.03
    //for e.g. 9% absorption at the bottom of the mirror, set MirrorAbsorptionBot = 0.03
    void
    Absorption(const utl::Photon& photonIn,
               const double mirrorAbsorptionFactorTop,
               const double mirrorAbsorptionFactorBot,
               const double verticalPosOnMirror,
               utl::Photon& photonOut,
               const double mirrorSize)
    {
      //con't linear interpolation between absorption factors at the top and bottom of the mirror
      double mirrorInterpolWeight = (verticalPosOnMirror + mirrorSize) / (2 * mirrorSize); //1 for top of the mirror, 0 for bottom
      double absorptionFactor = (mirrorInterpolWeight * mirrorAbsorptionFactorTop) + ((1 - mirrorInterpolWeight) * mirrorAbsorptionFactorBot);
      const double weightIn = photonIn.GetWeight();
      const double weightOut = weightIn * (1.0 - absorptionFactor);
      photonOut = photonIn;
      photonOut.SetWeight(weightOut);
    }

    double 
    Plane(const utl::Point& point,
          const utl::Vector& normal,
          const utl::Photon& photonIn,
          utl::Photon& photonOut) 
    {      
      photonOut = photonIn;
      const Vector& nIn = photonIn.GetDirection();
      const Point& pIn = photonIn.GetPosition();

      double b = normal * nIn;
      double a = (point - pIn) * normal;

      if (b==0)
        return (0);

      a /= b;

      photonOut.SetPosition(pIn + a * nIn);
      return (a);
    }

    
    int
    Refraction(const double n12, 
               const utl::Photon& photonIn, 
               const Vector& normal, 
               PhotonList& photonsOut) 
    {

      const Vector& nIn = photonIn.GetDirection();


      // Snell's law ------------------------------------
      const double cosTheta_in = nIn * normal;    // projection on normal

      Vector projOnSurface = nIn - cosTheta_in * normal;  // projected direction on surface

      projOnSurface *= n12;                             // take refraction into account -> outgoing

      double sin2theta_out = projOnSurface.GetMag2();

      // the reflected photon
      Photon photonReflected(photonIn);
      photonReflected.SetDirection(nIn - 2. * cosTheta_in * normal); // -> new direction


      if (sin2theta_out >= 1) {                            // TOTAL INTERNAL REFLECTION
        photonsOut.push_back(photonReflected);
        return 1;
      }

      // check direction of incidence
      double cosTheta_out;
      if (cosTheta_in > 0)
        cosTheta_out = +sqrt(1. - sin2theta_out);
      else
        cosTheta_out = -sqrt(1. - sin2theta_out);  // ray coming out of media ! should not happen


      // Fresnel's law ----------------------------------

      // T = t^2   transmisivity
      // R = r^2   reflectivity
      // R_unpolarized = 1/2 * (R_perp + R_par)
      // T_unpolarized = 1/2 * (T_perp + T_par)
      //
      // because R and T gives the ratio of the energy flux per unit area (intensity),
      // one has to take into acount the difference between theta1 and theta2 !
      // 1. = R +  [n2/n1 * cos(theta2)/cos(theta1)] * T

      double fact = 1.;
      if (cosTheta_in)
        fact = 1. / n12 * cosTheta_out / cosTheta_in;

      const double T_par = pow(2. * n12 * cosTheta_in / (n12 * cosTheta_out + cosTheta_in), 2);
      const double T_perp = pow(2. * n12 * cosTheta_in / (n12 * cosTheta_in + cosTheta_out), 2);
      const double T =  fact * .5 * (T_par + T_perp); // unpolarized light

      const double weightIn = photonIn.GetWeight();

      Photon photonRefracted(photonIn);
      photonRefracted.SetDirection(projOnSurface + cosTheta_out*normal);
      photonRefracted.SetWeight(weightIn * T);
      photonsOut.push_back(photonRefracted);

      photonReflected.SetWeight(weightIn * (1 - T));
      photonsOut.push_back(photonReflected);

      return 2;
    }



    /* calculate the point and the normal vector on a sphere, where
       the photon photonIn hits it.

       -----
       origin:    center of sphere
       radius:    radius
       photonIn:  incoming photon
       photonOut: photon hitting the sphere
       normal:    surface normal on sphere at photonOut.GetPosition()
       *****
  */

    bool
    Sphere(const Point& origin,
           const double radius,
           const utl::Photon& photonIn,
           IntersectionList& intersection) 
    {
      intersection.clear();

      const Vector& nIn = photonIn.GetDirection();
      const Point& pIn = photonIn.GetPosition();

      const Vector pointToOrigin = pIn - origin;
      const double dist2 = pointToOrigin.GetMag2();

      const double projection = pointToOrigin * nIn; // project nIn on radial Vector pointToOrigin
      double delta = projection*projection + radius*radius - dist2;

      if (delta < 0)
        return false;                   // no solution

      if (delta == 0) {
        Photon photonOut(photonIn);
        photonOut.SetPosition(pIn - projection * nIn);
        Vector normal = origin - photonOut.GetPosition();
        normal.Normalize();
        intersection.push_back(pair<Photon,Vector>(photonOut, normal));
      }
      else {

        delta = sqrt(delta);
        
        const double move1 = -projection - delta;
        const double move2 = -projection + delta;
        
        Photon photonOut1(photonIn);
        photonOut1.SetPosition(pIn + min(move1, move2) * nIn);
        Vector normal1 = origin - photonOut1.GetPosition();
        normal1.Normalize();
        intersection.push_back(pair<Photon,Vector>(photonOut1, normal1));
        
        Photon photonOut2(photonIn);
        photonOut2.SetPosition(pIn + max(move1, move2) * nIn);
        Vector normal2 = origin - photonOut2.GetPosition();
        normal2.Normalize();
        intersection.push_back(pair<Photon,Vector>(photonOut2, normal2));
      }

      return true;
    }



    Vector
    TorusNormal(const utl::Point& origin,
                const utl::Vector& inwards,
                const utl::Point& point,
                const double torusRadius,
                const double tubeRadius)
    {
      const double R = torusRadius;
      const double r = tubeRadius;
      const Vector p = point - origin;
      // p = (x,y,z), inwards = (0,0,-1)
      // (df/dx,df/dy,df/dz) = 4 * p * (p*p + R*R - r*r) - 8 * R*R * (x,y,0)
      const Vector n = 4 * p * (p * p + R * R - r * r) - 8 * R * R * (p - inwards * (inwards * p));
      // normalize and point into the torus (out of the telescope)
      return -n / n.GetMag();
    }

    

    
    //bool 
    void
    Torus(const utl::Point& origin,
          const utl::Vector& inwards,
          const double torusRadius,
          const double tubeRadius,
          const utl::Photon& photonIn, 
          IntersectionList& intersections)
    {
      intersections.clear();


     /*  formulas calculated analog to pdf by Max Wagner mwagner@digipen.edu
      * but with torus equation
      * R*R - 2 * R * sqrt(x*x + y*y) + x*x + y*y + z*z - r*r = 0
      * eliminating the square root leads to:
      * f(x,y,z) = pow(R*R - r*r + x*x + y*y + z*z , 2) - 4 * R*R * (x*x + y*y) = 0
      * intersections with photon (x,y,z) = p + t * d
      * leads to the quartic equation
      * a4 * t*t*t*t + a3 * t*t*t + a2 * t*t + a1 * t + a0 = 0
      * solved below
      */
      
      const double R = torusRadius;
      const double r = tubeRadius;

      const Vector& d = photonIn.GetDirection();
      const Point& inPos = photonIn.GetPosition(); 
      const Vector p = inPos - origin;

      const double pd = p * d;
      const double dd = 1.; // = d*d, and d is Normalized
      const double pp = p * p;

      const double dz = d * inwards; // z-component of d
      const double dx2dy2 = 1. - dz * dz;
      const double pz = p * inwards; // z-component of p
      const double px2py2 = p.GetMag2() - pz * pz;    

      const double r2 = r * r;
      const double R2 = R * R;
      const double pprrRR = pp - r2 + R2;

      const double a4 = 1; // = dd * dd;
      const double a3 = 4 * dd * pd; // dd = 1
      const double a2 = 4 * pd * pd + 2 * dd * pprrRR - 4 * R2 * dx2dy2;
      const double a1 = 4 * pd * pprrRR - 8 * R2 * (pd - dz * pz);
      const double a0 = pprrRR * pprrRR - 4 * R2 * px2py2;
      
      vector<double> roots;
      Quartic(a0, a1, a2, a3, a4, roots);
      
      for (vector<double>::const_iterator iRoot = roots.begin();
           iRoot != roots.end();
           ++iRoot) {
        const Point pInt = inPos + (*iRoot) * d;
        Photon photonOut(photonIn);
        photonOut.SetPosition(pInt);
        intersections.push_back(pair<Photon,Vector>(photonOut, TorusNormal(origin, inwards, pInt, torusRadius, tubeRadius)));
      }

    } // ----------------
    
    
    void
    Cubic(const double a0,
          const double a1,
          const double a2,
          const double a3,
          vector<double>& roots)
    {
      if (a3 == 0) { // only check the extreme case
        // x^2 + a1/a2 * x + a0/a2 = 0
        const double Phalve = a1 / (2 * a2);
        const double Phalve2_Q = Phalve * Phalve - a0 / a2;
        if (Phalve2_Q < 0)
          return;
        const double sqrtPhalve2_Q = sqrt(Phalve2_Q);
        roots.push_back(-Phalve + sqrtPhalve2_Q);
        if (Phalve2_Q > 0)
          roots.push_back(-Phalve - sqrtPhalve2_Q);
        return;
      }

      //if (a2==0) {
      // quadratic eq after substiting  u = x^2
      //exit(1);
      //}
      

      // depressed cubic
      // v^3 + P v + Q  = 0
      
      
      const double P = (3 * a3 * a1 - a2 * a2) / (3 * a3 * a3);
      const double Q = (2 * a2 * a2 * a2 - 9 * a3 * a2 * a1 + 27 * a3 * a3 * a0) / (27 * a3 * a3 * a3);
      

      const double trans = - a2 / (3 * a3);
      
      if (fabs(P) < 1e-11) {
        
        roots.push_back(trans - cbrt(Q));
        
      } else if (fabs(Q) < 1e-11) {
        
        // one solution is v=0, take it and do nothing
        roots.push_back(trans);
        
      } else {

        const double PPP4QQ27 = 4 * P * P * P + 27 * Q * Q;
        if (fabs(PPP4QQ27) < 1e-11) { 

          roots.push_back(trans + 3 * Q / P);
          roots.push_back(trans - 3 * Q / (2 * P));
          
        } else {
                  
          // Cardano's method
          
          const double PPP27 = P * P * P / 27;
          const double QQ4PPP27 = Q * Q / 4 + PPP27;
        
          if ( QQ4PPP27 >= 0) {
            
            const double sqrtQQ4PPP27 = sqrt(QQ4PPP27);
            
            roots.push_back(trans + cbrt(-Q / 2 + sqrtQQ4PPP27) + cbrt(-Q / 2 - sqrtQQ4PPP27)); 

          } else {
            
            if (PPP27 >= 0) {
              cout << " THIS CANNOT HAPPEN " << endl;
              exit(1);
            }
            
            const double sqrt_PPP27 = sqrt( -PPP27 ); // length
            const double mag = cbrt(sqrt_PPP27);

            const double sqrt_QQ4PPP27 = sqrt(-QQ4PPP27);           
            const double phi = atan2(sqrt_QQ4PPP27, -Q / 2);

            // cout << " THIS CANNOT HAPPEND   PHI=" << phi << " sqrt_PPP27=" << sqrt_PPP27 << " sqrt_QQ4PPP27=" << sqrt_QQ4PPP27 << " -Q/2=" << -Q/2 <<  endl;
            
                    
            roots.push_back(trans + 2 * mag * cos(phi / 3));
            roots.push_back(trans + 2 * mag * cos(phi / 3 + kPi * 2 / 3));
            roots.push_back(trans + 2 * mag * cos(phi / 3 - kPi * 2 / 3));
            
          }       
        }       
      }
    }


    void
    Quartic(const double a0,
            const double a1,
            const double a2,
            const double a3,
            const double a4,
            vector<double>& roots)
    {      
      /* solve quartic function according to: http://en.wikipedia.org/wiki/Quartic_function 
         using
      
         a4 is A
         a3 is B
         a2 is C
         a1 is D
         a0 is E 

         convert to depressed quartic with x = u - a3/(4*a4)

         const double alpha = -3./8. * a3*a3 / (a4*a4) + a2 / a4;
         const double beta = a3*a3*a3 / (8 * a4*a4*a4) - a3*a2 / (2 * a4*a4) + a1/a4;
         const double gamma = -3 * a3*a3*a3*a3 / (256 * a4*a4*a4*a4) + a2*a3*a3 / (16 * a4*a4*a4) - a3*a1 / (4 * a4*a4) + a0 / a4;

         
         where:  x = u - a3 / (4*a4)
                 u^4 + alpha u^2 + beta u + gamma = 0

        [[ a4 = 1 leads to simplified coefficients ]]

      */
      
      if (a4 == 0) {
        Cubic(a0, a1, a2, a3, roots);
        return;
      }
      
      const double alpha = -3. / 8. * a3 * a3 / (a4 * a4) + a2 / a4;
      const double beta  = a3 * a3 * a3 / (8 * a4 * a4 * a4) - a3 * a2 / (2 * a4 * a4) + a1 / a4;
      const double gamma = -3 * a3 * a3 * a3 * a3 / (256 * a4 * a4 * a4 * a4) + a2 * a3 * a3 / (16 * a4 * a4 * a4) - a3 * a1 / (4 * a4 * a4) + a0 / a4;

      const double subst = - a3 / (4 * a4);      
      
      if (fabs(beta) < 1.e-11) { 
        
        // boils down to quadratic equation, after: t = u^2
        //
        //   t^2 + alpha t + gamma = 0
        //
        //  with x: = sqrt(t) - a3 / (4a4)
        // 
        
        const double tmp = alpha * alpha / 4 - gamma;
        
        if (tmp < 0) 
          return;
        
        const double sTmp = sqrt(tmp);
        const double tmp2 = -alpha / 2 + sTmp;
        const double tmp3 = -alpha / 2 - sTmp;
        
        if (tmp2 >= 0) {
          
          const double sTmp2 = sqrt(tmp2);
          const double x1 = subst + sTmp2;
          roots.push_back(x1);

          if (sTmp2 > 0) {
            const double x2 = subst - sTmp2;
            roots.push_back(x2);
          }       
        }
        
        if (tmp3 >= 0 && sTmp != 0) {

          const double sTmp3 = sqrt(tmp3);
          const double x3 = subst + sTmp3;
          roots.push_back(x3);

          if (sTmp3 > 0) {
            const double x4 = subst - sTmp3;
            roots.push_back(x4);
          }       
        }
        
        return;
        
      } // if beta==0   

      
      // solve the cubic equation in y
      const double b0 = alpha * alpha * alpha / 2. - alpha * gamma / 2. - beta * beta / 8.;
      const double b1 = 2 * alpha * alpha - gamma;
      const double b2 = alpha * 5./2;
      const double b3 = 1.;

      vector<double> rootsC;
      Cubic(b0, b1, b2, b3, rootsC);

      if (rootsC.size() == 0)
        return;
      
      const double y = rootsC[0]; // any solution is fine, take first one            
      const double alpha2y = alpha + 2 * y;
      
      if (alpha2y <- 1) {
        cout << "THIS SHOULD NOT HAPPEN!!!!!!!!!!!!" << endl;
        exit(5);
      }

      const double W = alpha2y < 0 ? 0 : sqrt(alpha2y);
      
      
      const double tmp1 = -(3 * alpha + 2 * y + 2 * beta / W);
      const double tmp2 = -(3 * alpha + 2 * y - 2 * beta / W);
      
      if (tmp1 < 0 && tmp2 < 0)
        return;
      

      if (tmp1 >= 0) {
        
        const double sTmp1 = sqrt(tmp1);

        const double x1 = subst + (W - sTmp1 ) / 2.;
        roots.push_back(x1);
        
        if (sTmp1 > 0) {
          const double x2 = subst + (W + sTmp1) / 2.;
          roots.push_back(x2);
        }       
      }

      if (tmp2 >= 0) {
        
        const double sTmp2 = sqrt(tmp2);

        const double x3 = subst + (-W - sTmp2) / 2.;
        roots.push_back(x3);
        
        if (sTmp2 > 0) {
          const double x4 = subst + (-W + sTmp2) / 2.;
          roots.push_back(x4);
        }       
      }
    }
    


    Vector
    RandomNormal(utl::RandomEngine& rndm, 
                 const Vector& normalIn, 
                 const double sigma_alpha)
    {
      if (sigma_alpha == 0) {
        return normalIn;    // nothing to do
      }

      // first create a reference system with normalIn as z normal
      const CoordinateSystemPtr internalCS = normalIn.GetCoordinateSystem();
      const double phiIn = normalIn.GetPhi (internalCS);
      const double thetaIn = normalIn.GetTheta (internalCS);
      const CoordinateSystemPtr tmpCS = CoordinateSystem::RotationZ (phiIn, internalCS);
      const CoordinateSystemPtr normalCS = CoordinateSystem::RotationY (thetaIn, tmpCS);

      const double phi =  RandFlat::shoot(&rndm.GetEngine(), 0., 2. * kPi);
      const double sinPhi = sin(phi);
      const double cosPhi = cos(phi);
      
      const double theta =  RandGauss::shoot(&rndm.GetEngine(), 0., sigma_alpha);
      const double sinTheta = sin(theta);
      const double cosTheta = cos(theta);

      const double unit_x = sinTheta * cosPhi;
      const double unit_y = sinTheta * sinPhi;
      const double unit_z = cosTheta;
      
      return Vector (unit_x, unit_y, unit_z, normalCS);
    } // end of RandomNormal

    /*

       Code equivalent to geant4 optical boundary process for unified model:

       return facet normal of rough surface, a ray rayIn will hit

       This function code alpha to a random value taken from the
       distribution p(alpha) = g(alpha; 0, sigma_alpha)*std::sin(alpha),
       for alpha > 0 and alpha < 90, where g(alpha; 0, sigma_alpha)
       is a gaussian distribution with mean 0 and standard deviation
       sigma_alpha.
    */

    Vector 
    RandomFacet(utl::RandomEngine& rndm,
                const Vector& normalIn, 
                const double sigma_alpha,
                const Vector& rayIn) {

      if (sigma_alpha == 0) {
        return normalIn;    // nothing to do
      }

      // first create a reference system with normalIn as z normal
      const CoordinateSystemPtr internalCS = normalIn.GetCoordinateSystem();
      const double phiIn = normalIn.GetPhi (internalCS);
      const double thetaIn = normalIn.GetTheta (internalCS);
      const CoordinateSystemPtr tmpCS = CoordinateSystem::RotationZ (phiIn, internalCS);
      const CoordinateSystemPtr normalCS = CoordinateSystem::RotationY (thetaIn, tmpCS);


      const double f_max = std::min(1., 4. * sigma_alpha);

      const double phi =  RandFlat::shoot (&rndm.GetEngine(), 0., 2. * kPi);
      const double sinPhi = sin(phi);
      const double cosPhi = cos(phi);

      Vector facetNormal;

      do {

        double alpha = 0;
        do {
          alpha =  RandGauss::shoot (&rndm.GetEngine(), 0., sigma_alpha);
        } while (RandFlat::shoot (&rndm.GetEngine(), 0., f_max) > sin(alpha) ||
                 alpha >= kPi / 2.);

        double sinAlpha = sin(alpha);
        double cosAlpha = cos(alpha);

        double unit_x = sinAlpha * cosPhi;
        double unit_y = sinAlpha * sinPhi;
        double unit_z = cosAlpha;

        facetNormal = Vector (unit_x, unit_y, unit_z, normalCS);

      } while (rayIn * facetNormal >= 0.0);

      return facetNormal;
    }//end of randomfacet


    /* put code for simulation of lambertian distribution here
     */

    Vector 
    LambertDiffusion(utl::RandomEngine& rndm,
                     const Vector& rayIn,
                     const Vector& normal,
                     const double amountOfStrayLight,
                     const double MaxTheta) {

      double strayOrNot = RandFlat::shoot(&rndm.GetEngine(), 0., 1.);

      if (strayOrNot > amountOfStrayLight) {
        return rayIn;   // nothing to do
      }
      else {

      // first create a reference system with normal as z normal
      const CoordinateSystemPtr internalCS = normal.GetCoordinateSystem();
      const double phiIn = normal.GetPhi (internalCS);
      const double thetaIn = normal.GetTheta (internalCS);
      const CoordinateSystemPtr tmpCS = CoordinateSystem::RotationZ (phiIn, internalCS);
      const CoordinateSystemPtr normalCS = CoordinateSystem::RotationY (thetaIn, tmpCS);

      // code taken from DrumPhotonGenerator.cc

      // values from DrumPhotonGenerator.xml
      const double MinTheta = 0. * degree;
      //const double MaxTheta = 25.*degree;
      const double MinPhi = 0. * degree;
      const double MaxPhi = 360. * degree;

      // Generating random direction on the diaphragm
      double phi = RandFlat::shoot(&rndm.GetEngine(),
                                   MinPhi, MaxPhi);
      double NCos = 1.; // = 1 is lambertian
                       // = 0 would be isotropic according to DrumPhotonGenerator
      double theta_1 = std::pow(cos(MinTheta), NCos + 1);
      double theta_2 = std::pow(cos(MaxTheta), NCos + 1);
      // th goes like cos^NCOS * d(cos)
      double tmp = (theta_2 - theta_1) * RandFlat::shoot(&rndm.GetEngine(), 0., 1.) + theta_1;

      double theta;
      if (NCos == 0)      theta = acos(tmp);
      else if (NCos == 1) theta = acos(sqrt(tmp));
      else                theta = acos(exp(log(tmp) / (1. + NCos) ) );

      Vector newDir (1., theta, phi, normalCS, Vector::kSpherical);
      return newDir;
      // end code from DrumPhotonGenerator.cc
      }

    }//end of lambertdiffusion

    Vector 
    MirrorDiffusion(utl::RandomEngine& rndm,
                    const Vector& specularDir,
                    const utl::TabulatedFunction* mirrorDiffusionTop,
                    const utl::TabulatedFunction* mirrorDiffusionBot,
                    const double verticalPosOnMirror) {

      // only do diffusion if diffusion function is given in xml aka the TabulatedFunction is not NULL
      utl::TabulatedFunction cumuCurveTop = *mirrorDiffusionTop;
      utl::TabulatedFunction cumuCurveBot = *mirrorDiffusionBot;
      
      double value = RandFlat::shoot(&rndm.GetEngine(), 0., 1.);
      //angles from specular direction (--> z-axis of spherical coordinate system)
      double theta = 0.;
      double phi = 0.;
      const double peripheralRegion = 18 * deg; //upper and lower periphery of mirror for which distribution was measured

      //photon hits central region of the mirror
      //con't linear interpolation btw. scattering angles (from spec. dir.) using cumulative distr.
      if ((verticalPosOnMirror < peripheralRegion) && (verticalPosOnMirror > -1 * peripheralRegion)) {
        float mirrorInterpolWeight = (verticalPosOnMirror+peripheralRegion) / (2 * peripheralRegion); //1 for top of the mirror, 0 for bottom
        theta = (mirrorInterpolWeight * cumuCurveTop.Y(value)) + ((1 - mirrorInterpolWeight) * cumuCurveBot.Y(value));
      }
      
      //upper part of the mirror, using the distr. measured for the upper mirror region
      else if (verticalPosOnMirror > peripheralRegion)
        theta = cumuCurveTop.Y(value);

      //lower part of the mirror, using the distr. measured for the lower mirror region
      else if (verticalPosOnMirror < -1 * peripheralRegion)
        theta = cumuCurveBot.Y(value);
           
      // calculate new photon direction

      // get coordinate system with specular direction as z axis
      const CoordinateSystemPtr internalCS = specularDir.GetCoordinateSystem();
      const double phiIn = specularDir.GetPhi (internalCS);
      const double thetaIn = specularDir.GetTheta (internalCS);
      const CoordinateSystemPtr tmpCS = CoordinateSystem::RotationZ (phiIn, internalCS);
      const CoordinateSystemPtr specularCS = CoordinateSystem::RotationY (thetaIn, tmpCS);

      //random phi direction
      phi = RandFlat::shoot(&rndm.GetEngine(), 0, 2. * kPi);

      Vector newDir (1., theta, phi, specularCS, Vector::kSpherical);
      return newDir;

    } //end of MirrorDiffusion

  } // namespace RTFunctions

} // namesapce TelescopeSimulatorKG2