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 <CLHEP/Random/Randomize.h>

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


using namespace utl;
using namespace std;

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

namespace TelescopeSimulatorKG {

  namespace RTFunctions {

    int Reflection(const utl::Photon& photonIn, const Vector& normal, utl::Photon& photonOut) {

      photonOut = photonIn;
      const Vector& nIn = photonIn.GetDirection();

      const double dot = nIn*normal;   // projection
      photonOut.SetDirection(nIn - 2.*dot*normal); // -> new direction

      return 1;
    }


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

      double cosTheta_in = nIn*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);
      const double dot = nIn*normal;   // projection
      photonReflected.SetDirection(nIn - 2.*dot*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 ! sould 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.0;
      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()
       *****

       result:
       0 : Error
       1 : Ok
  */

    int 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 0;                   // 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 1;

    }







    /*

       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 RandomNormal (utl::RandomEngine& rndm, const Vector& normalIn, 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
      CoordinateSystemPtr inCS = normalIn.GetCoordinateSystem();
      double phiIn = normalIn.GetPhi (inCS);
      double thetaIn = normalIn.GetTheta (inCS);
      CoordinateSystemPtr tmpCS = CoordinateSystem::RotationZ (-phiIn, inCS);
      CoordinateSystemPtr normalCS = CoordinateSystem::RotationY (-thetaIn, tmpCS);


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

      double phi =  RandFlat::shoot (&rndm.GetEngine(), 0., 2.*kPi);
      double sinPhi = sin(phi);
      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;
    }

  } // namespace RTFunctions

} // namesapce TelescopeSimulatorKG