theia.optics.optic

24 ''' 25 26 Optic class. 27 28 This class is a base class for optics which may interact with Gaussian 29 beams and return transmitted and reflected beams (mirrors, lenses, etc.) 30 31 The way an optic interacts with a beam (if it adds to the order of a beam 32 upon reflection or transmission on HR or AR etc) is specified by the 33 action integers RonHR, TonHR, RonAR, TonAR of the optic. A beam which 34 reflects on HR will have its order increased by RonHR, etc. 35 36 All the optics which transmit or reflect beams (beam splitters, mirrors, 37 thin and thick lenses and special optics) inherit from this class. A 38 particular type of optic is characterized by its action integers and by the 39 inputs provided to the constuctors by the users. Everything else of the 40 optics follow the shape of this Optic class. 41 42 43 *=== Attributes ===* 44 SetupCount (inherited): class attribute, counts all setup components. 45 [integer] 46 OptCount: class attribute, counts optical components. [integer] 47 Name: class attribute. [string] 48 HRCenter (inherited): center of the 'chord' of the HR surface. [3D vector] 49 ARCenter (inherited): center of the 'chord' of the AR surface. [3D vector] 50 HRNorm (inherited): unitary normal to the 'chord' of the HR (always pointing 51 towards the outside of the component). [3D vector] 52 Thick (inherited): thickness of the optic, counted in opposite direction to 53 HRNorm. [float] 54 Dia (inherited): diameter of the component. [float] 55 Ref (inherited): reference string (for keeping track with the lab). [string] 56 ARNorm: unitary normal to the 'chord' of the AR (always pointing 57 towards the outside of the component). [3D vector] 58 N: refraction index of the material. [float] 59 HRK, ARK: curvature of the HR, AR surfaces. [float] 60 HRr, HRt, ARr, ARt: power reflectance and transmission coefficients of 61 the HR and AR surfaces. [float] 62 KeepI: whether of not to keep data of rays for interference calculations 63 on the HR. [boolean] 64 Wedge: wedge angle on the optic. [float] 65 Alpha: angle of the rotation to describe the orientation of the wedge. 66 See the documentation for details on this angle. [float] 67 TonHR, RonHR, TonAR, RonAR: amount by which the orders of the beams will 68 be increased upon relfection or transmission on AR or HR surfaces. 69 These are the principal parameters which distinguish mirrors and lenses 70 and beamsplitters, etc. 71 72 **Note**: the curvature of any surface is positive for a concave surface 73 (coating inside the sphere). 74 Thus kurv*HRNorm/|kurv| always points to the center 75 of the sphere of the surface, as is the convention for the lineSurfInter of 76 geometry module. Same for AR. 77 78 ******* HRK > 0 and ARK > 0 ******* HRK > 0 and ARK < 0 79 ***** ******** and |ARK| > |HRK| 80 H***A H*********A 81 ***** ******** 82 ******* ******* 83 84 ''' 85 86 OptCount = 0 #counts the setup components 87 Name = "Optic" 88

89 - def __init__(self, ARCenter, HRCenter, HRNorm, ARNorm, N, HRK, ARK, ARr, 90 ARt, HRr, HRt, KeepI, Thickness, Diameter, Wedge, Alpha, 91 RonHR, TonHR, RonAR, TonAR, Ref):

92 '''Optic base initializer. 93 94 Parameters are the attributes of the object to construct. 95 96 Returns an Optic. 97 98 ''' 99 # allow empty initializer 100 if Ref is None: 101 Ref = "Opt" + str(Optic.OptCount) 102 103 # initialize with input data 104 self.ARNorm = ARNorm 105 self.N = N 106 self.HRK = HRK 107 self.ARK = ARK 108 self.HRr = HRr 109 self.HRt = HRt 110 self.ARr = ARr 111 self.ARt = ARt 112 self.KeepI = KeepI 113 self.Wedge = Wedge 114 self.Alpha = Alpha 115 self.RonHR = RonHR 116 self.TonHR = TonHR 117 self.RonAR = RonAR 118 self.TonAR = TonAR 119 120 #call mother initializer 121 super(Optic, self).__init__(HRCenter = HRCenter, HRNorm = HRNorm, 122 Ref = Ref, Thickness = Thickness, 123 Diameter = Diameter, ARCenter = ARCenter) 124 Optic.OptCount = Optic.OptCount + 1

125

126 - def isHitDics(self, beam):

127 '''Determine the dictionaries to evaluate if a beam hits the optic. 128 129 Uses the line***Inter functions from the geometry module to calculate 130 the dictionaries of data on interaction with HR and AR and side of 131 optics. 132 133 Returns a tuple of 3 dictionaries with keys: 134 'isHit': whether the optic is hit by the beam. [bool] 135 'intersection point': 3D point where the beam impacts. 136 [3D np-array] 137 'distance': distance from beam origin to interaction point. [float] 138 ''' 139 140 # get impact parameters on HR, AR and side: 141 if np.abs(self.HRK) > 0.: 142 HRDict = geometry.lineSurfInter(beam.Pos, 143 beam.Dir, self.HRCenter, 144 self.HRK*self.HRNorm/np.abs(self.HRK), 145 np.abs(self.HRK), 146 self.Dia) 147 else: 148 HRDict = geometry.linePlaneInter(beam.Pos, beam.Dir, self.HRCenter, 149 self.HRNorm, self.Dia) 150 151 if np.abs(self.ARK) > 0.: 152 ARDict = geometry.lineSurfInter(beam.Pos, 153 beam.Dir, self.ARCenter, 154 self.ARK*self.ARNorm/np.abs(self.ARK), 155 np.abs(self.ARK), 156 self.Dia/np.cos(self.Wedge)) 157 else: 158 ARDict = geometry.linePlaneInter(beam.Pos, beam.Dir, self.ARCenter, 159 self.ARNorm, 160 self.Dia/np.cos(self.Wedge)) 161 162 SideDict = geometry.lineCylInter(beam.Pos, beam.Dir, 163 self.HRCenter, self.HRNorm, 164 self.Thick, self.Dia) 165 166 return HRDict, ARDict, SideDict

167

168 - def isHit(self, beam):

169 '''Determine if a beam hits the Optic. 170 171 This is a function uses the dictionaries provided by isHitDics to 172 find the closest face hit by the beam. 173 174 beam: incoming beam. [GaussianBeam] 175 176 Returns a dictionary with keys: 177 'isHit': whether the beam hits the optic. [boolean] 178 'intersection point': point in space where it is first hit. 179 [3D vector] 180 'face': to indicate which face is first hit, can be 'HR', 'AR' or 181 'Side'. [string] 182 'distance': geometrical distance from beam origin to impact. [float] 183 184 ''' 185 186 # Get isHit dictionaries 187 HRDict, ARDict, SideDict = self.isHitDics(beam) 188 189 # face tags 190 HRDict['face'] = 'HR' 191 ARDict['face'] = 'AR' 192 SideDict['face'] = 'Side' 193 194 # determine first hit 195 hitFaces = filter(lambda dic: dic['isHit'], [HRDict, ARDict, SideDict]) 196 197 if len(hitFaces) == 0: 198 return {'isHit': False, 199 'intersection point': np.array([0., 0., 0.], 200 dtype=np.float64), 201 'face': None, 202 'distance': 0.} 203 204 ans = hitFaces[0] 205 for dic in hitFaces: 206 if dic['distance'] < ans['distance']: 207 ans = dic 208 209 return {'isHit': True, 210 'intersection point': ans['intersection point'], 211 'face': ans['face'], 212 'distance': ans['distance'] 213 }

214

215 - def hit(self, beam, order, threshold):

216 '''Compute the refracted and reflected beams after interaction. 217 218 The beams returned are those selected after the order and threshold 219 criterion. 220 221 beam: incident beam. [GaussianBeam] 222 order: maximum strayness of daughter beams, which are not returned if 223 their strayness is over this order. [integer] 224 threshold: idem for the power of the daughter beams. [float] 225 226 Returns a dictionary of beams with keys: 227 't': refracted beam. [GaussianBeam] 228 'r': reflected beam. [GaussianBeam] 229 230 ''' 231 # get impact parameters and update beam 232 dic = self.isHit(beam) 233 beam.Length = dic['distance'] 234 beam.OptDist = beam.N * beam.Length 235 beam.TargetOptic = self.Ref 236 beam.TargetFace = dic['face'] 237 endSize = beam.width(beam.Length) 238 beam.TWx = endSize[0] 239 beam.TWy = endSize[1] 240 241 return { 242 'HR': lambda beam: self.hitHR(beam, 243 dic['intersection point'], order, threshold), 244 'AR': lambda beam: self.hitAR(beam, 245 dic['intersection point'], order, threshold), 246 'Side': lambda beam: self.hitSide(beam) 247 }[dic['face']](beam)

248

249 - def hitHR(self, beam, point, order, threshold):

250 '''Compute the daughter beams after interaction on HR at point. 251 252 beam: incident beam. [GaussianBeam] 253 point: point in space of interaction. [3D vector] 254 order: maximum strayness of daughter beams, which are not returned if 255 their strayness is over this order. [integer] 256 threshold: idem for the power of the daughter beams. [float] 257 258 Returns a dictionary of beams with keys: 259 't': refracted beam. [GaussianBeam] 260 'r': reflected beam. [GaussianBeam] 261 262 ''' 263 ans = dict() 264 265 # Reference to beam acccording to settings 266 BRef = beam.Ref if not settings.short else shortRef(beam.Ref) 267 268 # Calculate the local normal in opposite direction 269 if self.HRK == 0.: 270 localNorm = self.HRNorm 271 else:#undertending angle 272 theta = np.arcsin(self.Dia * self.HRK/2.)\ 273 if np.abs(self.Dia * self.HRK/2.) < 1. else np.pi/2. 274 275 sphereC = self.HRCenter + np.cos(theta)*self.HRNorm/self.HRK 276 localNorm = (sphereC - point)/np.linalg.norm(sphereC - point) 277 278 if np.dot(beam.Dir, localNorm) > 0.: 279 localNorm = - localNorm 280 K = np.abs(self.HRK) 281 else: 282 K = -np.abs(self.HRK) 283 284 # determine whether we're entering or exiting the substrate 285 if np.dot(beam.Dir, self.HRNorm) < 0.: 286 #entering 287 n1 = beam.N 288 n2 = self.N 289 else: 290 #exiting 291 n1 = self.N 292 n2 = 1. 293 294 # daughter directions 295 dir2 = geometry.newDir(beam.Dir, localNorm, n1, n2) 296 297 #warn on total reflection 298 if dir2['TR'] and settings.info \ 299 and (beam.N == 1. or not settings.short): 300 print "theia: Info: Total reflection of %s on HR of %s." \ 301 % (BRef, self.Ref) 302 303 # if there is no refracted 304 if beam.P * self.HRt < threshold\ 305 or beam.StrayOrder + self.TonHR > order or dir2['t'] is None: 306 ans['t'] = None 307 308 # if there is no reflected 309 if beam.P * self.HRr < threshold\ 310 or beam.StrayOrder + self.RonHR > order: 311 ans['r'] = None 312 313 # we're done if there are two Nones 314 if len(ans) == 2: 315 #warn if its due to power 316 if settings.info\ 317 and (beam.P * self.HRt < threshold\ 318 or beam.P * self.HRr < threshold)\ 319 and (beam.N == 1. or not settings.short): 320 print ("theia: Info: Reached power threshold by interaction"\ 321 +" (%s on %s, HR).") %(BRef, self.Ref) 322 return ans 323 324 # Calculate new basis 325 if not 'r' in ans: # for reflected 326 Uxr, Uyr = geometry.basis(dir2['r']) 327 Uzr = dir2['r'] 328 329 if not 't' in ans: # for refracted 330 Uxt, Uyt = geometry.basis(dir2['t']) 331 Uzt = dir2['t'] 332 333 Lx, Ly = geometry.basis(localNorm) 334 335 # Calculate daughter curv tensors 336 d = np.linalg.norm(point - beam.Pos) 337 C = np.array([[K, 0.], [0, K]]) 338 Ki = np.array([[np.dot(beam.U[0], Lx), np.dot(beam.U[0], Ly)], 339 [np.dot(beam.U[1], Lx), np.dot(beam.U[1], Ly)]]) 340 Qi = beam.Q(d) 341 Kit = np.transpose(Ki) 342 Xi = np.dot(np.dot(Kit, Qi), Ki) 343 344 if not 't' in ans: 345 Kt = np.array([[np.dot(Uxt, Lx), np.dot(Uxt, Ly)], 346 [np.dot(Uyt, Lx), np.dot(Uyt, Ly)]]) 347 Ktt = np.transpose(Kt) 348 Ktinv = np.linalg.inv(Kt) 349 Kttinv = np.linalg.inv(Ktt) 350 Xt = (np.dot(localNorm, beam.Dir) -n2*np.dot(localNorm, Uzt)/n1) * C 351 Qt = n1*np.dot(np.dot(Kttinv, Xi - Xt), Ktinv)/n2 352 353 if not 'r' in ans: 354 Kr = np.array([[np.dot(Uxr, Lx), np.dot(Uxr, Ly)], 355 [np.dot(Uyr, Lx), np.dot(Uyr, Ly)]]) 356 Krt = np.transpose(Kr) 357 Krinv = np.linalg.inv(Kr) 358 Krtinv = np.linalg.inv(Krt) 359 Xr = (np.dot(localNorm, beam.Dir) - np.dot(localNorm, Uzr)) * C 360 Qr = np.dot(np.dot(Krtinv, Xi - Xr), Krinv) 361 362 # Create new beams 363 if not 'r' in ans: 364 ans['r'] = GaussianBeam(Q = Qr, 365 Pos = point, Dir = Uzr, Ux = Uxr, Uy = Uyr, 366 N = n1, Wl = beam.Wl, P = beam.P * self.HRr, 367 StrayOrder = beam.StrayOrder + self.RonHR, 368 Ref = beam.Ref + 'r', 369 Optic = self.Ref, Face = 'HR', 370 Length = 0., OptDist = 0.) 371 372 if not 't' in ans: 373 ans['t'] = GaussianBeam(Q = Qt, Pos = point, 374 Dir = Uzt, Ux = Uxt, Uy = Uyt, N = n2, Wl = beam.Wl, 375 P = beam.P * self.HRt, 376 StrayOrder = beam.StrayOrder + self.TonHR, 377 Ref = beam.Ref + 't', Optic = self.Ref, Face = 'HR', 378 Length = 0., OptDist = 0.) 379 380 return ans

381

382 - def hitAR(self, beam, point, order, threshold):

383 '''Compute the daughter beams after interaction on AR at point. 384 385 beam: incident beam. [GaussianBeam] 386 point: point in space of interaction. [3D vector] 387 order: maximum strayness of daughter beams, which are not returned if 388 their strayness is over this order. [integer] 389 threshold: idem for the power of the daughter beams. [float] 390 391 Returns a dictionary of beams with keys: 392 't': refracted beam. [GaussianBeam] 393 'r': reflected beam. [GaussianBeam] 394 395 ''' 396 397 ans = dict() 398 399 # Reference according to settings 400 BRef = beam.Ref if not settings.short else shortRef(beam.Ref) 401 402 # Calculate the local normal 403 if self.ARK == 0.: 404 localNorm = self.ARNorm 405 else: #undertending angle 406 theta = np.arcsin(self.Dia * self.ARK/2.)\ 407 if np.abs(self.Dia * self.ARK/2.) < 1. else pi/2. 408 409 sphereC = self.ARCenter + np.cos(theta)*self.ARNorm/self.ARK 410 localNorm = (sphereC - point)/np.linalg.norm(sphereC - point) 411 412 if np.dot(beam.Dir, localNorm) > 0.: 413 localNorm = - localNorm 414 K = np.abs(self.ARK) 415 else: 416 K = - np.abs(self.ARK) 417 # determine whether we're entering or exiting the substrate 418 if np.dot(beam.Dir, self.ARNorm) < 0.: 419 #entering 420 n1 = beam.N 421 n2 = self.N 422 else: 423 #exiting 424 n1 = self.N 425 n2 = 1. 426 427 # daughter directions 428 dir2 = geometry.newDir(beam.Dir, localNorm, n1, n2) 429 430 #warn on total reflection 431 if dir2['TR'] and settings.info\ 432 and (beam.N == 1. or not settings.short): 433 print "theia: Info: Total reflection of %s on AR of %s." \ 434 % (BRef, self.Ref) 435 436 # if there is no refracted 437 if beam.P * self.ARt < threshold\ 438 or beam.StrayOrder + self.TonAR > order or dir2['t'] is None: 439 ans['t'] = None 440 441 # if there is no reflected 442 if beam.P * self.ARr < threshold\ 443 or beam.StrayOrder + self.RonAR > order: 444 ans['r'] = None 445 446 # we're done if there are two Nones 447 if len(ans) == 2: 448 if settings.info and (beam.N == 1. or not settings.short): 449 print ("theia: Info: Reached leaf of tree by interaction"\ 450 +" (%s on %s, HR).") % (BRef, self.Ref) 451 return ans 452 453 # Calculate new basis 454 if not 'r' in ans: # for reflected 455 Uxr, Uyr = geometry.basis(dir2['r']) 456 Uzr = dir2['r'] 457 458 if not 't' in ans: # for refracted 459 Uxt, Uyt = geometry.basis(dir2['t']) 460 Uzt = dir2['t'] 461 462 Lx, Ly = geometry.basis(localNorm) 463 464 # Calculate daughter curv tensors 465 d = np.linalg.norm(point - beam.Pos) 466 C = np.array([[K, 0.], [0, K]]) 467 Ki = np.array([[np.dot(beam.U[0], Lx), np.dot(beam.U[0], Ly)], 468 [np.dot(beam.U[1], Lx), np.dot(beam.U[1], Ly)]]) 469 Qi = beam.Q(d) 470 Kit = np.transpose(Ki) 471 Xi = np.dot(np.dot(Kit, Qi), Ki) 472 473 if not 't' in ans: 474 Kt = np.array([[np.dot(Uxt, Lx), np.dot(Uxt, Ly)], 475 [np.dot(Uyt, Lx), np.dot(Uyt, Ly)]]) 476 Ktt = np.transpose(Kt) 477 Ktinv = np.linalg.inv(Kt) 478 Kttinv = np.linalg.inv(Ktt) 479 Xt = (np.dot(localNorm, beam.Dir) -n2*np.dot(localNorm, Uzt)/n1) * C 480 Qt = n1*np.dot(np.dot(Kttinv, Xi - Xt), Ktinv)/n2 481 482 if not 'r' in ans: 483 Kr = np.array([[np.dot(Uxr, Lx), np.dot(Uxr, Ly)], 484 [np.dot(Uyr, Lx), np.dot(Uyr, Ly)]]) 485 Krt = np.transpose(Kr) 486 Krinv = np.linalg.inv(Kr) 487 Krtinv = np.linalg.inv(Krt) 488 Xr = (np.dot(localNorm, beam.Dir) - np.dot(localNorm, Uzr)) * C 489 Qr = np.dot(np.dot(Krtinv, Xi - Xr), Krinv) 490 491 # Create new beams 492 if not 'r' in ans: 493 ans['r'] = GaussianBeam(Q = Qr, 494 Pos = point, Dir = Uzr, Ux = Uxr, Uy = Uyr, 495 N = n1, Wl = beam.Wl, P = beam.P * self.ARr, 496 StrayOrder = beam.StrayOrder + self.RonAR, 497 Ref = beam.Ref + 'r', 498 Optic = self.Ref, Face = 'AR', 499 Length = 0., OptDist = 0.) 500 501 if not 't' in ans: 502 ans['t'] = GaussianBeam(Q = Qt, Pos = point, 503 Dir = Uzt, Ux = Uxt, Uy = Uyt, N = n2, Wl = beam.Wl, 504 P = beam.P * self.ARt, 505 StrayOrder = beam.StrayOrder + self.TonAR, 506 Ref = beam.Ref + 't', Optic = self.Ref, Face = 'AR', 507 Length = 0., OptDist = 0.) 508 509 return ans

510

511 - def hitSide(self, beam):

512 '''Compute the daughter beams after interaction on Side at point. 513 514 Generic function: all sides stop beams. 515 516 beam: incident beam. [GaussianBeam] 517 518 Returns {'t': None, 'r': None} 519 520 ''' 521 if settings.info and (beam.N == 1. or not settings.short): 522 print ("theia: Info: Reached leaf of tree by interaction "\ 523 +"(%s on %s, Side).") %(beam.Ref\ 524 if not settings.short else shortRef(beam.Ref), self.Ref) 525 return {'t': None, 'r': None}

526

527 - def apexes(self):

528 '''Returns the positions of the apexes of HR and AR as a tuple.''' 529 if self.HRK == 0.: 530 apex1 = self.HRCenter 531 else: 532 theta1 = np.arcsin(self.Dia * self.ARK/2.)\ 533 if np.abs(self.Dia * self.ARK/2.) < 1. else np.pi/2. 534 apex1 = self.HRCenter - (1-np.cos(theta1))*self.HRNorm/self.HRK 535 536 if self.ARK == 0.: 537 apex2 = self.ARCenter 538 else: 539 theta2 = np.arcsin(self.Dia * self.ARK/(2.*np.cos(self.Wedge)))\ 540 if np.abs(self.Dia * self.ARK/(2.*np.cos(self.Wedge))) < 1.\ 541 else np.pi/2. 542 apex2 = self.ARCenter - (1-np.cos(theta2))*self.ARNorm/self.ARK 543 544 return apex1, apex2

545

546 - def collision(self):

547 '''Determine whether the HR and AR surfaces intersect. 548 549 Returns True if there is an intersection, False if not. 550 551 ''' 552 553 if self.ARK <= 0. and self.HRK <= 0.: 554 # no collision if negative curvatures 555 return False 556 557 apex = self.apexes() 558 559 #vector from apex1 to apex2 560 vec = apex[1] - apex[0] 561 562 return np.dot(vec, self.HRNorm) > 0.

563

564 - def geoCheck(self, word):

565 '''Makes geometrical checks on surfaces and warns when necessary.''' 566 567 if self.HRt + self.HRr > 1.: 568 print "theia: Warning: In %s %s on HR, R + T > 1." %(word, self.Ref) 569 570 if self.ARt + self.ARr > 1.: 571 print "theia: Warning: In %s %s on AR, R + T > 1." %(word, self.Ref) 572 573 if self.N < 1.: 574 print "theia: Warning: In %s %s, optical index < 1."\ 575 % (word, self.Ref) 576 577 if self.HRK != 0. and np.abs(1./self.HRK) < self.Dia/2.: 578 print ("theia: Warning: In %s, the diameter of the %s exceeds the"\ 579 + " diameter of the HR surface.") % (self.Ref, word) 580 581 if self.ARK != 0. and np.abs(1./self.ARK) < self.Dia/2.: 582 print ("theia: Warning: In %s, the diameter of the %s exceeds the"\ 583 + " diameter of the AR surface.") % (self.Ref, word) 584 585 if self.collision(): 586 print "theia: Warning: In %s, HR and AR surfaces intersect."\ 587 %self.Ref

588

589 - def translate(self, X = 0., Y = 0., Z = 0.):

590 '''Move the optic to (current position + (X, Y, Z)). 591 592 This version takes care of HRcenter and ARCenter and overwrites the 593 SetupComponent version. 594 595 X, Y, Z: components of the translation vector. 596 597 No return value. 598 ''' 599 self.HRCenter = self.HRCenter + np.array([X, Y, Z], dtype = np.float64) 600 self.ARCenter = self.ARCenter + np.array([X, Y, Z], dtype = np.float64)

Source Code for Module theia.optics.optic