Index: /trunk/Mars/msimreflector/MOptics.cc =================================================================== --- /trunk/Mars/msimreflector/MOptics.cc (revision 19628) +++ /trunk/Mars/msimreflector/MOptics.cc (revision 19629) @@ -30,4 +30,9 @@ #include "MOptics.h" +#include + +#include "MQuaternion.h" +#include "MReflection.h" + ClassImp(MOptics); @@ -43,2 +48,292 @@ fTitle = title ? title : "Optics base class"; } + +// -------------------------------------------------------------------------- +// +// Critical angle for total reflection asin(n2/n1), or 90deg of 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 n1n2, 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()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()=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; +} Index: /trunk/Mars/msimreflector/MOptics.h =================================================================== --- /trunk/Mars/msimreflector/MOptics.h (revision 19628) +++ /trunk/Mars/msimreflector/MOptics.h (revision 19629) @@ -6,4 +6,5 @@ #endif +class TVector3; class MQuaternion; @@ -21,4 +22,20 @@ 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) };