Index: trunk/Mars/msimreflector/MFresnelLens.cc
===================================================================
--- trunk/Mars/msimreflector/MFresnelLens.cc	(revision 19596)
+++ trunk/Mars/msimreflector/MFresnelLens.cc	(revision 19597)
@@ -193,7 +193,9 @@
     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);
+    const double rc = r*c;
+
+    const double R = 1 - r*r + rc*rc; // = 1 - r*r *(1-c*c);
+
+    const double v = rc + sqrt(R<0 ? 0 : R);
 
     u = r*u - n*v;
@@ -202,28 +204,54 @@
 Int_t MFresnelLens::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);
+    // Corsika Coordinates are in cm!
+
+    const double D = 54.92;  // [cm] Diameter
+    const double F = 50.21;  // [cm] focal length (see also MGeomCamFAMOUS!)
+
+    const double R = D/2;    // [cm] radius of lens
+    const double w = 0.0100; // [cm] Width of a single groove
+    const double H = 0.18;   // [cm] Thickness
+
+    // Focal length is given at this wavelength
+    const double n0 = CalcRefractiveIndex(546);
+
+    const double n = CalcRefractiveIndex(wavelength==0 ? 546 : wavelength);
+
+    // Ray has missed the lens
+    if (p.R()>R)
+        return -1;
+
+    const double phi = atan(p.R()/F);
 
     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
+    const double r0 =  nth     *w; // Lower radius of nth groove
+    const double r1 = (nth+1)  *w; // Upper radius of nth groove
+    const double rc = (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));
+    // Angle of grooves
+    // https://shodhganga.inflibnet.ac.in/bitstream/10603/131007/13/09_chapter%202.pdf
+
+    // This might be a better reference (Figure 1 Variant 1)
+    // https://link.springer.com/article/10.3103/S0003701X13010064
+
+    // I suggest we just do the calculation ourselves (but we need to
+    // find out how our lens looks like)
+    // Not 100% sure, but I think we have SC943
+    // https://www.orafol.com/en/europe/products/optic-solutions/productlines#pl1
+
+    // sin(omega) = R / sqrt(R^2+f^2)
+    // tan(alpha) = sin(omega) / [ 1 - sqrt(n^2-sin(omega)^2) ]
+
+    const double so    = rc / sqrt(rc*rc + F*F);
+    const double theta = -atan(so / (1-sqrt(n0*n0 - so*so)));
+
+    //             theta                peak_z
+    const Groove g(TMath::Pi()/2 - theta, r0*tan(theta));
 
     // Calculate the insersection between the ray and the cone
-    float dz = CalcIntersection(p, u, g, -3);
+    float dz = CalcIntersection(p, u, g, -H);
 
     // No groove was hit
@@ -244,5 +272,5 @@
     // FIXME: Surface roughness?
     TVector3 norm;
-    norm.SetMagThetaPhi(1, g.fTheta+TMath::Pi()/2, p.XYvector().Phi());
+    norm.SetMagThetaPhi(1, g.fTheta-TMath::Pi()/2, p.XYvector().Phi());
 
     // Apply refraction at lens entrance (change directional vector)
@@ -253,20 +281,20 @@
     const double v = u.fRealPart; // c changes in the medium 
     u.fRealPart = n*v;
-    p.PropagateZ(u, -3);
+    p.PropagateZ(u, -H);
     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)
+    // Normal surface (bottom of lens)
     // FIXME: Surace roughness
-    ApplyRefraction(u, norm, n, 1);
-
-    // To make this consistent with a mirror system, we now change our coordinate system
+    ApplyRefraction(u, TVector3(0, 0, 1), 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);
+    // In the datasheet, it looks as if F is calculated
+    // towards the center of the lens.
+    p.fVectorPart.SetZ(-H/2);
 
     return nth; // Keep track of groove index