Changeset 19632


Ignore:
Timestamp:
Sep 10, 2019, 8:32:15 PM (9 days ago)
Author:
tbretz
Message:
Some updates ... not yet final.
Location:
trunk/Mars/msimreflector
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/Mars/msimreflector/MFresnelLens.cc

    r19616 r19632  
    1818!   Author(s): Thomas Bretz,  6/2019 <mailto:tbretz@physik.rwth-aachen.de>
    1919!
    20 !   Copyright: CheObs Software Development, 2000-2009
     20!   Copyright: CheObs Software Development, 2000-2019
    2121!
    2222!
     
    2626//
    2727//  MFresnelLens
     28//
     29//  For some details on definitions please refer to
     30//  https://application.wiley-vch.de/berlin/journals/op/07-04/OP0704_S52_S55.pdf
    2831//
    2932//////////////////////////////////////////////////////////////////////////////
     
    3336#include <errno.h>
    3437
     38#include <TRandom.h>
     39
    3540#include "MQuaternion.h"
     41#include "MReflection.h"
    3642
    3743#include "MLog.h"
     
    5258
    5359    fMaxR = 55/2.;
     60}
     61
     62// --------------------------------------------------------------------------
     63//
     64// Refractive Index of PMMA, according to
     65// https://refractiveindex.info/?shelf=organic&book=poly(methyl_methacrylate)&page=Szczurowski
     66//
     67// n^2-1=\frac{0.99654 l^2}{l^2-0.00787}+\frac{0.18964 l^2}{l^2-0.02191}+\frac{0.00411 l^2}{l^2-3.85727}
     68//
     69// Returns the refractive index n as a function of wavelength (in nanometers)
     70//
     71double MFresnelLens::RefractiveIndex(double lambda)
     72{
     73    const double l2 = lambda*lambda;
     74
     75    const double c0 = 0.99654/(1-0.00787e6/l2);
     76    const double c1 = 0.18964/(1-0.02191e6/l2);
     77    const double c2 = 0.00411/(1-3.85727e6/l2);
     78
     79    return sqrt(1+c0+c1+c2);
     80}
     81
     82// --------------------------------------------------------------------------
     83//
     84// Ideal angle of the Fresnel surfaces at a distance r from the center
     85// to achieve a focal distance F for a positive Fresnel lens made
     86// from a material with a refractive index n.
     87// A positive Fresnel lens is one which focuses light from infinity
     88// (the side with the grooves) to a point (the flat side of the lens).
     89//
     90// The calculation follows
     91// https://shodhganga.inflibnet.ac.in/bitstream/10603/131007/13/09_chapter%202.pdf
     92// Here, a thin lens is assumed
     93//
     94// The HAWC's Eye lens is an Orafol SC943
     95// https://www.orafol.com/en/europe/products/optic-solutions/productlines#pl1
     96//
     97// sin(omega) = r / sqrt(r^2+F^2)
     98// tan(alpha) = sin(omega) / [ 1 - sqrt(n^2-sin(omega)^2) ]
     99//
     100// Return alpha [rad] as a function of the radial distance r, the
     101// focal length F and the refractive index n. r and F have to have
     102// the same units. The Slope angle is defined with respect to the plane
     103// of the lens. (0 at the center, decreasing with increasing radial
     104// distance)
     105//
     106double MFresnelLens::SlopeAngle(double r, double F, double n)
     107{
     108    double so = r / sqrt(r*r + F*F);
     109    return atan(so / (1-sqrt(n*n - so*so))); // alpha<0, Range [0deg; -50deg]
     110}
     111
     112
     113//
     114// Draft angle of the Orafol SC943 According to the thesis of Eichler
     115// and NiggemannTim Niggemann:
     116//
     117// The surface of the lens follows the shape of a parabolic lens to compensate spherical aberration
     118// Draft angle: psi(r) = 3deg + r * 0.0473deg/mm
     119//
     120// The draft angle is returned in radians and is defined w.r.t. to the
     121// normal of the lens surface. (almost 90deg at the center,
     122// decreasing with increasing radial distance)
     123//
     124double MFresnelLens::DraftAngle(double r)
     125{
     126    return (3 + r*0.473)*TMath::DegToRad(); // Range [0deg; 15deg]
    54127}
    55128
     
    97170double CalcIntersection(const MQuaternion &p, const MQuaternion &u, const Groove &g, const double &fThickness)
    98171{
     172/*
     173    Z: position of peak
     174
     175    r_cone(z) = (Z-z)*tan(theta)
     176
     177    (x)   (p.x)          (u.x)
     178    (y) = (p.y)  + dz *  (u.y)
     179    (z)   (p.z)          (u.z)
     180
     181    r_line(z) = sqrt(x^2 + y^2) = sqrt( (p.x+dz*u.x)^2 + (p.y+dz*u.y)^2)
     182
     183    Solved with wmmaxima:
     184
     185    solve((H-z)^2*t^2 = (px+(z-pz)/uz*ux)^2+(py+(z-pz)/uz*uy)^2, z);
     186
     187     [
     188      z=-(uz*sqrt((py^2+px^2)*t^2*uz^2+((2*H*py-2*py*pz)*t^2*uy+(2*H*px-2*px*pz)*t^2*ux)*uz+((pz^2-2*H*pz+H^2)*t^2-px^2)*uy^2+2*px*py*ux*uy+((pz^2-2*H*pz+H^2)*t^2-py^2)*ux^2)-H*t^2*uz^2+(-py*uy-px*ux)*uz+pz*uy^2+pz*ux^2)/(t^2*uz^2-uy^2-ux^2)
     189      z= (uz*sqrt((py^2+px^2)*t^2*uz^2+((2*H*py-2*py*pz)*t^2*uy+(2*H*px-2*px*pz)*t^2*ux)*uz+((pz^2-2*H*pz+H^2)*t^2-px^2)*uy^2+2*px*py*ux*uy+((pz^2-2*H*pz+H^2)*t^2-py^2)*ux^2)+H*t^2*uz^2+( py*uy+px*ux)*uz-pz*uy^2-pz*ux^2)/(t^2*uz^2-uy^2-ux^2)
     190     ]
     191*/
     192
     193    const double Ux = u.X()/u.Z();
     194    const double Uy = u.Y()/u.Z();
     195
     196    const double px = p.X();
     197    const double py = p.Y();
     198    const double pz = p.Z();
     199
     200    const double H   = g.fPeakZ;
     201
     202    const double t   = g.fTanTheta;
     203    const double t2  = t*t;
     204
     205    const double Ur2 = Ux*Ux + Uy*Uy;
     206    const double pr2 = px*px + py*py;
     207    const double Up2 = Ux*px + Uy*py;
     208
     209    const double cr2 = Ux*py - Uy*px;
     210
     211    const double a   = t2 - Ur2;
     212    const double b   = Ur2*pz - Up2 - H*t2;
     213
     214    const double h   = H-pz;
     215    const double h2  = h*h;
     216
     217    // [ -b +-sqrt(b^2 - 4 ac) ] / [ 2a ]
     218
     219    const double radix = (Ur2*h2 + 2*Up2*h + pr2)*t2 - cr2*cr2;
     220
     221    if (radix<0)
     222        return 0;
     223
     224    double dz[2] =
     225    {
     226        (-b+sqrt(radix))/a,
     227        (-b-sqrt(radix))/a
     228    };
     229
     230     if (dz[0]<0 && dz[0]>fThickness)
     231        return dz[0];
     232
     233    if (dz[1]<0 && dz[1]>fThickness)
     234        return dz[1];
     235
     236    return 0;
     237}
     238
     239/*
     240double CalcIntersection(const MQuaternion &p, const MQuaternion &u, const Groove &g, const double &fThickness)
     241{
    99242   // p is the position  in the plane  of the lens (px, py, 0,  t)
    100243   // u is the direction of the photon             (ux, uy, dz, dt)
     
    103246
    104247   // Radius at distance z+dz from the tip and at distance r from the axis of the cone
    105    // r = (z+dz) * tan(theta)
     248   // r = (z+dz*u.z) * tan(theta)
    106249
    107250   // Defining
     
    116259   // (X)        (u.x)   (p.x)
    117260   // (Y) = dz * (u.y) + (p.y)
    118    // (Z)        (u.z)   ( 0 )
     261   // (Z)        (u.z)   (p.z)
    119262
    120263   // Equality
     
    151294        return 0;
    152295
     296    const int sign = g.fPeakZ>=0 ? 1 : -1;
     297
    153298    const double dz[2] =
    154299    {
    155         (sqrtBac+B) / a,
    156         (sqrtBac-B) / a
     300        (sqrtBac+B) / a * sign,
     301        (sqrtBac-B) / a * sign
    157302    };
    158303
     
    167312    return 0;
    168313}
    169 
    170 double CalcRefractiveIndex(double lambda)
    171 {
    172     // https://refractiveindex.info/?shelf=organic&book=poly(methyl_methacrylate)&page=Szczurowski
    173 
    174     // n^2-1=\frac{0.99654λ^2}{λ^2-0.00787}+\frac{0.18964λ^2}{λ^2-0.02191}+\frac{0.00411λ^2}{λ^2-3.85727}
    175 
    176     const double l2 = lambda*lambda;
    177 
    178     const double c0 = 0.99654/(1-0.00787e6/l2);
    179     const double c1 = 0.18964/(1-0.02191e6/l2);
    180     const double c2 = 0.00411/(1-3.85727e6/l2);
    181 
    182     return sqrt(1+c0+c1+c2);
    183 }
    184 
    185 bool ApplyRefraction(MQuaternion &u, const TVector3 &n, const double &n1, const double &n2)
    186 {
    187     // From: https://stackoverflow.com/questions/29758545/how-to-find-refraction-vector-from-incoming-vector-and-surface-normal
    188     // Note that as a default u is supposed to be normalized such that z=1!
    189 
    190     // u: incoming direction (normalized!)
    191     // n: normal vector of surface (pointing back into the old medium?)
    192     // n1: refractive index of old medium
    193     // n2: refractive index of new medium
    194 
    195     // V_refraction = r*V_incedence + (rc - sqrt(1- r^2 (1-c^2)))n
    196     // r = n1/n2
    197     // c = -n dot V_incedence.
    198 
    199     // The vector should be normalized already
    200     // u.NormalizeVector();
    201 
    202     const double r = n2/n1;
    203     const double c = n*u.fVectorPart;
    204 
    205     const double rc = r*c;
    206 
    207     const double R = 1 - r*r + rc*rc; // = 1 - (r*r*(1-c*c));
    208 
    209     // What is this? Total reflection?
    210     if (R<0)
    211         return false;
    212 
    213     const double v = rc + sqrt(R);
    214 
    215     u.fVectorPart = r*u.fVectorPart - v*n;
    216 
    217     return true;
    218 }
     314*/
    219315
    220316Int_t MFresnelLens::ExecuteOptics(MQuaternion &p, MQuaternion &u, const Short_t &wavelength) const
     
    226322
    227323    const double R = D/2;    // [cm] radius of lens
    228     const double w = 0.0100; // [cm] Width of a single groove
     324    const double w = 0.01 // [cm] Width of a single groove
    229325    const double H = 0.25;   // [cm] Thickness (from pdf)
    230326
    231327    // Focal length is given at this wavelength
    232     const double n0 = CalcRefractiveIndex(546);
    233 
    234     const double n = CalcRefractiveIndex(wavelength==0 ? 546 : wavelength);
     328    const double n0 = RefractiveIndex(546);
     329    const double n  = RefractiveIndex(wavelength==0 ? 546 : wavelength);
     330
     331    const double r = p.R();
    235332
    236333    // Ray has missed the lens
    237     if (p.R()>R)
     334    if (r>R)
    238335        return -1;
    239336
    240     const int nth = TMath::FloorNint(p.R()/w); // Ray will hit the nth groove
    241 
    242     const double r0 =  nth     *w; // Lower radius of nth groove
    243     const double r1 = (nth+1)  *w; // Upper radius of nth groove
    244     const double rc = (nth+0.5)*w; // Center of nth groove
    245 
    246     // FIXME: Do we have to check the nth-1 and nth+1 groove as well?
    247 
    248     // Angle of grooves
    249     // https://shodhganga.inflibnet.ac.in/bitstream/10603/131007/13/09_chapter%202.pdf
    250 
    251     // This might be a better reference (Figure 1 Variant 1)
    252     // https://link.springer.com/article/10.3103/S0003701X13010064
    253 
    254     // I suggest we just do the calculation ourselves (but we need to
    255     // find out how our lens looks like)
    256     // Not 100% sure, but I think we have SC943
    257     // https://www.orafol.com/en/europe/products/optic-solutions/productlines#pl1
    258 
    259     // From the web-page:
    260     // Positive Fresnel Lenses ... are usually corrected for spherical aberration.
    261 
    262     // sin(omega) = R / sqrt(R^2+f^2)
    263     // tan(alpha) = sin(omega) / [ 1 - sqrt(n^2-sin(omega)^2) ]
    264 
    265     const double so    = rc / sqrt(rc*rc + F*F);
    266     const double alpha = atan(so / (1-sqrt(n0*n0 - so*so))); // alpha<0, Range [0deg; -50deg]
    267 
    268     // Tim Niggemann:
    269     // The surface of the lens follows the shape of a parabolic lens to compensate spherical aberration
    270     // Draft angle: psi(r) = 3deg + r * 0.0473deg/mm
    271 
    272     const double psi = (3 + r0*4.73e-3)*TMath::DegToRad(); // Range [0deg; 15deg]
     337    // Calculate the ordinal number of the groove which is hit
     338    const unsigned int nth = TMath::FloorNint(r/w);
     339
     340    // Calculate its lower and upper radius and its center
     341    const double r0 = nth*w;       // Lower radius of nth groove
     342    const double rc = nth*w + w/2; // Upper radius of nth groove
     343    const double r1 = nth*w + w;   // Center of nth groove
     344
     345    // Calculate the slope and draft angles
     346    const double alpha = -SlopeAngle(rc, F, n0);
     347    const double psi   =  DraftAngle(r1);
     348
     349    // Define the corresponding surfaces
     350    const Groove slope(TMath::Pi()/2-alpha,  r0*tan(alpha));
     351    const Groove draft(-psi,                -r1/tan(psi));
     352
     353    // Calculate the insersection of the ray between the two surfaces
     354    // with z between 0 and -H
     355    const double dz_slope = CalcIntersection(p, u, slope, -H);
     356    const double dz_draft = CalcIntersection(p, u, draft, -H);
     357
     358    // valid means that the photon has hit the fresnel surface,
     359    // otherwise the draft surface
     360    bool valid = dz_slope>dz_draft || dz_draft>=0;
     361
     362    // Theta angle of the normal vector of the surface that was hit
     363    TVector3 norm;
     364
     365    // Propagate to the hit point. To get the correct normal vector of
     366    // the surface, the hit point must be knwone
     367    p.PropagateZ(u, valid ? dz_slope : dz_draft);
     368    norm.SetMagThetaPhi(1, valid ? alpha : psi-TMath::Pi()/2, p.fVectorPart.Phi());
     369
     370    // Estimate reflectivity for this particular hit (n1<n2 => check before)
     371    // Probability that the ray gets reflected
     372    if (gRandom->Uniform()<SchlickReflectivity(u.fVectorPart, norm, 1, n))
     373    {
     374        // Reflect at the surface
     375        u *= MReflection(norm);
     376
     377        // Upgoing rays are lost
     378        if (u.Z()>0)
     379            return valid ? -3 : -4;
     380
     381        // Propgate back to z=0 (FIXME: CalcIntersection requires z=0)
     382        p.PropagateZ(u, 0);
     383
     384        // Calc new intersection the other of the two surfaces
     385        valid = !valid;
     386
     387        const double dz = CalcIntersection(p, u, valid ? slope : draft, -H);
     388
     389        // Propagate to the hit point at z=dz (dz<0) and get the new normal vector
     390        p.PropagateZ(u, dz);
     391        norm.SetMagThetaPhi(1, valid ? alpha : psi-TMath::Pi()/2, p.fVectorPart.Phi());
     392    }
     393
     394    const double check_r = p.R();
     395
     396    // Some sanity checks (we loose a few permille due to numerical uncertainties)
     397    // If the ray has not hit within the right radii.. reject
     398    if (check_r<r0)
     399        return valid ? -5 : -6;
     400
     401    if (check_r>r1)
     402        return valid ? -7 : -8;
    273403
    274404    // Find dw value of common z-value
    275405    //
    276     // tan(90-psi)  =  z/dw
    277     // tan(alpha)   = -z/(w-dw)
     406    // gz*tan(psi) = w - gw
     407    // gz = dw*tan(alpha)
    278408    //
    279     // -> dw = w/(1-tan(90-psi)/tan(alpha))
    280 
    281     // In a simplified world, all photons which hit the draft surface get lost
    282     // FIXME: This needs proper calculation in the same manner than for the main surface
    283     const double dw = w/(1-tan(TMath::Pi()/2-psi)/tan(alpha));
    284     if (p.R()>r1-dw)
    285         return -1;
    286 
    287     //             theta                   peak_z
    288     const Groove g(TMath::Pi()/2 + alpha, -r0*tan(alpha));
    289 
    290     // Calculate the insersection between the ray and the cone and z between 0 and -H
    291     const double dz = CalcIntersection(p, u, g, -H);
    292 
    293     // No groove was hit
    294     if (dz>=0)
    295         return -1;
    296 
    297     // Propagate to the hit point at z=dz (dz<0)
    298     p.PropagateZ(u, dz);
    299 
    300     // Check where the ray has hit
    301     // If the ray has not hit within the right radius.. reject
    302     if (p.R()<=r0 || p.R()>r1)
    303         return -1;
    304 
    305     // Normal vector of the surface at the hit point
    306     // FIXME: Surface roughness?
    307     TVector3 norm;;
    308     norm.SetMagThetaPhi(1, alpha, p.XYvector().Phi());
    309 
    310     // Apply refraction at lens entrance (change directional vector)
     409    const double Gw =  w/(tan(alpha)*tan(psi)+1);
     410    const double Gz = -Gw*tan(alpha);
     411
     412    if (valid && check_r>r0+Gw)
     413        return -9;
     414    if (!valid && check_r<r0+Gw)
     415        return -10;
     416
     417    // Apply refraction at the lens entrance surface (changes directional vector)
    311418    // FIXME: Surace roughness
    312419    if (!ApplyRefraction(u, norm, 1, n))
    313         return -1;
     420        return valid ? -11 : -12;
    314421
    315422    // Propagate ray until bottom of lens
    316     const double v = u.fRealPart; // c changes in the medium
    317     u.fRealPart = n*v;
    318423    p.PropagateZ(u, -H);
    319     u.fRealPart = v;  // Set back to c
    320 
    321     // Apply refraction at lens exit (change directional vector)
    322     // Normal surface (bottom of lens)
     424
     425    // The lens exit surface is a plane in X/Y at Z=-H
     426    norm = TVector3(0, 0, 1);
     427
     428    // Apply refraction at the lens exit surface (changes directional vector)
    323429    // FIXME: Surace roughness
    324     if (!ApplyRefraction(u, TVector3(0, 0, 1), n, 1))
    325         return -1;
     430    if (!ApplyRefraction(u, norm, n, 1))
     431        return valid ? -13 : -14;
     432
     433    // Estimate reflectivity for this particular hit (n1>n2 => check after)
     434    // Probability that the ray gets reflected
     435    // (total internal reflection -> we won't track that further)
     436    if (gRandom->Uniform()<SchlickReflectivity(u.fVectorPart, norm, n, 1))
     437        return valid ? -15 : -16;
    326438
    327439    // To make this consistent with a mirror system,
     
    331443
    332444    // In the datasheet, it looks as if F is calculated
    333     // towards the center of the lens.
    334     p.fVectorPart.SetZ(-H/2);
     445    // towards the center of the lens
     446    // (Propagating to F means not propagating a distance of F-H/2)
     447    p.fVectorPart.SetZ(H/2);
    335448
    336449    return nth; // Keep track of groove index
  • trunk/Mars/msimreflector/MFresnelLens.h

    r19595 r19632  
    66#endif
    77
     8class TVector3;
    89class MQuaternion;
    910
     
    2728    Bool_t IsValid() const { return kTRUE; }
    2829
     30    // -----------------------------------------------------------
     31
     32    static double RefractiveIndex(double lambda);
     33    static double SlopeAngle(double r, double F, double n);
     34    static double DraftAngle(double r);
     35
    2936    ClassDef(MFresnelLens, 1) // Parameter container storing the description of a lens
    3037};
Note: See TracChangeset for help on using the changeset viewer.