Changeset 19597


Ignore:
Timestamp:
Sep 2, 2019, 8:12:27 PM (2 weeks ago)
Author:
tbretz
Message:
Improvements and fixed, adapted to the Orafol lens characteristics.
File:
1 edited

Legend:

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

    r19595 r19597  
    193193    const double c = n*u.fVectorPart;
    194194
    195     const double R = 1 - r*r*(1-c*c);
    196 
    197     const auto v = r*c + sqrt(R<0 ? 0 : R);
     195    const double rc = r*c;
     196
     197    const double R = 1 - r*r + rc*rc; // = 1 - r*r *(1-c*c);
     198
     199    const double v = rc + sqrt(R<0 ? 0 : R);
    198200
    199201    u = r*u - n*v;
     
    202204Int_t MFresnelLens::ExecuteOptics(MQuaternion &p, MQuaternion &u, const Short_t &wavelength) const
    203205{
    204     const double F = 50;    // [mm] focal length
    205     const double R = 25;    // [mm] radius of lens
    206     const int    N = 25;    // [cnt] Nuber of grooves within radius
    207 
    208     const double w = R/N;   // [mm] Width of a single groove
    209 
    210     const double n = CalcRefractiveIndex(wavelength==0 ? 450 : wavelength);
     206    // Corsika Coordinates are in cm!
     207
     208    const double D = 54.92;  // [cm] Diameter
     209    const double F = 50.21;  // [cm] focal length (see also MGeomCamFAMOUS!)
     210
     211    const double R = D/2;    // [cm] radius of lens
     212    const double w = 0.0100; // [cm] Width of a single groove
     213    const double H = 0.18;   // [cm] Thickness
     214
     215    // Focal length is given at this wavelength
     216    const double n0 = CalcRefractiveIndex(546);
     217
     218    const double n = CalcRefractiveIndex(wavelength==0 ? 546 : wavelength);
     219
     220    // Ray has missed the lens
     221    if (p.R()>R)
     222        return -1;
     223
     224    const double phi = atan(p.R()/F);
    211225
    212226    const int nth = TMath::FloorNint(p.R()/w); // Ray will hit the nth groove
    213227
    214     const double r0  =  nth*w;      // Lower radius of nth groove
    215     const double r1  = (nth+1)*w;  // Upper radius of nth groove
    216     const double cen = (nth+0.5)*w; // Center of nth groove
     228    const double r0 =  nth     *w; // Lower radius of nth groove
     229    const double r1 = (nth+1)  *w; // Upper radius of nth groove
     230    const double rc = (nth+0.5)*w; // Center of nth groove
    217231
    218232    // FIXME: Do we have to check the nth-1 and nth+1 groove as well?
    219233
    220     // Angle to the center of the groove
    221     const double phi = acos(cen / (F/1.125));
    222 
    223     //        theta                peak_z
    224     const Groove g( TMath::Pi()/2 - phi, r0*tan(phi));
     234    // Angle of grooves
     235    // https://shodhganga.inflibnet.ac.in/bitstream/10603/131007/13/09_chapter%202.pdf
     236
     237    // This might be a better reference (Figure 1 Variant 1)
     238    // https://link.springer.com/article/10.3103/S0003701X13010064
     239
     240    // I suggest we just do the calculation ourselves (but we need to
     241    // find out how our lens looks like)
     242    // Not 100% sure, but I think we have SC943
     243    // https://www.orafol.com/en/europe/products/optic-solutions/productlines#pl1
     244
     245    // sin(omega) = R / sqrt(R^2+f^2)
     246    // tan(alpha) = sin(omega) / [ 1 - sqrt(n^2-sin(omega)^2) ]
     247
     248    const double so    = rc / sqrt(rc*rc + F*F);
     249    const double theta = -atan(so / (1-sqrt(n0*n0 - so*so)));
     250
     251    //             theta                peak_z
     252    const Groove g(TMath::Pi()/2 - theta, r0*tan(theta));
    225253
    226254    // Calculate the insersection between the ray and the cone
    227     float dz = CalcIntersection(p, u, g, -3);
     255    float dz = CalcIntersection(p, u, g, -H);
    228256
    229257    // No groove was hit
     
    244272    // FIXME: Surface roughness?
    245273    TVector3 norm;
    246     norm.SetMagThetaPhi(1, g.fTheta+TMath::Pi()/2, p.XYvector().Phi());
     274    norm.SetMagThetaPhi(1, g.fTheta-TMath::Pi()/2, p.XYvector().Phi());
    247275
    248276    // Apply refraction at lens entrance (change directional vector)
     
    253281    const double v = u.fRealPart; // c changes in the medium
    254282    u.fRealPart = n*v;
    255     p.PropagateZ(u, -3);
     283    p.PropagateZ(u, -H);
    256284    u.fRealPart = v;  // Set back to c
    257285
    258     // New normal surface (bottom of lens)
    259     norm.SetMagThetaPhi(1, 0, 0);
    260 
    261286    // Apply refraction at lens exit (change directional vector)
     287    // Normal surface (bottom of lens)
    262288    // FIXME: Surace roughness
    263     ApplyRefraction(u, norm, n, 1);
    264 
    265     // To make this consistent with a mirror system, we now change our coordinate system
     289    ApplyRefraction(u, TVector3(0, 0, 1), n, 1);
     290
     291    // To make this consistent with a mirror system,
     292    // we now change our coordinate system
    266293    // Rays from the lens to the camera are up-going (positive sign)
    267294    u.fVectorPart.SetZ(-u.Z());
    268295
    269     // This is only correct if the focal distance is calculated from the inner lens surface!
    270     p.fVectorPart.SetZ(0);
     296    // In the datasheet, it looks as if F is calculated
     297    // towards the center of the lens.
     298    p.fVectorPart.SetZ(-H/2);
    271299
    272300    return nth; // Keep track of groove index
Note: See TracChangeset for help on using the changeset viewer.