Lens.cc

Go to the documentation of this file.

#include "Lens.h"
#include "RTFunctions.h"

#include <TPolyLine3D.h>
#include <TObjArray.h>

#include <utl/MathConstants.h>
#include <utl/AugerUnits.h>

#include <utl/Photon.h>
#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/RandomEngine.h>
#include <utl/CoordinateSystemPtr.h>

#include <fdet/Telescope.h>
#include <fdet/Corrector.h>

#include <det/Detector.h>

#include <atm/Atmosphere.h>
#include <atm/ProfileResult.h>

#include <iostream>
#include <sstream>
#include <utility>

using namespace TelescopeSimulatorKG;
using namespace utl;
using namespace std;


Lens::Lens (utl::RandomEngine& rndm, const fdet::Telescope& tel) {

  fRandom = &rndm;

  const fdet::Corrector& corrector = tel.GetCorrector();

  fTelCS = tel.GetTelescopeCoordinateSystem();
  fOrigin = Point (0,0,0, fTelCS);

  fR1 = corrector.GetInnerRadius();
  fR2 = corrector.GetOuterRadius();

  // const double Tm = corrector.GetMeanLensThickness();      // average lens thickness
  fSigRho = corrector.GetSigmaNormal();
  fRefractiveIndex = &(corrector.GetRefractiveIndex());
  fTransmittance = &(corrector.GetTransmittance());

  // air refractive index
  const atm::ProfileResult & refractiveIndexVsHeight =
    det::Detector::GetInstance().GetAtmosphere().EvaluateRefractionIndexVsHeight();
  fInd1 = refractiveIndexVsHeight.Y(refractiveIndexVsHeight.MinX());
  
} // End of Lens::Lens


RTResult Lens::Trace(const utl::Photon& photonIn, utl::Photon& photonOut) {

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

  if (nIn.GetZ(fTelCS)>0) {
    // the light is going out of the telescope
    photonOut = photonIn;
    return eInvalidLens;
  }
  
  const double x = pIn.GetX(fTelCS);
  const double y = pIn.GetY(fTelCS);
  const double r2 = x*x + y*y;
  const double r = sqrt(r2);
  if (r > fR2) {
    photonOut = photonIn;
    return eMissedDiaphragm;     //the light don't pass through the aperture window
  }
  
  if (r < fR1) {
    photonOut = photonIn;
    // the light went through the aperture window,
    return eOK;
  }

  const Vector normal1 = Curvature(x, y);
  const double z = Profile(r);
  
  // calculate the direction of the light inside the corrector lens
  const double r_index_12 = fInd1 / fRefractiveIndex->Y(lambda);
  //double r_index_12 = 1.0/1.5;
  RTFunctions::PhotonList photonsRefracted;
  RTFunctions::Refraction(r_index_12, photonIn, normal1, photonsRefracted);
  if (photonsRefracted.size() != 2) {
    return eReflectedByLens;
  }
  Photon photonRefracted = photonsRefracted[0];
    
  const double cosThetaPhoton = - photonRefracted.GetDirection().GetCosTheta(fTelCS);
  const double thickness = z / cosThetaPhoton;

  if (fTransmittance->Y(lambda) == 0) 
    return eAbsorbedByLens; 

  const double internal_absorption = exp(thickness/(1.*cm)*log(fTransmittance->Y(lambda)));
  
  // introduces imperfections at corrector ring surface
  // fRt->Normal(n1, fSigRho, &fISeed, n);               // NOT DONE



  // ------------- SECOND SURFACE of corrector ring ------------
  // calculate the position of light at the second surface of the lens.
  photonRefracted.SetPosition(photonRefracted.GetPosition() + thickness * photonRefracted.GetDirection());
  
  // fRt->Normal(n1, fSigRho, &fISeed, n);          // NOT DONE YET
  Vector normal2 (0., 0., 1., fTelCS); //second surface is plane.
  Vector nOut;
  RTFunctions::PhotonList photonsOut;
  RTFunctions::Refraction(1 / r_index_12, photonRefracted, normal2, photonsOut);
  if (photonsOut.size() != 2) {
    return eReflectedByLens; // total internal reflection inside lens
  }
  
  photonOut = photonsOut[0];
  const double weight = photonOut.GetWeight();
  photonOut.SetWeight(weight*internal_absorption);

  return eOK;
} // End of Lens::Trace


TObjArray* Lens::Draw() {

  TObjArray* objs = new TObjArray();

  int color = 1;
  int width = 2;

  int nSeg = 15;
  for (int i=0; i<nSeg; i++) {

    double theta_1 = 360.*deg/nSeg * i;
    double theta_2 = 360.*deg/nSeg * (i+1);

    TPolyLine3D *l_in = new TPolyLine3D (1);
    l_in->SetPoint (0, cos(theta_1)*fR1, sin(theta_1)*fR1, 0);
    l_in->SetPoint (1, cos(theta_2)*fR1, sin(theta_2)*fR1, 0);
    l_in->SetLineColor(color);
    l_in->SetLineWidth(width);
    l_in->Draw();
    objs->AddLast(l_in);

    TPolyLine3D *l_out = new TPolyLine3D (1);
    l_out->SetPoint (0, cos(theta_1)*fR2, sin(theta_1)*fR2, 0);
    l_out->SetPoint (1, cos(theta_2)*fR2, sin(theta_2)*fR2, 0);
    l_out->SetLineColor(color);
    l_out->SetLineWidth(width);
    l_out->Draw();
    objs->AddLast(l_out);

    TPolyLine3D *l_seg = new TPolyLine3D (1);
    l_seg->SetPoint (0, cos(theta_1)*fR1, sin(theta_1)*fR1, 0);
    l_seg->SetPoint (1, cos(theta_1)*fR2, sin(theta_1)*fR2, 0);
    l_seg->SetLineColor(color);
    l_seg->SetLineWidth(width);
    l_seg->Draw();
    objs->AddLast(l_seg);
  }
  return objs;
}


// thikness profile as a function of radial distance
Vector Lens::Curvature(const double x, const double y) const {
  
  const double r = sqrt(x*x + y*y);

  static const double f = 3.4*m - 1.743*m;
  static const double Rd = .85*m;
  //static const double z0 = 10.*mm;
  static const double A = 3./2. * Rd*Rd;
  static const double n = 1.5;

  static const double denominator = 32.*(n-1.)*f*f*f;
  static const double a1 = 1./denominator; // r^4
  static const double a2 = A/denominator;  // r^2

  double dTdr = 4.*a1*r*r*r - 2.*a2*r;
  // dTdr *= .9;                     ///// this would be a bit better

  /*
  // Tilos KA prototyp
  double dTdr = 6.*4.e-18 * mm * pow(r,5) * pow (mm,-6);
  */

  
  // new KG default from GAP ...

  double nr = dTdr;
  double nz = 1.;
  double norm = sqrt (nr*nr + nz*nz);

  nr /= norm;
  nz /= norm;

  double phi = atan2(y,x);

  return Vector (-nr*cos(phi), -nr*sin(phi), nz, fTelCS);
  

  /*
  // old FDSim code:
  double r2 = r*r;
  double aux = 2*(-1.2461221E-02 + 2*8.6476346E-03*r2 + 3*2.0361284E-03*r2*r2);

  Vector norm = Vector (-x*aux, -y*aux, 1.0, fTelCS);
  norm.Normalize();

  return norm;
  */

}


// thikness profile as a function of radial distance
double Lens::Profile(const double r) const {

  double r2 = r*r;

  /*
  // Tilos KA prototyp
  r2 /= (mm*mm);
  double z = 4.e-18*mm * pow (r2,3) + 0.91*mm;
  */

  
  // new KG default from GAP ...
  static const double f = 3.4*m - 1.743*m;
  static const double Rd = .85*m;
  static const double A = 3./2. * Rd*Rd;
  static const double n = 1.5;

  static const double denominator = 32.*(n-1.)*f*f*f;
  static const double a1 = 1./denominator; // r^4
  static const double a2 = A/denominator;  // r^2

  static const double z0 = 10.*mm;

  double z = z0 + a1*r2*r2 - a2*r2;
  

  /*
  // old FDSim code:
  double aux = 2*(-1.2461221E-02 + 2*8.6476346E-03*r2 + 3*2.0361284E-03*r2*r2);

  Vector norm = Vector (-x*aux, -y*aux, 1.0, fTelCS);
  norm.Normalize();

  double z = sqrt (pow(norm.GetX (fTelCS),2) + pow(norm.GetY (fTelCS),2));
  */

  return z;

} // End of Lens::Profile



// Configure (x)emacs for this file ...
// Local Variables:
// mode:c++
// compile-command: "make -C .. -k"
// End: