general_optics.py

Go to the documentation of this file.
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. #m/sec
00029 deg=pi/180.0
00030 
00031 Infinity="Infinity"
00032 
00033 ambient_index=1.0003 #air
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         #print "\nplanesolve: ", x0, u, x1, v, a, a*u+x0, "\n"
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         #if we are passed a q-tensor for which 2b/(d-a) isn't real, this is wrong!
00133         y=(2*b/(d-a))
00134         qtheta=math.atan(y.real)/2.0 #principal axis direction of q
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: #we are given a real inverse-q tensor, just use it
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)) #one is negative, use it
00213                 else:
00214                         return -max((qxr, qyr)) #the least negative one is the first waist, or no waist if both positive
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)) # +" qi=("+str(qi)+")"
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 #make daisy-chaining easy
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 #accumulated travel distance
00305                 self.total_drift_time=0.0 #accumulated distance*index
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() #keep a copy of the localized q for markers
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] #make sure unit vectors stay that way!
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) #copy q a little deeper
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 #to make daisy-chains easy
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]=[] #empty list for this optic
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 : #z is not even close to perpendicular to up!
00453                         return 0, Numeric.identity(3).astype(numeric_float) #punt, return false and identity
00454                 else:
00455                         x=x/xmag
00456                         y=cross(z,x) #may be tilted if the beam isn't completely level
00457                         direction=tr(ar((x,y,z)))
00458                         return 1, dot(tr(direction), dir) #coordinates are rational, return true and transform
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 : #z is not even close to perpendicular to up!
00466                         theta, qix, qiy=qi.qi_moments()
00467                         return qix, qiy #so just use principal moments  
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) #predefine this as a default abcd matrix (in 2 dimensions)
00482                 self.driftlength=0.0
00483                 
00484                 #copy any extra information directly into the class dictionary
00485                 for i in extras.keys():
00486                         setattr(self, i, extras[i])
00487 
00488                 if len(center)==2: #2-d optics live in the x-z plane
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 #in case anyone needs to know (e.g. dielectric interfaces)
00557                 beam.localize(self)
00558                 self.check_hit_optic()
00559                 self.local_transform()
00560                 beam.globalize(self)
00561                 del self.beam, self.backwards #get it out of our dictionary so it is an error if not localized
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)  #it better be a small angle, so sin(theta)=theta for this to work
00570                 dx, dy, dz=tuple(self.beam.local_x0) #dz should always be zero in these coordinates
00571                 x1,xp1, y1, yp1 =self.abcd_transform((dx, xp, dy, yp))
00572                 sx=-(xp1-xp) #this is sort of a sine of a rotation angle of x about y, with the right sign
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) #this is sort of a sine of a rotation angle of y about x, with the right sign
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: #convert conventional angle to rotation about y only
00610                         self.euler=euler(angle, 0, 0)
00611                         theta=angle
00612                 else: #this should be a triple angles theta (yaw), eta (roll around absolute z), phi(roll around local z)                       
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                                         # a fraction of 1 top justifies, 0 center justifies, and -1 bottom justifies
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()] #for simple optics, the list has one item
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) #assume is is a tuple or array
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) #assume is is a tuple or array
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 #make daisy-chaining easy!
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) #assume is is a tuple or array
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) #assume is is a tuple or array
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]))  #cross product perpendicular to plane of incidence, for reference
00766                 return self #make daisy-chaining easy
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) #it's easy!
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): #short format string
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 #make daisy-chaining easy
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                 #default lens draws as 1" diameter x 5 mm thick
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                 #raise "Warning!", "dielectric_interface is not ready to use! Do not try this."
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                 #must rotate by the appropriate refraction angle _in the plane of incidence_ 
00962                 sin_theta_in=math.sqrt(kx*kx+ky*ky) #incoming k had better be a unit vector!
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 #we're going backwards through the interface
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                 #this matrix is a rotation by by phi to put the incident beam on x, rotate about x-z by dtheta, then undo the phi
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     #always diagonal, for scalar media, so we don't have to localize        
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: #going the wrong way!
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                 #default lens draws as 1" diameter x 5 mm thick
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) #pick up our abcd transform for the lens
01048                 base_reflector.local_transform(self) #and do the reflection
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) #it's easy!
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                 #default lens draws as 1" diameter x 5 mm thick
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                 # t=-(a,b,c).(x1,y1,z1) +- sqrt( ( (a,b,c).(x1,y1,z1) )^2 - (x1,y1,z1)^2 + r^2 ) from Mathematica for line (a,b,c)*t+(x1,y1,z1) and sphere of radius r at origin
01077                 #the center of curvature of the mirror is derived from the apex position (our location) and the rotation matrix
01078                                 
01079                 radius=self.f*2.0 #radius of curvature
01080                 
01081                 centerpos=self.center-self.matrix_to_global[:,2]*radius #position of center of sphere           
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] #select only forward solutions
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] #pick solution with surface pointing the right way
01100                 #there will always only be one of the solutions which passes this test, at most
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  #unit surface normal in general direction of beam, useful later
01109                 return soln #this is the only possible solution
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 #in case anyone needs to know (e.g. dielectric interfaces)
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) #do optics with exact local matrix
01118                 self.matrix_to_global=Numeric.transpose(self.matrix_to_local) #do optics with exact local matrix
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 #get it out of our dictionary so it is an error if not localized
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) #use our real yhat x our zhat to generate xhat
01133                 xmag=math.sqrt(Numeric.dot(xhat, xhat))
01134                 if xmag < 0.707: #bad guess for basis... mirror normal is very close to our yhat, use our xhat instead (they can't BOTH be bad!)
01135                         yhat=cross( zhat, self.matrix_to_global[:,0]) #use our zhat x our real xhat to generate yhat (swapped to keep handedness)
01136                         ymag=math.sqrt(Numeric.dot(yhat, yhat))
01137                         yhat/=ymag
01138                         xhat=cross(yhat, zhat) #xhat is automatically a unit vector
01139                 else:
01140                         xhat /= xmag
01141                         yhat=cross(zhat, xhat) #yhat is automatically a unit vector
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))))) #reflect along local z axis
01147                 self.beam.q.focus(self.strength) #it's easy!
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                 #need to find basis set for rotation which is perp. to beam and grating rulings (yhat)
01185                 kx, ky, kz = v2 = self.beam.local_direction
01186                 v3 = cross((1,0,0), v2) #this is rotation axis (!) in grating frame 
01187                 v3 /= vec_mag(v3)
01188                 v1 = cross(v2, v3) #this is first basis vector for rotation 
01189                 v1 /= vec_mag(v1)
01190                 coord=Numeric.array((v1,v2,v3)) #matrix to convert coordinates in grating system into diffraction rotation system
01191                 cinv=tr(coord) #matrix to convert diffraction basis vectors to grating basis i.e. cinv.(0,0,1) is rotation axis in grating system
01192                 
01193                 #print "**grating**"
01194                 #print Numeric.array_str(coord, precision=4)
01195                 
01196                 theta=math.atan2(kx, math.sqrt(ky*ky+kz*kz))
01197                 littrow, out=self.angles(theta, self.beam.get_lambda())
01198                 #print littrow/deg, theta/deg, out/deg, (out+theta)/deg
01199                 
01200                 #remember, outgoing angle is _not_ outgoing direction... beam has changed from incoming to outgoing, too
01201                 dtheta=out+theta
01202                 s,c=-math.sin(dtheta), math.cos(dtheta)
01203                 
01204                 #this funny matrix is a rotation by dtheta, followed by a reflection in z
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                 #print Numeric.array_str(rot1, precision=4)
01207                 self.beam.transform(self.globalize_transform(rot1))
01208 
01209         def degree_angles(self, theta, lam):
01210                 #return angle info for incoming and outgoing angles in degrees (for user convenience)
01211                 litt, beta = self.angles(theta*deg, lam)
01212                 return litt/deg, beta/deg
01213                 
01214         def angles(self, theta, lam):
01215                 #return angle info for incoming and outgoing angles in degrees (for user convenience)
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) #always wrap to front of grating!
01223                 beta0=asin(gratparm-sin(phi))
01224                 self.parafocus_scale=math.cos(phi)/math.cos(beta0) #angular magnification along x
01225                 #print "grating angles: ", lam, theta/deg, littrow/deg, beta0/deg
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: #refer all coordinates to first optic in list by default
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: #only pull out keys which are used, and let python eliminate duplicates
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) #update our top-level matrices, etc.
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): #only mark primitive optics
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) #'backwards' must toggle the reversed flag, so backwards-backwards is forwards
01310                         else:
01311                                 optic.transform(beam, back)
01312                                 
01313                         if not isinstance(optic, composite_optic): #only mark primitive optics
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                         #use additions semantics to flatten nesting of lists
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) #assume is is a tuple or array
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: #prefer signed theta and unsigned eta
01359                         eta=180-eta
01360                         theta=-theta
01361                 self.update_coordinates(self.center, self.center, euler(theta, eta, 0).mat)     
01362                 return self #make daisy-chaining easy
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]] #in case the tag had our object identifier attached
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                 #copy any extra information directly into the class dictionary
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 #just return an incomplete trace if something goes wrong!
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                 #if beam is going forward, shift (+), if backwards, shift (-)
01449                 #where forward/backwards comes from the local z-component of the beam direction
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                 #specify reference as (0,0,0) to make our coordinates referenced to (0,0,0)
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

Generated on Wed Nov 21 10:18:31 2007 for analysis by  doxygen 1.5.4