source: trunk/Mars/msimreflector/MOptics.cc @ 19629

Last change on this file since 19629 was 19629, checked in by tbretz, 7 months ago
Some ray-tracing functions.
File size: 10.3 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of CheObs, the Modular Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appears in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18!   Author(s): Thomas Bretz, 6/2019 <mailto:tbretz@physik.rwth-aachen.de>
19!
20!   Copyright: CheObs Software Development, 2000-2019
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27//  MOptics
28//
29//////////////////////////////////////////////////////////////////////////////
30#include "MOptics.h"
31
32#include <TRandom.h>
33
34#include "MQuaternion.h"
35#include "MReflection.h"
36
37ClassImp(MOptics);
38
39using namespace std;
40
41// --------------------------------------------------------------------------
42//
43// Default constructor
44//
45MOptics::MOptics(const char *name, const char *title)
46{
47    fName  = name  ? name  : "MOptics";
48    fTitle = title ? title : "Optics base class";
49}
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}
Note: See TracBrowser for help on using the repository browser.