'''
Define optical components for gtrace.

'''
#{{{ Import modules
import numpy as np
pi = np.pi
array = np.array
sqrt = np.lib.scimath.sqrt
from numpy.linalg import norm

from traits.api import HasTraits, Int, Float, CFloat, CArray, List, Str

import gtrace.optics as optics
import gtrace.optics.geometric
from unit import *
import copy
import gtrace.draw as draw

#}}}

#{{{ Author and License Infomation

#Copyright (c) 2011-2012, Yoichi Aso
# All rights reserved.

# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions
# are met:
#
# * Redistributions of source code must retain the above copyright
#   notice, this list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright
#   notice, this list of conditions and the following disclaimer in the
#   documentation and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.

__author__ = "Yoichi Aso"
__copyright__ = "Copyright 2011-2012, Yoichi Aso"
__credits__ = ["Yoichi Aso"]
__license__ = "BSD"
__version__ = "0.2.0"
__maintainer__ = "Yoichi Aso"
__email__ = "aso@granite.phys.s.u-tokyo.ac.jp"
__status__ = "Beta"

#}}}

#{{{ Generic Optics Class

class Optics(HasTraits):
    '''
    A general optics class from which other specific
    optics classes are derived.

    === Attributes ===
    name: Name of the optics. [string]
    center: Center position of the optics. [(2,) float array]
    rotationAngle: This angle defines the orientation of the optics.                   
    '''
    name = Str()
    center = CArray(dtype=np.float64, shape=(2,))
    rotationAngle = CFloat(0.0) #in rad

#{{{ isHit(beam)

    def isHit(beam):
        '''
        A function to see if a beam hits this optics or not.

        **Input parameters**
        
          ``beam``:
          A GaussianBeam object to be interacted by the optics.

        **Returned value**
        
          The return value is a dictionary with the following keys:
          ``isHit, position, distance, face``

          ``isHit``:
          This is a boolean to answer whether the beam hit the optics
          or not.

          ``position``:
          A numpy array containing the coordinate values of the intersection
          point between the beam and the optics. If isHit is False, this parameter
          does not mean anything.

          ``distance``
          The distance between the beam origin and the intersection point.

          ``face``:
          An optional string identifying which face of the optics was hit.
          For example, ``face`` can be either "HR" or "AR" for a mirror.
          ``face`` can also be "side", meaning that the beam hits a side
          of the optics, which is not meant to be used, e.g. the side of a mirror.
          In this case, the beam have reached a dead end.
        '''
        #This is an abstract function
        return {'isHit': False, 'position': np.array((0,0)),
                'distance': 0.0, 'face':''}

#}}}

#{{{ hit(beam, order=0, threshold=0.0):

    def hit(beam, order=0, threshold=0.0):
        '''
        A function to hit the optics with a beam.

        This function attempts to hit the optics with the source beam, ``beam``.

        **Input parameters**:

          ``beam``:

          ``order``:
          An integer to specify how many times the internal reflections
          are computed.

          ``threshold``
          The power threshold for internal reflection calculation.
          If the power of an auxiliary beam falls below this threshold,
          further propagation of this beam will not be performed.

        **Return values**
        ``(isHit, beamDict, face)``

        ``isHit``
        This is a boolean to answer whether the beam hit the optics
        or not.

        ``beamDict``
        A dictionary containing resultant beams.
        
        ``face``:
          For a mirror, ``face`` is any of "HR", "AR" or "side".
          
        '''
        return {False, {}, "side"}

#}}}

#{{{  _isHitSurface_()

    def _isHitSurface_(self, beam, surface_center, normal_vector,
                       surface_size=1.0, inv_ROC=0.0):
        '''
        Determine if a beam hit a surface

        *Input parameters*

        ``beam``
        A GaussianBeam object

        *Return value*
        A dictionary with the following keys:
        "isHit": A boolean value whether the beam hit the surface or not.
        "Intersection Point": numpy array of the coordinates of the intersection point.
        "distance": Distance between the origin of the beam and the intersection point.
        "localNormVect": A numpy array representing the normal vector of the surface
                                    at the intersection point.
        "localNormAngle": The angle of the localNormVect.
        '''
        if np.abs(inv_ROC) < 1e-5:

            ans = optics.geometric.line_plane_intersection(pos=beam.pos, dirVect=beam.dirVect,
                                          plane_center=surface_center, normalVector=normal_vector,
                                          diameter=surface_size)
            localNormVect = normal_vector
            localNormAngle = np.mod(np.arctan2(localNormVect[1],
                                               localNormVect[0]), 2*pi)
            ans['localNormVect'] = localNormVect
            ans['localNormAngle'] = localNormAngle
            return ans
        else:
            ans = optics.geometric.line_arc_intersection(pos=beam.pos, dirVect=beam.dirVect,
                                                         chord_center=surface_center,
                                                         chordNormVect=normal_vector,
                                                         invROC=inv_ROC,
                                                         diameter=surface_size)

            return ans
        
#}}}

#}}}

#{{{ Mirror Class

class Mirror(Optics):
    '''
    Representing a partial reflective mirror.

    HRcenter: The position of the center of the arc of the HR surface.
                          [(2,) float array]    
              
    HRcenterC: The position of the center of the chord of the HR surface.
                        [(2,) float array]    

    normVectHR: Normal vector of the HR surface.
    
    normAngleHR: Angle of the HR normal vector.

    ARcenter: The position of the center of the AR surface.
              [(2,) float array]

    normVectAR: Normal vector of the HR surface.
    
    normAngleAR: Angle of the HR normal vector.

    HRtransmissive: A boolean value defaults to False. If True, this mirror
                              is supposed to transmit beams on the HR surface. Therefore,
                              for the first encounter of a beam on the HR surface of this mirror
                              will not increase the stray_order. This flag should be set to True for
                              beam splitters and input test masses.

    term_on_HR: If this is True, a beam with stray_order=0 will be terminated when
                          it hits on HR. This is to avoid the inifinite loop of non-sequencial
                          trace by forming a cavity.
    '''

#{{{ Traits definitions

    HRcenter = CArray(dtype=np.float64, shape=(2,))
    HRcenterC = CArray(dtype=np.float64, shape=(2,))
    sagHR = CFloat()
    normVectHR = CArray(dtype=np.float64, shape=(2,))
    normAngleHR = CFloat()
    
    ARcenter = CArray(dtype=np.float64, shape=(2,))
    normVectAR = CArray(dtype=np.float64, shape=(2,))
    normAngleAR = CFloat()
    
    diameter = CFloat(25.0*cm) #
    ARdiameter = CFloat()
    thickness = CFloat(15.0*cm) # 
    wedgeAngle = CFloat(0.25*pi/180) # in rad
    n = CFloat(1.45) #Index of refraction

    inv_ROC_HR = CFloat(1.0/7000.0) #Inverse of the ROC of the HR surface.
    inv_ROC_AR = CFloat(0.0) #Inverse of the ROC of the AR surface.

    Refl_HR = CFloat(99.0) #Power reflectivity of the HR side.
    Trans_HR = CFloat(1.0) #Power transmittance of the HR side.

    Refl_AR = CFloat(0.01) #Power reflectivity of the AR side.
    Trans_AR = CFloat(99.99) #Power transmittance of the HR side.

#}}}

#{{{ __init__

    def __init__(self, HRcenter=[0.0,0.0], normAngleHR=0.0,
                 normVectHR=None, diameter=25.0*cm, thickness=15.0*cm,
                 wedgeAngle=0.25*pi/180., inv_ROC_HR=1.0/7000.0, inv_ROC_AR=0.0,
                 Refl_HR=0.99, Trans_HR=0.01, Refl_AR=0.01, Trans_AR=0.99, n=1.45,
                 name="Mirror", HRtransmissive=False, term_on_HR=False):
        '''
        Create a mirror object.

        === Arguments ===
        HRcenter: Position of the center of the HR surface. [(2,) array]
        
        normAngleHR: Direction angle of the normal vector of the HR surface. [rad]
        
        normVectHR: Normal vector of the HR surface. [(2,) array]
        
        diameter: Diameter of the mirror.
        
        thickness: Thickness of the mirror.
        
        wedgeAngle: Wedge angle between the HR and AR surfaces.

        inv_ROC_HR: 1/ROC of the HR surface.

        inv_ROC_AR: 1/ROC of the AR surface.

        Refl_HR: Power reflectivity of the HR surface.

        Trans_HR: Power transmissivity of the HR surface.
        
        Refl_AR: Power reflectivity of the AR surface.

        Trans_AR: Power transmissivity of the AR surface.

        n: Index of refraction.

        name: Name of the mirror.
        '''
        self.diameter = diameter

        #Compute the sag.
        #Sag is positive for convex mirror.
        if np.abs(inv_ROC_HR) > 1./(10*km):
            R = 1./inv_ROC_HR
            r = self.diameter/2
            self.sagHR =  - np.sign(R)*(np.abs(R) - np.sqrt(R**2 - r**2))
        else:
            self.sagHR = 0.0;

        #Convert rotationAngle to normVectHR or vice versa.
        if normVectHR is not None:
            self.normVectHR = normVectHR
        else:
            self.normAngleHR = normAngleHR

        self.HRcenter = HRcenter
        self._HRcenter_changed(0,0)
        
        self.thickness = thickness
        self.wedgeAngle = wedgeAngle
        self.ARdiameter = self.diameter/np.cos(self.wedgeAngle)
        self.inv_ROC_HR = inv_ROC_HR
        self.inv_ROC_AR = inv_ROC_AR
        self.Refl_HR = Refl_HR
        self.Trans_HR = Trans_HR
        self.Refl_AR = Refl_AR
        self.Trans_AR = Trans_AR
        self.n = n
        self._normAngleHR_changed(0,0)
        self.name = name
        self.HRtransmissive = HRtransmissive
        self.term_on_HR = term_on_HR

#}}}

#{{{ copy

    def copy(self):
        return Mirror(HRcenter=self.HRcenter, normAngleHR=self.normAngleHR,
                      diameter=self.diameter, thickness=self.thickness,
                      wedgeAngle=self.wedgeAngle, inv_ROC_HR=self.inv_ROC_HR,
                      inv_ROC_AR=self.inv_ROC_AR, Refl_HR=self.Refl_HR,
                      Trans_HR=self.Trans_HR, Refl_AR=self.Refl_AR, Trans_AR=self.Trans_AR,
                      n=self.n, name=self.name, HRtransmissive=self.HRtransmissive,
                      term_on_HR=self.term_on_HR)

#}}}

#{{{ rotate

    def rotate(self, angle, center=False):
        '''
        Rotate the mirror. If center is not specified, the center of rotation is
        HRcenter. If center is given (as a vector), the center of rotation is
        center. center is a position vector in the global coordinates.
        '''

        if center:
            center = np.array(center)
            pointer = self.HRcenter - center
            pointer = optics.geometric.vector_rotation_2D(pointer, angle)
            self.HRcenter = center + pointer
            
        self.normAngleHR = self.normAngleHR + angle    
#}}}

#{{{ Translate

    def translate(self, trVect):
        trVect = np.array(trVect)
        self.center = self.center + trVect
        
#}}}

#{{{ Draw

    def draw(self, cv, drawName=False):
        '''
        Draw itself
        '''
        
        plVect = optics.geometric.vector_rotation_2D(self.normVectHR, pi/2)
        p1 = self.HRcenterC + plVect * self.diameter/2
        p2 = p1 - plVect * self.diameter
        p3 = p2 - self.normVectHR * (self.thickness - np.tan(self.wedgeAngle)*self.diameter/2)
        p4 = p1 - self.normVectHR * (self.thickness + np.tan(self.wedgeAngle)*self.diameter/2) 

        cv.add_shape(draw.Line(p2,p3), layername="Mirrors")
        cv.add_shape(draw.Line(p4,p1), layername="Mirrors")        
        
        d = self.thickness/10
        l1 = p1 - self.normVectHR * d
        l2 = p2 - self.normVectHR * d
        cv.add_shape(draw.Line(l1,l2), layername="Mirrors")                

        #Draw Curved surface

        #HR
        if np.abs(self.inv_ROC_HR) > 1.0/1e5:
            R = 1/self.inv_ROC_HR
            theta = np.arcsin(self.diameter/2/R)
            sag = R*(1-np.cos(theta))
            x = np.linspace(0, self.diameter/2, 30)
            y = R*(1.0 - np.sqrt(1.0 - x**2/(R**2))) -sag
            x2 = -np.flipud(x)
            y2 = np.flipud(y)
            x = np.hstack((x2,x))
            y = np.hstack((y2,y))
            v = np.vstack((x,y))
            v = optics.geometric.vector_rotation_2D(v, self.normAngleHR - pi/2)
            v = v.T + self.HRcenterC
            cv.add_shape(draw.PolyLine(x=v[:,0], y=v[:,1]), layername="Mirrors")                            
            #dxf.append(sdxf.LwPolyLine(points=list(v), layer="Mirrors"))
        else:
            cv.add_shape(draw.Line(p1,p2), layername="Mirrors")
            #dxf.append(sdxf.Line(points=[p1,p2], layer="Mirrors"))

        #AR
        if np.abs(self.inv_ROC_AR) > 1.0/1e5:
            diameter = self.diameter/np.cos(self.wedgeAngle)
            
            R = 1/self.inv_ROC_AR
            theta = np.arcsin(diameter/2/R)
            x = np.linspace(0, diameter/2, 30)
            v = optics.geometric.vector_rotation_2D(v, self.normAngleAR - pi/2)
            v = v.T + self.ARcenter
            cv.add_shape(draw.PolyLine(x=v[:,0], y=v[:,1]), layername="Mirrors")
        else:
            cv.add_shape(draw.Line(p3,p4), layername="Mirrors")
            #dxf.append(sdxf.Line(points=[p3,p4], layer="Mirrors"))

        if drawName:
            center = (p1+p2+p3+p4)/4.
            height = self.thickness/4.
            width = height*len(self.name)
            center = center - np.array([width/2, height/2])
            cv.add_shape(draw.Text(text=self.name, point=center,height=height),
                         layername="text")
            # dxf.append(sdxf.Text(text=self.name, point=center, # 
            #                      height=height, layer='text'))
            
#}}}

#{{{ isHit

    def isHit(self, beam):
        '''

          ``beam``:

          or not.

          ``face``:
        '''

        HRsurface = {'center': self.HRcenterC, 'normal_vector': self.normVectHR,
                     'size': self.diameter, 'inv_ROC': self.inv_ROC_HR, 'name': 'HR'}
        ARsurface = {'center': self.ARcenter, 'normal_vector': self.normVectAR,
                     'size': self.diameter, 'inv_ROC': 0.0, 'name': 'AR'}

        #The vector parallel to the HR surface, pointing left.
        v1 = np.array((-self.normVectHR[1], self.normVectHR[0]))
        #Left corner of the HR surface
        c1 = self.HRcenterC + self.diameter/2 * v1
        #Right corner of the HR surface
        c2 = self.HRcenterC - self.diameter/2 * v1
        #Center of the side 1
        side_center_1 = c1 + self.thickness/2 * (- self.normVectHR)
        #Center of the side 2
        side_center_2 = c2 + self.thickness/2 * (- self.normVectHR)
        
        Side1 = {'center': side_center_1, 'normal_vector': v1,
                     'size': self.thickness, 'inv_ROC': 0.0, 'name': 'side'}
        Side2 = {'center': side_center_2, 'normal_vector': -v1,
        
        faceList = [HRsurface, ARsurface, Side1, Side2]

        min_dist = 1e16
        final_answer = None
        for face in faceList:
            ans = self._isHitSurface_(beam, surface_center=face['center'],
                                normal_vector=face['normal_vector'],
                                surface_size=face['size'], inv_ROC=face['inv_ROC'])
            if ans['isHit']:
                if min_dist > ans['distance']:
                    min_dist = ans['distance']
                    final_answer = ans
                    face_name = face['name']

        if final_answer is None:
            return {'isHit': False, 'position': np.array((0,0)),
        else:
            return {'isHit': True, 'position': final_answer['Intersection Point'],
                'distance': min_dist, 'face': face_name}
                
#}}}

#{{{ hit()
    def hit(self, beam, order=0, threshold=0.0, face=False):
        '''

          ``beam``:

        ``isHit``
        or not.

        ``face``:
          
        '''

        #If an optional argument ``face`` is specified
        if face:
            if face == 'HR':
                beams = self.hitFromHR(beam, order=order, threshold=threshold) 
            elif face == 'AR':
                beams = self.hitFromAR(beam, order=order, threshold=threshold)
            else:
                print('Wrong face %s is specified'%face)
                return (False, {}, "")
        #If face is not specified
        else:
            #Check if the beam hit the mirror
            ans = self.isHit(beam)
            face = ans['face']
                if face == 'HR':
                    beams = self.hitFromHR(beam, order=order, threshold=threshold) 
                elif face == 'AR':
                    beams = self.hitFromAR(beam, order=order, threshold=threshold)
                else:
                    #The beam hit a side of the mirror
                    inputBeam = beam.copy()
                    inputBeam.length=ans['distance']
                    return (True, {inputBeam}, "side")
            else:
                #The beam did not hit the mirror
                    
        return (True, beams, face) 
            
#}}}

#{{{ hitFromHR

    def hitFromHR(self, beam, order=0, threshold=0.0, verbose=False):
        '''
        Compute the reflected and deflected beams when
        an input beam hit the HR surface.

        The internal reflections are computed as long as the number
        of internal reflections are below the ``order`` and the power
        of the reflected beams is over the threshold.
        '''

        #A dictionary to hold beams
        beams={}
        
        #Get the intersection point
        ans = optics.geometric.line_arc_intersection(pos=beam.pos, dirVect=beam.dirVect,
                                                     chord_center=self.HRcenterC,
                                                     chordNormVect=self.normVectHR,
                                                     invROC=self.inv_ROC_HR,
                                                     diameter=self.diameter)
        if not ans['isHit']:
            #The input beam does not hit the mirror.
            if verbose:
                print(self.name + ': The beam does not hit the mirror')
            return beams

        #Local normal angle
        localNormAngle = ans['localNormAngle']

        beam_in = beam.copy() #Make a copy
        beam_in.length = ans['distance']
        beam_in.incSurfAngle = localNormAngle
        beam_in.incSurfInvROC = self.inv_ROC_HR
        beams['input']= beam_in

        #Propagate the input beam to the intersection point
        beam_on_HR = beam_in.copy()
        beam_on_HR.propagate(ans['distance'])

        #Calculate reflection and deflection angles along with the ABCD matrices
        #for reflection and deflection.
        (reflAngle, deflAngle, Mrx, Mry, Mtx, Mty) = \
                    optics.geometric.refl_defl_angle(beam_on_HR.dirAngle,
                                               localNormAngle,
                                               1.0, self.n, invROC=self.inv_ROC_HR)        
        #Reflected beam
        beam_r1 = beam_on_HR.copy()
        beam_r1.P = beam_r1.P * self.Refl_HR
        if beam_r1.P > threshold:
            beam_r1.dirAngle = reflAngle
            beam_r1.ABCDTrans(Mrx, Mry)
            beam_r1.departSurfAngle = localNormAngle
            beam_r1.departSurfInvROC = self.inv_ROC_HR        
            beam_r1.incSurfAngle = None
            beam_r1.incSurfInvROC = None
            beam_r1.name = self.name+':r1'
            beams['r1'] = beam_r1
        
        #Transmitted beam
        beam_s1 = beam_on_HR.copy()
        beam_s1.P = beam_s1.P * self.Trans_HR
        if not self.HRtransmissive:
            beam_s1.stray_order = beam_s1.stray_order+1
        if beam_s1.P < threshold or beam_s1.stray_order > order:
        beam_s1.dirAngle = deflAngle
        beam_s1.n = self.n
        beam_s1.ABCDTrans(Mtx, Mty)
        beam_s1.departSurfAngle = np.mod(localNormAngle+pi, 2*pi)
        beam_s1.departSurfInvROC = -self.inv_ROC_HR

        #Hit AR from back
        ans = optics.geometric.line_arc_intersection(pos=beam_s1.pos, dirVect=beam_s1.dirVect,
                                                     chord_center=self.ARcenter,
                                                     chordNormVect=-self.normVectAR,
                                                     invROC=-self.inv_ROC_AR,
                                                     diameter=self.ARdiameter)

        beam_s1.length = ans['distance']
        beam_s1.incSurfAngle = localNormAngle
        beam_s1.incSurfInvROC = -self.inv_ROC_AR
        beam_s1.name = self.name+':s1'
        beams['s1'] = beam_s1

        #Propagate the beam to the AR surface
        beam_on_AR = beam_s1.copy()
        beam_on_AR.propagate(ans['distance'])

                    optics.geometric.refl_defl_angle(beam_on_AR.dirAngle,
                                               self.n, 1.0, invROC=-self.inv_ROC_AR)

        beam_t1 = beam_on_AR.copy()
        beam_t1.P = beam_on_AR.P * self.Trans_AR
        if beam_t1.P > threshold:
            beam_t1.dirAngle = deflAngle
            beam_t1.n = 1.0
            beam_t1.ABCDTrans(Mtx, Mty)        
            beam_t1.departSurfAngle = np.mod(localNormAngle+pi, 2*pi)
            beam_t1.departSurfInvROC = self.inv_ROC_AR
            beam_t1.incSurfAngle = None
            beam_t1.incSurfInvROC = None
            beam_t1.name = self.name+':t%d'%(1)
            beams['t1'] = beam_t1

        beam_sr = beam_on_AR.copy()
        beam_sr.P = beam_sr.P * self.Refl_AR
        beam_sr.stray_order = beam_sr.stray_order+1
        if beam_sr.P < threshold or beam_sr.stray_order > order:
        beam_sr.dirAngle = reflAngle
        beam_sr.ABCDTrans(Mrx, Mry)
        beam_sr.departSurfAngle = localNormAngle
        beam_sr.departSurfInvROC = -self.inv_ROC_AR
        
        #Calculate higher order reflections

        ii = 1
        while ii <= 10*order:

            #Hit the HR from the back
            
            #Get the intersection point
            ans = optics.geometric.line_arc_intersection(pos=beam_sr.pos, dirVect=beam_sr.dirVect,
                                                         chord_center=self.HRcenterC,
                                                         chordNormVect=-self.normVectHR,
                                                         invROC=-self.inv_ROC_HR,
                                                         diameter=self.diameter)

            if not ans['isHit']:
                break

            #Local normal angle
            localNormAngle = ans['localNormAngle']            

            beam_sr.length = ans['distance']
            beam_sr.layer = 'aux_beam'
            beam_sr.incSurfAngle = localNormAngle
            beam_sr.incSurfInvROC = -self.inv_ROC_HR
            beam_sr.name = self.name+':s%d'%(2*ii)              
            beams['s'+str(2*ii)]= beam_sr

            #Propagate the input beam to the intersection point
            beam_on_HR = beam_sr.copy()
            beam_on_HR.propagate(ans['distance'])

            #Calculate reflection and deflection angles along with the ABCD matrices
            #for reflection and deflection.
            (reflAngle, deflAngle, Mrx, Mry, Mtx, Mty) = \
                        optics.geometric.refl_defl_angle(beam_on_HR.dirAngle,
                                                         localNormAngle,
                                                         self.n, 1.0, invROC=-self.inv_ROC_HR)

            #Transmitted through HR
            beam_r1 = beam_on_HR.copy()
            beam_r1.P = beam_r1.P * self.Trans_HR
            beam_r1.stray_order = beam_r1.stray_order+1            
            if beam_r1.P > threshold and beam_r1.stray_order <= order:
                beam_r1.dirAngle = deflAngle
                beam_r1.n = 1.0
                beam_r1.ABCDTrans(Mtx, Mty)
                beam_r1.departSurfAngle = np.mod(localNormAngle+pi, 2*pi)
                beam_r1.departSurfInvROC = self.inv_ROC_HR
                beam_r1.incSurfAngle = None
                beam_r1.incSurfInvROC = None
                beam_r1.name = self.name+':r%d'%(ii+1)                  
                beams['r'+str(ii+1)] = beam_r1

            #Reflected by HR
            beam_s1 = beam_on_HR.copy()
            beam_s1.P = beam_s1.P * self.Refl_HR
            if beam_s1.P < threshold:
                break
            beam_s1.dirAngle = reflAngle
            beam_s1.ABCDTrans(Mrx, Mry)
            beam_s1.departSurfAngle = localNormAngle
            beam_s1.departSurfInvROC = -self.inv_ROC_HR

            #Hit AR from back
            ans = optics.geometric.line_arc_intersection(pos=beam_s1.pos, dirVect=beam_s1.dirVect,
                                                         chord_center=self.ARcenter,
                                                         chordNormVect=-self.normVectAR,
                                                         invROC=-self.inv_ROC_AR,
                                                         diameter=self.ARdiameter)
            
                break

            localNormAngle = ans['localNormAngle']

            beam_s1.incSurfAngle = localNormAngle
            beam_s1.incSurfInvROC = -self.inv_ROC_AR
            beam_s1.length = ans['distance']
            beam_s1.name = self.name+':s%d'%(2*ii+1)              
            beams['s'+str(2*ii+1)] = beam_s1

            #Propagate the beam to the AR surface
            beam_on_AR = beam_s1.copy()
            beam_on_AR.propagate(ans['distance'])

                        optics.geometric.refl_defl_angle(beam_on_AR.dirAngle,
                                                         self.n, 1.0, invROC=-self.inv_ROC_AR)
            #Transmitted beam
            beam_t1 = beam_on_AR.copy()
            beam_t1.P = beam_on_AR.P * self.Trans_AR
            if beam_t1.P > threshold:
                beam_t1.dirAngle = deflAngle
                beam_t1.n = 1.0
                beam_t1.ABCDTrans(Mtx, Mty)
                beam_t1.departSurfAngle = np.mod(localNormAngle+pi, 2*pi)
                beam_t1.departSurfInvROC = self.inv_ROC_AR
                beam_t1.incSurfAngle = None
                beam_t1.incSurfInvROC = None
                beam_t1.name = self.name+':t%d'%(ii+1)  
                beams['t'+str(ii+1)] = beam_t1

            #Reflected beam
            beam_sr = beam_on_AR.copy()
            beam_sr.P = beam_sr.P * self.Refl_AR
            beam_sr.stray_order = beam_sr.stray_order+1
            if beam_sr.P < threshold or beam_sr.stray_order > order:
                break
            beam_sr.dirAngle = reflAngle
            beam_sr.ABCDTrans(Mrx, Mry)
            beam_sr.departSurfAngle = localNormAngle
            beam_sr.departSurfInvROC = -self.inv_ROC_AR

            ii=ii+1

        return beams
#}}}

#{{{ hitFromAR

    def hitFromAR(self, beam, order=0, threshold=0.0, verbose=False):
        '''
        an input beam hit the AR surface.

        '''

        beams={}
        
                                                     chordNormVect=self.normVectAR,
                                                     invROC=self.inv_ROC_AR,

        localNormAngle = ans['localNormAngle']        

        beam_in.incSurfInvROC = self.inv_ROC_AR 

        beam_on_AR = beam_in.copy()

                                               1.0, self.n, invROC=self.inv_ROC_AR)        
        beam_r1 = beam_on_AR.copy()
        beam_r1.P = beam_r1.P * self.Refl_AR
        beam_r1.stray_order = beam_r1.stray_order+1        
        if beam_r1.P > threshold and beam_r1.stray_order <= order:
            beam_r1.departSurfInvROC = self.inv_ROC_AR        
            beam_r1.name = self.name+':r1'            
        
        beam_s1 = beam_on_AR.copy()
        beam_s1.P = beam_s1.P * self.Trans_AR
        if beam_s1.P < threshold:
        beam_s1.departSurfInvROC = -self.inv_ROC_AR

        #Hit HR from back
                                                     chordNormVect=-self.normVectHR,
                                                     invROC=-self.inv_ROC_HR,

        localNormAngle = ans['localNormAngle']                
        beam_s1.incSurfInvROC = -self.inv_ROC_HR
        beam_s1.name = self.name+':s1'        

        #Propagate the beam to the HR surface
        beam_on_HR = beam_s1.copy()

                                               self.n, 1.0, invROC=-self.inv_ROC_HR)

        beam_t1 = beam_on_HR.copy()
        beam_t1.P = beam_on_HR.P * self.Trans_HR
            beam_t1.stray_order = beam_t1.stray_order+1
        if beam_t1.P > threshold and beam_t1.stray_order <= order:
            beam_t1.ABCDTrans(Mtx, Mty)
            beam_t1.departSurfInvROC = self.inv_ROC_HR
            beam_t1.name = self.name+':t1'            
                                               
        beam_sr = beam_on_HR.copy()
        beam_sr.P = beam_sr.P * self.Refl_HR
        if beam_sr.P < threshold:
        beam_sr.departSurfInvROC = -self.inv_ROC_HR
        
        ii = 1

                break

            beam_sr.incSurfInvROC = -self.inv_ROC_AR
            beam_sr.name = self.name+':s%d'%(2*ii)            

            beam_on_AR = beam_sr.copy()

            #Transmitted through AR
            beam_r1 = beam_on_AR.copy()
            beam_r1.P = beam_r1.P * self.Trans_AR
            if beam_r1.P > threshold:
                beam_r1.departSurfInvROC = self.inv_ROC_AR
                beam_r1.name = self.name+':r%d'%(ii+1)                              

            #Reflected by AR
            beam_s1 = beam_on_AR.copy()
            beam_s1.P = beam_s1.P * self.Refl_AR
            beam_s1.stray_order = beam_s1.stray_order+1            
            if beam_s1.P < threshold or beam_s1.stray_order > order:
                break
            beam_s1.departSurfInvROC = -self.inv_ROC_AR

            #Hit HR from back
            
                break
            
           #Local normal angle
            beam_s1.incSurfInvROC = -self.inv_ROC_HR
            beam_s1.name = self.name+':s%d'%(2*ii+1)                          

            #Propagate the beam to the HR surface
            beam_on_HR = beam_s1.copy()

            beam_t1 = beam_on_HR.copy()
            beam_t1.P = beam_t1.P * self.Trans_HR
            beam_t1.stray_order = beam_t1.stray_order+1            
            if beam_t1.P > threshold and beam_t1.stray_order <= order:
                beam_t1.departSurfInvROC = self.inv_ROC_HR
                beam_t1.name = self.name+':t%d'%(ii+1)                

            beam_sr = beam_on_HR.copy()
            beam_sr.P = beam_sr.P * self.Refl_HR
            if beam_sr.P < threshold:
                break
            beam_sr.departSurfInvROC = -self.inv_ROC_HR

            ii=ii+1

#}}}

#{{{ Notification handlers

    def _normAngleHR_changed(self, old, new):
        self.set(trait_change_notify=False,
                 normVectHR = array([np.cos(self.normAngleHR), np.sin(self.normAngleHR)]))
        self.set(trait_change_notify=False,        
                 normAngleHR = np.mod(self.normAngleHR, 2*pi))
        
        self.normVectAR = optics.geometric.vector_rotation_2D(self.normVectHR, pi+self.wedgeAngle)
        self.normAngleAR = np.mod(self.normAngleHR + pi + self.wedgeAngle, 2*pi)
        self.HRcenterC = self.HRcenter - self.normVectHR * self.sagHR
        
    def _normVectHR_changed(self, old, new):
        #Normalize
                 normVectHR = self.normVectHR/np.linalg.norm(array(self.normVectHR)))
        #Update dirAngle accordingly
                 normAngleHR = np.mod(np.arctan2(self.normVectHR[1],
                                                   self.normVectHR[0]), 2*pi))

    def _HRcenterC_changed(self, old, new):
                 ARcenter = self.HRcenterC - self.normVectHR * self.thickness)
                 center = (self.HRcenterC + self.ARcenter)/2.0)
                 HRcenter = self.HRcenterC + self.sagHR*self.normVectHR)

    def _HRcenter_changed(self, old, new):
                 HRcenterC = self.HRcenter - self.sagHR*self.normVectHR)

    def _center_changed(self, old, new):
                 HRcenterC = self.center + self.normVectHR * self.thickness/2.0)

    def _wedgeAngle_changed(self, old, new):
                 normAngleAR = np.mod(self.normAngleHR + pi + self.wedgeAngle, 2*pi))
                 normVectAR = optics.geometric.vector_rotation_2D(self.normVectHR, pi+self.wedgeAngle))

    def _inv_ROC_HR_changed(self, old, new):
        #First update the sag
        if np.abs(self.inv_ROC_HR) > 1./(10*km):
            R = 1./self.inv_ROC_HR
        else:
        #Update the HRcenterC
        
#}}}        

#}}}

'''
Renderer module for gtrace.draw
'''

import gtrace.draw.dxf as dxf
from gtrace.unit import *

#}}}

# are met:
#
#
#

#}}}

#{{{ DXF renderer

class UnknownShapeError(BaseException):
    def __initi__(self, message):
        self.message = message

def renderDXF(canvas, filename):
    '''
    Render a canvas into a DXF file
    '''

    if canvas.unit == 'm':
        scale_factor = 1000.0
    elif canvas.unit == 'km':
        scale_factor = 1e6
    else:
        scale_factor = 1.0
        
    d = dxf.DXF(filename)

    for ly in canvas.layers.values():
        d.add_layer(ly.name, color=dxf.color_encode(ly.color))
        for s in ly.shapes:
            if isinstance(s, draw.Line):
                d.add_entity(dxf.Line(np.array(s.start)*scale_factor, np.array(s.stop)*scale_factor), ly.name)
            elif isinstance(s, draw.PolyLine):
                x = np.array(s.x)*scale_factor
                y = np.array(s.y)*scale_factor
                d.add_entity(dxf.LwPolyLine(x, y), ly.name)
            elif isinstance(s, draw.Rectangle):
                d.add_entity(dxf.Rectangle(np.array(s.point)*scale_factor,
                                           s.width*scale_factor, s.height*scale_factor), ly.name)
            elif isinstance(s, draw.Circle):
                d.add_entity(dxf.Circle(np.array(s.center)*scale_factor, s.radius*scale_factor), ly.name)
            elif isinstance(s, draw.Arc):
                d.add_entity(dxf.Arc(np.array(s.center)*scale_factor, s.radius*scale_factor,
                                     rad2deg(s.startangle), rad2deg(s.stopangle)), ly.name)
            elif isinstance(s, draw.Text):
                d.add_entity(dxf.Text(s.text, np.array(s.point)*scale_factor, s.height*scale_factor, rad2deg(s.rotation)), ly.name)
            else:
                raise UnknownShapeError('Shape not supported')

    d.save_to_file()
                
#}}}

'''
dxf.py - a DXF export library for python

Sample code:

import dxf
d = dxf.DXF('test.dxf')
d.add_layer('ABC', color=5)
d.add_entity(dxf.Line((1.5,5), (56,-89)), 'ABC')
d.save_to_file()

'''

import os
import math
pi = math.pi

#}}}

#{{{ Classes

#{{{ DXF class

class DXF(object):
    '''
    A DXF file class.
    '''

    def __init__(self, filename='drawing.dxf'):
        self.filename = filename
        self.layers = {}  #A dictionary to hold the layers
        self.handseed = 120

    def add_layer(self, name, color=1, ltype='Continuous'):
        self.layers[name] = Layer(name, hex(self.handseed)[2:], color=color, ltype=ltype)
        self.handseed += 1

    def add_entity(self, entity, layername):
        if not self.layers.has_key(layername):
            self.add_layer(layername)

        entity.set_handle(hex(self.handseed)[2:])
        self.layers[layername].add_entity(entity)

#{{{ def save_to_file(self):

    def save_to_file(self):
        '''
        Save the DXF file
        '''

        #Open a file
        dxffile = open(self.filename, 'w')

#{{{ Header

        #Compute the lower left and the upper left corners of the drawing
        xmin = 0.0
        xmax = 0.0
        ymin = 0.0
        ymax = 0.0

        for ly in self.layers.values():
            for ent in ly.entities:
                a = ent.report_min_max()
                xmin = a[0][0] if a[0][0] < xmin else xmin
                xmax = a[1][0] if a[1][0] > xmax else xmax                
                ymin = a[0][1] if a[0][1] < ymin else ymin
                ymax = a[1][1] if a[1][1] > ymax else ymax                
        
        dxffile.write(header_template1)

        dxffile.write('  9\n$EXTMIN\n 10\n')
        dxffile.write(str(xmin)+'\n 20\n'+str(ymin)+'\n 30\n0.0\n')        
        dxffile.write('  9\n$EXTMAX\n 10\n')
        dxffile.write(str(xmax)+'\n 20\n'+str(ymax)+'\n 30\n0.0\n')        
        dxffile.write('  9\n$LIMMIN\n 10\n')
        dxffile.write(str(xmin)+'\n 20\n'+str(ymin)+'\n')        
        dxffile.write('  9\n$LIMMAX\n 10\n')
        dxffile.write(str(xmax)+'\n 20\n'+str(ymax)+'\n')        

        dxffile.write(header_template3)

#}}}

        dxffile.write(classes_template)
#}}}

#{{{ Tables
        dxffile.write(tables_template1)
        dxffile.write(vport_template)    
        dxffile.write(linetype_template)

        #Layers
        dxffile.write(layer_template1)
        #Loop for layers
            dxffile.write('  0\nLAYER\n  5\n')
            dxffile.write(ly.handle+'\n')
            dxffile.write('330\n2\n100\nAcDbSymbolTableRecord\n')
            dxffile.write('100\nAcDbLayerTableRecord\n  2\n')
            dxffile.write(ly.name+'\n')
            dxffile.write(' 70\n0\n 62\n')
            dxffile.write(str(ly.color)+'\n')
            dxffile.write('  6\n')
            dxffile.write(ly.ltype+'\n')
            dxffile.write('370\n    -3\n390\nF\n347\n3E\n')

        dxffile.write('  0\nENDTAB\n')

        dxffile.write(style_template)
        dxffile.write(view_template)
        dxffile.write(misctable_template)    
        dxffile.write(blockrecords_template)
        dxffile.write(tables_template2)

#}}}

#{{{ Blocks
        dxffile.write(blocks_template)
#}}}

#{{{ Entities

        dxffile.write(entities_template1)
                dxffile.write(ent.draw(ly.name))

        dxffile.write('  0\nENDSEC\n')
#}}}

#{{{ Objects
        dxffile.write(objects_template)
#}}}

        dxffile.write('0\nEOF\n')
        dxffile.close()

#}}}

#}}}

#{{{ Layer
class Layer(object):
    '''
    Layer class
    '''

    def __init__(self, name, handle, color=1, ltype='Continuous'):
        self.handle = handle
        self.color = color
        self.ltype = ltype
        self.entities = [] #A list of entities, such as lines, circles, etc.

    def add_entity(self, entity):
        self.entities.append(entity)

#}}}

#{{{ Entity

class Entity(object):
    '''
    A graphic entity
    '''

    def __init__(self):
        self.handle = ''

    def set_handle(self, handle):

    def report_min_max(self):
        '''
        Return the coordinates of the lower left and the upper right corners
        of the drawing.
        Return value: ((xmin, ymin), (xmax, ymax))
        '''
        return ((0.0,0.0), (1.0,1.0))

    def draw(self):
        return ''

#}}}

#{{{ Line

class Line(Entity):
    '''
    A line entity
    '''

    def __init__(self, start, stop, thickness=0):
        super(Line, self).__init__()
        self.start = start
        self.stop = stop

        '''
        '''
        xmin = min((self.start[0], self.stop[0]))
        xmax = max((self.start[0], self.stop[0]))
        ymin = min((self.start[1], self.stop[1]))
        ymax = max((self.start[1], self.stop[1]))
        
        return ((xmin, ymin), (xmax, ymax))

    def draw(self, layername):
        output ='''  0\nLINE\n  5\n'''
        output += self.handle+'\n'
        output +='330\n1F\n100\nAcDbEntity\n  8\n'
        output += layername+'\n'
        output += '100\nAcDbLine\n 10\n'
        output += str(self.start[0])+'\n 20\n'
        output += str(self.start[1])+'\n 30\n0.0\n  11\n'
        output += str(self.stop[0])+'\n 21\n'
        output += str(self.stop[1])+'\n 31\n0.0\n'        

        return output

#}}}    

#{{{ LWPOLYLINE

class NumberOfElementError(BaseException):
        
class LwPolyLine(Entity):
    '''
    A light weight poly-line
    '''

    def __init__(self, x, y, thickness=0):
        '''
        = Arguments =
        x: x coordinates of the vertices
        y: y coordinates of the vertices        
        '''
        super(LwPolyLine, self).__init__()
        self.x = x
        self.y = y
        if len(x) != len(y):
            raise NumberOfElementError('The numbers of elements of x and y do not match.')
        self.numpoints = len(x)

        '''
        '''
        xmin = min(self.x)
        xmax = max(self.x)
        ymin = min(self.y)
        ymax = max(self.y)
        
        output ='''  0\nLWPOLYLINE\n  5\n'''
        output += '100\nAcDbPolyline\n 90\n'
        output += str(self.numpoints)+'\n 70\n0\n43\n0.0\n'
        for ii in range(self.numpoints):
            output += ' 10\n'
            output += str(self.x[ii])+'\n 20\n'
            output += str(self.y[ii])+'\n'
            
#}}}

#{{{ Rectangle

class Rectangle(Entity):
    '''
    A rectangle
    '''

    def __init__(self, point, width, height, thickness=0):
        '''

        '''
        super(Rectangle, self).__init__()
        self.point = point
        self.width = width
        self.height = height
        x0 = point[0]
        y0 = point[1]        
        self.x = [x0, x0+width, x0+width, x0, x0]
        self.y = [y0, y0, y0+height, y0+height, y0]        

        '''
        '''
        xmin = min(self.point[0], self.point[0]+self.width)
        xmax = max(self.point[0], self.point[0]+self.width)        
        ymin = min(self.point[1], self.point[1]+self.height)
        ymax = max(self.point[1], self.point[1]+self.height)
        
        output += str(5)+'\n 70\n0\n43\n0.0\n'
        for ii in range(5):
            
#}}}

#{{{ Circle

class Circle(Entity):
    '''
    A circle entity
    '''

    def __init__(self, center, radius, thickness=0):
        super(Circle, self).__init__()
        self.center = center
        self.radius = radius

        '''
        '''
        xmin = self.center[0]-self.radius
        xmax = self.center[0]+self.radius
        ymin = self.center[1]-self.radius
        ymax = self.center[1]+self.radius
        
        output ='''  0\nCIRCLE\n  5\n'''
        output += '100\nAcDbCircle\n 10\n'
        output += str(self.center[0])+'\n 20\n'
        output += str(self.center[1])+'\n 30\n0.0\n  40\n'
        output += str(self.radius)+'\n'

#}}}

#{{{ Arc

class Arc(Entity):
    '''
    An arc entity
    '''

    def __init__(self, center, radius, startangle, stopangle, thickness=0):
        super(Arc, self).__init__()
        self.startangle = startangle
        self.stopangle = stopangle

        '''
        '''
        if self.startangle > self.stopangle:
            stopangle = self.stopangle + 360.0
        else:
            stopangle = self.stopangle

        startangle = self.startangle

        dx1 = self.radius*math.cos(pi*startangle/180)
        dy1 = self.radius*math.sin(pi*startangle/180)        
        dx2 = self.radius*math.cos(pi*stopangle/180)
        dy2 = self.radius*math.sin(pi*stopangle/180)

        #Assign temporary values for xmin etc.
        xmin = min(dx1, dx2)+self.center[0]
        xmax = max(dx1, dx2)+self.center[0]            
        ymin = min(dy1, dy2)+self.center[1]
        ymax = max(dy1, dy2)+self.center[1]            

        #Quadrature of the start point
        start_quad = int(math.floor(startangle/90.0))
        #Quadrature of the end point
        stop_quad = int(math.floor(stopangle/90.0))
        #List of poles in between the start and the stop quadratures
        poles = range(start_quad, stop_quad)

        for n in poles:
            a = n%4
            if a == 0:
                ymax = self.center[1]+self.radius
            elif a == 1:
                xmin = self.center[0]-self.radius
            elif a == 2:
                ymin = self.center[1]-self.radius
            elif a == 3:
                xmax = self.center[0]+self.radius
        
        output ='''  0\nARC\n  5\n'''
        output += str(self.radius)+'\n100\nAcDbArc\n 50\n'
        output += str(self.startangle)+'\n 51\n'
        output += str(self.stopangle)+'\n'
        
#}}}

#{{{ Text

class Text(Entity):
    '''
    Text
    '''

    def __init__(self, text, point, height=1.0, rotation=0.0):
        super(Text, self).__init__()
        self.text = text
        self.rotation = rotation

        '''
        '''
        length = len(self.text)*self.height
        dx = length*math.cos(pi*self.rotation/180)
        dy = length*math.sin(pi*self.rotation/180)        
        xmin = min(self.point[0], self.point[0]+dx)
        xmax = max(self.point[0], self.point[0]+dx)
        ymin = min(self.point[1], self.point[1]+dy)
        ymax = max(self.point[1], self.point[1]+dy)
        
        output ='''  0\nTEXT\n  5\n'''
        output += '100\nAcDbText\n 10\n'
        output += str(self.point[0])+'\n 20\n'
        output += str(self.point[1])+'\n 30\n0.0\n  40\n'
        output += str(self.height)+'\n  1\n'
        output += self.text+'\n 50\n'
        output += str(self.rotation)+'\n 7\nStandard\n100\nAcDbText\n'
        
#}}}

#}}}

#{{{ Utilities

def color_encode(color):
    '''
    Given a set of RGB values for a color, find the closest matching
    one from the pre-defined colors in the DXF specification.
    Then return its color code (an integer in 1 to 255).

    = Input =
    color: A tuple of three numbers in the range of 0-255, i.e. (R,G,B)

    = Return =
    best_color_num: integer
    '''

    min_dist = 1e10
    best_color_num = 1
    for c in color_table.keys():
        candidate = color_table[c]
        dist = (candidate[0]-color[0])**2 + (candidate[1]-color[1])**2+(candidate[2]-color[2])**2
        if dist < min_dist:
            min_dist = dist
            best_color_num = c

    return best_color_num

#}}}

#{{{ Color Table

color_table = {1:(255,0,0),
2:(255,255,0),
3:(0,255,0),
4:(0,255,255),
5:(0,0,255),
6:(255,0,255),
7:(255,255,255),
8:(128,128,128),
9:(192,192,192),
10:(255,0,0),
11:(255,127,127),
12:(165,0,0),
13:(165,82,82),
14:(127,0,0),
15:(127,63,63),
16:(76,0,0),
17:(76,38,38),
18:(38,0,0),
19:(38,19,19),
20:(255,63,0),
21:(255,159,127),
22:(165,41,0),
23:(165,103,82),
24:(127,31,0),
25:(127,79,63),
26:(76,19,0),
27:(76,47,38),
28:(38,9,0),
29:(38,23,19),
30:(255,127,0),
31:(255,191,127),
32:(165,82,0),
33:(165,124,82),
34:(127,63,0),
35:(127,95,63),
36:(76,38,0),
37:(76,57,38),
38:(38,19,0),
39:(38,28,19),
40:(255,191,0),
41:(255,223,127),
42:(165,124,0),
43:(165,145,82),
44:(127,95,0),
45:(127,111,63),
46:(76,57,0),
47:(76,66,38),
48:(38,28,0),
49:(38,33,19),
50:(255,255,0),
51:(255,255,127),
52:(165,165,0),
53:(165,165,82),
54:(127,127,0),
55:(127,127,63),
56:(76,76,0),
57:(76,76,38),
58:(38,38,0),
59:(38,38,19),
60:(191,255,0),
61:(223,255,127),
62:(124,165,0),
63:(145,165,82),
64:(95,127,0),
65:(111,127,63),
66:(57,76,0),
67:(66,76,38),
68:(28,38,0),
69:(33,38,19),
70:(127,255,0),
71:(191,255,127),
72:(82,165,0),
73:(124,165,82),
74:(63,127,0),
75:(95,127,63),
76:(38,76,0),
77:(57,76,38),
78:(19,38,0),
79:(28,38,19),
80:(63,255,0),
81:(159,255,127),
82:(41,165,0),
83:(103,165,82),
84:(31,127,0),
85:(79,127,63),
86:(19,76,0),
87:(47,76,38),
88:(9,38,0),
89:(23,38,19),
90:(0,255,0),
91:(127,255,127),
92:(0,165,0),
93:(82,165,82),
94:(0,127,0),
95:(63,127,63),
96:(0,76,0),
97:(38,76,38),
98:(0,38,0),
99:(19,38,19),
100:(0,255,63),
101:(127,255,159),
102:(0,165,41),
103:(82,165,103),
104:(0,127,31),
105:(63,127,79),
106:(0,76,19),
107:(38,76,47),
108:(0,38,9),
109:(19,38,23),
110:(0,255,127),
111:(127,255,191),
112:(0,165,82),
113:(82,165,124),
114:(0,127,63),
115:(63,127,95),
116:(0,76,38),
117:(38,76,57),
118:(0,38,19),
119:(19,38,28),
120:(0,255,191),
121:(127,255,223),
122:(0,165,124),
123:(82,165,145),
124:(0,127,95),
125:(63,127,111),
126:(0,76,57),
127:(38,76,66),
128:(0,38,28),
129:(19,38,33),
130:(0,255,255),
131:(127,255,255),
132:(0,165,165),
133:(82,165,165),
134:(0,127,127),
135:(63,127,127),
136:(0,76,76),
137:(38,76,76),
138:(0,38,38),
139:(19,38,38),
140:(0,191,255),
141:(127,223,255),
142:(0,124,165),
143:(82,145,165),
144:(0,95,127),
145:(63,111,127),
146:(0,57,76),
147:(38,66,76),
148:(0,28,38),
149:(19,33,38),
150:(0,127,255),
151:(127,191,255),
152:(0,82,165),
153:(82,124,165),
154:(0,63,127),
155:(63,95,127),
156:(0,38,76),
157:(38,57,76),
158:(0,19,38),
159:(19,28,38),
160:(0,63,255),
161:(127,159,255),
162:(0,41,165),
163:(82,103,165),
164:(0,31,127),
165:(63,79,127),
166:(0,19,76),
167:(38,47,76),
168:(0,9,38),
169:(19,23,38),
170:(0,0,255),
171:(127,127,255),
172:(0,0,165),
173:(82,82,165),
174:(0,0,127),
175:(63,63,127),
176:(0,0,76),
177:(38,38,76),
178:(0,0,38),
179:(19,19,38),
180:(63,0,255),
181:(159,127,255),
182:(41,0,165),
183:(103,82,165),
184:(31,0,127),
185:(79,63,127),
186:(19,0,76),
187:(47,38,76),
188:(9,0,38),
189:(23,19,38),
190:(127,0,255),
191:(191,127,255),
192:(82,0,165),
193:(124,82,165),
194:(63,0,127),
195:(95,63,127),
196:(38,0,76),
197:(57,38,76),
198:(19,0,38),
199:(28,19,38),
200:(191,0,255),
201:(223,127,255),
202:(124,0,165),
203:(145,82,165),
204:(95,0,127),
205:(111,63,127),
206:(57,0,76),
207:(66,38,76),
208:(28,0,38),
209:(33,19,38),
210:(255,0,255),
211:(255,127,255),
212:(165,0,165),
213:(165,82,165),
214:(127,0,127),
215:(127,63,127),
216:(76,0,76),
217:(76,38,76),
218:(38,0,38),
219:(38,19,38),
220:(255,0,191),
221:(255,127,223),
222:(165,0,124),
223:(165,82,145),
224:(127,0,95),
225:(127,63,111),
226:(76,0,57),
227:(76,38,66),
228:(38,0,28),
229:(38,19,33),
230:(255,0,127),
231:(255,127,191),
232:(165,0,82),
233:(165,82,124),
234:(127,0,63),
235:(127,63,95),
236:(76,0,38),
237:(76,38,57),
238:(38,0,19),
239:(38,19,28),
240:(255,0,63),
241:(255,127,159),
242:(165,0,41),
243:(165,82,103),
244:(127,0,31),
245:(127,63,79),
246:(76,0,19),
247:(76,38,47),
248:(38,0,9),
249:(38,19,23),
250:(84,84,84),
251:(118,118,118),
252:(160,160,160),
253:(192,192,192),
254:(224,224,224),
255:(0,0,0)}

#}}}

#{{{ Templates

header_template1 = \
'''  0
SECTION
  2
HEADER
  9
$ACADVER
  1
AC1024
  9
$ACADMAINTVER
 70
     6
  9
$DWGCODEPAGE
  3
ANSI_1252
  9
$INSBASE
 10
0.0
 20
0.0
 30
0.0
'''

header_template2 =\
''' 9
$EXTMIN
 10
1.0
 20
0.0
 30
0.0
  9
$EXTMAX
 10
73.0
 20
41.0
 30
0.0
  9
$LIMMIN
 10
0.0
 20
0.0
  9
$LIMMAX
 10
1000.0
 20
1000.0
'''

header_template3 =\
'''  9
$ORTHOMODE
 70
     0
  9
$REGENMODE
 70
     1
  9
$FILLMODE
 70
     1
  9
$QTEXTMODE
 70
     0
  9
$MIRRTEXT
 70
     1
  9
$LTSCALE
 40
1.0
  9
$ATTMODE
 70
     1
  9
$TEXTSIZE
 40
2.5
  9
$TRACEWID
 40
1.0
  9
$TEXTSTYLE
  7
Standard
  9
$CLAYER
  8
Outline
  9
$CELTYPE
  6
ByLayer
  9
$CECOLOR
 62
   256
  9
$CELTSCALE
 40
1.0
  9
$DISPSILH
 70
     0
  9
$DIMSCALE
 40
1.0
  9
$DIMASZ
 40
2.5
  9
$DIMEXO
 40
0.625
  9
$DIMDLI
 40
3.75
  9
$DIMRND
 40
0.0
  9
$DIMDLE
 40
0.0
  9
$DIMEXE
 40
1.25
  9
$DIMTP
 40
0.0
  9
$DIMTM
 40
0.0
  9
$DIMTXT
 40
2.5
  9
$DIMCEN
 40
2.5
  9
$DIMTSZ
 40
0.0
  9
$DIMTOL
 70
     0
  9
$DIMLIM
 70
     0
  9
$DIMTIH
 70
     0
  9
$DIMTOH
 70
     0
  9
$DIMSE1
 70
     0
  9
$DIMSE2
 70
     0
  9
$DIMTAD
 70
     1
  9
$DIMZIN
 70
     8
  9
$DIMBLK
  1

  9
$DIMASO
 70
     1
  9
$DIMSHO
 70
     1
  9
$DIMPOST
  1

  9
$DIMAPOST
  1

  9
$DIMALT
 70
     0
  9
$DIMALTD
 70
     3
  9
$DIMALTF
 40
0.03937007874016
  9
$DIMLFAC
 40
1.0
  9
$DIMTOFL
 70
     1
  9
$DIMTVP
 40
0.0
  9
$DIMTIX
 70
     0
  9
$DIMSOXD
 70
     0
  9
$DIMSAH
 70
     0
  9
$DIMBLK1
  1

  9
$DIMBLK2
  1

  9
$DIMSTYLE
  2
ISO-25
  9
$DIMCLRD
 70
     0
  9
$DIMCLRE
 70
     0
  9
$DIMCLRT
 70
     0
  9
$DIMTFAC
 40
1.0
  9
$DIMGAP
 40
0.625
  9
$DIMJUST
 70
     0
  9
$DIMSD1
 70
     0
  9
$DIMSD2
 70
     0
  9
$DIMTOLJ
 70
     0
  9
$DIMTZIN
 70
     8
  9
$DIMALTZ
 70
     0
  9
$DIMALTTZ
 70
     0
  9
$DIMUPT
 70
     0
  9
$DIMDEC
 70
     2
  9
$DIMTDEC
 70
     2
  9
$DIMALTU
 70
     2
  9
$DIMALTTD
 70
     3
  9
$DIMTXSTY
  7
Standard
  9
$DIMAUNIT
 70
     0
  9
$DIMADEC
 70
     0
  9
$DIMALTRND
 40
0.0
  9
$DIMAZIN
 70
     0
  9
$DIMDSEP
 70
    44
  9
$DIMATFIT
 70
     3
  9
$DIMFRAC
 70
     0
  9
$DIMLDRBLK
  1

  9
$DIMLUNIT
 70
     2
  9
$DIMLWD
 70
    -2
  9
$DIMLWE
 70
    -2
  9
$DIMTMOVE
 70
     0
  9
$DIMFXL
 40
1.0
  9
$DIMFXLON
 70
     0
  9
$DIMJOGANG
 40
0.7853981633974483
  9
$DIMTFILL
 70
     0
  9
$DIMTFILLCLR
 70
     0
  9
$DIMARCSYM
 70
     0
  9
$DIMLTYPE
  6

  9
$DIMLTEX1
  6

  9
$DIMLTEX2
  6

  9
$DIMTXTDIRECTION
 70
     0
  9
$LUNITS
 70
     2
  9
$LUPREC
 70
     4
  9
$SKETCHINC
 40
1.0
  9
$FILLETRAD
 40
10.0
  9
$AUNITS
 70
     0
  9
$AUPREC
 70
     0
  9
$MENU
  1
.
  9
$ELEVATION
 40
0.0
  9
$PELEVATION
 40
0.0
  9
$THICKNESS
 40
0.0
  9
$LIMCHECK
 70
     0
  9
$CHAMFERA
 40
10.0
  9
$CHAMFERB
 40
10.0
  9
$CHAMFERC
 40
20.0
  9
$CHAMFERD
 40
0.0
  9
$SKPOLY
 70
     0
  9
$TDCREATE
 40
2456231.748506945
  9
$TDUCREATE
 40
2456231.373506945
  9
$TDUPDATE
 40
  9
$TDUUPDATE
 40
  9
$TDINDWG
 40
0.0000000116
  9
$TDUSRTIMER
 40
  9
$USRTIMER
 70
     1
  9
$ANGBASE
 50
0.0
  9
$ANGDIR
 70
     0
  9
$PDMODE
 70
     0
  9
$PDSIZE
 40
0.0
  9
$PLINEWID
 40
0.0
  9
$SPLFRAME
 70
     0
  9
$SPLINETYPE
 70
     6
  9
$SPLINESEGS
 70
     8
  9
$HANDSEED
  5
6F
  9
$SURFTAB1
 70
     6
  9
$SURFTAB2
 70
     6
  9
$SURFTYPE
 70
     6
  9
$SURFU
 70
     6
  9
$SURFV
 70
     6
  9
$UCSBASE
  2

  9
$UCSNAME
  2

  9
$UCSORG
 10
0.0
 20
0.0
 30
0.0
  9
$UCSXDIR
 10
1.0
 20
0.0
 30
0.0
  9
$UCSYDIR
 10
0.0
 20
1.0
 30
0.0
  9
$UCSORTHOREF
  2

  9
$UCSORTHOVIEW
 70
     0
  9
$UCSORGTOP
 10
0.0
 20
0.0
 30
0.0
  9
$UCSORGBOTTOM
 10
0.0
 20
0.0
 30
0.0
  9
$UCSORGLEFT
 10
0.0
 20
0.0
 30
0.0
  9
$UCSORGRIGHT
 10
0.0
 20
0.0
 30
0.0
  9
$UCSORGFRONT
 10
0.0
 20
0.0
 30
0.0
  9
$UCSORGBACK
 10
0.0
 20
0.0
 30
0.0
  9
$PUCSBASE
  2

  9
$PUCSNAME
  2

  9
$PUCSORG
 10
0.0
 20
0.0
 30
0.0
  9
$PUCSXDIR
 10
1.0
 20
0.0
 30
0.0
  9
$PUCSYDIR
 10
0.0
 20
1.0
 30
0.0
  9
$PUCSORTHOREF
  2

  9
$PUCSORTHOVIEW
 70
     0
  9
$PUCSORGTOP
 10
0.0
 20
0.0
 30
0.0
  9
$PUCSORGBOTTOM
 10
0.0
 20
0.0
 30
0.0
  9
$PUCSORGLEFT
 10
0.0
 20
0.0
 30
0.0
  9
$PUCSORGRIGHT
 10
0.0
 20
0.0
 30
0.0
  9
$PUCSORGFRONT
 10
0.0
 20
0.0
 30
0.0
  9
$PUCSORGBACK
 10
0.0
 20
0.0
 30
0.0
  9
$USERI1
 70
     0
  9
$USERI2
 70
     0
  9
$USERI3
 70
     0
  9
$USERI4
 70
     0
  9
$USERI5
 70
     0
  9
$USERR1
 40
0.0
  9
$USERR2
 40
0.0
  9
$USERR3
 40
0.0
  9
$USERR4
 40
0.0
  9
$USERR5
 40
0.0
  9
$WORLDVIEW
 70
     1
  9
$SHADEDGE
 70
     3
  9
$SHADEDIF
 70
    70
  9
$TILEMODE
 70
     1
  9
$MAXACTVP
 70
    64
  9
$PINSBASE
 10
0.0
 20
0.0
 30
0.0
  9
$PLIMCHECK
 70
     0
  9
$PEXTMIN
 10
1.000000000000000E+20
 20
 30
  9
$PEXTMAX
 10
-1.000000000000000E+20
 20
 30
  9
$PLIMMIN
 10
0.0
 20
0.0
  9
$PLIMMAX
 10
420.0
 20
297.0
  9
$UNITMODE
 70
     0
  9
$VISRETAIN
 70
     1
  9
$PLINEGEN
 70
     0
  9
$PSLTSCALE
 70
     1
  9
$TREEDEPTH
 70
  3020
  9
$CMLSTYLE
  2
Standard
  9
$CMLJUST
 70
     0
  9
$CMLSCALE
 40
20.0
  9
$PROXYGRAPHICS
 70
     1
  9
$MEASUREMENT
 70
     1
  9
$CELWEIGHT
370
    -1
  9
$ENDCAPS
280
     0
  9
$JOINSTYLE
280
     0
  9
$LWDISPLAY
290
     0
  9
$INSUNITS
 70
     4
  9
$HYPERLINKBASE
  1

  9
$STYLESHEET
  1

  9
$XEDIT
290
     1
  9
$CEPSNTYPE
380
     0
  9
$PSTYLEMODE
290
     1
  9
$FINGERPRINTGUID
  2
{92905EEA-4798-D3AF-0738-BBE039DCC74E}
  9
$VERSIONGUID
  2
{FAEB1C32-E019-11D5-929B-00C0DF256EC4}
  9
$EXTNAMES
290
     1
  9
$PSVPSCALE
 40
0.0
  9
$OLESTARTUP
290
     0
  9
$SORTENTS
280
   127
  9
$INDEXCTL
280
     0
  9
$HIDETEXT
280
     1
  9
$XCLIPFRAME
280
     2
  9
$HALOGAP
280
     0
  9
$OBSCOLOR
 70
   257
  9
$OBSLTYPE
280
     0
  9
$INTERSECTIONDISPLAY
280
     0
  9
$INTERSECTIONCOLOR
 70
   257
  9
$DIMASSOC
280
     2
  9
$PROJECTNAME
  1

  9
$CAMERADISPLAY
290
     0
  9
$LENSLENGTH
 40
50.0
  9
$CAMERAHEIGHT
 40
0.0
  9
$STEPSPERSEC
 40
2.0
  9
$STEPSIZE
 40
6.0
  9
$3DDWFPREC
 40
2.0
  9
$PSOLWIDTH
 40
5.0
  9
$PSOLHEIGHT
 40
80.0
  9
$LOFTANG1
 40
1.570796326794897
  9
$LOFTANG2
 40
  9
$LOFTMAG1
 40
0.0
  9
$LOFTMAG2
 40
0.0
  9
$LOFTPARAM
 70
     7
  9
$LOFTNORMALS
280
     1
  9
$LATITUDE
 40
37.795
  9
$LONGITUDE
 40
-122.394
  9
$NORTHDIRECTION
 40
0.0
  9
$TIMEZONE
 70
 -8000
  9
$LIGHTGLYPHDISPLAY
280
     1
  9
$TILEMODELIGHTSYNCH
280
     1
  9
$CMATERIAL
347
3C
  9
$SOLIDHIST
280
     1
  9
$SHOWHIST
280
     1
  9
$DWFFRAME
280
     2
  9
$DGNFRAME
280
     2
  9
$REALWORLDSCALE
290
     1
  9
$INTERFERECOLOR
 62
   256
  9
$CSHADOW
280
     0
  9
$SHADOWPLANELOCATION
 40
0.0
  0
ENDSEC
'''
#}}}

classes_template =\
 '''  0
SECTION
  2
CLASSES
  0
CLASS
  1
ACDBDICTIONARYWDFLT
  2
AcDbDictionaryWithDefault
  3
ObjectDBX Classes
 90
        0
 91
        4
280
     0
281
     0
  0
CLASS
  1
VISUALSTYLE
  2
AcDbVisualStyle
  3
 90
     4095
 91
        4
280
     0
281
     0
  0
CLASS
  1
MATERIAL
  2
AcDbMaterial
  3
 90
     1153
 91
        4
280
     0
281
     0
  0
CLASS
  1
SUN
  2
AcDbSun
  3
SCENEOE
 90
     1024
 91
        4
280
     0
281
     0
  0
ENDSEC
'''

#}}}

tables_template1 =\
'''  0
SECTION
  2
TABLES
'''
#{{{ vport

vport_template =\
'''  0
TABLE
  2
VPORT
  5
8
330
0
100
AcDbSymbolTable
 70
     1
  0
VPORT
  5
29
330
8
100
AcDbSymbolTableRecord
100
AcDbViewportTableRecord
  2
*Active
 70
     0
 10
0.0
 20
0.0
 11
1.0
 21
1.0
 12
360.5069142256372
 22
77.37740091485402
 13
0.0
 23
0.0
 14
10.0
 24
10.0
 15
10.0
 25
10.0
 16
0.0
 26
0.0
 36
1.0
 17
0.0
 27
0.0
 37
0.0
 40
643.00113180099
 41
2.116173120728929
 42
50.0
 43
0.0
 44
0.0
 50
0.0
 51
0.0
 71
     0
 72
   100
 73
     1
 74
     3
 75
     0
 76
     0
 77
     0
 78
     0
281
     0
 65
     1
110
0.0
120
0.0
130
0.0
111
1.0
121
0.0
131
0.0
112
0.0
122
1.0
132
0.0
 79
     0
146
0.0
 60
     3
 61
     5
292
     1
282
     1
141
0.0
142
0.0
 63
   250
361
3F
  0
ENDTAB
'''

#}}}

#{{{ linetype

linetype_template =\
'''  0
TABLE
  2
LTYPE
  5
5
330
0
100
 70
     7
  0
LTYPE
  5
14
330
5
100
100
AcDbLinetypeTableRecord
  2
ByBlock
 70
     0
  3

 72
    65
 73
     0
 40
0.0
  0
LTYPE
  5
15
330
5
100
100
  2
ByLayer
 70
     0
  3

 72
    65
 73
     0
 40
0.0
  0
LTYPE
  5
16
330
5
100
100
  2
Continuous
 70
     0
  3
Solid line
 72
    65
 73
     0
 40
0.0
  0
LTYPE
  5
40
330
5
100
100
  2
DOT
 70
     0
  3
Dot . . . . . . . . . . . . . . . . . . . . . . . .
 72
    65
 73
     2
 40
2.032
 49
1.016
 74
     0
 49
-1.016
 74
     0
  0
LTYPE
  5
41
330
5
100
100
  2
CENTER
 70
     0
  3
Center ____ _ ____ _ ____ _ ____ _ ____ _ ____
 72
    65
 73
     4
 40
13.208
 49
8.128
 74
     0
 49
-1.524
 74
     0
 49
2.032
 74
     0
 49
-1.524
 74
     0
  0
LTYPE
  5
42
330
5
100
100
  2
HIDDEN2
 70
     0
  3
Hidden (.5x) _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
 72
    65
 73
     2
 40
4.762499999999999
 49
3.175
 74
     0
 49
-1.5875
 74
     0
  0
LTYPE
  5
43
330
5
100
100
  2
DASHED
 70
     0
  3
Dashed __ __ __ __ __ __ __ __ __ __ __ __ __ _
 72
    65
 73
     2
 40
9.652
 49
8.128
 74
     0
 49
-1.524
 74
     0
  0
LTYPE
  5
44
330
5
100
100
  2
PHANTOM2
 70
     0
  3
Phantom (.5x) ___ _ _ ___ _ _ ___ _ _ ___ _ _
 72
    65
 73
     6
 40
14.732
 49
8.128
 74
     0
 49
-1.524
 74
     0
 49
1.016
 74
     0
 49
-1.524
 74
     0
 49
1.016
 74
     0
 49
-1.524
 74
     0
  0
LTYPE
  5
45
330
5
100
100
  2
DASHDOT
 70
     0
  3
Dash dot __ . __ . __ . __ . __ . __ . __ . __
 72
    65
 73
     4
 40
11.176
 49
8.128
 74
     0
 49
-1.524
 74
     0
 49
0.0
 74
     0
 49
-1.524
 74
     0
  0
ENDTAB
'''
#}}}

#{{{ layer

layer_template1 =\
'''  0
TABLE
  2
LAYER
  5
2
330
0
100
 70
     7
  0
LAYER
  5
10
330
2
100
100
AcDbLayerTableRecord
  2
0
 70
     0
 62
     5
  6
370
    -3
390
F
347
3E
  0
LAYER
  5
46
330
2
100
100
  2
Default
 70
     0
 62
     255
  6
370
    -3
390
F
347
3E
'''
layer_template2 =\
'''  0
ENDTAB
'''
#}}}

#{{{ style

style_template =\
'''  0
TABLE
  2
STYLE
  5
3
330
0
100
 70
    33
  0
STYLE
  5
11
330
3
100
100
AcDbTextStyleTableRecord
  2
Standard
 70
     0
 40
0.0
 41
1.0
 50
0.0
 71
     0
 42
2.5
  3
txt
  4

  0
ENDTAB
'''

#}}}

#{{{ view

view_template =\
'''  0
TABLE
  2
VIEW
  5
6
330
0
100
 70
     0
  0
ENDTAB
'''

#}}}

#{{{ UCS, APPID, DIMSTYLE

misctable_template =\
'''  0
TABLE
  2
UCS
  5
7
330
0
100
 70
     0
  0
ENDTAB
  0
TABLE
  2
APPID
  5
9
330
0
100
 70
     1
  0
APPID
  5
12
330
9
100
100
AcDbRegAppTableRecord
  2
ACAD
 70
     0
  0
ENDTAB
  0
TABLE
  2
DIMSTYLE
  5
A
330
0
100
 70
     1
100
AcDbDimStyleTable
  0
DIMSTYLE
105
27
330
A
100
100
AcDbDimStyleTableRecord
  2
ISO-25
 70
     0
 41
2.5
 42
0.625
 43
3.75
 44
1.25
 73
     0
 74
     0
 77
     1
 78
     8
140
2.5
141
2.5
143
147
0.625
171
     3
172
     1
178
     0
271
     2
272
     2
274
     3
278
    44
283
     0
284
     8
340
11
  0
ENDTAB
'''

#}}}

#{{{ blockrecords

blockrecords_template =\
'''  0
TABLE
  2
BLOCK_RECORD
  5
1
330
0
100
 70
     1
  0
  5
1F
330
1
100
100
AcDbBlockTableRecord
  2
*Model_Space
340
22
 70
     0
280
     1
281
     0
  0
  5
1B
330
1
100
100
  2
*Paper_Space
340
1E
 70
     0
280
     1
281
     0
  0
  5
23
330
1
100
100
  2
*Paper_Space0
340
26
 70
     0
280
     1
281
     0
  0
ENDTAB
'''

#}}}

tables_template2 =\
'''    0
ENDSEC
'''

#}}}

blocks_template =\
'''  0
SECTION
  2
BLOCKS
  0
BLOCK
  5
20
330
1F
100
AcDbEntity
  8
0
100
AcDbBlockBegin
  2
 70
     0
 10
0.0
 20
0.0
 30
0.0
  3
  1

  0
ENDBLK
  5
21
330
1F
100
  8
0
100
AcDbBlockEnd
  0
BLOCK
  5
1C
330
1B
100
 67
     1
  8
0
100
  2
 70
     0
 10
0.0
 20
0.0
 30
0.0
  3
  1

  0
ENDBLK
  5
1D
330
1B
100
 67
     1
  8
0
100
  0
BLOCK
  5
24
330
23
100
  8
0
100
  2
 70
     0
 10
0.0
 20
0.0
 30
0.0
  3
  1

  0
ENDBLK
  5
25
330
23
100
  8
0
100
  0
ENDSEC
'''
#}}}

entities_template1 =\
'''  0
SECTION
  2
ENTITIES
'''
entities_template0 =\
'''  0
LINE
  5
7C
330
1F
100
  8
Default
100
AcDbLine
 10
1.0
 20
0.0
 30
0.0
 11
73.0
 21
41.0
 31
0.0
  0
ENDSEC
'''

#}}}

objects_template =\
'''  0
SECTION
  2
OBJECTS
  0
DICTIONARY
  5
C
330
0
100
AcDbDictionary
281
     1
  3
ACAD_GROUP
350
D
  3
ACAD_LAYOUT
350
1A
  3
ACAD_MATERIAL
350
3B
  3
ACAD_MLEADERSTYLE
350
6E
  3
ACAD_MLINESTYLE
350
17
  3
ACAD_PLOTSETTINGS
350
19
  3
ACAD_PLOTSTYLENAME
350
E
  3
ACAD_TABLESTYLE
350
6D
  3
ACAD_VISUALSTYLE
350
2A
  0
SUN
  5
3F
330
29
100
AcDbSun
 90
        1
290
     0
 63
     7
421
 16777215
 40
1.0
291
     1
 91
  2455826
 92
 54000000
292
     0
 70
     2
 71
   256
280
     1
  0
  5
D
102
{ACAD_REACTORS
330
C
102
}
330
C
100
281
     1
  0
  5
1A
102
330
C
102
}
330
C
100
281
     1
  3
Layout1
350
1E
  3
Layout2
350
26
  3
Model
350
22
  0
  5
3B
102
330
C
102
}
330
C
100
281
     1
  3
ByBlock
350
3D
  3
ByLayer
350
3C
  3
Global
350
3E
  0
  5
6E
102
330
C
102
}
330
C
100
281
     1
  0
  5
17
102
330
C
102
}
330
C
100
281
     1
  3
Standard
350
18
  0
  5
19
102
330
C
102
}
330
C
100
281
     1
  0
  5
E
102
330
C
102
}
330
C
100
281
     1
  3
Normal
350
F
100
340
F
  0
  5
6D
102
330
C
102
}
330
C
100
281
     1
  0
  5
2A
102
330
C
102
}
330
C
100
281
     1
  3
2dWireframe
350
2F
  3
3D Hidden
350
31
  3
3dWireframe
350
30
  3
Basic
350
32
  3
Brighten
350
36
  3
ColorChange
350
3A
  3
Conceptual
350
34
  3
Dim
350
35
  3
Facepattern
350
39
  3
Flat
350
2B
  3
FlatWithEdges
350
2C
  3
Gouraud
350
2D
  3
GouraudWithEdges
350
2E
  3
Linepattern
350
38
  3
Realistic
350
33
  3
Thicken
350
37
  0
LAYOUT
  5
1E
102
330
1A
102
}
330
1A
100
AcDbPlotSettings
  1

  2
none_device
  4

  6

 40
0.0
 41
0.0
 42
0.0
 43
0.0
 44
0.0
 45
0.0
 46
0.0
 47
0.0
 48
0.0
 49
0.0
140
0.0
141
0.0
142
1.0
143
1.0
 70
   688
 72
     0
 73
     0
 74
     5
  7

 75
    16
 76
     0
 77
     2
 78
   300
147
1.0
148
0.0
149
0.0
100
AcDbLayout
  1
Layout1
 70
     1
 71
     1
 10
0.0
 20
0.0
 11
420.0
 21
297.0
 12
0.0
 22
0.0
 32
0.0
 14
 24
 34
 15
 25
 35
146
0.0
 13
0.0
 23
0.0
 33
0.0
 16
1.0
 26
0.0
 36
0.0
 17
0.0
 27
1.0
 37
0.0
 76
     0
330
1B
  0
LAYOUT
  5
26
102
330
1A
102
}
330
1A
100
  1

  2
  4

  6

 40
0.0
 41
0.0
 42
0.0
 43
0.0
 44
0.0
 45
0.0
 46
0.0
 47
0.0
 48
0.0
 49
0.0
140
0.0
141
0.0
142
1.0
143
1.0
 70
   688
 72
     0
 73
     0
 74
     5
  7

 75
    16
 76
     0
 77
     2
 78
   300
147
1.0
148
0.0
149
0.0
100
  1
Layout2
 70
     1
 71
     2
 10
0.0
 20
0.0
 11
0.0
 21
0.0
 12
0.0
 22
0.0
 32
0.0
 14
0.0
 24
0.0
 34
0.0
 15
0.0
 25
0.0
 35
0.0
146
0.0
 13
0.0
 23
0.0
 33
0.0
 16
1.0
 26
0.0
 36
0.0
 17
0.0
 27
1.0
 37
0.0
 76
     0
330
23
  0
LAYOUT
  5
22
102
330
1A
102
}
330
1A
100
  1

  2
  4

  6

 40
0.0
 41
0.0
 42
0.0
 43
0.0
 44
0.0
 45
0.0
 46
0.0
 47
0.0
 48
0.0
 49
0.0
140
0.0
141
0.0
142
1.0
143
1.0
 70
  1712
 72
     0
 73
     0
 74
     0
  7

 75
     0
 76
     0
 77
     2
 78
   300
147
1.0
148
0.0
149
0.0
100
  1
Model
 70
     1
 71
     0
 10
0.0
 20
0.0
 11
420.0
 21
297.0
 12
0.0
 22
0.0
 32
0.0
 14
1.0
 24
0.0
 34
0.0
 15
73.0
 25
41.0
 35
0.0
146
0.0
 13
0.0
 23
0.0
 33
0.0
 16
1.0
 26
0.0
 36
0.0
 17
0.0
 27
1.0
 37
0.0
 76
     0
330
1F
331
29
  0
MATERIAL
  5
3D
102
330
3B
102
}
330
3B
100
  1
ByBlock
 72
     1
 94
      127
  0
MATERIAL
  5
3C
102
330
3B
102
}
330
3B
100
  1
ByLayer
 72
     1
 94
      127
  0
MATERIAL
  5
3E
102
330
3B
102
}
330
3B
100
  1
Global
 72
     1
 94
      127
  0
MLINESTYLE
  5
18
102
330
17
102
}
330
17
100
AcDbMlineStyle
  2
Standard
 70
     0
  3

 62
   256
 51
90.0
 52
90.0
 71
     2
 49
0.5
 62
   256
  6
BYLAYER
 49
-0.5
 62
   256
  6
BYLAYER
  0
ACDBPLACEHOLDER
  5
F
102
330
E
102
}
330
E
  0
  5
2F
102
330
2A
102
}
330
2A
100
  2
 70
     4
177
     2
291
     0
 71
     0
176
     1
 72
     2
176
     1
 73
     0
176
     1
 90
        0
176
     1
 40
-0.6
176
     1
 41
-30.0
176
     1
 62
     5
 63
     7
421
 16777215
176
     1
 74
     1
176
     1
 91
        4
176
     1
 64
     7
176
     1
 65
   257
176
     1
 75
     1
176
     1
175
     1
176
     1
 42
1.0
176
     1
 92
        0
176
     1
 66
   257
176
     1
 43
1.0
176
     1
 76
     1
176
     1
 77
     6
176
     1
 78
     2
176
     1
 67
     7
176
     1
 79
     5
176
     1
170
     0
176
     1
171
     0
176
     1
290
     0
176
     1
 93
        1
176
     1
 44
0.0
176
     1
173
     0
176
     1
  0
  5
31
102
330
2A
102
}
330
2A
100
  2
3D Hidden
 70
     6
177
     2
291
     0
 71
     1
176
     1
 72
     2
176
     1
 73
     2
176
     1
 90
        0
176
     1
 40
-0.6
176
     1
 41
-30.0
176
     1
 62
     5
 63
     7
421
 16777215
176
     1
 74
     2
176
     1
 91
        2
176
     1
 64
     7
176
     1
 65
   257
176
     1
 75
     2
176
     1
175
     1
176
     1
 42
40.0
176
     1
 92
        0
176
     1
 66
   257
176
     1
 43
1.0
176
     1
 76
     1
176
     1
 77
     6
176
     1
 78
     2
176
     1
 67
     7
176
     1
 79
     3
176
     1
170
     0
176
     1
171
     0
176
     1
290
     0
176
     1
 93
        1
176
     1
 44
0.0
176
     1
173
     0
176
     1
  0
  5
30
102
330
2A
102
}
330
2A
100
  2
 70
     5
177
     2
291
     0
 71
     0
176
     1
 72
     2
176
     1
 73
     0
176
     1
 90
        0
176
     1
 40
-0.6
176
     1
 41
-30.0
176
     1
 62
     5
 63
     7
421
 16777215
176
     1
 74
     1
176
     1
 91
        4
176
     1
 64
     7
176
     1
 65
   257
176
     1
 75
     1
176
     1
175
     1
176
     1
 42
1.0
176
     1
 92
        0
176
     1
 66
   257
176
     1
 43
1.0
176
     1
 76
     1
176
     1
 77
     6
176
     1
 78
     2
176
     1
 67
     7
176
     1
 79
     5
176
     1
170
     0
176
     1
171
     0
176
     1
290
     0
176
     1
 93
        1
176
     1
 44
0.0
176
     1
173
     0
176
     1
  0
  5
32
102
330
2A
102
}
330
2A
100
  2
Basic
 70
     7
177
     2
291
     1
 71
     1
176
     1
 72
     0
176
     1
 73
     1
176
     1
 90
        0
176
     1
 40
-0.6
176
     1
 41
-30.0
176
     1
 62
     5
 63
     7
421
 16777215
176
     1
 74
     0
176
     1
 91
        4
176
     1
 64
     7
176
     1
 65
   257
176
     1
 75
     1
176
     1
175
     1
176
     1
 42
1.0
176
     1
 92
        8
176
     1
 66
     7
176
     1
 43
1.0
176
     1
 76
     1
176
     1
 77
     6
176
     1
 78
     2
176
     1
 67
     7
176
     1
 79
     5
176
     1
170
     0
176
     1
171
     0
176
     1
290
     0
176
     1
 93
        1
176
     1
 44
0.0
176
     1
173
     0
176
     1
  0
  5
36
102
330
2A
102
}
330
2A
100
  2
Brighten
 70
    12
177
     2
291
     1
 71
     2
176
     1
 72
     2
176
     1
 73
     0
176
     1
 90
        0
176
     1
 40
-0.6
176
     1
 41
-30.0
176
     1
 62
     5
 63
     7
421
 16777215
176
     1
 74
     1
176
     1
 91
        4
176
     1
 64
     7
176
     1
 65
   257
176
     1
 75
     1
176
     1
175
     1
176
     1
 42
1.0
176
     1
 92
        8
176
     1
 66
     7
176
     1
 43
1.0
176
     1
 76
     1
176
     1
 77
     6
176
     1
 78
     2
176
     1
 67
     7
176
     1
 79
     5
176
     1
170
     0
176
     1
171
     0
176
     1
290
     0
176
     1
 93
        1
176
     1
 44
50.0
176
     1
173
     0
176
     1
  0
  5
3A
102
330
2A
102
}
330
2A
100
  2
 70
    16
177
     2
291
     1
 71
     2
176
     1
 72
     2
176
     1
 73
     3
176
     1
 90
        0
176
     1
 40
-0.6
176
     1
 41
-30.0
176
     1
 62
     5
 63
     8
421
  8421504
176
     1
 74
     1
176
     1
 91
        4
176
     1
 64
     7
176
     1
 65
   257
176
     1
 75
     1
176
     1
175
     1
176
     1
 42
1.0
176
     1
 92
        8
176
     1
 66
     8
424
  8421504
176
     1
 43
1.0
176
     1
 76
     1
176
     1
 77
     6
176
     1
 78
     2
176
     1
 67
     7
176
     1
 79
     5
176
     1
170
     0
176
     1
171
     0
176
     1
290
     0
176
     1
 93
        1
176
     1
 44
0.0
176
     1
173
     0
176
     1
  0
  5
34
102
330
2A
102
}
330
2A
100
  2
 70
     9
177
     2
291
     0
 71
     3
176
     1
 72
     2
176
     1
 73
     0
176
     1
 90
        0
176
     1
 40
-0.6
176
     1
 41
-30.0
176
     1
 62
     5
 63
     7
421
 16777215
176
     1
 74
     2
176
     1
 91
        2
176
     1
 64
     7
176
     1
 65
   257
176
     1
 75
     1
176
     1
175
     1
176
     1
 42
40.0
176
     1
 92
        8
176
     1
 66
     7
176
     1
 43
1.0
176
     1
 76
     1
176
     1
 77
     6
176
     1
 78
     2
176
     1
 67
     7
176
     1
 79
     3
176
     1
170
     0
176
     1
171
     0
176
     1
290
     0
176
     1
 93
        1
176
     1
 44
0.0
176
     1
173
     0
176
     1
  0
  5
35
102
330
2A
102
}
330
2A
100
  2
Dim
 70
    11
177
     2
291
     1
 71
     2
176
     1
 72
     2
176
     1
 73
     0
176
     1
 90
        0
176
     1
 40
-0.6
176
     1
 41
-30.0
176
     1
 62
     5
 63
     7
421
 16777215
176
     1
 74
     1
176
     1
 91
        4
176
     1
 64
     7
176
     1
 65
   257
176
     1
 75
     1
176
     1
175
     1
176
     1
 42
1.0
176
     1
 92
        8
176
     1
 66
     7
176
     1
 43
1.0
176
     1
 76
     1
176
     1
 77
     6
176
     1
 78
     2
176
     1
 67
     7
176
     1
 79
     5
176
     1
170
     0
176
     1
171
     0
176
     1
290
     0
176
     1
 93
        1
176
     1
 44
-50.0
176
     1
173
     0
176
     1
  0
  5
39
102
330
2A
102
}
330
2A
100
  2
 70
    15
177
     2
291
     1
 71
     2
176
     1
 72
     2
176
     1
 73
     0
176
     1
 90
        0
176
     1
 40
-0.6
176
     1
 41
-30.0
176
     1
 62
     5
 63
     7
421
 16777215
176
     1
 74
     1
176
     1
 91
        4
176
     1
 64
     7
176
     1
 65
   257
176
     1
 75
     1
176
     1
175
     1
176
     1
 42
1.0
176
     1
 92
        8
176
     1
 66
     7
176
     1
 43
1.0
176
     1
 76
     1
176
     1
 77
     6
176
     1
 78
     2
176
     1
 67
     7
176
     1
 79
     5
176
     1
170
     0
176
     1
171
     0
176
     1
290
     0
176
     1
 93
        1
176
     1
 44
0.0
176
     1
173
     0
176
     1
  0
  5
2B
102
330
2A
102
}
330
2A
100
  2
Flat
 70
     0
177
     2
291
     1
 71
     2
176
     1
 72
     1
176
     1
 73
     1
176
     1
 90
        2
176
     1
 40
-0.6
176
     1
 41
30.0
176
     1
 62
     5
 63
     7
421
 16777215
176
     1
 74
     0
176
     1
 91
        4
176
     1
 64
     7
176
     1
 65
   257
176
     1
 75
     1
176
     1
175
     1
176
     1
 42
1.0
176
     1
 92
        8
176
     1
 66
     7
176
     1
 43
1.0
176
     1
 76
     1
176
     1
 77
     6
176
     1
 78
     2
176
     1
 67
     7
176
     1
 79
     5
176
     1
170
     0
176
     1
171
     0
176
     1
290
     0
176
     1
 93
       13
176
     1
 44
0.0
176
     1
173
     0
176
     1
  0
  5
2C
102
330
2A
102
}
330
2A
100
  2
 70
     1
177
     2
291
     1
 71
     2
176
     1
 72
     1
176
     1
 73
     1
176
     1
 90
        2
176
     1
 40
-0.6
176
     1
 41
30.0
176
     1
 62
     5
 63
     7
421
 16777215
176
     1
 74
     1
176
     1
 91
        4
176
     1
 64
     7
176
     1
 65
   257
176
     1
 75
     1
176
     1
175
     1
176
     1
 42
1.0
176
     1
 92
        0
176
     1
 66
   257
176
     1
 43
1.0
176
     1
 76
     1
176
     1
 77
     6
176
     1
 78
     2
176
     1
 67
     7
176
     1
 79
     5
176
     1
170
     0
176
     1
171
     0
176
     1
290
     0
176
     1
 93
       13
176
     1
 44
0.0
176
     1
173
     0
176
     1
  0
  5
2D
102
330
2A
102
}
330
2A
100
  2
Gouraud
 70
     2
177
     2
291
     1
 71
     2
176
     1
 72
     2
176
     1
 73
     1
176
     1
 90
        2
176
     1
 40
-0.6
176
     1
 41
30.0
176
     1
 62
     5
 63
     7
421
 16777215
176
     1
 74
     0
176
     1
 91
        4
176
     1
 64
     7
176
     1
 65
   257
176
     1
 75
     1
176
     1
175
     1
176
     1
 42
1.0
176
     1
 92
        0
176
     1
 66
     7
176
     1
 43
1.0
176
     1
 76
     1
176
     1
 77
     6
176
     1
 78
     2
176
     1
 67
     7
176
     1
 79
     5
176
     1
170
     0
176
     1
171
     0
176
     1
290
     0
176
     1
 93
       13
176
     1
 44
0.0
176
     1
173
     0
176
     1
  0
  5
2E
102
330
2A
102
}
330
2A
100
  2
 70
     3
177
     2
291
     1
 71
     2
176
     1
 72
     2
176
     1
 73
     1
176
     1
 90
        2
176
     1
 40
-0.6
176
     1
 41
30.0
176
     1
 62
     5
 63
     7
421
 16777215
176
     1
 74
     1
176
     1
 91
        4
176
     1
 64
     7
176
     1
 65
   257
176
     1
 75
     1
176
     1
175
     1
176
     1
 42
1.0
176
     1
 92
        0
176
     1
 66
   257
176
     1
 43
1.0
176
     1
 76
     1
176
     1
 77
     6
176
     1
 78
     2
176
     1
 67
     7
176
     1
 79
     5
176
     1
170
     0
176
     1
171
     0
176
     1
290
     0
176
     1
 93
       13
176
     1
 44
0.0
176
     1
173
     0
176
     1
  0
  5
38
102
330
2A
102
}
330
2A
100
  2
 70
    14
177
     2
291
     1
 71
     2
176
     1
 72
     2
176
     1
 73
     0
176
     1
 90
        0
176
     1
 40
-0.6
176
     1
 41
-30.0
176
     1
 62
     5
 63
     7
421
 16777215
176
     1
 74
     1
176
     1
 91
        4
176
     1
 64
     7
176
     1
 65
   257
176
     1
 75
     7
176
     1
175
     7
176
     1
 42
1.0
176
     1
 92
        8
176
     1
 66
     7
176
     1
 43
1.0
176
     1
 76
     1
176
     1
 77
     6
176
     1
 78
     2
176
     1
 67
     7
176
     1
 79
     5
176
     1
170
     0
176
     1
171
     0
176
     1
290
     0
176
     1
 93
        1
176
     1
 44
0.0
176
     1
173
     0
176
     1
  0
  5
33
102
330
2A
102
}
330
2A
100
  2
Realistic
 70
     8
177
     2
291
     0
 71
     2
176
     1
 72
     2
176
     1
 73
     0
176
     1
 90
        0
176
     1
 40
-0.6
176
     1
 41
-30.0
176
     1
 62
     5
 63
     7
421
 16777215
176
     1
 74
     1
176
     1
 91
        0
176
     1
 64
     7
176
     1
 65
   257
176
     1
 75
     1
176
     1
175
     1
176
     1
 42
1.0
176
     1
 92
        8
176
     1
 66
     8
424
  7895160
176
     1
 43
1.0
176
     1
 76
     1
176
     1
 77
     6
176
     1
 78
     2
176
     1
 67
     7
176
     1
 79
     5
176
     1
170
     0
176
     1
171
     0
176
     1
290
     0
176
     1
 93
       13
176
     1
 44
0.0
176
     1
173
     0
176
     1
  0
  5
37
102
330
2A
102
}
330
2A
100
  2
Thicken
 70
    13
177
     2
291
     1
 71
     2
176
     1
 72
     2
176
     1
 73
     0
176
     1
 90
        0
176
     1
 40
-0.6
176
     1
 41
-30.0
176
     1
 62
     5
 63
     7
421
 16777215
176
     1
 74
     1
176
     1
 91
        4
176
     1
 64
     7
176
     1
 65
   257
176
     1
 75
     1
176
     1
175
     1
176
     1
 42
1.0
176
     1
 92
       12
176
     1
 66
     7
176
     1
 43
1.0
176
     1
 76
     1
176
     1
 77
     6
176
     1
 78
     2
176
     1
 67
     7
176
     1
 79
     5
176
     1
170
     0
176
     1
171
     0
176
     1
290
     0
176
     1
 93
        1
176
     1
 44
0.0
176
     1
173
     0
176
     1
  0
ENDSEC
'''

#}}}

#}}}

#{{{ Test

def test_func():
    d = DXF('test.dxf')
    d.add_layer('ABC', color=5)
    d.add_entity(Line((15,5), (56,-89)), 'ABC')
    x = (2.9, 5, 9, 8.0, 20)
    y = (-5,-12.0, -3, 6, 18)
    d.add_entity(LwPolyLine(x, y), 'ABC')
    d.add_entity(Rectangle((-50,50),200, 100), 'ABC')    
    d.add_entity(Circle((-10,10), 40), 'ABC')
    d.add_entity(Arc((100,50), 40, 30, 200), 'ABC')
    d.add_entity(Text('This is my test DXF file.',(-50,-50), 5.0, -30), 'Text')    
    
#}}}

'''
gtrace.nonsequential

A module to perform non-sequential trace of a beam
in an optical system.
'''

# are met:
#
#
#

#}}}

#{{{ non_seq_trace

def non_seq_trace(optList, src_beam, order=10, power_threshold=0.1):
    '''
    Perform non-sequential trace of the source beam, src_beam,
    through the optical system represented by a collection of optics,
    optList.

    The return value of this function is a list of beams.
    '''

    #Loop over all the optics to see if the source beam hit them.
    #Then select the closest optics being hit.
    min_dist = 1e15
    final_answer = None
    hit_optics = None
    for opt in optList:
        #See if the beam hit the optics
        ans = opt.isHit(src_beam)
        #If the beam hit the optics
        if ans['isHit']:
            #If the intersection point is closest one so far
            if min_dist > ans['distance']:
                    hit_optics = opt

    if final_answer is None:
        #The beam does not hit any optics
        return [src_beam]

    if final_answer['face'] == 'side':
        #The beam is terminated on a side of an optics
        src_beam.length = final_answer['distance']

    #The beam hits an actual optical surface

    #Avoid forming a cavity
    if hit_optics.term_on_HR  and final_answer['face'] == 'HR' and src_beam.stray_order == 0:
    
    ans = hit_optics.hit(src_beam, order=order, threshold=power_threshold,
                         face=final_answer['face'])

    terminated_beam_list = [b for b in ans[1].values() if b.incSurfAngle is not None]
    open_beam_list = [b for b in ans[1].values() if b.incSurfAngle is None]

    #For each open beam, perform the non sequential trace
    for b in open_beam_list:
        beams = non_seq_trace(optList=optList, src_beam=b.copy(), order=order,
                              power_threshold=power_threshold)
        terminated_beam_list.extend(beams)

    return terminated_beam_list
        
#}}}

#{{{ Imports

pi = np.pi

#}}}

#{{{ Length units
km = 1e3
m = 1.0
cm = 0.01
mm = 1e-3
um = 1e-6
nm = 1e-9
#}}}

#{{{ Power units

kW = 1e3
W = 1.0
mW = 1e-3
uW = 1e-6

#}}}

#{{{ Frequencies

THz = 1e12
GHz = 1e9
MHz = 1e6
kHz = 1e3
Hz = 1.0
mHz = 1e-3
uHz = 1e-6

#}}}

#{{{ Other

ppm = 1e-6

#}}}

#{{{ Conversion Functions

def rad2deg(rad):
    return rad * 180.0 / pi

def deg2rad(deg):
    return deg * pi /180.0

#}}}

pi = np.pi

#}}}

#Copyright (c) 2011, Yoichi Aso

# are met:
#
#
#

__copyright__ = "Copyright 2011, Yoichi Aso"

#}}}

km = 1e3
m = 1.0
cm = 0.01
mm = 1e-3
um = 1e-6
nm = 1e-9
#}}}

kW = 1e3
W = 1.0
mW = 1e-3
uW = 1e-6

#}}}

THz = 1e12
GHz = 1e9
MHz = 1e6
kHz = 1e3
Hz = 1.0
mHz = 1e-3
uHz = 1e-6

#}}}

#{{{ Other

ppm = 1e-6

#}}}

#}}}

'''
Drawing classes for gtrace
'''

pi = np.pi
#}}}

# are met:
#
#
#

#}}}

#{{{ Canvas class

class Canvas(object):
    '''
    Canvas class
    '''

    def __init__(self, unit='m'):
        self.layers = {}
        self.unit = unit

    def add_layer(self, name, color=(0,0,0)):
        self.layers[name] = Layer(name, color=color)

    def add_shape(self, shape, layername):

        self.layers[layername].add_shape(shape)

#}}}

#{{{ Layer class

    '''
    '''

    def __init__(self, name, color=(0,0,0)):
        self.shapes = []

    def add_shape(self, shape):
        self.shapes.append(shape)
        
#}}}

#{{{ Shape class

class Shape(object):
    '''
    Shape class
    '''

        pass
        
#}}}

#{{{ Line class

class Line(Shape):
    '''
    Line class
    '''

#}}}

#{{{ PolyLine

class PolyLine(Shape):
    '''
    '''
    
        '''
        '''
        super(PolyLine, self).__init__()
        self.x = x
        self.y = y

#}}}

class Rectangle(Shape):
    '''
    '''
    
        '''
        point: lower left corner of the rectangle
        width:
        height:
        '''

#}}}

class Circle(Shape):
    '''
    A circle
    '''

#}}}

#{{{ Arc

class Arc(Shape):
    '''
    An arc

    Note that angles are stored in rad.
    '''

    def __init__(self, center, radius, startangle, stopangle, thickness=0, angle_in_rad=True):
        if angle_in_rad:
            self.startangle = startangle
            self.stopangle = stopangle
        else:
            self.startangle = pi*startangle/180.0
            self.stopangle = pi*stopangle/180.0
            
#}}}

#{{{ Text

class Text(Shape):
    '''
    Text

    Note that angles are stored in rad.    
    '''

    def __init__(self, text, point, height=1.0, rotation=0.0, angle_in_rad=True):
        if angle_in_rad:        
            self.rotation = rotation
        else:
            self.rotation = pi*rotation/180.0

#}}}

"""
gaussian - Gaussian Optics Module

This module contains several utility functions for gaussian optics.
"""

pi = np.pi
import scipy.special as spf

#}}}

#{{{ modeSpacing

def modeSpacing(g1, g2):
    return np.arccos(np.sqrt(g1*g2))/pi

#}}}

#{{{ q-parameter related functions

def q2zr(q):
    '''
    Convert a q-parameter to Rayleigh range.
    '''
    zr = np.float(np.imag(q))
    return zr

def q2w(q, wl=1064*nm):
    '''
    Convert a q-parameter to the beam size
    '''
    S = -1.0/np.imag(1.0/q)
    w = np.sqrt(wl*S/pi)
    return w

def q2R(q):
    '''
    Convert a q-parameter to the ROC
    '''
    return 1.0/np.real(1.0/q)

def Rw2q(ROC=1.0, w=1.0, wl=1064e-9):
    '''
    Get the q-parameter from the ROC and w.
    '''
    k = 2.0*pi/wl
    S = w**2 * k/2

    return 1.0/(1.0/ROC + 1.0/(1j*S))

def zr2w0(zr, wl=1064*nm):
    '''
    Convert Rayleigh range to the waist size
    '''
    return np.sqrt(2*zr*wl/(2*pi))

def w02zr(w0, wl=1064*nm):
    '''
    '''
    return (2*pi/wl)*(w0**2)/2

def modeMatching(q1, q2x, q2y=False):
    '''
    Mode matching between two beams with different q-parameters.
    The axes of the two beams are assumed to be matched.

    q1: q-parameter of the first beam. This beam is assumed to be circular.

    q2x: q-parameter of the second beam in x-direction. If the second beam
         is also circular, omit the next argument.

    q2y: q-parameter of the second beam in y-direction. Specify this parameter
         if the second beam is eliptic.
    '''
    
    zr1 = np.imag(q1)
    d1 = np.real(q1)

    if q2y:
        #Eliptic beam
        zrx = np.imag(q2x)
        dx = np.real(q2x)
        zry = np.imag(q2y)
        dy = np.real(q2y)
        return np.abs(2*sqrt(zr1*sqrt(zrx*zry)/((zr1+zrx+1j*(dx-d1))*(zr1+zry+1j*(dy-d1)))))**2
    else:
        #Circular beam
        zr2 = np.imag(q2x)
        d2 = np.real(q2x)
        ec = zr2 - zr1
        az = d2 -d1

        return np.abs(2*zr1*sqrt(1+ec/zr1)/(2*zr1+ec+1j*az))**2

def optimalMatching(q1, q2):
    '''
    Returns a mode (q-parameter) which best matches the given
    two q-parameters, q1 and q2.

    Returned values: (q, match)

    q: The best matching q-parameter

    match: Mode matching rate
    '''

    zr2 = np.imag(q2)
    d2 = np.real(q2)

    zr = np.sqrt(zr1*zr2)*np.sqrt(1+((d2-d1)/(zr1+zr2))**2)
    d = (zr1*d2 + zr2*d1)/(zr1+zr2)

    q = d +1j*zr
    match = modeMatching(q, q1, q2)

    return (q, match)

#{{{ For compatibility

def qToRadius(q, wl=1064e-9):
    k = 2*pi/wl
    return sqrt(-2.0/(k*np.imag(1.0/q)))

def qToROC(q):

def ROCandWtoQ(ROC=1.0, w=1.0, wl=1064e-9):

#}}}
    
#}}}

#{{{ Beam clip

def beamClip(a = 1.0,  w = 3.0):
    return (1+spf.erf(np.sqrt(2)*a/w))/2

def appertureCut(r = 1.0,  w = 3.0):
    return 1-np.exp(-2*(r/w)**2)

#}}}

import optparse

class hello():
    def greetings(self):
        print('hello')

    def __call__(self):
        self.greetings()

def saysomething():
    parser = optparse.OptionParser(usage='-v')
    parser.add_option("-v", "--verbose", dest="verbose",
                  action="store_true", default=False)
    (options, args) = parser.parse_args()
    hello()()
    if options.verbose:
        print('Yes')

#}}}

#{{{ Index of refractions

def sellmeier(wl, B1, B2, B3, C1, C2, C3):
    '''
    Calculate index of refraction using Sellmeiers equation

    n^2 = 1+B1*wl^2/(wl^2 - C1) + B2*wl^2/(wl^2 - C2) + B3*wl^2/(wl^2 - C3)

    See below for the coefficients for specific materials.
    http://www.cvimellesgriot.com/products/Documents/Catalog/Dispersion_Equations.pdf
    
    '''

    n = np.sqrt(1+B1*wl**2/(wl**2 - C1) + B2*wl**2/(wl**2 - C2) + B3*wl**2/(wl**2 - C3))
    return n
    
def n_fused_silica(wl):
    '''
    Calculate the index of refraction of fused silica for a given wavelength.
    '''
    #Convert the wavelength to the unit of um.
    wl = wl*1e6
    
    B1 = 6.961663e-1
    B2 = 4.079426e-1
    B3 = 8.974794e-1
    C1 = 4.67914826e-3
    C2 = 1.35120631e-2
    C3 = 9.79340025e1
    
    return sellmeier(wl, B1, B2, B3, C1, C2, C3)

def n_sapphire_ordinary(wl):
    '''
    Calculate the index of refraction of Sapphire ordinary axis for a given wavelength.
    '''
    
    wl = wl*1e6
    
    B1 = 1.4313493
    B2 = 6.5054713e-1
    B3 = 5.3414021
    C1 = 5.2799261e-3
    C2 = 1.42382647e-2
    C3 = 3.25017834e2
    
def n_sapphire_extraordinary(wl):
    '''
    Calculate the index of refraction of Sapphire extraordinary axis for a given wavelength.
    '''
    
    wl = wl*1e6
    
    B1 = 1.5039759
    B2 = 5.5069141e-1
    B3 = 6.5927379
    C1 = 5.48041129e-3
    C2 = 1.47994281e-2
    C3 = 4.0289514e2
    
n_fused_silica_532nm = 1.46071
n_fused_silica_1064nm = 1.45
n_fused_silica_1550nm = 1.444

no_sapphire_532nm = 1.77170
ne_sapphire_532nm = 1.76355

no_sapphire_1064nm = 1.75449
ne_sapphire_1064nm = 1.74663

no_sapphire_1550nm = 1.74618
ne_sapphire_1550nm = 1.73838

#}}}

"""
cavity.py - A Cavity class and related functions for representing a Fabry-Perot cavity
"""

pi = np.pi
c = 299792458.0
import gtrace.optics.gaussian as gauss
from enthought.traits.api import CFloat, HasTraits
#}}}

#{{{ Class Cavity

class Cavity(HasTraits):
    """
    A class to represent a Fabry-Perot cavity.

    *Attributes*
    
    **r1:** Input mirror reflectivity (amplitude)
    
    **r2:** End mirror reflectivity (amplitude)
    
    **rp1:** Input mirror reflectivity (power)
    
    **rp2:** End mirror reflectivity (power)
    
    **L:** Length

    **R1:** ROC of the input mirror (positive when concave to incident light, i.e.
        convex seen from inside the cavity)

    **R2:** ROC of the end mirror (positive when concave to incident light, i.e.
        concave seen from inside the cavity)

    **wl:** Wavelength
    """
    r1 = CFloat(0.9)
    r2 = CFloat(0.99)
    rp1 = CFloat()
    rp2 = CFloat()
    tp1 = CFloat()
    tp2 = CFloat()
    t1 = CFloat()
    t2 = CFloat()

    L = CFloat(1.0)
    R1 = CFloat(-1.5)
    R2 = CFloat(1.5)
    g1 = CFloat()
    g2 = CFloat()
    wl=CFloat()
    
    def __init__(self, r1=0.9, r2=0.99, L=1.0, R1=-1.5, R2=1.5, wl=1064e-9, power=False):
        """
        **r1:** Input mirror reflectivity 

        **r2:** End mirror reflectivity 

        **L:** Length

        **R1:** ROC of the input mirror (positive when concave to incident light, i.e.
            convex seen from inside the cavity)

        **R2:** ROC of the end mirror (positive when concave to incident light, i.e.
            concave seen from inside the cavity)

        **wl:** Wavelength

        **power:** If True, r1 and r2 are treated as power reflectivities.
               Otherwise, r1 and r2 are regarded as amplitude reflectivities.
        """
        if power:
            self.r1 = np.sqrt(r1)
            self.r2 = np.sqrt(r2)
        else:
            self.r1 = r1
            self.r2 = r2

        self.L = L
        self.R1 = R1
        self.R2 = R2
        self.wl = wl

    #{{{ Trait change handlers
    def _R1_changed(self, old, new):
        self.g1 = 1 + self.L/new

    def _R2_changed(self, old, new):
        self.g2 = 1 - self.L/new

    def _L_changed(self, old, new):
        self.g1 = 1 + new/self.R1
        self.g2 = 1 - new/self.R2

    def _r1_changed(self, old, r1):
        self.rp1 = r1**2
        self.tp1 = 1 - self.rp1
        self.t1 = np.sqrt(self.tp1)

    def _r2_changed(self, old, r2):
        self.rp2 = r2**2
        self.tp2 = 1 - self.rp2
        self.t2 = np.sqrt(self.tp2)

    #}}}

    def finesse(self):
        '''
        Returns the finesse of the cavity.
        '''
        return finesse(self.r1, self.r2)

    def storageTime(self):
        '''
        Storage time
        '''

        return 2*self.L*self.finesse()/(2*pi*c)

    def pole(self):
        '''
        Cavity pole frequency [Hz]
        '''

        return 1/(2*pi*self.storageTime())

    def Nbounce(self):
        '''
        Bounce number
        '''

        return 2*self.finesse()/pi

    def powerGain(self):
        '''
        Ratio of the intra-cavity power to the input power.
        '''
        return (self.t1/(1-self.r1*self.r2))**2

    def FSR(self):
        '''
        Returns the free spectral range of the cavity.
        '''
        return c/(2*self.L)

    def modeSpacing(self):
        '''
        Return the transverse mode spacing of the cavity (commonly called gamma).
        It is a fractional number defined by gamma = (mode spacing frequency)/FSR.
        
        '''
        return gauss.modeSpacing(self.g1, self.g2)

    def waist(self, size=False):
        """
        Return the q-parameter or the radius of the beam at the cavity waist.

        *Input arguments*
        
        **size:** (optional) if set to true, the first element of the returned tuple will be the waist size, rather than the q-parameter.

        *==Returned parameters ==*
        
        **(q0, d):**  This function returns a tuple with two elements.
                       The first element is the q-parameter of the cavity mode at
                       the cavity waist. If size=True is given, it becomes the waist
                       size (1/e^2 radius).
                       The second element is the distance of the cavity waist from
                       the input mirror.
        """
        # q0 = 1j*np.sqrt(self.L)*np.sqrt(-(self.L+self.R1)*(self.L-self.R2)\
        #                                 *(self.L+self.R1-self.R2))\
        #                                 /(2*self.L+self.R1-self.R2)

        q0 = 1j*sqrt(-self.L*(self.L+self.R1)**2*(self.L-self.R2)*\
                         (self.L+self.R1-self.R2)/((-2*self.L-self.R1+self.R2)**2))\
                         /sqrt(self.L+self.R1)
        
        d = self.L*(self.L-self.R2)/(2*self.L+self.R1-self.R2)

        if size:
            return (gauss.q2w(q0, wl=self.wl), d)
        else:
            return (q0, d)
        
    def spotSize(self):
        """
        Returns the beam spot sizes on the input and end mirrors
        as a tuple (w1,w2).
        """
        (q0,d) = self.waist()
        w1=gauss.q2w(q0-d, wl=self.wl)
        w2=gauss.q2w(q0+self.L-d, wl=self.wl)
        return (w1, w2)
        
    def trans(self, f=0, d=0):
        """
        Returns the amplitude transmissivity of the cavity.
        It assumes the cavity was locked to the incident light first. Then computes the
        amplitude transmissivity for the light with a frequency shift f from the original light
        with the cavity length changed by d from the initial state.

        *== Input arguments ==*
        
        **f:** Frequency shift of the light in Hz.

        **d:** Cavity length detuning in m.

        *== Returned parameter ==*
        
        The amplitude transmissivity of the cavity (a complex number).
        
        """
        #One way phase change
        phi=2*pi*(self.L*f/c + d/self.wl + d*f/c)
        
        return self.t1 * self.t2 * np.exp(-1j*phi)/\
                   (1-self.r1*self.r2*np.exp(-2j*phi))
    
    def refl(self, f=0, d=0):
        """
        Returns the amplitude reflectivity of the cavity.
        amplitude reflectivity for the light with a frequency shift f from the original light

        The amplitude reflectivity of the cavity (a complex number).
        
        """

        return -self.r1+self.t1**2 * self.r2 *np.exp(-2j*phi)/\
    
    def intra(self, f=0, d=0):
        """
        Returns the intra cavity field amplitude.
        intra-cavity field amplitude for the light with a frequency shift f from the original light

        The intra-cavity field amplitude at the input mirror surface (a complex number).
        
        """
        phi=2*pi*((self.L+d)*f/c + d/self.wl)
        
        return self.t1/(1-self.r1*self.r2*np.exp(-2j*phi))

#}}}

#{{{ Often used parameters

def finesse(r1, r2, power=False):
    '''
    Returns the finesse of a cavity
    '''
    if power:
        r1 = np.sqrt(r1)
        r2 = np.sqrt(r2)
        
    return pi*np.sqrt(r1*r2)/(1-r1*r2)

#}}}

#{{{ Import
pi = np.pi
#}}}

#{{{ Snell's Law

def deflection_angle(theta, n1, n2, deg=True):
    if deg:
        factor = pi/180.0
    else:
        factor = 1.0
        
    return np.arcsin(n1*np.sin(theta*factor)/n2)/factor

#}}}

#{{{ Geometry utilities

#{{{ line_plane_intersection

def line_plane_intersection(pos,
                            dirVect,
                            plane_center,
                            normalVector,
                            diameter):
    '''
    Compute the intersection point between a line
    and a plane

    A line is specified by its origin (pos) and the
    direction vector (dirVect).
    A plane is specfied by its center coordinates (plane_center)
    and the normal vector (normalVector).
    The plane has its size (diameter).

    The returned value is a dictionary of with the following keys:
    "Intersection Point": numpy array of the coordinates of the intersection point.
    "isHit": A boolean value of whether the line intersects with the plane or not.
    "distance": Distance between the origin of the line and the intersection point.
    "distance from center": Distance between the center of the plane and the intersection point.
    '''

    #Make sure the inputs are ndarrays
    pos = np.array(pos, dtype=np.float64)
    dirVect = np.array(dirVect, dtype=np.float64)
    plane_center = np.array(plane_center, dtype=np.float64)
    normalVector = np.array(normalVector, dtype=np.float64)
    diameter = float(diameter)
    
    #Get a normalized vector along the plane
    plVect = np.array([-normalVector[1], normalVector[0]])
    plVect = plVect/np.linalg.norm(plVect)

    #Normalize
    dirVect = dirVect/np.linalg.norm(dirVect)

    #Make sure that the plVect and dirVect are not parallel
    if np.abs(np.dot(dirVect, plVect)) > 1 - 1e-10:
        return {'Intersection Point': np.array((0.,0.)), 'isHit': False,
            'distance': 0.0,
            'distance from center': 0.0}

    #Solve line equations to get the intersection point
    M = np.vstack((dirVect, -plVect)).T
    ans = np.linalg.solve(M, plane_center - pos)    
    intersection_point = pos + ans[0]*dirVect

    #How far the intersection point is from the center
    #of the plane
    dist_from_center = np.abs(ans[1])
    if dist_from_center > diameter/2.0\
           or ans[0] < 0.\
           or np.dot(dirVect, normalVector) > 0.:
        
        hit = False
    else:
        hit = True
        
    return {'Intersection Point': intersection_point, 'isHit': hit,
            'distance': np.abs(ans[0]),
            'distance from center': ans[1]}

#}}}

#{{{ line_arc_intersection

def line_arc_intersection(pos,
                          dirVect,
                          chord_center,
                          chordNormVect,
                          invROC,
                          diameter,
                          verbose=False):
    '''
    and an arc.
    
    pos: Origin of the line
    dirVect: Direction of the line
    chord_center: The center of the chord made by the arc.
    chordNormVect: Normal vector of the chord.
    invROC: Inverse of the ROC of the arc. Positive for concave surface.
    diameter: Length of the chord.
    '''
    chord_center = np.array(chord_center, dtype=np.float64)
    chordNormVect = np.array(chordNormVect, dtype=np.float64)
    invROC = float(invROC)

    chordNormVect = chordNormVect/np.linalg.norm(chordNormVect)

    #Check if the ROC is too large.
    if np.abs(invROC) < 1e-5:
        #It is almost a plane
        ans = line_plane_intersection(pos, dirVect, chord_center, chordNormVect, diameter)
        localNormVect = chordNormVect
        localNormAngle = np.mod(np.arctan2(localNormVect[1],
                                           localNormVect[0]), 2*pi)

        ans['localNormVect'] = localNormVect
        ans['localNormAngle'] = localNormAngle

        return ans

    ROC = 1/invROC
    
    #Compute the center of the arc
    theta = np.arcsin(diameter/(2*ROC))
    l = ROC*np.cos(theta)
    arc_center = chord_center + chordNormVect*l

    #For convex surface, pos has to be outside the circle.
    if ROC < 0 and np.linalg.norm(pos - arc_center) < np.abs(ROC):
        if verbose:
            print('The line does not hit the arc.')
        return {'isHit': False}

    #First, decompose the vector connecting from the arc_center
    #to pos into the components parallel to the line and orthogonal to it.
    # s is the component in the orthogonal direction and t is the one along
    #the line.
    #A vector orthogonal to the line
    k = np.array([-dirVect[1], dirVect[0]])
    #Solve the equation to decompose the vector pos-arc_center
    M = np.vstack((k, -dirVect)).T
    ans = np.linalg.solve(M, pos - arc_center)
    s = ans[0]
    t = ans[1]
    
    if np.abs(s) > np.abs(ROC):

    #Compute two cross points
    #Work with the chord formed by the line and the circle.
    #d is half the length of the chord.
    d = np.sqrt(ROC**2 - s**2)
    if ROC > 0:
        intersection_point = k*s+arc_center + d*dirVect
        localNormVect = arc_center - intersection_point
    else:
        intersection_point = k*s+arc_center - d*dirVect
        localNormVect = intersection_point - arc_center

    #Check if dirVect and the vector connecting from pos to intersection_point
    #are pointing the same direction.
    if np.dot(dirVect, intersection_point - pos) < 0:
        
    localNormVect = localNormVect/np.linalg.norm(localNormVect)
    localNormAngle = np.mod(np.arctan2(localNormVect[1],
                                       localNormVect[0]), 2*pi)

    #Check if the intersection point is within the
    #diameter
    v0 = - np.sign(ROC) * chordNormVect*(1-1e-16)  #(1-1e-16) is necessary to avoid rounding error
    v1 = intersection_point - arc_center
    v1 = v1/np.linalg.norm(v1)*(1-1e-16)
    if np.arccos(np.dot(v0,v1)) > np.abs(theta):

    distance = np.linalg.norm(intersection_point - pos)

    return {'Intersection Point': intersection_point, 'isHit': True,
            'distance': distance, 'localNormVect': localNormVect,
            'localNormAngle': localNormAngle}

#}}}

#{{{ vector_rotation_2D

def vector_rotation_2D(vect, angle):
    vect = np.array(vect)
    angle = float(angle)
    
    M = np.array([[np.cos(angle), -np.sin(angle)],
                  [np.sin(angle),np.cos(angle)]])
    return np.dot(M, vect)

#}}}

#{{{ normSpheric

def normSpheric(normAngle, invROC, dist_from_center):
    '''
    Returns the local normal angle of a spheric mirror
    at a distance from the center.

    normAngle: The angle formed by the normal vector of the mirror
               at the center and the x-axis.

    invROC: 1/R, where R is the ROC of the mirror.

    dist_from_center: The distance from the center of the point where
                      the local normal is requested.
                      This is a signed value.
                      For a mirror facing +x (the normal vector points
                      towards positive x direction), this distance
                      is positive for points with positive y coordinate,
                      and negative for points with negative y coordinate.
    '''

    normAngle = np.mod(normAngle, 2*pi)
    return np.mod(np.arcsin(- dist_from_center * invROC) + normAngle, 2*pi)

#}}}

#{{{ reflection and deflection angle

def refl_defl_angle(beamAngle, normAngle, n1, n2, invROC=None):
    '''
    Returns a tuples of reflection and deflection angles.
    
    beamAngle: The angle formed by the propagation direction vector
               of the incident beam and the x-axis.

    normAngle: The angle formed by the normal vector of the surface
               and the x-axis.

    n1: Index of refraction of the incident side medium.

    n2: Index of refraction of the transmission side medium.
    '''

    beamAngle = np.mod(beamAngle, 2*pi)
    
    incidentAngle = np.mod(beamAngle - normAngle, 2*pi) - pi

    reflAngle = np.mod(normAngle - incidentAngle, 2*pi)
    
    deflAngle = np.arcsin(n1*np.sin(incidentAngle)/n2)
    deflAngle = np.mod(deflAngle + pi + normAngle, 2*pi)

    if not invROC == None:
        #Calculate ABCD matrices

        #Absolute value of the incident angle
        theta1 = np.abs(incidentAngle)

        #For reflection
        Mrx = np.array([[1., 0.], [-2*n1*invROC/np.cos(theta1), 1.]])
        Mry = np.array([[1., 0.], [-2*n1*invROC*np.cos(theta1), 1.]])

        #For transmission
        theta2 = np.arcsin(n1*np.sin(theta1)/n2)

        nex = (n2*np.cos(theta2)-n1*np.cos(theta1))/(np.cos(theta1)*np.cos(theta2))
        Mtx = np.array([[np.cos(theta2)/np.cos(theta1), 0.],
                        [nex*invROC, np.cos(theta1)/np.cos(theta2)]])

        ney = n2*np.cos(theta2)-n1*np.cos(theta1)
        Mty = np.array([[1., 0.],[ney*invROC, 1.]])

        return (reflAngle, deflAngle, Mrx, Mry, Mtx, Mty)

    else:
        return (reflAngle, deflAngle)
    
#}}}

#}}}

#{{{ VariCAD utility functions

def vc_deflect(theta, theta1, n1, n2):
    '''
    Deflection angle helper function for VariCAD.
    theta is the angle of the surface measured from right.
    theta1 is the angle of the incident beam measured from right.
    It returns an angle of the deflected beam measured from right.
    
    '''
    #Combert theta and theta1 to 0-360 format
    if theta < 0:
        theta = 360.0 + theta

    if theta > 180:
        theta = theta -180.0
        
    if theta1 < 0:
        theta1 = 360.0 + theta1
    
    #Determine the incident angle
    phi = abs(theta - theta1)
    phi1 = 90.0-np.arcsin(np.abs(np.sin(pi*phi/180.0)))*180.0/pi

    #Calculate deflection angle
    phi2 = deflection_angle(phi1, n1, n2)

    #Convert to the 0-360 angle
    s1 = np.sign(np.sin(pi*(theta1 - theta)/180.0))
    s2 = -np.sign(np.cos(pi*(theta1 - theta)/180.0))
    phi2 = theta + s1*90 + s1*s2*phi2
    return phi2 

def vc_reflect(theta, theta1):

    return theta - (theta1 - theta)

#}}}

import gtrace.beam as beam
import gtrace.optcomp as opt
import gtrace.draw.renderer as renderer

#}}}

# are met:
#
#
#

#}}}

#{{{ Draw optical system

def drawOptSys(optList, beamList, filename, fontSize=False):

    d = draw.Canvas()
    d.unit = 'm'

    d.add_layer("main_beam", color=(255,0,0))
    d.add_layer("main_beam_width", color=(255,0,255))    
    d.add_layer("stray_beam", color=(0,255,0))
    d.add_layer("stray_beam_width", color=(0,255,255))    

    for b in beamList:
        if b.stray_order > 0:
            b.layer = 'stray_beam'
            sigma = 1.0
            drawWidth=False
        else:
            b.layer = 'main_beam'
            sigma = 3.0
            drawWidth=True
            
        b.draw(d, sigma=sigma, drawWidth=drawWidth, drawPower=True, drawName=True, fontSize=fontSize)
    
    drawAllOptics(d, optList, drawName=True)

    render.renderDXF(d, filename)
    
#}}}

#{{{ Draw all beams

def drawAllBeams(d, beamList, sigma=3.0, drawWidth=True, drawPower=False,
                 drawROC=False, drawGouy=False, drawOptDist=False, layer=None, mode='x',
                    fontSize=0.01):
    
    for ii in range(len(beamList)):
        if layer is not None:
            beamList[ii].layer = layer
            
        beamList[ii].draw(d, sigma=sigma, mode=mode, drawWidth=drawWidth, drawPower=drawPower,
                    drawROC=drawROC, drawGouy=drawGouy, drawOptDist=drawOptDist,
                    fontSize=fontSize)

#}}}

#{{{ Draw all optics
def drawAllOptics(d, opticsList, drawName=True, layer=None):
    for ii in range(len(opticsList)):
            opticsList[ii].layer = layer

        opticsList[ii].draw(d, drawName=drawName)

#}}}

#{{{ Translate all

def transAll(objList, transVect):
    for ii in range(len(objList)):
        objList[ii].translate(transVect)
#}}}

#{{{ Rotata all

def rotateAll(objList, angle, center):
        objList[ii].rotate(angle, center)
#}}}

'''
gtrace.beam

A module to define GaussianBeam class.

'''

pi = np.pi

from traits.api import HasTraits, Int, Float, CFloat, CComplex, CArray, List, Str

from gtrace.optics.gaussian import q2zr, q2w, q2R, optimalMatching
#import gtrace.sdxf as sdxf
import scipy.optimize as sopt

#}}}

# are met:
#
#
#

#}}}

#{{{ GaussianBeam Class

class GaussianBeam(HasTraits):
    '''
        
    This is a class to represent a Gaussian beam.
    A GaussianBeam object has its origin (pos) and a propagation
    direction (dirVect or dirAngle).
    A GaussianBeam is characterized by q-parameter(s) at its origin.
    The beam can be either circular or elliptic. In order to deal with
    elliptic beams, some parameters are stored in pairs like (q0x, q0y).
    x and y denote the axes of the cross section of the beam. x-axis is
    parallel to the paper and the y-axis is perpendicular to the paper.

    A beam object can be propagated through a free space or made
    to interact with an optics.

    As a beam propagate through optical system, optical distance and
    Gouy phase are accumerated.

    *=== Attributes ===*
    
    **q:** q-parameter of the beam. If the beam is eliptic, q is the q-parameter
        of the best matching circular mode.
        
    **qx:** q-parameter of the beam in the x-direction. [complex float]
    
    **qy:** q-parameter of the beam in the y-direction. [complex float]
    
    **pos:** Position of the beam origin. [(2,) float array]
    
    **dirVect:** Propagation direction vector. [(2,) float array]
    
    **dirAngle:** Propagation direction angle measured from the positive x-axis. [float]
              
    **length:** Length of the beam (used for DXF export) [float]
    
    **layer:** Layer name of the beam when exported to a DXF file.

    **name:** Name of the beam
              
    **wl:** Wavelength in vacuum. Not the wavelength in the medium. [float]

    **n:** Index of refraction of the medium the beam is passing through. [float]

    **P:** Power [float]
    
    **wx:** Beamwidth in x-direction [float]

    **wy:** Beamwidth in y-direction [float]

    **optDist:** Accumulated optical distance.

    **Gouyx:** Accumulated Gouy phase in x-direction.

    **Gouyy:** Accumulated Gouy phase in y-direction.

    **Mx:** ABCD matrix in x-direction.
           This is a 2x2 matrix representing the product
           of ABCD transformations applied to this beam.
           It defaults to an identity matrix.
           Whenever a beam experience an ABCD matrix
           transformation, such as propagation in the space
           or reflection by a curved mirror, the applied ABCD
           matrix is multiplied to this matrix, so that we can
           keep track of what kind of transformations were
           made during beam propagation.

    **My:** ABCD matrix in y-direction. The meaning is the same as Mx.

    **departSurfAngle:** The angle formed by x-axis and the normal vector of
                                  the surface from which the beam is departing.
                                  Default is None. Used by the drawing routine.

    **departSurfInvROC:** Inverse of the ROC of the surface from which the beam is departing.
                                     The ROC is positive for a concave surface seen from
                                     the beam side.
                                     Default is None. Used by the drawing routine.                                    
    
    **incSurfAngle:** The angle formed by the x-arm and the normal vector of
                            the surface to which the beam is incident.
                            Default is None. Used by the drawing routine.    

    **incSurfInvROC:** Inverse of the ROC of the surface to which the beam is incident.
                               The ROC is positive for a concave surface seen from
                               the beam side.
                               Default is None. Used by the drawing routine.

    **stray_order:** An integer indicating if this beam is a stray light or not.
                        The default value is 0. Every time a beam is reflected by an AR surface
                        or transmits an HR surface, this couter is increased by 1.

    '''

#{{{ Traits Definitions

    wl = CFloat(1064.0*nm)  #Wavelength
    P = CFloat(1*W)  #Power
    q = CComplex()  #q-parameter at the origin (best matching circular mode)
    qx = CComplex(1j)  #q-parameter at the origin (x-direction)
    qy = CComplex(1j)  #q-parameter at the origin (y-direction)
    qrx = CComplex()  #Reduced q-parameter at the origin (x-direction)
                                  # qrx = qx/n
    qry = CComplex()  #Reduced q-parameter at the origin (x-direction)
                                  # qry = qy/n
    Gouyx = CFloat(0.0) #Accumurated Gouy phase
    Gouyy = CFloat(0.0) #Accumurated Gouy phase
    wx = CFloat()
    wy = CFloat()
    n = CFloat(1.0)

    pos = CArray(dtype=np.float64, shape=(2,))
    length = CFloat(1.0)
    layer = Str()
    dirVect = CArray(dtype=np.float64, shape=(2,))
    dirAngle = CFloat()
    optDist = CFloat(0.0)

    Mx = CArray(value=[[1,0],[0,1]],dtype=np.float64, shape=(2,2))
    My = CArray(value=[[1,0],[0,1]],dtype=np.float64, shape=(2,2))
    
#}}}

    def __init__(self, q0=1j*2*pi/(1064*nm)*1e-6/2, q0x=False, q0y=False,
                 pos=[0.0,0.0], length=1.0, dirAngle=0.0, dirVect=None,
                 wl=1064*nm, P=1*W, n=1.0, name="Beam", layer='main_beam'):
        '''

        q0: q-parameter. If q0x or q0y is not given, this parameter is used
            as both qx and qy. [complex float]
            
            q0x: q-parameter in x-direction [complex float]
            q0y: q-parameter in y-direction [complex float]
            
            pos: Position of the origin of the beam [(2,) float array]
            length: Length of the beam (used for DXF export) [float]
            dirAngle: Propagation direction angle measured from the positive x-axis.
                      [float]    
                      dirVect: Propagation direction vector. [(2,) float array]
                      wl: wavelength [float]
                      P: Power [float]
                      n: Index of refraction of the medium the beam is passing through. [float]
                      name: Name of the beam
                      layer: Layer name of the beam when exported to a DXF file.
        '''

        self.P = P
        self.pos = pos
        self.length = length
        self.layer = layer
        self.n = n
        
        if q0x:
            self.qx = q0x
        else:
            self.qx = q0
            
        if q0y:
            self.qy = q0y
        else:
            self.qy = q0

        if dirVect != None:
            self.dirVect = dirVect
        else:
            self.dirAngle = dirAngle
            self._dirAngle_changed(0,0)

        self.optDist = 0.0

        self.departSurfAngle = None
        self.departSurfInvROC = None        
        self.incSurfAngle = None
        self.incSurfInvROC = None
        self.stray_order = 0

#}}}

#{{{ copy

        '''
        Make a deep copy.
        '''
        b = copy.deepcopy(self)
        b.qrx = self.qrx
        b.qry = self.qry
        return b

#}}}

#{{{ propagate

    def propagate(self, d):
        '''
        Propagate the beam by a distance d from the current position.
        self.n is used as the index of refraction.
        During this process, the optical distance traveled is added
        to self.optDist.
        self.Goux and self.Gouyy are also updated to record the Gouy
        phase change.
        '''

        qx0 = self.qx
        qy0 = self.qy

        ABCD = np.array([[1.0, d/self.n],[0.0, 1.0]])
        #ABCD = np.array([[1.0, d/self.n**2],[0.0, 1.0]])
        self.ABCDTrans(ABCD)
        self.pos = self.pos + self.dirVect*d

        #Increase the optical distance
        self.optDist = self.optDist + self.n*d

        #Increase the Gouy phase
        self.Gouyx = self.Gouyx + np.arctan(np.real(self.qx)/np.imag(self.qx))\
                     - np.arctan(np.real(qx0)/np.imag(qx0))

        self.Gouyy = self.Gouyy + np.arctan(np.real(self.qy)/np.imag(self.qy))\
                     - np.arctan(np.real(qy0)/np.imag(qy0))

#}}}

#{{{ ABCD Trans
    def ABCDTrans(self, ABCDx, ABCDy=None):
        '''
        Apply ABCD transformation to the beam.

        **ABCDx:** ABCD matrix for x-direction

        **ABCDy:** ABCD matrix for y-direction
        '''

        if ABCDy == None:
            ABCDy = ABCDx
            
        #Update q-parameters
        self.qrx = (ABCDx[0,0]*self.qrx + ABCDx[0,1])/(ABCDx[1,0]*self.qrx + ABCDx[1,1])
        self.qry = (ABCDy[0,0]*self.qry + ABCDy[0,1])/(ABCDy[1,0]*self.qry + ABCDy[1,1])

        #Update Mx and My
        self.Mx = np.dot(ABCDx, self.Mx)
        self.My = np.dot(ABCDy, self.My)        
        
#}}}

        '''
        Rotate the beam around 'center'.
        If center is not given, the beam is rotated
        around self.pos.

        *Arguments*

        **angle:** Rotation angle in radians.

        **center:** Center for rotation. [(2,) array of float]
        '''
        
        if center:
            pointer = self.pos - center
            self.pos = center + pointer
            
        self.dirAngle = self.dirAngle + angle

#}}}

        '''
        Translate the beam by the direction and the distance
        specified by a vector.

       *Arguments*

        **trVect:** A vector to specify the translation direction and
                distance. [(2,) float array]
        '''
        
        self.pos = self.pos + trVect
        
#}}}

#{{{ Flip
    def flip(self, flipDirVect=True):
        '''
        Change the propagation direction of the beam
        by 180 degrees.
        This is equivalent to the reflection of the beam
        by a spherical mirror with the same ROC as the beam.

        If optional argument flipDirVect is set to False,
        the propagation direction of the beam is not changed.
        '''
        self.qx = -np.real(self.qx)+1j*np.imag(self.qx)
        self.qy = -np.real(self.qy)+1j*np.imag(self.qy)
        if flipDirVect:
            self.dirVect = - self.dirVect

#}}}

#{{{ width

    def width(self, dist):
        '''
        Returns the beam width at a distance dist
        from the origin of the beam.
        The width is the readius where the light power becomes 1/e^2.
        '''

        dist = np.array(dist)
        k = 2*pi/(self.wl/self.n)
        qx = self.qx + dist
        qy = self.qy + dist

        return (np.sqrt(-2.0/(k*np.imag(1.0/qx))), np.sqrt(-2.0/(k*np.imag(1.0/qy))))

#}}}

#{{{ R

    def R(self, dist=0.0):
        '''
        Returns the beam ROC at a distance dist
        '''

        k = 2*pi/self.wl
        qx = self.qx + dist/self.n
        qy = self.qy + dist/self.n

        return (q2R(qx), q2R(qy))

#}}}

#{{{ Waist

    def waist(self):
        '''
        Return the tuples of waist size and distance
        '''

        #Waist positions
        dx = -np.real(self.qx)
        dy = -np.real(self.qy)

        #Waist sizes
        wx = self.width(dx)[0]
        wy = self.width(dy)[1]

        return {"Waist Size":(wx, wy), "Waist Position":(dx, dy)}
    
#}}}

#{{{ draw

#{{{ Main Function

    def draw(self, cv, sigma=3., mode='x', drawWidth=True,
             fontSize=False, drawPower=False, drawROC=False, drawGouy=False,
             drawOptDist=False, drawName=False, debug=False):
        '''
        Draw the beam into a DXF object.

        **cv:** draw.Canvas object.
        
        **sigma:** The width of the beam drawn is sigma * (1/e^2 radius of the beam).
               The default is sigma = 3. sigma = 2.7 gives 1ppm diffraction loss.
               [float]
               
        **mode:** 'avg', 'x', or 'y'. A beam can have different widths for x- and y-
              directions. If 'avg' is specified, the average of them are drawn.
              'x' and 'y' specifies to show the width of the respective directions.
        
        **fontSize:** Size of the font used to show supplemental informations. [float]

        **drawWidth:** Whether to draw width or not. [Boolean]

        **drawPower:** Whether to show the beam power. [Boolean]

        **drawROC:** Whether to show the ROC or not. [Boolean]

        **drawGouy:** Whether to show the Gouy phase or not. [Boolean]

        **drawOptDist:** Whether to show the accumulated optical distance or not. [Boolean]

        **drawName:** Whether draw the name of the beam or not. [Boolean]
        '''
        if not fontSize:
            fontSize = self.width(self.length/2)[0]*sigma/5

        start = tuple(self.pos)
        stop = tuple(self.pos + self.dirVect * self.length)

        #Location to put texts
        #mid = self.pos + self.dirVect * self.length/2
        #side = mid+fontSize*1.2*k        
        k = np.array((-self.dirVect[1], self.dirVect[0]))
        text_location = tuple(self.pos + (self.length - 10*fontSize)*self.dirVect + k*fontSize*1)

        #Draw the center line
        cv.add_shape(draw.Line(start, stop), self.layer)

        if drawWidth:
            self.drawWidth(cv, sigma, mode)
            
        annotation = ''
            annotation = annotation+'%s '%self.name
            
        if drawPower:
            annotation = annotation+'P=%.2E '%self.P

        if drawROC:
            annotation = annotation+'ROCx=%.2E '%q2R(self.qx) 
            annotation = annotation+'ROCy=%.2E '%q2R(self.qy)             

        if drawGouy:
            annotation = annotation+'Gouyx=%.2E '%self.Gouyx
            annotation = annotation+'Gouyy=%.2E '%self.Gouyy

        if drawOptDist:
            annotation = annotation+'Optical distance=%.2E '%self.optDist

        cv.add_shape(draw.Text(text=annotation, point=text_location, 
                             height=fontSize), layername='text')

        #Indicate the beam direction
        # text_location = tuple(self.pos + 10*fontSize*self.dirVect +k*fontSize*1)
        # dxf.append(sdxf.Text(text=self.name+':origin P=%.2E'%self.P, point=text_location,
        #                      height=fontSize, layer='text'))

        # text_location = tuple(self.pos + (self.length - 10*fontSize)*self.dirVect + k*fontSize*1)
        # dxf.append(sdxf.Text(text=self.name+':end P=%.2E'%self.P, point=text_location,
        
#}}}

#{{{ drawWidth

    def drawWidth(self, cv, sigma, mode):
        #Draw the width

        #Determine the number of line segments to draw
        zr = q2zr(self.qx)
        resolution = zr/10.0
        if resolution > self.length/10.0:
            resolution = self.length/10.0

        numSegments = self.length/resolution

        if numSegments > 100:
            numSegments = 100

        #q0 to be used
        if mode == 'x':
            q0 = self.qx
        elif mode == 'y':
            q0 = self.qy
        else:
            q0 = (self.qx + self.qy)/2.0

        #Determine the start z coordinates
        if self.departSurfAngle is not None:
            theta = self.departSurfAngle - self.dirAngle
            k = 2*pi/(self.wl/self.n)
            if np.abs(self.departSurfInvROC) > 1.0/100:
                R = 1.0/self.departSurfInvROC
                (z1u, z1d) = optimStartPointR(theta, R, q0, k, sigma)
            else:
                (z1u, z1d) = optimCrossPointFlat(theta, q0, k, sigma)
        else:
            z1u = 0.0
            z1d = 0.0

        #Determine the end z coordinates
        if self.incSurfAngle is not None:
            theta = np.mod(self.incSurfAngle+pi, 2*pi) - self.dirAngle
            k = 2*pi/(self.wl/self.n)            
            if np.abs(self.incSurfInvROC) > 1.0/100:
                R = 1.0/self.incSurfInvROC
                (z2u, z2d) = optimEndPointR(theta, R, q0+self.length, k, sigma)
            else:
                (z2u, z2d) = optimCrossPointFlat(theta, q0+self.length, k, sigma)
        else:
            z2u = 0.0
            z2d = 0.0
            
        du = np.linspace(z1u, self.length+z2u, numSegments)
        dd = np.linspace(z1d, self.length+z2d, numSegments)        
        au = self.width(du)
        ad = self.width(dd)
            wu = au[0]*sigma
            wd = ad[0]*sigma            
            wu = au[1]*sigma
            wd = ad[1]*sigma            
        else:
            wu = sigma*(au[0]+au[1])/2
            wd = sigma*(ad[0]+ad[1])/2            

        v = np.vstack((du,wu))
        v = optics.geometric.vector_rotation_2D(v, self.dirAngle)
        v = v + np.array([self.pos]).T
        cv.add_shape(draw.PolyLine(x=v[0,:], y=v[1,:]), layername=self.layer+"_width")

        v = np.vstack((dd,-wd))
        cv.add_shape(draw.PolyLine(x=v[0,:], y=v[1,:]), layername=self.layer+"_width")        

#}}}        

#{{{ drawWidthOld

    def drawWidthOld(self, dxf, sigma, mode):

        d = np.linspace(0,self.length, numSegments)
        a = self.width(d)
            w = a[0]*sigma
            w = a[1]*sigma
        else:
            w = sigma*(a[0]+a[1])/2

        v = np.vstack((d,w))
        dxf.append(sdxf.LwPolyLine(points=list(v.T), layer=self.layer+"_width"))

        v = np.vstack((d,-w))

#}}}

#}}}

#{{{ Notification Handlers

    def _dirAngle_changed(self, old, new):
                 dirVect=array([np.cos(self.dirAngle), np.sin(self.dirAngle)]))
                 dirAngle = np.mod(self.dirAngle, 2*pi))
#        self.dirVect = array([np.cos(self.dirAngle), np.sin(self.dirAngle)])
#        self.dirAngle = np.mod(self.dirAngle, 2*pi)
            
    def _dirVect_changed(self, old, new):
                 dirVect = self.dirVect/np.linalg.norm(array(self.dirVect)))
                 dirAngle = np.mod(np.arctan2(self.dirVect[1],
                                              self.dirVect[0]), 2*pi))

    def _qx_changed(self, old, new):
        self.wx = q2w(self.qx, wl=self.wl/self.n)
                 qrx = self.qx/self.n)
        self.q = optimalMatching(self.qx, self.qy)[0]

    def _qy_changed(self, old, new):
        self.wy = q2w(self.qy, wl=self.wl/self.n)
                 qry = self.qy/self.n)
        
    def _qrx_changed(self, old, new):
        self.qx = self.qrx*self.n
        
    def _qry_changed(self, old, new):
        self.qy = self.qry*self.n

    def _n_changed(self, old, new):
                 qx = self.qrx*self.n)
                 qy = self.qry*self.n)

#}}}        

#}}}

#{{{ Functions to determine the cross points

def optFunForStartPointR(phi, Mrot, R, q0, k, sigma, side):
    '''
    A function to return the distance between the point on
    the spherical surface at an angle phi and the beam width
    at the same z.
    '''
    cp = np.cos(phi)
    sp = np.sin(phi)
    v = R*np.array([[1-cp], [-sp]])

    a = np.dot(Mrot, v)
    z = a[0]
    w1 = a[1]*side
    w2 = sigma*np.sqrt(-2.0/(k*np.imag(1/(q0+z))))
    return (w1-w2)[0]

def optimStartPointR(theta, R, q0, k, sigma):
    '''
    '''
    ct = np.cos(theta)
    st = np.sin(theta)
    Mrot = np.array([[ct, -st], [st, ct]])
    
    phi1 = sopt.newton(optFunForStartPointR, 0, args=(Mrot, R, q0, k, sigma, 1))
    cp = np.cos(phi1)
    sp = np.sin(phi1)
    z1 = np.dot(Mrot, v)[0][0]

    phi2 = sopt.newton(optFunForStartPointR, 0, args=(Mrot, R, q0, k, sigma, -1))
    cp = np.cos(phi2)
    sp = np.sin(phi2)
    z2 = np.dot(Mrot, v)[0][0]
    
    return (z1,z2) 

def optFunForEndPointR(phi, Mrot, R, q0, k, sigma, side):
    '''
    '''
    v = -R*np.array([[1 - cp], [-sp]])

    z = a[0]

def optimEndPointR(theta, R, q0, k, sigma):
    '''
    '''
    
    phi1 = sopt.newton(optFunForEndPointR, 0, args=(Mrot, R, q0, k, sigma, 1))
    v = R*np.array([[cp - 1], [sp]])

    phi2 = sopt.newton(optFunForEndPointR, 0, args=(Mrot, R, q0, k, sigma, -1))
    
def optFunForFlat(a, Mrot, q0, k, sigma, side):
    '''
    '''
    v = np.array([[0],[a]])

    b = np.dot(Mrot, v)
    z = b[0]
    w1 = b[1]*side

def optimCrossPointFlat(theta, q0, k, sigma):
    '''
    '''
    
    a1 = sopt.newton(optFunForFlat, 0, args=(Mrot, q0, k, sigma, 1))
    v = np.array([[0],[a1]])

    a2 = sopt.newton(optFunForFlat, 0, args=(Mrot, q0, k, sigma, -1))
    v = np.array([[0],[a2]])
    
#}}}