Changeset 19632
 Timestamp:
 Sep 10, 2019, 8:32:15 PM (9 days ago)
 Location:
 trunk/Mars/msimreflector
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Mars/msimreflector/MFresnelLens.cc
r19616 r19632 18 18 ! Author(s): Thomas Bretz, 6/2019 <mailto:tbretz@physik.rwthaachen.de> 19 19 ! 20 ! Copyright: CheObs Software Development, 200020 0920 ! Copyright: CheObs Software Development, 20002019 21 21 ! 22 22 ! … … 26 26 // 27 27 // MFresnelLens 28 // 29 // For some details on definitions please refer to 30 // https://application.wileyvch.de/berlin/journals/op/0704/OP0704_S52_S55.pdf 28 31 // 29 32 ////////////////////////////////////////////////////////////////////////////// … … 33 36 #include <errno.h> 34 37 38 #include <TRandom.h> 39 35 40 #include "MQuaternion.h" 41 #include "MReflection.h" 36 42 37 43 #include "MLog.h" … … 52 58 53 59 fMaxR = 55/2.; 60 } 61 62 //  63 // 64 // Refractive Index of PMMA, according to 65 // https://refractiveindex.info/?shelf=organic&book=poly(methyl_methacrylate)&page=Szczurowski 66 // 67 // n^21=\frac{0.99654 l^2}{l^20.00787}+\frac{0.18964 l^2}{l^20.02191}+\frac{0.00411 l^2}{l^23.85727} 68 // 69 // Returns the refractive index n as a function of wavelength (in nanometers) 70 // 71 double MFresnelLens::RefractiveIndex(double lambda) 72 { 73 const double l2 = lambda*lambda; 74 75 const double c0 = 0.99654/(10.00787e6/l2); 76 const double c1 = 0.18964/(10.02191e6/l2); 77 const double c2 = 0.00411/(13.85727e6/l2); 78 79 return sqrt(1+c0+c1+c2); 80 } 81 82 //  83 // 84 // Ideal angle of the Fresnel surfaces at a distance r from the center 85 // to achieve a focal distance F for a positive Fresnel lens made 86 // from a material with a refractive index n. 87 // A positive Fresnel lens is one which focuses light from infinity 88 // (the side with the grooves) to a point (the flat side of the lens). 89 // 90 // The calculation follows 91 // https://shodhganga.inflibnet.ac.in/bitstream/10603/131007/13/09_chapter%202.pdf 92 // Here, a thin lens is assumed 93 // 94 // The HAWC's Eye lens is an Orafol SC943 95 // https://www.orafol.com/en/europe/products/opticsolutions/productlines#pl1 96 // 97 // sin(omega) = r / sqrt(r^2+F^2) 98 // tan(alpha) = sin(omega) / [ 1  sqrt(n^2sin(omega)^2) ] 99 // 100 // Return alpha [rad] as a function of the radial distance r, the 101 // focal length F and the refractive index n. r and F have to have 102 // the same units. The Slope angle is defined with respect to the plane 103 // of the lens. (0 at the center, decreasing with increasing radial 104 // distance) 105 // 106 double MFresnelLens::SlopeAngle(double r, double F, double n) 107 { 108 double so = r / sqrt(r*r + F*F); 109 return atan(so / (1sqrt(n*n  so*so))); // alpha<0, Range [0deg; 50deg] 110 } 111 112 113 // 114 // Draft angle of the Orafol SC943 According to the thesis of Eichler 115 // and NiggemannTim Niggemann: 116 // 117 // The surface of the lens follows the shape of a parabolic lens to compensate spherical aberration 118 // Draft angle: psi(r) = 3deg + r * 0.0473deg/mm 119 // 120 // The draft angle is returned in radians and is defined w.r.t. to the 121 // normal of the lens surface. (almost 90deg at the center, 122 // decreasing with increasing radial distance) 123 // 124 double MFresnelLens::DraftAngle(double r) 125 { 126 return (3 + r*0.473)*TMath::DegToRad(); // Range [0deg; 15deg] 54 127 } 55 128 … … 97 170 double CalcIntersection(const MQuaternion &p, const MQuaternion &u, const Groove &g, const double &fThickness) 98 171 { 172 /* 173 Z: position of peak 174 175 r_cone(z) = (Zz)*tan(theta) 176 177 (x) (p.x) (u.x) 178 (y) = (p.y) + dz * (u.y) 179 (z) (p.z) (u.z) 180 181 r_line(z) = sqrt(x^2 + y^2) = sqrt( (p.x+dz*u.x)^2 + (p.y+dz*u.y)^2) 182 183 Solved with wmmaxima: 184 185 solve((Hz)^2*t^2 = (px+(zpz)/uz*ux)^2+(py+(zpz)/uz*uy)^2, z); 186 187 [ 188 z=(uz*sqrt((py^2+px^2)*t^2*uz^2+((2*H*py2*py*pz)*t^2*uy+(2*H*px2*px*pz)*t^2*ux)*uz+((pz^22*H*pz+H^2)*t^2px^2)*uy^2+2*px*py*ux*uy+((pz^22*H*pz+H^2)*t^2py^2)*ux^2)H*t^2*uz^2+(py*uypx*ux)*uz+pz*uy^2+pz*ux^2)/(t^2*uz^2uy^2ux^2) 189 z= (uz*sqrt((py^2+px^2)*t^2*uz^2+((2*H*py2*py*pz)*t^2*uy+(2*H*px2*px*pz)*t^2*ux)*uz+((pz^22*H*pz+H^2)*t^2px^2)*uy^2+2*px*py*ux*uy+((pz^22*H*pz+H^2)*t^2py^2)*ux^2)+H*t^2*uz^2+( py*uy+px*ux)*uzpz*uy^2pz*ux^2)/(t^2*uz^2uy^2ux^2) 190 ] 191 */ 192 193 const double Ux = u.X()/u.Z(); 194 const double Uy = u.Y()/u.Z(); 195 196 const double px = p.X(); 197 const double py = p.Y(); 198 const double pz = p.Z(); 199 200 const double H = g.fPeakZ; 201 202 const double t = g.fTanTheta; 203 const double t2 = t*t; 204 205 const double Ur2 = Ux*Ux + Uy*Uy; 206 const double pr2 = px*px + py*py; 207 const double Up2 = Ux*px + Uy*py; 208 209 const double cr2 = Ux*py  Uy*px; 210 211 const double a = t2  Ur2; 212 const double b = Ur2*pz  Up2  H*t2; 213 214 const double h = Hpz; 215 const double h2 = h*h; 216 217 // [ b +sqrt(b^2  4 ac) ] / [ 2a ] 218 219 const double radix = (Ur2*h2 + 2*Up2*h + pr2)*t2  cr2*cr2; 220 221 if (radix<0) 222 return 0; 223 224 double dz[2] = 225 { 226 (b+sqrt(radix))/a, 227 (bsqrt(radix))/a 228 }; 229 230 if (dz[0]<0 && dz[0]>fThickness) 231 return dz[0]; 232 233 if (dz[1]<0 && dz[1]>fThickness) 234 return dz[1]; 235 236 return 0; 237 } 238 239 /* 240 double CalcIntersection(const MQuaternion &p, const MQuaternion &u, const Groove &g, const double &fThickness) 241 { 99 242 // p is the position in the plane of the lens (px, py, 0, t) 100 243 // u is the direction of the photon (ux, uy, dz, dt) … … 103 246 104 247 // Radius at distance z+dz from the tip and at distance r from the axis of the cone 105 // r = (z+dz ) * tan(theta)248 // r = (z+dz*u.z) * tan(theta) 106 249 107 250 // Defining … … 116 259 // (X) (u.x) (p.x) 117 260 // (Y) = dz * (u.y) + (p.y) 118 // (Z) (u.z) ( 0)261 // (Z) (u.z) (p.z) 119 262 120 263 // Equality … … 151 294 return 0; 152 295 296 const int sign = g.fPeakZ>=0 ? 1 : 1; 297 153 298 const double dz[2] = 154 299 { 155 (sqrtBac+B) / a ,156 (sqrtBacB) / a 300 (sqrtBac+B) / a * sign, 301 (sqrtBacB) / a * sign 157 302 }; 158 303 … … 167 312 return 0; 168 313 } 169 170 double CalcRefractiveIndex(double lambda) 171 { 172 // https://refractiveindex.info/?shelf=organic&book=poly(methyl_methacrylate)&page=Szczurowski 173 174 // n^21=\frac{0.99654λ^2}{λ^20.00787}+\frac{0.18964λ^2}{λ^20.02191}+\frac{0.00411λ^2}{λ^23.85727} 175 176 const double l2 = lambda*lambda; 177 178 const double c0 = 0.99654/(10.00787e6/l2); 179 const double c1 = 0.18964/(10.02191e6/l2); 180 const double c2 = 0.00411/(13.85727e6/l2); 181 182 return sqrt(1+c0+c1+c2); 183 } 184 185 bool ApplyRefraction(MQuaternion &u, const TVector3 &n, const double &n1, const double &n2) 186 { 187 // From: https://stackoverflow.com/questions/29758545/howtofindrefractionvectorfromincomingvectorandsurfacenormal 188 // Note that as a default u is supposed to be normalized such that z=1! 189 190 // u: incoming direction (normalized!) 191 // n: normal vector of surface (pointing back into the old medium?) 192 // n1: refractive index of old medium 193 // n2: refractive index of new medium 194 195 // V_refraction = r*V_incedence + (rc  sqrt(1 r^2 (1c^2)))n 196 // r = n1/n2 197 // c = n dot V_incedence. 198 199 // The vector should be normalized already 200 // u.NormalizeVector(); 201 202 const double r = n2/n1; 203 const double c = n*u.fVectorPart; 204 205 const double rc = r*c; 206 207 const double R = 1  r*r + rc*rc; // = 1  (r*r*(1c*c)); 208 209 // What is this? Total reflection? 210 if (R<0) 211 return false; 212 213 const double v = rc + sqrt(R); 214 215 u.fVectorPart = r*u.fVectorPart  v*n; 216 217 return true; 218 } 314 */ 219 315 220 316 Int_t MFresnelLens::ExecuteOptics(MQuaternion &p, MQuaternion &u, const Short_t &wavelength) const … … 226 322 227 323 const double R = D/2; // [cm] radius of lens 228 const double w = 0.01 00;// [cm] Width of a single groove324 const double w = 0.01; // [cm] Width of a single groove 229 325 const double H = 0.25; // [cm] Thickness (from pdf) 230 326 231 327 // Focal length is given at this wavelength 232 const double n0 = CalcRefractiveIndex(546); 233 234 const double n = CalcRefractiveIndex(wavelength==0 ? 546 : wavelength); 328 const double n0 = RefractiveIndex(546); 329 const double n = RefractiveIndex(wavelength==0 ? 546 : wavelength); 330 331 const double r = p.R(); 235 332 236 333 // Ray has missed the lens 237 if ( p.R()>R)334 if (r>R) 238 335 return 1; 239 336 240 const int nth = TMath::FloorNint(p.R()/w); // Ray will hit the nth groove 241 242 const double r0 = nth *w; // Lower radius of nth groove 243 const double r1 = (nth+1) *w; // Upper radius of nth groove 244 const double rc = (nth+0.5)*w; // Center of nth groove 245 246 // FIXME: Do we have to check the nth1 and nth+1 groove as well? 247 248 // Angle of grooves 249 // https://shodhganga.inflibnet.ac.in/bitstream/10603/131007/13/09_chapter%202.pdf 250 251 // This might be a better reference (Figure 1 Variant 1) 252 // https://link.springer.com/article/10.3103/S0003701X13010064 253 254 // I suggest we just do the calculation ourselves (but we need to 255 // find out how our lens looks like) 256 // Not 100% sure, but I think we have SC943 257 // https://www.orafol.com/en/europe/products/opticsolutions/productlines#pl1 258 259 // From the webpage: 260 // Positive Fresnel Lenses ... are usually corrected for spherical aberration. 261 262 // sin(omega) = R / sqrt(R^2+f^2) 263 // tan(alpha) = sin(omega) / [ 1  sqrt(n^2sin(omega)^2) ] 264 265 const double so = rc / sqrt(rc*rc + F*F); 266 const double alpha = atan(so / (1sqrt(n0*n0  so*so))); // alpha<0, Range [0deg; 50deg] 267 268 // Tim Niggemann: 269 // The surface of the lens follows the shape of a parabolic lens to compensate spherical aberration 270 // Draft angle: psi(r) = 3deg + r * 0.0473deg/mm 271 272 const double psi = (3 + r0*4.73e3)*TMath::DegToRad(); // Range [0deg; 15deg] 337 // Calculate the ordinal number of the groove which is hit 338 const unsigned int nth = TMath::FloorNint(r/w); 339 340 // Calculate its lower and upper radius and its center 341 const double r0 = nth*w; // Lower radius of nth groove 342 const double rc = nth*w + w/2; // Upper radius of nth groove 343 const double r1 = nth*w + w; // Center of nth groove 344 345 // Calculate the slope and draft angles 346 const double alpha = SlopeAngle(rc, F, n0); 347 const double psi = DraftAngle(r1); 348 349 // Define the corresponding surfaces 350 const Groove slope(TMath::Pi()/2alpha, r0*tan(alpha)); 351 const Groove draft(psi, r1/tan(psi)); 352 353 // Calculate the insersection of the ray between the two surfaces 354 // with z between 0 and H 355 const double dz_slope = CalcIntersection(p, u, slope, H); 356 const double dz_draft = CalcIntersection(p, u, draft, H); 357 358 // valid means that the photon has hit the fresnel surface, 359 // otherwise the draft surface 360 bool valid = dz_slope>dz_draft  dz_draft>=0; 361 362 // Theta angle of the normal vector of the surface that was hit 363 TVector3 norm; 364 365 // Propagate to the hit point. To get the correct normal vector of 366 // the surface, the hit point must be knwone 367 p.PropagateZ(u, valid ? dz_slope : dz_draft); 368 norm.SetMagThetaPhi(1, valid ? alpha : psiTMath::Pi()/2, p.fVectorPart.Phi()); 369 370 // Estimate reflectivity for this particular hit (n1<n2 => check before) 371 // Probability that the ray gets reflected 372 if (gRandom>Uniform()<SchlickReflectivity(u.fVectorPart, norm, 1, n)) 373 { 374 // Reflect at the surface 375 u *= MReflection(norm); 376 377 // Upgoing rays are lost 378 if (u.Z()>0) 379 return valid ? 3 : 4; 380 381 // Propgate back to z=0 (FIXME: CalcIntersection requires z=0) 382 p.PropagateZ(u, 0); 383 384 // Calc new intersection the other of the two surfaces 385 valid = !valid; 386 387 const double dz = CalcIntersection(p, u, valid ? slope : draft, H); 388 389 // Propagate to the hit point at z=dz (dz<0) and get the new normal vector 390 p.PropagateZ(u, dz); 391 norm.SetMagThetaPhi(1, valid ? alpha : psiTMath::Pi()/2, p.fVectorPart.Phi()); 392 } 393 394 const double check_r = p.R(); 395 396 // Some sanity checks (we loose a few permille due to numerical uncertainties) 397 // If the ray has not hit within the right radii.. reject 398 if (check_r<r0) 399 return valid ? 5 : 6; 400 401 if (check_r>r1) 402 return valid ? 7 : 8; 273 403 274 404 // Find dw value of common zvalue 275 405 // 276 // tan(90psi) = z/dw277 // tan(alpha) = z/(wdw)406 // gz*tan(psi) = w  gw 407 // gz = dw*tan(alpha) 278 408 // 279 // > dw = w/(1tan(90psi)/tan(alpha)) 280 281 // In a simplified world, all photons which hit the draft surface get lost 282 // FIXME: This needs proper calculation in the same manner than for the main surface 283 const double dw = w/(1tan(TMath::Pi()/2psi)/tan(alpha)); 284 if (p.R()>r1dw) 285 return 1; 286 287 // theta peak_z 288 const Groove g(TMath::Pi()/2 + alpha, r0*tan(alpha)); 289 290 // Calculate the insersection between the ray and the cone and z between 0 and H 291 const double dz = CalcIntersection(p, u, g, H); 292 293 // No groove was hit 294 if (dz>=0) 295 return 1; 296 297 // Propagate to the hit point at z=dz (dz<0) 298 p.PropagateZ(u, dz); 299 300 // Check where the ray has hit 301 // If the ray has not hit within the right radius.. reject 302 if (p.R()<=r0  p.R()>r1) 303 return 1; 304 305 // Normal vector of the surface at the hit point 306 // FIXME: Surface roughness? 307 TVector3 norm;; 308 norm.SetMagThetaPhi(1, alpha, p.XYvector().Phi()); 309 310 // Apply refraction at lens entrance (change directional vector) 409 const double Gw = w/(tan(alpha)*tan(psi)+1); 410 const double Gz = Gw*tan(alpha); 411 412 if (valid && check_r>r0+Gw) 413 return 9; 414 if (!valid && check_r<r0+Gw) 415 return 10; 416 417 // Apply refraction at the lens entrance surface (changes directional vector) 311 418 // FIXME: Surace roughness 312 419 if (!ApplyRefraction(u, norm, 1, n)) 313 return 1;420 return valid ? 11 : 12; 314 421 315 422 // Propagate ray until bottom of lens 316 const double v = u.fRealPart; // c changes in the medium317 u.fRealPart = n*v;318 423 p.PropagateZ(u, H); 319 u.fRealPart = v; // Set back to c 320 321 // Apply refraction at lens exit (change directional vector) 322 // Normal surface (bottom of lens) 424 425 // The lens exit surface is a plane in X/Y at Z=H 426 norm = TVector3(0, 0, 1); 427 428 // Apply refraction at the lens exit surface (changes directional vector) 323 429 // FIXME: Surace roughness 324 if (!ApplyRefraction(u, TVector3(0, 0, 1), n, 1)) 325 return 1; 430 if (!ApplyRefraction(u, norm, n, 1)) 431 return valid ? 13 : 14; 432 433 // Estimate reflectivity for this particular hit (n1>n2 => check after) 434 // Probability that the ray gets reflected 435 // (total internal reflection > we won't track that further) 436 if (gRandom>Uniform()<SchlickReflectivity(u.fVectorPart, norm, n, 1)) 437 return valid ? 15 : 16; 326 438 327 439 // To make this consistent with a mirror system, … … 331 443 332 444 // In the datasheet, it looks as if F is calculated 333 // towards the center of the lens. 334 p.fVectorPart.SetZ(H/2); 445 // towards the center of the lens 446 // (Propagating to F means not propagating a distance of FH/2) 447 p.fVectorPart.SetZ(H/2); 335 448 336 449 return nth; // Keep track of groove index 
trunk/Mars/msimreflector/MFresnelLens.h
r19595 r19632 6 6 #endif 7 7 8 class TVector3; 8 9 class MQuaternion; 9 10 … … 27 28 Bool_t IsValid() const { return kTRUE; } 28 29 30 //  31 32 static double RefractiveIndex(double lambda); 33 static double SlopeAngle(double r, double F, double n); 34 static double DraftAngle(double r); 35 29 36 ClassDef(MFresnelLens, 1) // Parameter container storing the description of a lens 30 37 };
Note: See TracChangeset
for help on using the changeset viewer.