Changeset 19629 for trunk/Mars/msimreflector

Timestamp:

09/08/19 23:05:09 (5 years ago)

Author:

tbretz

Message:

Some ray-tracing functions.

Location:

trunk/Mars/msimreflector

Files:

• MOptics.cc (modified) (2 diffs)
• MOptics.h (modified) (2 diffs)

trunk/Mars/msimreflector/MOptics.cc

@@ -30 +30 @@
 #include "MOptics.h"
 
+#include <TRandom.h>
+
+#include "MQuaternion.h"
+#include "MReflection.h"
+
 ClassImp(MOptics);
 
@@ -43 +48 @@
     fTitle = title ? title : "Optics base class";
 }
+
+// --------------------------------------------------------------------------
+//
+// Critical angle for total reflection asin(n2/n1), or 90deg of n1<n2
+// https://graphics.stanford.edu/courses/cs148-10-summer/docs/2006--degreve--reflection_refraction.pdf
+//
+double MOptics::CriticalAngle(double n1, double n2)
+{
+    return n1>=n2 ? asin(n2/n1) : TMath::Pi()/2;
+}
+
+// --------------------------------------------------------------------------
+//
+// This returns an approximation for the reflectivity for a ray
+// passing from a medium with refractive index n1 to a medium with
+// refractive index n2. This approximation is only valid for similar
+// values of n1 and n2 (e.g. 1 and 1.5 but not 1 and 5).   /
+//
+// For n1<n2 the function has to be called with theta being the
+// angle between the surface-normal and the incoming ray, for
+// n1>n2, alpha is the angle between the outgoing ray and the
+// surface-normal.
+//
+// It is not valid to call the function for n1>n2 when alpha exceeds
+// the critical angle. Alpha is defined in the range [0deg, 90deg]
+// For Alpha [90;180], The angle is assumed  the absolute value of the cosine is used.
+//
+// alpha must be given in radians and defined between [0deg and 90deg]
+//
+// Taken from:
+// https://graphics.stanford.edu/courses/cs148-10-summer/docs/2006--degreve--reflection_refraction.pdf
+//
+double MOptics::SchlickReflectivity(double alpha, double n1, double n2)
+{
+    const double r0 = (n1-n2)/(n1+n2);
+    const double r2 = r0 * r0;
+
+    const double p  = 1-cos(alpha);
+
+    return r2 + (1-r2) * p*p*p*p*p;
+}
+
+// --------------------------------------------------------------------------
+//
+// Same as SchlickReflectivity(double,double,double) but takes the direction
+// vector of the incoming or outgoing ray as u and the normal vector of the
+// surface. The normal vector points from n2 to n1.
+//
+double MOptics::SchlickReflectivity(const TVector3 &u, const TVector3 &n, double n1, double n2)
+{
+    const double r0 = (n1-n2)/(n1+n2);
+    const double r2 = r0 * r0;
+
+    const double c = -n*u;
+
+    const double p  = 1-c;
+
+    return r2 + (1-r2) * p*p*p*p*p;
+}
+
+// --------------------------------------------------------------------------
+//
+// Applies refraction to the direction vector u on a surface with the normal
+// vector n for a ray coming from a medium with refractive index n1
+// passing into a medium with refractive index n2. Note that the normal
+// vector is defined pointing from medium n2 to medium n1 (so opposite
+// of u).
+//
+// Note that u and n must be normalized before calling the function!
+//
+// Solution accoridng to (Section 6)
+// https://graphics.stanford.edu/courses/cs148-10-summer/docs/2006--degreve--reflection_refraction.pdf
+//
+// u_out = n1/n2 * u_in + (n1/n2 * cos(theta_in) - sqrt(1-sin^2(theta_out)) * n
+// sin^2(theta_out) = (n1/n2)^2 * sin^2(theta_in) = (n1/n2)^2 * (1-cos^2(theta_in))
+// cos(theta_in) = +- u_in * n
+//
+// In case of success, the vector u is altered and true is returned. In case
+// of failure (incident angle above critical angle for total internal reflection)
+// vector u stays untouched and false is returned.
+//
+bool MOptics::ApplyRefraction(TVector3 &u, const TVector3 &n, double n1, double n2)
+{
+    // The vector should be normalized already
+    // u.NormalizeVector();
+    // n.NormalizeVector();
+
+    // Usually: Theta [0;90]
+    // Here:    Theta [90;180] => c=|n*u| < 0
+
+    const double r = n1/n2;
+    const double c = -n*u;
+
+    const double rc = r*c;
+
+    const double s2 = r*r - rc*rc; // sin(theta_out)^2 = r^2*(1-cos(theta_in)^2)
+
+    if (s2>1)
+        return false;
+
+    const double v = rc - sqrt(1-s2);
+
+    u = r*u + v*n;
+
+    return true;
+}
+
+// --------------------------------------------------------------------------
+//
+// In addition to calculating ApplyRefraction(u.fVectorPart, n, n1, n2),
+// the fourth component (inverse of the speed) is multiplied with n1/n2
+// of tramission was successfull (ApplyRefraction retruned true)
+//
+// Note that .fVectorPart and n must be normalized before calling the function!
+//
+bool MOptics::ApplyRefraction(MQuaternion &u, const TVector3 &n, double n1, double n2)
+{
+    if (!ApplyRefraction(u.fVectorPart, n, n1, n2))
+        return false;
+
+    u.fRealPart *= n1/n2;
+    return true;
+}
+
+// --------------------------------------------------------------------------
+//
+// Applies a transition from one medium (n1) to another medium (n2)
+// n1 is always the medium where the photon comes from.
+//
+// Uses Schlick's approximation for reflection
+//
+// The direction of the normal vector does not matter, it is automatically
+// aligned (opposite) of the direction of the incoming photon.
+//
+// Based on
+// https://graphics.stanford.edu/courses/cs148-10-summer/docs/2006--degreve--reflection_refraction.pdf
+//
+// Returns:
+//
+//  0 error (total internal reflection from optical thin to optical thick medium)
+//
+//  1 total internal reflection applied
+//  2 Monte Carlo    reflection applied
+//
+//  3 nothing done (n1==n2, transmission)
+//  4 refraction applied from optical thick to optically thinner medium
+//  5 refraction applied from optical thin  to optically thicker medium
+//
+int MOptics::ApplyTransitionFast(TVector3 &dir, TVector3 n, double n1, double n2)
+{
+    if (n1==n2)
+        return 3;
+
+    // The normal vector must point in the same direction as
+    // the photon is moving and thus cos(theta)=[90;180]
+    if (dir*n>0)
+        n *= -1;
+
+    TVector3 u(dir);
+
+    // From optical thick to optical thin medium
+    if (n1>n2)
+    {
+        // Check for refraction...
+        // ... if possible:     use exit angle for calculating reflectivity
+        // ... if not possible: reflect ray
+        if (!ApplyRefraction(u, n, n1, n2))
+        {
+            // Total Internal Reflection
+            dir *= MReflection(n);
+            return 1;
+        }
+    }
+
+    const double reflectivity = SchlickReflectivity(u, n, n1, n2);
+
+    // ----- Case of reflection ----
+
+    if (gRandom->Uniform()<reflectivity)
+    {
+        dir *= MReflection(n);
+        return 2;
+    }
+
+    // ----- Case of refraction ----
+
+    // From optical thick to optical thin medium
+    if (n1>n2)
+    {
+        // Now we know that refraction was correctly applied
+        dir = u;
+        return 4;
+    }
+
+    // From optical thin to optical thick medium
+    // Still need to apply refraction
+
+    if (ApplyRefraction(dir, n, n1, n2))
+        return 5;
+
+    // ERROR => Total Internal Reflection not possible
+    // This should never happen
+    return 0;
+}
+
+// --------------------------------------------------------------------------
+//
+// Applies a transition from one medium (n1) to another medium (n2)
+// n1 is always the medium where the photon comes from.
+//
+// Uses Fresnel's equation for calculating reflection
+//
+// The direction of the normal vector does not matter, it is automatically
+// aligned (opposite) of the direction of the incoming photon.
+//
+// Based on
+// https://graphics.stanford.edu/courses/cs148-10-summer/docs/2006--degreve--reflection_refraction.pdf
+//
+// Returns:
+//
+//  0 error (total internal reflection from optical thin to optical thick medium)
+//
+//  1 total internal reflection applied
+//  2 Monte Carlo    reflection applied
+//
+//  3 nothing done (n1==n2, transmission)
+//  4 refraction applied
+//
+int MOptics::ApplyTransition(TVector3 &u, TVector3 n, double n1, double n2)
+{
+    if (n1==n2)
+        return 3;
+
+    // The normal vector must point in the same direction as
+    // the photon is moving and thus cos(theta)=[90;180]
+    if (u*n>0)
+        n *= -1;
+
+    // Calculate refraction outgoing direction (Snell's law)
+    const double r   = n1/n2;
+    const double ci  = -n*u;           // cos(theta_in)
+    const double rci =  r*ci;          // n1/n2 * cos(theta_in)
+    const double st2 = r*r - rci*rci;  // sin(theta_out)^2 = r^2*(1-cos(theta_in)^2)
+
+    if (st2>1)
+    {
+        // Total Internal Reflection
+        u *= MReflection(n);
+        return 1;
+
+    }
+
+    const double ct  = sqrt(1-st2);    // cos(theta_out) = sqrt(1 - sin(theta_out)^2)
+    const double rct = r*ct;           // n1/n2 * cos(theta_out)
+
+    // Calculate reflectivity for none polarized rays (Fresnel's equation)
+    const double Rt = (rci - ct)/(rci + ct);
+    const double Rp = (rct - ci)/(rct + ci);
+
+    const double reflectivity = (Rt*Rt + Rp*Rp)/2;
+
+    if (gRandom->Uniform()<reflectivity)
+    {
+        // ----- Case of reflection ----
+        u *= MReflection(n);
+        return 2;
+    }
+    else
+    {
+        // ----- Case of refraction ----
+        u = r*u + (rci-ct)*n;
+        return 4;
+    }
+}
+
+int MOptics::ApplyTransitionFast(MQuaternion &u, const TVector3 &n, double n1, double n2)
+{
+    const int rc = ApplyTransitionFast(u.fVectorPart, n, n1, n2);
+    if (rc>=3)
+        u.fRealPart *= n1/n2;
+    return rc;
+}
+
+int MOptics::ApplyTransition(MQuaternion &u, const TVector3 &n, double n1, double n2)
+{
+    const int rc = ApplyTransition(u.fVectorPart, n, n1, n2);
+    if (rc>=3)
+        u.fRealPart *= n1/n2;
+    return rc;
+}

trunk/Mars/msimreflector/MOptics.h

@@ -6 +6 @@
 #endif
 
+class TVector3;
 class MQuaternion;
 
@@ -21 +22 @@
     virtual Bool_t IsValid() const = 0;
 
+    // -----------------------------------------------------------
+
+    static double CriticalAngle(double n1, double n2);
+
+    static double SchlickReflectivity(double alpha, double n1, double n2);
+    static double SchlickReflectivity(const TVector3 &u, const TVector3 &n, double n1, double n2);
+
+    static bool ApplyRefraction(TVector3 &u, const TVector3 &n, double n1, double n2);
+    static bool ApplyRefraction(MQuaternion &u, const TVector3 &n, double n1, double n2);
+
+    static int ApplyTransitionFast(TVector3 &u, TVector3 n, double n1, double n2);
+    static int ApplyTransitionFast(MQuaternion &u, const TVector3 &n, double n1, double n2);
+
+    static int ApplyTransition(TVector3 &u, TVector3 n, double n1, double n2);
+    static int ApplyTransition(MQuaternion &u, const TVector3 &n, double n1, double n2);
+
     ClassDef(MOptics, 1) // Base class for optics (reflector, lens)
 };