source: trunk/Mars/msimreflector/MFresnelLens.cc@ 19594

/* ======================================================================== *\
!
! *
! * This file is part of CheObs, the Modular Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appears in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Thomas Bretz,  6/2019 <mailto:tbretz@physik.rwth-aachen.de>
!
!   Copyright: CheObs Software Development, 2000-2009
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//  MLens
//
//////////////////////////////////////////////////////////////////////////////
#include "MLens.h"

#include <fstream>
#include <errno.h>

#include "MQuaternion.h"

#include "MLog.h"
#include "MLogManip.h"

ClassImp(MLens);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor
//
MLens::MLens(const char *name, const char *title)
{
    fName  = name  ? name  : "MLens";
    fTitle = title ? title : "Parameter container storing a collection of several mirrors (reflector)";

    fMaxR = 25;
}

// --------------------------------------------------------------------------
//
// Return the total Area of all mirrors. Note, that it is recalculated
// with any call.
//
Double_t MLens::GetA() const
{
    return fMaxR*fMaxR*TMath::Pi();
}

// --------------------------------------------------------------------------
//
// Check with a rough estimate whether a photon can hit the reflector.
//
Bool_t MLens::CanHit(const MQuaternion &p) const
{
    // p is given in the reflectory coordinate frame. This is meant as a
    // fast check without lengthy calculations to omit all photons which
    // cannot hit the reflector at all
    return p.R2()<fMaxR*fMaxR;
}

struct Groove
{
    double fTheta;
    double fPeakZ;

    double fTanTheta;
    double fTanTheta2[3];

    Groove() { }
    Groove(double theta, double z) : fTheta(theta), fPeakZ(z)
    {
        fTanTheta = tan(theta);

        fTanTheta2[0] = fTanTheta*fTanTheta;
        fTanTheta2[1] = fTanTheta*fTanTheta * fPeakZ;
        fTanTheta2[2] = fTanTheta*fTanTheta * fPeakZ*fPeakZ;
    }
};

float CalcIntersection(const MQuaternion &p, const MQuaternion &u, const Groove &g, const double &fThickness)
{
   // p is the position  in the plane  of the lens (px, py, 0,  t)
   // u is the direction of the photon             (ux, uy, dz, dt)

   // Assuming a cone with peak at z and opening angle theta

   // Radius at distance z+dz from the tip and at distance r from the axis of the cone
   // r = (z-dz) * tan(theta)

   // Defining
   // zc := z+dz

   // Surface of the cone
   // (X)         ( tan(theta) * sin(phi) )
   // (Y)  = zc * ( tan(theta) * cos(phi) )
   // (Z)         (          1            )

   // Line
   // (X)        (u.x)   (p.x)
   // (Y) = dz * (u.y) + (p.y)
   // (Z)        (u.z)   ( 0 )

   // Equality
   //
   // X^2+Y^2 = r^2
   //
   // (dz*u.x+p.x)^2 + (dz*u.y+p.y)^2 = (z-dz)^2 * tan(theta)^2
   //
   // dz^2*ux^2 + px^2 + 2*dz*ux*px + dz^2*uy^2 + py^2 + 2*dz*uy*py = (z^2*tan(theta)^2 + dz^2*tan(theta)^2 - 2*z*dz*tan(theta)^2
   //
   // dz^2*(ux^2 + uy^2 - tan(theta)^2)  +  2*dz*(ux*px + uy*py + z*tan(theta)^2)  +  (px^2+py^2 - z^2*tan(theta)^2)  =  0

   // t   := tan(theta)
   // a   := ux^2    + uy^2    -     t^2
   // b/2 := ux*px   + uy*py   - z  *t^2 := B
   // c   :=    px^2 +    py^2 - z^2*t^2

   // dz = [ -b +- sqrt(b^2 - 4*a*c) ] / [2*a]
   // dz = [ -B +- sqrt(B^2 -   a*c) ] /    a

    const double a = u.R2()                    - g.fTanTheta2[0];
    const double B = u.XYvector()*p.XYvector() + g.fTanTheta2[1];
    const double c = p.R2()                    - g.fTanTheta2[2];

    const double B2 = B*B;
    const double ac = a*c;

    if (B2<ac)
        return 0;

    const double sqrtBac = sqrt(B2 - a*c);

    if (a==0)
        return 0;

    const double dz[2] =
    {
        (sqrtBac-B) / a,
        (sqrtBac+B) / a
    };

    // Only one dz[i] can be within the lens (due to geometrical reasons

    if (dz[0]<0 && dz[0]>fThickness)
        return dz[0];

    if (dz[1]<0 && dz[1]>fThickness)
        return dz[1];

    return 0;
}

double CalcRefractiveIndex(const double &lambda)
{
    // https://refractiveindex.info/?shelf=organic&book=poly(methyl_methacrylate)&page=Szczurowski

    // n^2-1=\frac{0.99654λ^2}{λ^2-0.00787}+\frac{0.18964λ^2}{λ^2-0.02191}+\frac{0.00411λ^2}{λ^2-3.85727}

    const double l2 = lambda*lambda;

    const double c0 = 0.99654/(1-0.00787/l2);
    const double c1 = 0.18964/(1-0.02191/l2);
    const double c2 = 0.00411/(1-3.85727/l2);

    return sqrt(1+c0+c1+c2);
}

void ApplyRefraction(MQuaternion &u, const TVector3 &n, const double &n1, const double &n2)
{
    // u: incoming direction
    // n: normal vector of surface
    // n2: refractive index of new(?) medium
    // n1: refractive index of old(?) medium

    const double r = n2/n1;
    const double c = n*u.fVectorPart;

    const double R = 1 - r*r*(1-c*c);

    const auto v = r*c + sqrt(R<0 ? 0 : R);

    u = r*u - n*v;
}

Int_t MLens::ExecuteOptics(MQuaternion &p, MQuaternion &u, const Short_t &wavelength) const
{
    const double F = 50;    // [mm] focal length
    const double R = 25;    // [mm] radius of lens
    const int    N = 25;    // [cnt] Nuber of grooves within radius

    const double w = R/N;   // [mm] Width of a single groove

    const double n = CalcRefractiveIndex(wavelength==0 ? 450 : wavelength);

    const int nth = TMath::FloorNint(p.R()/w); // Ray will hit the nth groove

    const double r0  =  nth*w;      // Lower radius of nth groove
    const double r1  = (nth+1)*w;   // Upper radius of nth groove
    const double cen = (nth+0.5)*w; // Center of nth groove

    // FIXME: Do we have to check the nth-1 and nth+1 groove as well?

    // Angle to the center of the groove
    const double phi = acos(cen / (F/1.125));

    //        theta                peak_z
    const Groove g( TMath::Pi()/2 - phi, r0*tan(phi));

    // Calculate the insersection between the ray and the cone
    float dz = CalcIntersection(p, u, g, -3);

    // No groove was hit
    if (dz>=0)
        return -1;

    // Propagate to the hit point
    p.PropagateZ(u, dz);

    // Check where the ray has hit
    const double pr = p.R();

    // If the ray has not hit within the right radius.. reject
    if (pr<=r0 || pr>r1)
        return -1;

    // Normal vector of the surface at the hit point
    // FIXME: Surface roughness?
    TVector3 norm;
    norm.SetMagThetaPhi(1, g.fTheta+TMath::Pi()/2, p.XYvector().Phi());

    // Apply refraction at lens entrance (change directional vector)
    // FIXME: Surace roughness
    ApplyRefraction(u, norm, 1, n);

    // Propagate ray until bottom of lens
    const double v = u.fRealPart; // c changes in the medium 
    u.fRealPart = n*v;
    p.PropagateZ(u, -3);
    u.fRealPart = v;  // Set back to c

    // New normal surface (bottom of lens)
    norm.SetMagThetaPhi(1, 0, 0);

    // Apply refraction at lens exit (change directional vector)
    // FIXME: Surace roughness
    ApplyRefraction(u, norm, n, 1);

    // To make this consistent with a mirror system, we now change our coordinate system
    // Rays from the lens to the camera are up-going (positive sign)
    u.fVectorPart.SetZ(-u.Z());

    // This is only correct if the focal distance is calculated from the inner lens surface!
    p.fVectorPart.SetZ(0);

    return nth; // Keep track of groove index
}