00001 """Compute the properties of a diffracting TEM00 Gaussian beam as it passes through an optical system.
00002 This provides utilities to enable the user to lay out a complete optical bench, with mirrors, lenses,
00003 diffraction gratings, etc., and run a laser beam through it.
00004 It correctly handles off-axis optics of most types (tilted lenses & mirrors, e.g.).
00005 It has been used to model a 10 Joule Nd:Glass CPA system at Vanderbilt University, for example
00006 """
00007 _rcsid="$Id: general_optics.py,v 1.19 2007/08/30 14:44:01 mendenhall Exp $"
00008
00009 from math import *
00010 import math
00011 try:
00012 import numpy as Numeric
00013 import numpy.linalg
00014 numeric_float=Numeric.float
00015 numeric_complex=Numeric.complex
00016 eigenvectors=numpy.linalg.eig
00017
00018 except:
00019 import Numeric
00020 import LinearAlgebra
00021 numeric_float=Numeric.Float
00022 numeric_complex=Numeric.Complex
00023 eigenvectors=LinearAlgebra.eigenvectors
00024
00025 import types
00026 import traceback
00027
00028 clight=299792458.
00029 deg=pi/180.0
00030
00031 Infinity="Infinity"
00032
00033 ambient_index=1.0003
00034
00035 def get_ambient_index(lambda0):
00036 if type(ambient_index) is types.FunctionType:
00037 return ambient_index(lambda0)
00038 else:
00039 return ambient_index
00040
00041 def reset():
00042 global ambient_index
00043 ambient_index=1.0003
00044
00045 def vec_mag(vec):
00046 return math.sqrt(Numeric.sum(vec*vec))
00047
00048 def planesolve(x0, u, x1, v):
00049 "planesolve(x0, u, x1, v) returns a, x such that x_=a*u_+x0_ solves (x_-x1_).v_=0"
00050 if vec_mag(x1-x0) < 1e-10:
00051 return 0, x1
00052 a=Numeric.dot((x1-x0), v) / Numeric.dot(u,v)
00053
00054 return a, a*u+x0
00055
00056 def cross(a,b):
00057 x1,y1,z1=tuple(a)
00058 x2,y2,z2=tuple(b)
00059 return Numeric.array((y1*z2-z1*y2, z1*x2-x1*z2, x1*y2-y1*x2))
00060
00061 def wrap_angle(theta):
00062 "wrap_angle always returns an angle in (-180,180)"
00063 return ( (theta + 180) % 360 ) -180
00064
00065 def sincosdeg(theta_deg):
00066 return math.sin(theta_deg*deg), math.cos(theta_deg*deg)
00067
00068 def normalize(x,y):
00069 rho=math.sqrt(x*x+y*y)
00070 if rho < 1e-10:
00071 return 1., 0.
00072 else:
00073 return x/rho, y/rho
00074
00075 def eulermat(theta, eta, phi):
00076 "eulermat(theta, eta, phi) returns r(eta)*r(theta)*r(phi) phi=orientation around normal, theta=yaw, eta=normal roll"
00077 pm=Numeric.identity(3).astype(numeric_float)
00078 em=Numeric.identity(3).astype(numeric_float)
00079 tm=Numeric.identity(3).astype(numeric_float)
00080 sp, cp=sincosdeg(phi)
00081 st, ct=sincosdeg(theta)
00082 se, ce=sincosdeg(eta)
00083 pm[0,0]=pm[1,1]=cp; pm[0,1]=-sp; pm[1,0]=sp
00084 em[0,0]=em[1,1]=ce; em[0,1]=-se; em[1,0]=se
00085 tm[0,0]=tm[2,2]=ct; tm[2,0]=-st; tm[0,2]=st
00086 return Numeric.dot(em, Numeric.dot(tm, pm))
00087
00088
00089 class euler:
00090 def __init__(self, theta, eta, phi):
00091 self.theta=theta
00092 self.eta=eta
00093 self.phi=phi
00094 self.mat=eulermat(theta, eta, phi)
00095
00096 def __str__(self):
00097 s=""
00098 if self.theta != 0:
00099 s=s+"theta = %.1f, " % self.theta
00100 if self.eta != 0:
00101 s=s+"eta = %.1f, " % self.eta
00102 if self.phi != 0:
00103 s=s+"phi = %.1f" % self.phi
00104 if s[-2:]==", ": s=s[:-2]
00105 return s
00106
00107 def general_simil_tens(transform, tens, reversed=0):
00108 "return transform^(-1)*tens*transform or the other way around if reversed=1"
00109 q,r,s,t=transform[0,0], transform[0,1], transform[1,0], transform[1,1]
00110 inverse=Numeric.array(((t,-r),(-s,q))) / (q*t-r*s)
00111 if reversed:
00112 inverse, transform=transform, inverse
00113 return Numeric.dot(inverse, Numeric.dot(tens, transform))
00114
00115 def simil_tens_cs(k, s, tens):
00116 "similarity transform a 2x2 tensor by angle k=cos(theta), s=sin(theta)."
00117 a,b,c,d=tens[0,0], tens[0,1], tens[1,0], tens[1,1]
00118 s2=s*s; k2=k*k; sk=s*k
00119 a1=a*k2+d*s2+sk*(b+c)
00120 b1=b*k2-c*s2+sk*(d-a)
00121 c1=c*k2-b*s2+sk*(d-a)
00122 d1=d*k2+a*s2-sk*(b+c)
00123 return Numeric.array(((a1,b1),(c1,d1)))
00124
00125 def simil_tens(theta, tens):
00126 "similarity transform a 2x2 tensor. If the tensor is deeper than 2nd rank, the inner indices are passed through"
00127 return simil_tens_cs(math.cos(theta), math.sin(theta), tens)
00128
00129 def principal_axis_angle(q22):
00130 a,b,c,d=tuple(Numeric.ravel(q22))
00131 if abs(d-a) < 1e-10: return 0
00132
00133 y=(2*b/(d-a))
00134 qtheta=math.atan(y.real)/2.0
00135 if abs(y.real) > 1e-20 and abs(y.imag)/abs(y.real) > 1e-6:
00136 raise "Bad tensor to diagonalize: y="+str(y)+" theta=: "+str(qtheta/deg)+"\n"+str(q22)
00137
00138 return qtheta
00139
00140 def expand_to_2x2tensor(object):
00141 "expand x to ((x,0),(0,x)), (x,y) to ((x,0),(0,y)), ((x,t),(z,t)) to itself"
00142 a=Numeric.array(object)
00143 if a.shape==():
00144 return a[0]*Numeric.identity(2)
00145 elif a.shape==(2,):
00146 return Numeric.array(((a[0],0),(0,a[1])))
00147 else:
00148 return a
00149
00150
00151 class qtens:
00152 "this is the heart of the general optics package: the 2x2 transverse tensor q parameter, and the book-keeping it needs"
00153
00154 def __init__(self, lambda_initial, q=None, w=None, r=None, name=None, qit=None, medium_index=1.0):
00155 """build a q tensor from
00156 a) another valid inverse q tensor qit,
00157 b) a scalar complex q parameter which creates a diagonal q inverse tensor ( (1/q, 0), (0,1/q) )
00158 c) a vector of two q parameters which creates a diagonal q inverse tensor ( (1/q[0], 0), (0, 1/q[1]) )
00159 d) a radius of curvature r, wavelength lambda_initial, and medium index, which creates a diagonal tensor.
00160 If r is general_optics.Infinity, beam is collimated
00161 """
00162
00163 self.name=name
00164 self.lambda0=lambda_initial*medium_index
00165 self.medium_index=medium_index
00166 if qit is not None:
00167 self.qit=qit
00168 else:
00169 if q is None:
00170 if r is Infinity:
00171 rinv=0
00172 else:
00173 rinv=1.0/r
00174 qi=complex(rinv, -self.lambda0 / (self.medium_index*(math.pi*w**2)) )
00175 self.qit=Numeric.array(((qi,0),(0,qi)))
00176 elif isinstance(q, complex):
00177 self.qit=Numeric.array(((1.0/q,0),(0,1.0/q)))
00178 else:
00179 self.qit=Numeric.array(((1.0/q[0],0),(0,1.0/q[1])))
00180
00181 def q_moments(self):
00182 "return the eigenvectors v and complex q eigenvalues for this q tensor"
00183 u,v=eigenvectors(self.qit)
00184 return v, 1.0/u[0], 1.0/u[1]
00185
00186 def qi_moments(self):
00187 "return the eigenvectors v and complex inverse q eigenvalues for this q tensor"
00188 u,v = eigenvectors(self.qit)
00189 return v, u[0], u[1]
00190
00191
00192 def set_medium_index(self, new_index):
00193 "shift the current q tensor through a perpendicular planar boundary with a different index"
00194 self.qit=self.qit*(self.medium_index/new_index)
00195 self.medium_index=new_index
00196
00197 def rw(self, qi):
00198 "get the radius of curvature r and the spot size w for a given inverse qi assuming our wavelength and index of refraction"
00199 w=math.sqrt(-self.lambda0/qi.imag/math.pi/self.medium_index)
00200 if abs(qi.real) < 1e-15:
00201 r=Infinity
00202 else:
00203 r=1.0/qi.real
00204 return r, w
00205
00206 def next_waist(self):
00207 "find the distance to the next waist on either of our principal axes. If the number is negative, we are past the waist."
00208 theta, qx, qy=self.q_moments()
00209 qxr=qx.real
00210 qyr=qy.real
00211 if qxr*qyr < 0:
00212 return -min((qxr,qyr))
00213 else:
00214 return -max((qxr, qyr))
00215
00216 def rwstr(self, qi):
00217 "format ourself as a nice string given r, w and the distance dz to the next waist"
00218 r,w=self.rw(qi)
00219 if r is not Infinity:
00220 r="%.5g" % r
00221 dz=(1.0/qi).real
00222 return ("r=%s w=%.5g dz=%.4f" % (r,w, dz))
00223
00224 def __str__(self):
00225 theta, q0, q1=self.qi_moments()
00226 if self.name is not None:
00227 s=self.name+": "
00228 else:
00229 s=""
00230
00231 if abs(q0-q1)/(abs(q0)+abs(q1)) < 1e-2:
00232 s=s+"q={"+self.rwstr(q0)+"}"
00233 else:
00234 s=s+"qxx={"+self.rwstr(q0)+"}, qyy={"+self.rwstr(q1)+"}"
00235
00236 return s
00237 def __repr__(self):
00238 return self.__str__()
00239
00240 def qw(self, element):
00241 "for a given optical element which provides a q_transform method, apply it to ourself"
00242 element.q_transform(self)
00243 return self
00244
00245 def abcd_transform(self, abcd):
00246 """for a given abcd matrix, apply it to ourself.
00247 If the matrix elements themselves are vectors, it is an on-axis transform with different x & y matrices.
00248 Such a matrox would look like ( ( (ax, ay), (by, by) ), ( (cx, cy), (dx, dy) ) ).
00249 If the matrix elements themselves are tensors, it is operating off axis.
00250 Such a matrox would look like ( ( ( ( axx, axy ) , ( ayx, ayy )), ( ( bxx, bxy ) , ( byx, byy )) ), ( ... ) )
00251 """
00252 dot=Numeric.dot
00253 ar=Numeric.array
00254 a, b, c, d = [expand_to_2x2tensor(i) for i in ar(abcd).flat]
00255 q,r,s,t=(a+dot(b, self.qit)).flat
00256 inverse=ar(((t,-r),(-s,q))) / (q*t-r*s)
00257 self.qit=dot(c+dot(d, self.qit), inverse)
00258
00259 def drift(self, length):
00260 "advance the beam forward a distance length"
00261 q,r,s,t=(Numeric.identity(2)+length*self.qit).flat
00262 inverse=Numeric.array(((t,-r),(-s,q))) *(1.0/ (q*t-r*s))
00263 self.qit=Numeric.dot(self.qit, inverse)
00264
00265 def focus(self, strength):
00266 "apply a pure focusing strength. It can be a scalar in diopters, a vector of x & y strengths, or a tensor for an off-axis cylindrical lens"
00267 self.qit=self.qit+expand_to_2x2tensor(strength)
00268
00269 def clone(self, newname=None):
00270 "make a complete clone of ourself"
00271 q=copy.copy(self)
00272 if newname is not None:
00273 q.name=newname
00274 return q
00275
00276 def transform(self, tensor):
00277 "transform ourself by tensor . qinv . transpose(conjugate(tensor)). Note that if tensor is not unitary, this isn't either."
00278 tr=Numeric.transpose
00279 dot=Numeric.dot
00280 self.qit=dot(tensor, dot(self.qit, Numeric.conjugate(tr(tensor))))
00281 return self
00282
00283 import exceptions
00284
00285 class OpticDirectionError(exceptions.AssertionError):
00286 "Light pointing the wrong way to hit this optic!"
00287 import copy
00288
00289 class beam:
00290 """ the beam class carries around information for a beam of light: starting point, polarization, wavelength... starting along +z.
00291 It also carries around a dictionary of markers which record snapshots of the beam at interesting pointds along its trajectory
00292 """
00293 def __init__(self, x0, q, lam=None, polarization=(1,0), medium_index=None):
00294
00295 self.x0=Numeric.array(x0,numeric_float)
00296 self.matrix_to_global=Numeric.identity(3).astype(numeric_float)
00297 self.polarization=Numeric.array(polarization, numeric_complex)
00298 if isinstance(q, qtens):
00299 self.q=q.clone()
00300 else:
00301 if medium_index is None:
00302 medium_index=get_ambient_index(lam)
00303 self.q=qtens(lam, q=q, medium_index=medium_index)
00304 self.total_drift=0.0
00305 self.total_drift_time=0.0
00306 self.marks={"order":[]}
00307
00308 def direction(self):
00309 "get the current propagation direction of the beam"
00310 return Numeric.array(self.matrix_to_global[:,2])
00311
00312 def set_medium_index(self, new_index):
00313 "adjust the beam for a new index of refraction"
00314 self.q.set_medium_index(new_index)
00315 return self
00316
00317 def get_medium_index(self):
00318 "get the index of the material in which the beam is currently propagating"
00319 return self.q.medium_index
00320
00321 def set_lambda(self, lambda0):
00322 "set the wavelength of the beam to a new value"
00323 self.q.lambda0=lambda0
00324 return self
00325
00326 def get_lambda(self):
00327 "get the wavelength that has been set"
00328 return self.q.lambda0
00329
00330 def free_drift(self, distance):
00331 "allow the beam to step forward a specified distance, and set some variables which can be recorded in a marker, if desired"
00332 self.q.drift(distance)
00333 self.total_drift=self.total_drift+distance
00334 self.x0=self.x0+self.direction()*distance
00335 self.total_drift_time=self.total_drift_time+distance*self.q.medium_index
00336 self.incoming_q=self.q.clone()
00337 self.footprint_q=self.q.clone()
00338 return self
00339
00340 def localize(self, optic):
00341 """transform the beam into the local coordinate system of an optic, taking care to handle intersecting an optic at any angle and position.
00342 This is typically used by optics to get the q parameter (etc.) of the beam as it appears in the optic's preferred coordinate system
00343 It also records information about the shape of the beam which can be captured by setting a marker.
00344 """
00345 ar=Numeric.array
00346 dot=Numeric.dot
00347 tr=Numeric.transpose
00348
00349 if hasattr(self,"local_x0"):
00350 raise "attempt to localize already localized beam at optic: "+optic.name
00351 self.incoming_direction=self.direction()
00352 cm=dot(optic.matrix_to_local, self.matrix_to_global)
00353 cm2=ar(cm[0:2,0:2])
00354 a,b,c,d=tuple(cm2.flat)
00355 cm2inv=ar(((d,-b),(-c,a)))/(a*d-b*c)
00356 self.incoming_q=self.q.clone()
00357 self.q.transform(cm2)
00358 self.footprint_q=self.q.clone()
00359 self.local_x0=dot(optic.matrix_to_local, self.x0-optic.center)
00360 self.local_direction=ar(cm[:,2])
00361 self.localize_tensor_transform=cm2inv
00362 self.local_polarization=dot(cm2, self.polarization)
00363
00364 def globalize(self, optic):
00365 "undo the localizing transform. Typically called by an optic after it has applied itself to us in its preferred coordinate system"
00366 ar=Numeric.array
00367 dot=Numeric.dot
00368 tr=Numeric.transpose
00369
00370 cm=dot(optic.matrix_to_local, self.matrix_to_global)
00371 cm2=ar(cm[0:2,0:2])
00372 a,b,c,d=tuple(cm2.flat)
00373 cm2inv=ar(((d,-b),(-c,a)))/(a*d-b*c)
00374 self.q.transform(cm2inv)
00375 self.polarization=dot(cm2inv, self.local_polarization)
00376
00377 del self.local_x0, self.local_direction, self.localize_tensor_transform, self.local_polarization
00378
00379 def transform(self, matrix):
00380 "transform the beam into a new coordinate system"
00381 ar=Numeric.array
00382 dot=Numeric.dot
00383 mat=dot(matrix, self.matrix_to_global)
00384 matl=[x/vec_mag(x) for x in mat]
00385 self.matrix_to_global=Numeric.array(matl)
00386
00387 def update_q(self, optic):
00388 "update our q parameter for a given optic. See qtens.qw for more info"
00389 self.q.qw(optic)
00390
00391 def __str__(self):
00392 return "beam: x0=%s, dir=%s, %s" % (Numeric.array_str(self.x0, precision=3, suppress_small=1),
00393 Numeric.array_str(self.direction(), precision=3, suppress_small=1),
00394 str(self.q))
00395
00396 def __repr__(self):
00397 return self.__str__()
00398
00399 def clone(self):
00400 "make a complete clone of the beam"
00401 b=copy.copy(self)
00402 b.q=copy.copy(self.q)
00403 return b
00404
00405 def clone_no_marks(self):
00406 "make a clone of the beam without any marks. This is used by the marks system to embed non-recursive snapshots of the beam into the marks array"
00407 marks=self.marks
00408 del self.marks
00409 q=self.clone()
00410 self.marks=marks
00411 return q
00412
00413 def shift_lambda(self, delta):
00414 """reset the beam to a wavelength shifted by delta from the current wavelength.
00415 Typically used on a cloned beam to prepare a set of beams with different wavelengths but otherwise identical for propagation
00416 """
00417 self.set_lambda(self.get_lambda()+delta)
00418 return self
00419
00420 def x(self, optic):
00421 """apply an optic's tranform method to us, and return self for daisy chaining. For example,
00422 b.x(mylens).x(mygrating).x(mydrift) etc. causes the beam to be passed theough mylens, mygrating, and mydrift in that order.
00423 """
00424 optic.transform(self)
00425 return self
00426
00427 def mark(self, label=None):
00428 "record an entry in the marks array as to what we look like here"
00429 if not self.marks.has_key(label):
00430 self.marks[label]=[]
00431
00432 self.marks[label].append(self.clone_no_marks())
00433 self.marks["order"].append((label, len(self.marks[label])-1))
00434
00435 def transformation_matrix_to_table(self, up):
00436 """find a best-effort transformation matrix which describes our current direction relative to the vector 'up'.
00437 This is used to find a matrix which transforms the beam into a coordinate system which looks natural on the optics table.
00438 Returns flag, matrix where flag is 0 if there is no obvious solution, and 1 if there is a pretty good solution.
00439 If our direction is close to perpendicular to 'up', the matrix puts z in our propagation direction, x perpendicular to 'up' and y perpendicular to x & z.
00440 If the beam is close to parallel to 'up', returns the identity matrix.
00441 The result when the flag is true is a nice coordinate system in which x lies in the plane of the table, y lies nearly in the cirection of 'up'
00442 and z is the beam direction.
00443 """
00444 ar=Numeric.array
00445 dot=Numeric.dot
00446 tr=Numeric.transpose
00447
00448 dir=self.matrix_to_global
00449 z=dir[:,2]
00450 x=cross(up,z)
00451 xmag=sqrt(dot(x,x))
00452 if xmag < 0.1 :
00453 return 0, Numeric.identity(3).astype(numeric_float)
00454 else:
00455 x=x/xmag
00456 y=cross(z,x)
00457 direction=tr(ar((x,y,z)))
00458 return 1, dot(tr(direction), dir)
00459
00460 def transform_q_to_table(self, qi, up=(0,1,0)):
00461 "return qix, qiy (our inverse q parameter diagonal components) in a coordinate system relative to the direction 'up', if possible"
00462 dot=Numeric.dot
00463 tr=Numeric.transpose
00464 ok, matrix=self.transformation_matrix_to_table(up)
00465 if not ok :
00466 theta, qix, qiy=qi.qi_moments()
00467 return qix, qiy
00468 else:
00469 xyt=matrix[0:2,0:2]
00470 qiprime=dot(xyt, dot(qi.qit, tr(xyt)))
00471 return qiprime[0,0], qiprime[1,1]
00472
00473 class general_optic(object):
00474 "general_optic a a class which provides coordinate-system and transport support for many 'atomic' optics types such as thin lenses, mirrors, and gratings"
00475 def __init__(self, name, center=(0,0,0), angle=0, **extras):
00476 "set up our name, center position, and angle, and any extra keyword, value pairs. See general_optic.reset_angle for the meaning of the angle"
00477 self.init(name, center, angle, extras)
00478
00479 def init(self, name, center, angle, extras):
00480
00481 self.abcd=Numeric.identity(2).astype(numeric_float)
00482 self.driftlength=0.0
00483
00484
00485 for i in extras.keys():
00486 setattr(self, i, extras[i])
00487
00488 if len(center)==2:
00489 center=(center[0], 0, center[1])
00490
00491 self.center=Numeric.array(center, numeric_float)
00492 self.name=name
00493
00494 self.reset_angle(angle)
00495 self.post_init()
00496
00497 def post_init(self):
00498 """override post_init in subclasses to do computations after the construction, usually based on extra keywords.
00499 This avoids the need for every general_optic to override the default constructor.
00500 """
00501 pass
00502
00503 def add_info(self, **keys):
00504 "copy extra keywrd-value pairs into our attributes"
00505 for i in keys.keys():
00506 setattr(self, i, keys[i])
00507
00508 def __str__(self):
00509 if self.name is not None:
00510 return (self.name+" "+Numeric.array_str(self.center, precision=3, suppress_small=1) +
00511 " "+Numeric.array_str(self.cosines, precision=3, suppress_small=1))
00512 else:
00513 return "Unnamed optic"
00514
00515 def entrance_center(self):
00516 """return the coordiinate at which this optic's front/entrance surface is centered. For 'atomic' optics, it is the center of the whole object'
00517 override this for thick optics, etc.
00518 """
00519 return self.center
00520 def exit_center(self):
00521 """return the coordiinate at which this optic's back/exit surface is centered. For 'atomic' optics, it is the center of the whole object'
00522 override this for thick optics, etc.
00523 """
00524 return self.center
00525
00526 def transform_into_local(self, direction):
00527 "take a vector in global coordinates, and express it in the local, preferred frame of this optic"
00528 return Numeric.dot(self.matrix_to_local, direction)
00529
00530 def transform_into_global(self, local_direction):
00531 """take a vector in the local frame of the optic, and express it globally. For example, transform_into_global((0,0,1))
00532 usually points in the global direction of forward propagation through the optic
00533 """
00534 return Numeric.dot(self.matrix_to_global, local_direction)
00535
00536 def check_hit_optic(self):
00537 "check to make sure local coordinates are within bounds of optic: default always true"
00538 return 1
00539
00540 def localize_tensor(self, tens):
00541 "compute what a tensor looks like in the coordinate system of self.beam "
00542 dot=Numeric.dot
00543 tr=Numeric.transpose
00544 return dot(self.beam.localize_tensor_transform, dot(tens, tr(self.beam.localize_tensor_transform)))
00545
00546 def globalize_transform(self, transform):
00547 """express a transform which is represented in our local frame as a global transform. Used typically to
00548 express a rotation about a specific axis in our preferred frame globally
00549 """
00550 dot=Numeric.dot
00551 return dot(self.matrix_to_global, dot(transform, self.matrix_to_local))
00552
00553 def transform(self, beam, backwards=0):
00554 "transform(beam, backwards) set up local coordinates, calls local_transform on the beam, and resets it to global coordinates. Returns self for chaining."
00555 self.beam=beam
00556 self.backwards=backwards
00557 beam.localize(self)
00558 self.check_hit_optic()
00559 self.local_transform()
00560 beam.globalize(self)
00561 del self.beam, self.backwards
00562 return self
00563
00564 def local_transform(self):
00565 "by default, do abcd_transform on self.beam, assuming small angles"
00566 ar=Numeric.array
00567 dot=Numeric.dot
00568 tr=Numeric.transpose
00569 xp, yp, zp=tuple(self.beam.local_direction)
00570 dx, dy, dz=tuple(self.beam.local_x0)
00571 x1,xp1, y1, yp1 =self.abcd_transform((dx, xp, dy, yp))
00572 sx=-(xp1-xp)
00573 if abs(sx) > 0.25:
00574 raise ValueError("x-angle is too big to by paraxial... %.3f" % sx)
00575
00576 cx=math.sqrt(1-sx*sx)
00577 sy=-(yp1-yp)
00578 if abs(sy) > 0.25:
00579 raise ValueError("y-angle is too big to by paraxial... %.3f" % sy)
00580 cy=math.sqrt(1-sy*sy)
00581 rot=dot(ar(((cx, 0, -sx),(0,1,0),(sx,0,cx))),ar(((1,0,0),(0,cy,-sy),(0,sy,cy))))
00582 self.beam.transform(self.globalize_transform(rot))
00583 self.q_transform()
00584
00585 def intersect(self, from_point, from_direction):
00586 "find the intersection of a beam coming from from_point, going in from_direction, with our center. Raise an exception if beam won't hit center going that way."
00587 a,x = planesolve(from_point, from_direction, self.center, self.cosines)
00588
00589 if a < -1e-9:
00590 raise OpticDirectionError, "Optic found on backwards side of beam: "+str(self)+" incoming (x,y,z)="+\
00591 Numeric.array_str(from_point, precision=3, suppress_small=1)+" cosines="+\
00592 Numeric.array_str(from_direction, precision=3, suppress_small=1) + \
00593 (" distance = %.3f, point=(%.3f %.3f %.3f)" % (a, x[0], x[1], x[2]))
00594 return a, x
00595
00596 def transport_to_here(self, beam):
00597 "transport beam from wherever it is to our center"
00598 distance, x=self.intersect(beam.x0, beam.direction())
00599 beam.free_drift(distance)
00600 return self
00601
00602 def format_geometry(self):
00603 return ("(x,y,z)=(%.3f, %.3f, %.3f) "%tuple(self.center))+str(self.euler)
00604
00605 def reset_angle(self, angle):
00606 """set the angle of the optic. If angle is a scalar, it is an absolute rotation about the y axis.
00607 Otherwise, angle should be a triple Euler angles theta (yaw), eta (roll around absolute z), phi(roll around local z)
00608 """
00609 if type(angle) != types.TupleType:
00610 self.euler=euler(angle, 0, 0)
00611 theta=angle
00612 else:
00613 self.euler=euler(angle[0], angle[1], angle[2])
00614 theta=angle[0]
00615
00616 self.matrix_to_global=self.euler.mat
00617 self.matrix_to_local=Numeric.transpose(self.euler.mat)
00618 self.cosines=self.transform_into_global(Numeric.array((0,0,1)))
00619
00620 def old_rotate_in_local_frame(self, matrix_to_global):
00621 self.matrix_to_global=Numeric.dot(self.matrix_to_global, matrix_to_global)
00622 self.matrix_to_local=Numeric.transpose(self.matrix_to_global)
00623 self.cosines=self.transform_into_global(Numeric.array((0,0,1)))
00624
00625 def rotate_in_local_frame(self, matrix_to_global):
00626 """Apply an incremental change to our global direction, expressed in local coordinates. For example, this makes it easy to
00627 roll an optic around its local z axis or turn it off center around its local y axis, even after it has already been set in its global frame
00628 """
00629 self.update_coordinates(self.center, self.center, self.globalize_transform(matrix_to_global))
00630
00631 def update_coordinates(self, parent_center=None, reference_coordinate=None, matrix_to_global=None):
00632 """do a complete transform of our angular coordinates and positions to a new frame.
00633 Our rotation is adjusted from its current value by matrix_to_global.
00634 Our center is moved to parent_center+rotated value of (self.center-reference_coordinate).
00635 This is used when we are embedded in a composite optic, and the whole composite optic is moved or rotated.
00636 """
00637 if parent_center is None:
00638 parent_center=self.center
00639 if reference_coordinate is None:
00640 reference_coordinate=self.center
00641 if matrix_to_global is not None:
00642 self.matrix_to_global=Numeric.dot(matrix_to_global, self.matrix_to_global)
00643 self.matrix_to_local=Numeric.transpose(self.matrix_to_global)
00644 self.center=parent_center+Numeric.dot(matrix_to_global, self.center-reference_coordinate)
00645 self.cosines=self.transform_into_global(Numeric.array((0,0,1)))
00646 self.euler=None
00647 return self
00648
00649 def polygon(self):
00650 "return a polygonal boundary for the object. By default it is just the bounding box"
00651 try:
00652
00653 face_centered=hasattr(self,"face_centered")
00654
00655 if hasattr(self,"justify"):
00656 justify=self.justify
00657 else:
00658 justify="center"
00659
00660 if justify=="left":
00661 left=0
00662 right=-self.width
00663 elif justify=="right":
00664 left=self.width
00665 right=0
00666 else:
00667 try:
00668
00669 fraction=(float(justify)+1)*0.5
00670 left=self.width*(1-fraction)
00671 right=-self.width*fraction
00672 except:
00673 left=self.width/2.0
00674 right=-left
00675
00676 if hasattr(self,"height"):
00677 height=self.height
00678 else:
00679 height=self.width
00680
00681 top=height/2.0
00682 bottom=-height/2.0
00683
00684 if(face_centered):
00685 front=0
00686 back=self.thickness
00687 else:
00688 front=-self.thickness/2.0
00689 back=-front
00690
00691 baserect=Numeric.array((
00692 (left, right, right, left, left, left, right, right, left, left, right, right),
00693 (top, top, bottom, bottom, top, top, top, bottom, bottom, top, bottom, bottom),
00694 (front, front, front, front, front, back, back, back, back, back, back, front)))
00695 return Numeric.transpose(Numeric.dot(self.matrix_to_global, baserect))+self.center
00696
00697 except:
00698 return None
00699
00700 def polygon_list(self):
00701 "return a list of polygons which might be used to draw this object. In this case, it is a single item."
00702 return [self.polygon()]
00703
00704 def place_between(self, from_obj, to_obj, distance):
00705 "place_between(from, to, distance, set_info) puts the optic between object or coordinate 'from' and 'to' and rotates it to the specified axis."
00706 if isinstance(from_obj, general_optic):
00707 fc=from_obj.exit_center()
00708 fn=from_obj.name
00709 else:
00710 fc=Numeric.array(from_obj)
00711 fn=Numeric.array_str(fc, precision=3, suppress_small=1)
00712 if isinstance(to_obj, general_optic):
00713 tc=to_obj.entrance_center()
00714 tn=to_obj.name
00715 else:
00716 tc=Numeric.array(to_obj)
00717 tn=Numeric.array_str(tc, precision=3, suppress_small=1)
00718
00719 dx1=tc-fc
00720 dx1hat=dx1/vec_mag(dx1)
00721 if distance > 0.0:
00722 self.center=fc+distance*dx1hat
00723 else:
00724 self.center=tc+distance*dx1hat
00725
00726 x,y,z=tuple(dx1.flat)
00727 rho=math.sqrt(x*x+y*y)
00728 theta=math.atan2(rho,z)/deg
00729 eta=math.atan2(y,x)/deg
00730 self.reset_angle((theta, eta, 0))
00731 return self
00732
00733 def set_direction(self, from_obj, to_obj):
00734 "set_direction(from, to, set_info) points a mirror from object or coordinate 'from' to 'to'"
00735 if isinstance(from_obj, general_optic):
00736 fc=from_obj.exit_center()
00737 fn=from_obj.name
00738 else:
00739 fc=Numeric.array(from_obj)
00740 fn=Numeric.array_str(fc, precision=3, suppress_small=1)
00741 if isinstance(to_obj, general_optic):
00742 tc=to_obj.entrance_center()
00743 tn=to_obj.name
00744 else:
00745 tc=Numeric.array(to_obj)
00746 tn=Numeric.array_str(tc, precision=3, suppress_small=1)
00747
00748 self.looking_from_name=fn
00749 self.looking_to_name=tn
00750 self.looking_from_obj=from_obj
00751 self.looking_to_obj=to_obj
00752
00753 dx1=fc-self.center
00754 dx2=tc-self.center
00755 dx1hat=dx1/vec_mag(dx1)
00756 dx2hat=dx2/vec_mag(dx2)
00757 halfway=(dx1hat+dx2hat)/vec_mag(dx1hat+dx2hat)
00758 x,y,z=tuple(halfway.flat)
00759 rho=math.sqrt(x*x+y*y)
00760 theta=math.atan2(rho,z)/deg
00761 eta=math.atan2(y,x)/deg
00762 self.reset_angle((180+theta, eta, 0))
00763 self.perp=Numeric.array((dx1hat[1]*dx2hat[2]-dx1hat[2]*dx2hat[1],
00764 dx1hat[2]*dx2hat[0]-dx1hat[0]*dx2hat[2],
00765 dx1hat[0]*dx2hat[1]-dx1hat[1]*dx2hat[0]))
00766 return self
00767
00768 def rotate_axis(self, axis_theta):
00769 self.rotate_in_local_frame(matrix_to_global=eulermat(0.0, 0.0, axis_theta))
00770 return self
00771
00772 def tilt_off_axis(self, angle):
00773 if type(angle) is not types.TupleType:
00774 theta=angle
00775 eta=0.0
00776 else:
00777 theta=angle[0]
00778 eta=angle[1]
00779 self.rotate_in_local_frame(matrix_to_global=eulermat(theta, eta, 0.0))
00780 return self
00781
00782 def clone(self, newname=None):
00783 a=copy.copy(self)
00784 if newname is not None:
00785 a.name=newname
00786 return a
00787
00788 def format_name(self):
00789 if self.name is None:
00790 return ""
00791 else:
00792 return self.name + ": "
00793
00794 class base_reflector:
00795 "reflector is a class for mirror-like objects... the default is a mirror"
00796 def local_transform(self):
00797 "a mirror just reverse the z-component of the motion of the beam"
00798 self.beam.transform(self.globalize_transform(Numeric.array(((1.0,0,0),(0,1.0,0),(0,0,-1.0)))))
00799
00800 def post_init(self):
00801 self.reflector_info=self.format_geometry()
00802
00803 class base_lens:
00804 """a base_lens handles simple thin lenses with either a scalar or vector (diagonal) focal length or strength specified. To put it off axis,
00805 apply a rotation after the object has been initialized.
00806 """
00807 def post_init(self):
00808 f=self.f
00809 if type(f)==types.TupleType:
00810 f1, f2 = f
00811 if f1 is not None:
00812 d1=-1.0/f1
00813 else:
00814 d1=0.0
00815 if f2 is not None:
00816 d2=-1.0/f2
00817 else:
00818 d2=0.0
00819 else:
00820 if f is not Infinity:
00821 d1=d2=-1.0/f
00822 else:
00823 d1=d2=0.0
00824
00825 self.d1=d1
00826 self.d2=d2
00827 self.strength=Numeric.array(((d1,0),(0,d2)))
00828 self.info_str()
00829
00830 def mat2(self):
00831 "get an abcd matrix, if we look like a spherical lens"
00832 if self.d1==self.d2:
00833 return Numeric.array(((1.0,0.),(-1.0/self.d1,1)))
00834 else:
00835 raise "Attempt to get 2x2 matrix of non-spherical lens"
00836
00837 def mat4(self):
00838 "get a full abcd tensor for this lens, as needed to transform q"
00839 mxx, mxy, myx, myy=[Numeric.array(((1.,0.),(d,1.))) for d in self.strength.flat]
00840 return Numeric.array(((mxx,mxy),(myx,myy)))
00841
00842 def rotate_axis(self, axis_theta):
00843 "rrotate our axis, and update our info string"
00844 general_optic.rotate_axis(self, axis_theta)
00845 self.info_str()
00846 return self
00847
00848 def q_transform(self):
00849 "apply our q-transform to a localized beam"
00850 self.beam.q.focus(self.strength)
00851
00852 def abcd_transform(self, vec):
00853 "apply our abcd matrix to a vector. Used for propagating the center of a beam"
00854 x, xp, y, yp = tuple(vec)
00855
00856 dxx,dxy,dyx,dyy=tuple(self.localize_tensor(self.strength).flat)
00857 return x, xp+dxx*x+dxy*y, y, yp+dyy*y+dyx*x
00858
00859 def info_str(self):
00860 if self.d1==self.d2:
00861 if self.d1 != 0:
00862 self.base_lens_info=("f=%.1f"%(-1.0/self.d1))
00863 else:
00864 self.base_lens_info=("f= Infinity")
00865
00866 else:
00867 if self.d1==0:
00868 f1=""
00869 else:
00870 f1="f1=%.1f " % (-1.0/self.d1)
00871 if self.d2==0:
00872 f2=""
00873 else:
00874 f2="f2=%.1f " % (-1.0/self.d2)
00875 self.base_lens_info=f1+f2+("axis theta=%.1f"%(self.axis_theta))
00876
00877 def __str__(self):
00878 return self.format_name()+self.base_lens_info
00879
00880 class reflector(general_optic, base_reflector):
00881 "reflector is a class for mirror-like objects... the default is a mirror. Like a base_reflector, but with size and drawing info set up"
00882
00883 def post_init(self):
00884 if not hasattr(self,"face_centered"):
00885 self.face_centered=1
00886 if not hasattr(self,"width"):
00887 self.width=0.0254
00888 if not hasattr(self,"height"):
00889 self.height=self.width
00890 if not hasattr(self,"thickness"):
00891 self.thickness=0.01
00892 base_reflector.post_init(self)
00893
00894 def __str__(self):
00895 return self.format_name()+self.reflector_info
00896
00897 def local_transform(self):
00898 base_reflector.local_transform(self)
00899
00900 def q_transform(self):
00901 pass
00902
00903 def set_direction(self, from_obj, to_obj, set_info=1):
00904 "set_direction(from, to, set_info) points the mirror from object or coordinate 'from' to 'to' and updates the info string"
00905 general_optic.set_direction(self, from_obj, to_obj)
00906
00907 if set_info:
00908 self.reflector_info=("looking from %s to %s " %
00909 (self.looking_from_name, self.looking_to_name))+str(self.euler)
00910 return self
00911
00912 class null_optic(general_optic):
00913 "null_optic is drawable, but has no effect on the beam except the beam will be transported to it... useful for markers in an optical system"
00914 def post_init(self):
00915 if not hasattr(self,"width"):
00916 self.width=0.0254
00917 if not hasattr(self,"thickness"):
00918 self.thickness=0.01
00919
00920 def abcd_transform(self, vec):
00921 return vec
00922
00923 def q_transform(self):
00924 pass
00925
00926 class marker_optic(null_optic):
00927 "marker_optic is drawable, but has no effect on the beam at all... useful for markers in an optical system"
00928 def transport_to_here(self, beam):
00929 "transport beam from wherever it is to our center"
00930 return self
00931
00932 class lens(general_optic, base_lens):
00933 "lens is like a base_lens, but has geometry information"
00934 def post_init(self):
00935 base_lens.post_init(self)
00936
00937 if not hasattr(self,'thickness'): self.thickness=0.005
00938 if not hasattr(self,'width'): self.width=0.0254
00939 if not hasattr(self,'height'): self.height=self.width
00940
00941 def __str__(self):
00942 return base_lens.__str__(self)+" "+self.format_geometry()
00943
00944 class dielectric_interface(general_optic, base_lens):
00945 """ warning... as of 9/24/2002, this is not even close to done. Do not use"""
00946
00947 def post_init(self):
00948
00949 self.last_lambda=0.0
00950 if not hasattr(self, "external_index"):
00951 self.external_index=get_ambient_index
00952
00953 def local_transform(self):
00954
00955 dot=Numeric.dot
00956
00957 self.update_strength()
00958
00959 kx, ky, kz=self.beam.local_direction
00960
00961
00962 sin_theta_in=math.sqrt(kx*kx+ky*ky)
00963 if abs(sin_theta_in) < 1e-6: return
00964 cphi, sphi = kx / sin_theta_in, ky/sin_theta_in
00965
00966 n1=self.ambient_index
00967 n2=self.medium_index
00968
00969 if kz < 0.0:
00970 n1, n2 = n2, n1
00971 sin_theta_in=-sin_theta_in
00972
00973 self.final_index=n2
00974
00975 sin_theta_out= n1/n2*sin_theta_in
00976 dtheta=math.asin(sin_theta_out)-math.asin(sin_theta_in)
00977 c,s = math.cos(dtheta), math.sin(dtheta)
00978 phimat= Numeric.array(((cphi,-sphi, 0),( sphi, cphi, 0),(0,0,1.)))
00979 phimati=Numeric.array(((cphi, sphi, 0),(-sphi, cphi, 0),(0,0,1.)))
00980
00981
00982 self.beam.transform(self.globalize_transform(dot(phimat,dot(Numeric.array(((c,0,s),(0,1,0),(-s,0,c))), phimati))))
00983 self.q_transform()
00984
00985
00986 def update_strength(self):
00987 oldlambda=self.beam.get_lambda()
00988 if oldlambda==self.last_lambda: return
00989
00990 if type(self.external_index) is types.FunctionType:
00991 self.ambient_index=self.external_index(oldlambda)
00992 else:
00993 self.ambient_index=self.external_index
00994 n1=self.ambient_index
00995
00996 if type(self.index) is types.FunctionType:
00997 self.medium_index=self.index(oldlambda)
00998 else:
00999 self.medium_index=self.index
01000
01001 n2=self.medium_index
01002
01003 if hasattr(self, "radius"):
01004 raise "Warning...", "Curved dielectric interfaces don't work yet! try again later"
01005 if type(self.radius) is not types.TupleType:
01006 rx=ry=self.radius
01007 else:
01008 rx,ry=self.radius
01009 if rx is not Infinity:
01010 cx=(n1-n2)/(rx*n2)
01011 else:
01012 cx=0.0
01013 if ry is not Infinity:
01014 cy=(n1-n2)/(ry*n2)
01015 else:
01016 cy=0.0
01017 self.d_scalar=n1/n2
01018 self.strength=Numeric.array(((cx,0),(0,cy)))
01019 else:
01020 self.strength=Numeric.array(((0.,0.),(0.,0.)))
01021
01022 self.last_lambda=oldlambda
01023
01024 def q_transform(self):
01025 self.beam.set_medium_index(self.final_index)
01026 kx, ky, kz=self.beam.local_direction
01027 if kz > 0.0:
01028 self.beam.q.focus(self.strength)
01029 else:
01030 self.beam.q.focus(self.strength*(-self.medium_index/self.ambient_index))
01031
01032
01033 class paraxial_spherical_mirror(reflector, base_lens):
01034 """spherical_mirror is a reflector which focuses.
01035 This is a paraxial approximation, and only makes sense for flairly 'flat' mirrors, for which the intersection can be computed by planar geometry.
01036 Since base_lens (q.v.) also handles optics with different x&y strengths, this can represent cylinders and ellipses, too.
01037 """
01038 def post_init(self):
01039 base_reflector.post_init(self)
01040 base_lens.post_init(self)
01041
01042 if not hasattr(self,'thickness'): self.thickness=0.005
01043 if not hasattr(self,'width'): self.width=0.0254
01044 if not hasattr(self,'height'): self.height=0.0254
01045
01046 def local_transform(self):
01047 general_optic.local_transform(self)
01048 base_reflector.local_transform(self)
01049
01050 def q_transform(self):
01051 "apply our q-transform to a localized beam, making sure this is found at the right level by general_optic"
01052 self.beam.q.focus(self.strength)
01053
01054 def __str__(self):
01055 return self.format_name()+self.reflector_info+" "+self.base_lens_info+" "+self.format_geometry()
01056
01057 class spherical_mirror(reflector, base_lens):
01058 """spherical_mirror is a reflector which focuses. This solves the exact intersection of a beam with the spherical surface, and is right even
01059 for tilted mirrors and strong curvature. For these mirrors, the width and height are assumed to be the same, and are important
01060 since the decision about where the beam hits the mirror depends on these. This logic is not fully implemented yet...
01061 This is also a useful base class for other highly curved surfaces, since transform() handles dynamic geometry.
01062 """
01063 def post_init(self):
01064 base_lens.post_init(self)
01065
01066 if not hasattr(self,'thickness'): self.thickness=0.005
01067 if not hasattr(self,'width'): self.width=0.0254
01068 if not hasattr(self,'height'): self.height=0.0254
01069
01070 def __str__(self):
01071 return self.format_name()+self.reflector_info+" "+self.base_lens_info+" "+self.format_geometry()
01072
01073 def intersect(self, from_point, from_direction):
01074 """find the intersection of a beam coming from from_point, going in from_direction, with our center. Raise an exception if beam won't hit center going that way."""
01075
01076
01077
01078
01079 radius=self.f*2.0
01080
01081 centerpos=self.center-self.matrix_to_global[:,2]*radius
01082 from0=from_point-centerpos
01083 abcxyz=Numeric.dot(from0, from_direction)
01084 disc=abcxyz*abcxyz - Numeric.dot(from0, from0) + radius*radius
01085 if disc < 0:
01086 raise OpticDirectionError, ("beam misses sphere: "+str(self)+" incoming (x,y,z)="+
01087 Numeric.array_str(from_point, precision=3, suppress_small=1)+" cosines="+
01088 Numeric.array_str(from_direction, precision=3, suppress_small=1) )
01089
01090 tvals =[ t for t in (-abcxyz - math.sqrt(disc), -abcxyz + math.sqrt(disc)) if t > 1e-9]
01091
01092 if not tvals:
01093 raise OpticDirectionError, ( "Optic found on backwards side of beam: "+str(self)+" incoming (x,y,z)="+
01094 Numeric.array_str(from_point, precision=3, suppress_small=1)+" cosines="+
01095 Numeric.array_str(from_direction, precision=3, suppress_small=1) )
01096
01097 solns = [ (t, from_point+from_direction*t) for t in tvals]
01098
01099 solns=[ (t,xx) for t,xx in solns if Numeric.dot(xx-centerpos, from_direction)*radius > 0]
01100
01101
01102 if not solns:
01103 raise OpticDirectionError, ( "Only interaction with spherical mirror is on wrong side: "+str(self)+" incoming (x,y,z)="+
01104 Numeric.array_str(from_point, precision=3, suppress_small=1)+" cosines="+
01105 Numeric.array_str(from_direction, precision=3, suppress_small=1) )
01106
01107 soln=solns[0]
01108 self.intersection_zhat=(soln[1]-centerpos)/radius
01109 return soln
01110
01111 def transform(self, beam, backwards=0):
01112 "transform(beam, backwards) set up local coordinates, calls local_transform on the beam, and resets it to global coordinates. Returns self for chaining."
01113 self.beam=beam
01114 self.backwards=backwards
01115 ml=self.matrix_to_local
01116 mg=self.matrix_to_global
01117 self.matrix_to_local=self.get_dynamic_matrix_to_local(beam.matrix_to_global)
01118 self.matrix_to_global=Numeric.transpose(self.matrix_to_local)
01119 beam.localize(self)
01120 self.check_hit_optic()
01121 self.local_transform()
01122 beam.globalize(self)
01123 self.matrix_to_local=ml
01124 self.matrix_to_global=mg
01125 del self.beam, self.backwards
01126 return self
01127
01128 def get_dynamic_matrix_to_local(self, matrix_to_global):
01129 """return the transformation to the exact coordinates for the intersection point last found by intersect() and optimizing based on beam's matrix.
01130 In this case, the beam matrix is not used, since our x,y are sensible choices since in the future, they are likely to be principal axes of an ellipsoid"""
01131 zhat=self.intersection_zhat
01132 xhat=cross(self.matrix_to_global[:,1], zhat)
01133 xmag=math.sqrt(Numeric.dot(xhat, xhat))
01134 if xmag < 0.707:
01135 yhat=cross( zhat, self.matrix_to_global[:,0])
01136 ymag=math.sqrt(Numeric.dot(yhat, yhat))
01137 yhat/=ymag
01138 xhat=cross(yhat, zhat)
01139 else:
01140 xhat /= xmag
01141 yhat=cross(zhat, xhat)
01142
01143 return Numeric.array((xhat,yhat,zhat))
01144
01145 def local_transform(self):
01146 self.beam.transform(self.globalize_transform(Numeric.array(((1.0,0,0),(0,1.0,0),(0,0,-1.0)))))
01147 self.beam.q.focus(self.strength)
01148
01149
01150 if 0:
01151 print "\n\n****start mirror tests"
01152 mir1=reflector("reflector1", center=(0,0,1))
01153 mir2=reflector("reflector2", center=(0,1,1))
01154 mir1.set_direction((0,0,0), mir2)
01155 mir2.set_direction(mir1, (1,1,1))
01156 mybeam=beam((0,0,0), qtens(1.054e-6, r=Infinity, w=.002))
01157 mybeam.q.qit[0,0]*=0.25
01158 for optic in (mir1, mir2):
01159 print optic.name
01160 optic.transport_to_here(mybeam)
01161 print mybeam
01162 optic.transform(mybeam)
01163 print mybeam
01164
01165 if 0:
01166 print "\n\n****start lens tests"
01167 optic=lens("reflector", center=(0,0,1), f=0.25)
01168 mybeam=beam((0,0,0), qtens(1.054e-6, r=Infinity, w=.002))
01169 print mybeam
01170 optic.transport_to_here(mybeam)
01171 print mybeam
01172 optic.transform(mybeam)
01173 print mybeam
01174
01175
01176 class grating(reflector):
01177 "grating is a reflector which diffracts. The ruling is along the y axis. It should be initialized with keywords order and pitch."
01178
01179 def local_transform(self):
01180 dot=Numeric.dot
01181 tr=Numeric.transpose
01182 deg=math.pi/180.0
01183
01184
01185 kx, ky, kz = v2 = self.beam.local_direction
01186 v3 = cross((1,0,0), v2)
01187 v3 /= vec_mag(v3)
01188 v1 = cross(v2, v3)
01189 v1 /= vec_mag(v1)
01190 coord=Numeric.array((v1,v2,v3))
01191 cinv=tr(coord)
01192
01193
01194
01195
01196 theta=math.atan2(kx, math.sqrt(ky*ky+kz*kz))
01197 littrow, out=self.angles(theta, self.beam.get_lambda())
01198
01199
01200
01201 dtheta=out+theta
01202 s,c=-math.sin(dtheta), math.cos(dtheta)
01203
01204
01205 rot1=dot(Numeric.array(((1,0,0),(0,1,0),(0,0,-1))), dot(cinv, dot(Numeric.array(((c,s,0),(-s,c,0),(0,0,1))), coord )))
01206
01207 self.beam.transform(self.globalize_transform(rot1))
01208
01209 def degree_angles(self, theta, lam):
01210
01211 litt, beta = self.angles(theta*deg, lam)
01212 return litt/deg, beta/deg
01213
01214 def angles(self, theta, lam):
01215
01216 if theta<0:
01217 sgn=-1
01218 else:
01219 sgn=1
01220 gratparm=self.order*lam*self.pitch*sgn
01221 littrow=asin(gratparm/2.0)
01222 phi=(((theta + math.pi/2.0) % math.pi ) - math.pi/2.0)
01223 beta0=asin(gratparm-sin(phi))
01224 self.parafocus_scale=math.cos(phi)/math.cos(beta0)
01225
01226 return (littrow, beta0)
01227
01228 def __str__(self):
01229 return reflector.__str__(self)+("pitch= %.3f" % self.pitch)
01230
01231 class key_tag:
01232 def __init__(self, key):
01233 self.key=key
01234
01235 class backwards(key_tag):
01236 pass
01237
01238 def get_tagged_key(object):
01239 if isinstance(object, key_tag):
01240 return object.key
01241 else:
01242 return object
01243
01244 class composite_optic(general_optic):
01245 "composite_optic is a container for a list of general_optics which can be rotated, etc. as a group"
01246
01247 def __init__(self, name, optics_dict, optics_order, reference_coordinate=None,
01248 center=(0,0,0), angle=0.0, **extras):
01249 """create a composite_optics. It is composed of a dictionary 'optics_dict' of named objects, and a list 'optics_order' of dictionary keys
01250 which defines the order in which elements of this optic should be traversed. A single element may appear multiple times in the list.
01251 The reference_coordinate is the coordinate around which this optic will rotate, it is the reference_coordinate that is placed at the position 'center'.
01252 If no reference_coordinate is specified, it is ssumed to be the entrance_center of the first optic in the optics_order list.
01253 """
01254 self.init( name, optics_dict, optics_order, reference_coordinate, center, angle, extras)
01255
01256 def init(self, name, optics_dict, optics_order, reference_coordinate, center, angle, extras):
01257
01258 if reference_coordinate is None:
01259 reference_coordinate=optics_dict[optics_order[0]].entrance_center()
01260
01261 self.optics_dict={}
01262 self.optics_order=list(optics_order)
01263
01264 for k in optics_order:
01265 self.optics_dict[get_tagged_key(k)]=optics_dict[get_tagged_key(k)]
01266
01267 general_optic.init(self, name, center=center, angle=angle, extras=extras)
01268 for k in self.optics_dict.keys():
01269 self.optics_dict[k].update_coordinates(center, reference_coordinate, self.matrix_to_global)
01270
01271 def exit_center(self):
01272 "return the position of the exit center of this optic"
01273 return self.optics_dict[self.exit_optics_tags()[1]].center
01274
01275 def entrance_center(self):
01276 "return the position of the entrance center of this optic"
01277 return self.optics_dict[self.entrance_optics_tags()[0]].center
01278
01279 def update_coordinates(self, new_center, parent_reference, matrix_to_global):
01280 "move us to a new positin and/or angle. See general_optic.update_coordinates"
01281 for k in self.optics_dict.keys():
01282 self.optics_dict[k].update_coordinates(new_center, parent_reference, matrix_to_global)
01283 general_optic.update_coordinates(self, new_center, parent_reference, matrix_to_global)
01284
01285 def mark_label(self, opticname):
01286 "provide a default label for automatic marking of the beam as it passes though this optic"
01287 return (self,opticname)
01288
01289 def transform(self, beam, back=0):
01290 "track a beam through this optics, forwards or backwards, and mark the beam whenever it strikes an atomic optic"
01291 if back:
01292 order=copy.copy(self.optics_order)
01293 order.reverse()
01294 else:
01295 order=self.optics_order
01296
01297 for opticname in order:
01298 optic=self.optics_dict[get_tagged_key(opticname)]
01299 if not isinstance(optic, composite_optic):
01300 waist=beam.q.next_waist()
01301 current_drift=beam.total_drift
01302 optic.transport_to_here(beam)
01303 distance=beam.total_drift-current_drift
01304 if waist > 0 and waist < distance:
01305 beam.free_drift(waist-distance)
01306 beam.mark("waist")
01307 optic.transport_to_here(beam)
01308 if isinstance(opticname, backwards):
01309 optic.transform(beam, not back)
01310 else:
01311 optic.transform(beam, back)
01312
01313 if not isinstance(optic, composite_optic):
01314 beam.mark(self.mark_label(get_tagged_key(opticname)))
01315
01316 def polygon_list(self):
01317 "return a list of polygons which will draw this optic"
01318 l=[]
01319 for o in self.optics_dict.keys():
01320
01321 l=l+self.optics_dict[o].polygon_list()
01322 return [poly for poly in l if poly is not None]
01323
01324 def exit_optics_tags(self):
01325 "return the last two elements of this optic, to allow one to get an exit axis direction"
01326 return self.optics_order[-2], self.optics_order[-1]
01327
01328 def entrance_optics_tags(self):
01329 "return the first two elements of this optic, to allow one to get an entrance axis direction"
01330 return self.optics_order[0], self.optics_order[1]
01331
01332 def set_exit_direction(self, look_to):
01333 "make the last optic of this look from the second-to-last optic to a final look_to point. Useful if the last optic is a mirror"
01334 l=self.optics_dict
01335 e1, e2 = self.exit_optics_tags()
01336 l[e2].set_direction(l[e1], look_to)
01337
01338 def set_entrance_direction(self, look_from):
01339 "make the first optic of this look from a look_from point to the second optic. Useful if the first optic is a mirror"
01340 l=self.optics_dict
01341 e1, e2 = self.entrance_optics_tags()
01342 l[e1].set_direction(look_from, l[e2])
01343
01344 def rotate_to_axis(self, from_obj):
01345 "rotate this object by an angle that makes it look at from_obj, assuming it started out looking along the z axis"
01346 if isinstance(from_obj, general_optic):
01347 fc=from_obj.exit_center()
01348 fn=from_obj.name
01349 else:
01350 fc=Numeric.array(from_obj)
01351 fn=Numeric.array_str(fc, precision=3, suppress_small=1)
01352
01353 dx1=self.center-fc
01354 x,y,z=tuple(dx1)
01355 rho=math.sqrt(x*x+y*y)
01356 theta=math.atan2(rho,z)/deg
01357 eta=math.atan2(y,x)/deg
01358 if abs(eta)>90:
01359 eta=180-eta
01360 theta=-theta
01361 self.update_coordinates(self.center, self.center, euler(theta, eta, 0).mat)
01362 return self
01363
01364 def __getitem__(self, tag):
01365 "get an optic from our list by its identifier. If a mark is passed in of the form (self, item), use the item"
01366 if type(tag) is types.TupleType and tag[0] is self:
01367 return self.optics_dict[tag[1]]
01368 else:
01369 return self.optics_dict[tag]
01370
01371 def __len__(self):
01372 return len(self.optics_order)
01373
01374 if 0:
01375 print "testing grating"
01376 optic=grating("grating", pitch=1.5e6, angle=38.0, order=1)
01377 mybeam=beam((0,0,0), qtens(1.054e-6, r=Infinity, w=.002))
01378 print mybeam
01379 optic.transport_to_here(mybeam)
01380 print mybeam
01381 optic.transform(mybeam)
01382 print mybeam
01383 print math.atan2(mybeam.direction()[0], mybeam.direction()[2])/(math.pi/180.0)
01384
01385 class optics_trace:
01386 def __init__(self, marks, color=None, **extras):
01387 self.color=color
01388 self.marks=marks
01389
01390 for i in extras.keys():
01391 setattr(self, i, extras[i])
01392
01393 def __getitem__(self, index):
01394 if type(index) is types.IntType:
01395 tag, n = self.marks["order"][index]
01396 return self.marks[tag][n]
01397 else:
01398 return self.marks[index[0]][index[1]]
01399
01400 def __len__(self):
01401 return len(self.marks["order"])
01402
01403 def trace_path(optics_system, beam):
01404 """tracks a beam through a composite_optic, and returns an optics_trace wrapper around the marks which were created.
01405 If an error is encountered which throws an exception, the exception is printed, and the partial trace returned.
01406 This makes it very easy to find a problem optic (e.g. a grating which is diffracting a beam away from the next item,
01407 or a mirror which isn't steered right).
01408 """
01409 beam.total_drift=0.0
01410 try:
01411 optics_system.transform(beam)
01412 except:
01413 traceback.print_exc()
01414 pass
01415 return optics_trace(beam.marks)
01416
01417 class phase_plate(null_optic):
01418 def post_init(self):
01419 theta=self.phase_angle*math.pi/180.0
01420 self.x_phase_shift=complex(math.cos(theta), math.sin(theta))
01421
01422 def local_transform(self):
01423 self.beam.local_polarization[1]*=self.x_phase_shift
01424
01425 def info_str(self):
01426 return "Phase plate: x shift = %.1f" % self.phase_angle
01427
01428 class halfwave_plate(phase_plate):
01429 def post_init(self):
01430 self.phase_angle=180.0
01431 phase_plate.post_init(self)
01432
01433 class quarterwave_plate(phase_plate):
01434 def post_init(self):
01435 self.phase_angle=90
01436 phase_plate.post_init(self)
01437
01438 class faraday_rotator(null_optic):
01439 def post_init(self):
01440 pass
01441
01442 def info_str(self):
01443 return "Faraday Rotator: rotation = %.1f" % self.rotation
01444
01445 def local_transform(self):
01446 theta=self.rotation*self.beam.local_direction[2]
01447 s,c=sincosdeg(theta)
01448
01449
01450 self.beam.local_polarization=Numeric.dot(Numeric.array(((c,-s),(s,c))), self.beam.local_polarization)
01451
01452 class dielectric_trapezoid(composite_optic):
01453 ENTRANCE='entrance'
01454 CENTER='center'
01455 EXIT='exit'
01456 def __init__(self, name, entrance_face_angle=0.0, exit_face_angle=0.0, thickness=0.05, width=0.025, index=None, center=(0.,0.,0.), angle=0.0, **extras):
01457
01458 if index is None:
01459 raise "No index of refraction (scalar or function) given as index=foo creating dielectric trapezoid"
01460
01461 my=dielectric_trapezoid
01462 self.thickness=thickness
01463 self.width=width
01464
01465 self.front_face_offset=math.tan(entrance_face_angle*deg)*width*0.5
01466 self.back_face_offset=math.tan(exit_face_angle*deg)*width*0.5
01467
01468 optics={}
01469 optics[my.ENTRANCE]=dielectric_interface("entrance face", angle=entrance_face_angle, index=index,
01470 center=(0,0,-0.5*thickness))
01471 optics[my.CENTER]=null_optic("trapezoid center", center=(0,0,0), width=width*0.05, thickness=width*0.05)
01472 optics[my.EXIT]=dielectric_interface("exit face", angle=exit_face_angle+180.0, index=index,
01473 center=(0,0,0.5*thickness) )
01474 order=[my.ENTRANCE, my.CENTER, my.EXIT]
01475
01476
01477 composite_optic.init(self, name, optics, order, (0,0,0), center, angle, extras=extras )
01478
01479 def polygon(self):
01480 if hasattr(self,"height"):
01481 height=self.height
01482 else:
01483 height=self.width
01484
01485 top=height/2.0
01486 bottom=-height/2.0
01487
01488 front=-self.thickness/2.0
01489 back=-front
01490
01491 fo=self.front_face_offset
01492 bo=self.back_face_offset
01493
01494 left=0.5*self.width
01495 right=-left
01496
01497 baserect=Numeric.array((
01498 (left, right, right, left, left, left, right, right, left, left, right, right),
01499 (top, top, bottom, bottom, top, top, top, bottom, bottom, top, bottom, bottom),
01500 (front-fo, front+fo, front+fo, front-fo, front-fo, back-bo, back+bo, back+bo, back-bo, back-bo, back+bo, front+fo)))
01501 return Numeric.transpose(Numeric.dot(self.matrix_to_global, baserect))+self.center
01502
01503 def polygon_list(self):
01504 my=dielectric_trapezoid
01505 return [self.polygon(), self.optics_dict[my.CENTER].polygon()]
01506
01507
01508 if 0:
01509 print "\n\n****start phase plate tests"
01510 optic=halfwave_plate("phase", center=(0,0,1)).rotate_axis(30)
01511 mybeam=beam((0,0,0), qtens(1.054e-6, r=Infinity, w=.002))
01512 print mybeam
01513 optic.transport_to_here(mybeam)
01514 print mybeam, mybeam.polarization
01515 optic.transform(mybeam)
01516 print mybeam, mybeam.polarization