Changeset 19629
 Timestamp:
 Sep 8, 2019, 11:05:09 PM (6 months ago)
 Location:
 trunk/Mars/msimreflector
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Mars/msimreflector/MOptics.cc
r19538 r19629 30 30 #include "MOptics.h" 31 31 32 #include <TRandom.h> 33 34 #include "MQuaternion.h" 35 #include "MReflection.h" 36 32 37 ClassImp(MOptics); 33 38 … … 43 48 fTitle = title ? title : "Optics base class"; 44 49 } 50 51 //  52 // 53 // Critical angle for total reflection asin(n2/n1), or 90deg of n1<n2 54 // https://graphics.stanford.edu/courses/cs14810summer/docs/2006degrevereflection_refraction.pdf 55 // 56 double 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 surfacenormal and the incoming ray, for 70 // n1>n2, alpha is the angle between the outgoing ray and the 71 // surfacenormal. 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/cs14810summer/docs/2006degrevereflection_refraction.pdf 81 // 82 double MOptics::SchlickReflectivity(double alpha, double n1, double n2) 83 { 84 const double r0 = (n1n2)/(n1+n2); 85 const double r2 = r0 * r0; 86 87 const double p = 1cos(alpha); 88 89 return r2 + (1r2) * 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 // 98 double MOptics::SchlickReflectivity(const TVector3 &u, const TVector3 &n, double n1, double n2) 99 { 100 const double r0 = (n1n2)/(n1+n2); 101 const double r2 = r0 * r0; 102 103 const double c = n*u; 104 105 const double p = 1c; 106 107 return r2 + (1r2) * 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/cs14810summer/docs/2006degrevereflection_refraction.pdf 122 // 123 // u_out = n1/n2 * u_in + (n1/n2 * cos(theta_in)  sqrt(1sin^2(theta_out)) * n 124 // sin^2(theta_out) = (n1/n2)^2 * sin^2(theta_in) = (n1/n2)^2 * (1cos^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 // 131 bool 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*(1cos(theta_in)^2) 146 147 if (s2>1) 148 return false; 149 150 const double v = rc  sqrt(1s2); 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 // 165 bool 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/cs14810summer/docs/2006degrevereflection_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 // 198 int 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/cs14810summer/docs/2006degrevereflection_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 // 278 int 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*(1cos(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(1st2); // 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 + (rcict)*n; 321 return 4; 322 } 323 } 324 325 int 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 333 int 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 6 6 #endif 7 7 8 class TVector3; 8 9 class MQuaternion; 9 10 … … 21 22 virtual Bool_t IsValid() const = 0; 22 23 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 23 40 ClassDef(MOptics, 1) // Base class for optics (reflector, lens) 24 41 };
Note: See TracChangeset
for help on using the changeset viewer.