Changeset 19629


Ignore:
Timestamp:
Sep 8, 2019, 11:05:09 PM (10 days ago)
Author:
tbretz
Message:
Some ray-tracing functions.
Location:
trunk/Mars/msimreflector
Files:
2 edited

Legend:

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

    r19538 r19629  
    3030#include "MOptics.h"
    3131
     32#include <TRandom.h>
     33
     34#include "MQuaternion.h"
     35#include "MReflection.h"
     36
    3237ClassImp(MOptics);
    3338
     
    4348    fTitle = title ? title : "Optics base class";
    4449}
     50
     51// --------------------------------------------------------------------------
     52//
     53// Critical angle for total reflection asin(n2/n1), or 90deg of n1<n2
     54// https://graphics.stanford.edu/courses/cs148-10-summer/docs/2006--degreve--reflection_refraction.pdf
     55//
     56double MOptics::CriticalAngle(double n1, double n2)
     57{
     58    return n1>=n2 ? asin(n2/n1) : TMath::Pi()/2;
     59}
     60
     61// --------------------------------------------------------------------------
     62//
     63// This returns an approximation for the reflectivity for a ray
     64// passing from a medium with refractive index n1 to a medium with
     65// refractive index n2. This approximation is only valid for similar
     66// values of n1 and n2 (e.g. 1 and 1.5 but not 1 and 5).   /
     67//
     68// For n1<n2 the function has to be called with theta being the
     69// angle between the surface-normal and the incoming ray, for
     70// n1>n2, alpha is the angle between the outgoing ray and the
     71// surface-normal.
     72//
     73// It is not valid to call the function for n1>n2 when alpha exceeds
     74// the critical angle. Alpha is defined in the range [0deg, 90deg]
     75// For Alpha [90;180], The angle is assumed  the absolute value of the cosine is used.
     76//
     77// alpha must be given in radians and defined between [0deg and 90deg]
     78//
     79// Taken from:
     80// https://graphics.stanford.edu/courses/cs148-10-summer/docs/2006--degreve--reflection_refraction.pdf
     81//
     82double MOptics::SchlickReflectivity(double alpha, double n1, double n2)
     83{
     84    const double r0 = (n1-n2)/(n1+n2);
     85    const double r2 = r0 * r0;
     86
     87    const double p  = 1-cos(alpha);
     88
     89    return r2 + (1-r2) * p*p*p*p*p;
     90}
     91
     92// --------------------------------------------------------------------------
     93//
     94// Same as SchlickReflectivity(double,double,double) but takes the direction
     95// vector of the incoming or outgoing ray as u and the normal vector of the
     96// surface. The normal vector points from n2 to n1.
     97//
     98double MOptics::SchlickReflectivity(const TVector3 &u, const TVector3 &n, double n1, double n2)
     99{
     100    const double r0 = (n1-n2)/(n1+n2);
     101    const double r2 = r0 * r0;
     102
     103    const double c = -n*u;
     104
     105    const double p  = 1-c;
     106
     107    return r2 + (1-r2) * p*p*p*p*p;
     108}
     109
     110// --------------------------------------------------------------------------
     111//
     112// Applies refraction to the direction vector u on a surface with the normal
     113// vector n for a ray coming from a medium with refractive index n1
     114// passing into a medium with refractive index n2. Note that the normal
     115// vector is defined pointing from medium n2 to medium n1 (so opposite
     116// of u).
     117//
     118// Note that u and n must be normalized before calling the function!
     119//
     120// Solution accoridng to (Section 6)
     121// https://graphics.stanford.edu/courses/cs148-10-summer/docs/2006--degreve--reflection_refraction.pdf
     122//
     123// u_out = n1/n2 * u_in + (n1/n2 * cos(theta_in) - sqrt(1-sin^2(theta_out)) * n
     124// sin^2(theta_out) = (n1/n2)^2 * sin^2(theta_in) = (n1/n2)^2 * (1-cos^2(theta_in))
     125// cos(theta_in) = +- u_in * n
     126//
     127// In case of success, the vector u is altered and true is returned. In case
     128// of failure (incident angle above critical angle for total internal reflection)
     129// vector u stays untouched and false is returned.
     130//
     131bool MOptics::ApplyRefraction(TVector3 &u, const TVector3 &n, double n1, double n2)
     132{
     133    // The vector should be normalized already
     134    // u.NormalizeVector();
     135    // n.NormalizeVector();
     136
     137    // Usually: Theta [0;90]
     138    // Here:    Theta [90;180] => c=|n*u| < 0
     139
     140    const double r = n1/n2;
     141    const double c = -n*u;
     142
     143    const double rc = r*c;
     144
     145    const double s2 = r*r - rc*rc; // sin(theta_out)^2 = r^2*(1-cos(theta_in)^2)
     146
     147    if (s2>1)
     148        return false;
     149
     150    const double v = rc - sqrt(1-s2);
     151
     152    u = r*u + v*n;
     153
     154    return true;
     155}
     156
     157// --------------------------------------------------------------------------
     158//
     159// In addition to calculating ApplyRefraction(u.fVectorPart, n, n1, n2),
     160// the fourth component (inverse of the speed) is multiplied with n1/n2
     161// of tramission was successfull (ApplyRefraction retruned true)
     162//
     163// Note that .fVectorPart and n must be normalized before calling the function!
     164//
     165bool MOptics::ApplyRefraction(MQuaternion &u, const TVector3 &n, double n1, double n2)
     166{
     167    if (!ApplyRefraction(u.fVectorPart, n, n1, n2))
     168        return false;
     169
     170    u.fRealPart *= n1/n2;
     171    return true;
     172}
     173
     174// --------------------------------------------------------------------------
     175//
     176// Applies a transition from one medium (n1) to another medium (n2)
     177// n1 is always the medium where the photon comes from.
     178//
     179// Uses Schlick's approximation for reflection
     180//
     181// The direction of the normal vector does not matter, it is automatically
     182// aligned (opposite) of the direction of the incoming photon.
     183//
     184// Based on
     185// https://graphics.stanford.edu/courses/cs148-10-summer/docs/2006--degreve--reflection_refraction.pdf
     186//
     187// Returns:
     188//
     189//  0 error (total internal reflection from optical thin to optical thick medium)
     190//
     191//  1 total internal reflection applied
     192//  2 Monte Carlo    reflection applied
     193//
     194//  3 nothing done (n1==n2, transmission)
     195//  4 refraction applied from optical thick to optically thinner medium
     196//  5 refraction applied from optical thin  to optically thicker medium
     197//
     198int MOptics::ApplyTransitionFast(TVector3 &dir, TVector3 n, double n1, double n2)
     199{
     200    if (n1==n2)
     201        return 3;
     202
     203    // The normal vector must point in the same direction as
     204    // the photon is moving and thus cos(theta)=[90;180]
     205    if (dir*n>0)
     206        n *= -1;
     207
     208    TVector3 u(dir);
     209
     210    // From optical thick to optical thin medium
     211    if (n1>n2)
     212    {
     213        // Check for refraction...
     214        // ... if possible:     use exit angle for calculating reflectivity
     215        // ... if not possible: reflect ray
     216        if (!ApplyRefraction(u, n, n1, n2))
     217        {
     218            // Total Internal Reflection
     219            dir *= MReflection(n);
     220            return 1;
     221        }
     222    }
     223
     224    const double reflectivity = SchlickReflectivity(u, n, n1, n2);
     225
     226    // ----- Case of reflection ----
     227
     228    if (gRandom->Uniform()<reflectivity)
     229    {
     230        dir *= MReflection(n);
     231        return 2;
     232    }
     233
     234    // ----- Case of refraction ----
     235
     236    // From optical thick to optical thin medium
     237    if (n1>n2)
     238    {
     239        // Now we know that refraction was correctly applied
     240        dir = u;
     241        return 4;
     242    }
     243
     244    // From optical thin to optical thick medium
     245    // Still need to apply refraction
     246
     247    if (ApplyRefraction(dir, n, n1, n2))
     248        return 5;
     249
     250    // ERROR => Total Internal Reflection not possible
     251    // This should never happen
     252    return 0;
     253}
     254
     255// --------------------------------------------------------------------------
     256//
     257// Applies a transition from one medium (n1) to another medium (n2)
     258// n1 is always the medium where the photon comes from.
     259//
     260// Uses Fresnel's equation for calculating reflection
     261//
     262// The direction of the normal vector does not matter, it is automatically
     263// aligned (opposite) of the direction of the incoming photon.
     264//
     265// Based on
     266// https://graphics.stanford.edu/courses/cs148-10-summer/docs/2006--degreve--reflection_refraction.pdf
     267//
     268// Returns:
     269//
     270//  0 error (total internal reflection from optical thin to optical thick medium)
     271//
     272//  1 total internal reflection applied
     273//  2 Monte Carlo    reflection applied
     274//
     275//  3 nothing done (n1==n2, transmission)
     276//  4 refraction applied
     277//
     278int MOptics::ApplyTransition(TVector3 &u, TVector3 n, double n1, double n2)
     279{
     280    if (n1==n2)
     281        return 3;
     282
     283    // The normal vector must point in the same direction as
     284    // the photon is moving and thus cos(theta)=[90;180]
     285    if (u*n>0)
     286        n *= -1;
     287
     288    // Calculate refraction outgoing direction (Snell's law)
     289    const double r   = n1/n2;
     290    const double ci  = -n*u;           // cos(theta_in)
     291    const double rci =  r*ci;          // n1/n2 * cos(theta_in)
     292    const double st2 = r*r - rci*rci;  // sin(theta_out)^2 = r^2*(1-cos(theta_in)^2)
     293
     294    if (st2>1)
     295    {
     296        // Total Internal Reflection
     297        u *= MReflection(n);
     298        return 1;
     299
     300    }
     301
     302    const double ct  = sqrt(1-st2);    // cos(theta_out) = sqrt(1 - sin(theta_out)^2)
     303    const double rct = r*ct;           // n1/n2 * cos(theta_out)
     304
     305    // Calculate reflectivity for none polarized rays (Fresnel's equation)
     306    const double Rt = (rci - ct)/(rci + ct);
     307    const double Rp = (rct - ci)/(rct + ci);
     308
     309    const double reflectivity = (Rt*Rt + Rp*Rp)/2;
     310
     311    if (gRandom->Uniform()<reflectivity)
     312    {
     313        // ----- Case of reflection ----
     314        u *= MReflection(n);
     315        return 2;
     316    }
     317    else
     318    {
     319        // ----- Case of refraction ----
     320        u = r*u + (rci-ct)*n;
     321        return 4;
     322    }
     323}
     324
     325int MOptics::ApplyTransitionFast(MQuaternion &u, const TVector3 &n, double n1, double n2)
     326{
     327    const int rc = ApplyTransitionFast(u.fVectorPart, n, n1, n2);
     328    if (rc>=3)
     329        u.fRealPart *= n1/n2;
     330    return rc;
     331}
     332
     333int MOptics::ApplyTransition(MQuaternion &u, const TVector3 &n, double n1, double n2)
     334{
     335    const int rc = ApplyTransition(u.fVectorPart, n, n1, n2);
     336    if (rc>=3)
     337        u.fRealPart *= n1/n2;
     338    return rc;
     339}
  • trunk/Mars/msimreflector/MOptics.h

    r19587 r19629  
    66#endif
    77
     8class TVector3;
    89class MQuaternion;
    910
     
    2122    virtual Bool_t IsValid() const = 0;
    2223
     24    // -----------------------------------------------------------
     25
     26    static double CriticalAngle(double n1, double n2);
     27
     28    static double SchlickReflectivity(double alpha, double n1, double n2);
     29    static double SchlickReflectivity(const TVector3 &u, const TVector3 &n, double n1, double n2);
     30
     31    static bool ApplyRefraction(TVector3 &u, const TVector3 &n, double n1, double n2);
     32    static bool ApplyRefraction(MQuaternion &u, const TVector3 &n, double n1, double n2);
     33
     34    static int ApplyTransitionFast(TVector3 &u, TVector3 n, double n1, double n2);
     35    static int ApplyTransitionFast(MQuaternion &u, const TVector3 &n, double n1, double n2);
     36
     37    static int ApplyTransition(TVector3 &u, TVector3 n, double n1, double n2);
     38    static int ApplyTransition(MQuaternion &u, const TVector3 &n, double n1, double n2);
     39
    2340    ClassDef(MOptics, 1) // Base class for optics (reflector, lens)
    2441};
Note: See TracChangeset for help on using the changeset viewer.