import math
import cmath
from string import *
from pdPythonLib import *
import os

class PCFShape:
    def __init__(self):
        self.shapeType = 1  #1=ellipse 2=regular polygon 3=rectangle
        self.NumHexPoints = 6
        self.IsRounding = 0
        self.polyRoundingDiameter = 1.0
        self.polyroundingResolution = 1   # 1 = "1 extra point"
        self.width = 1
        self.height = 1
        self.rix = 1
        self.poly_xlist = []        # Only for use if the shape is a reg. polygon
        self.poly_ylist = []        # Only for use if the shape is a reg. polygon
        
    def InsertShapeInMWG(self,appClass,wgname,xpos,ypos):
        f = appClass
        if (self.shapeType == 1):
            # shape is a circle
            f.AddCmd(wgname+".insertshape(1,ellipse,"+repr(xpos)+","+repr(ypos)+")")
            f.AddCmd(wgname+".shapes[1].width="+repr(self.width))
            f.AddCmd(wgname+".shapes[1].height="+repr(self.height))
            
        elif (self.shapeType == 2):
            # shape is a regular polygon
            f.AddCmd(wgname+".insertshape(1,polygon,"+repr(xpos)+","+repr(ypos)+")")
            # default polygon has 4 points. delete them.
            f.AddCmd(wgname+".shapes[1].removepoint(1)")
            f.AddCmd(wgname+".shapes[1].removepoint(1)")
            f.AddCmd(wgname+".shapes[1].removepoint(1)")
            f.AddCmd(wgname+".shapes[1].removepoint(1)")
            
            # Now insert shape using points list
            l = len(self.poly_xlist)
            for i in range (0,l):
                x = self.poly_xlist[i]
                y = self.poly_ylist[i]
                f.AddCmd(wgname+".shapes[1].insertpoint(1)")
                f.AddCmd(wgname+".shapes[1].plist[1].x="+repr(x))
                f.AddCmd(wgname+".shapes[1].plist[1].y="+repr(y))
                
        elif (self.shapeType == 3):
            # shape is a rectangle
            f.AddCmd(wgname+".insertshape(1,rectangle,"+repr(xpos)+","+repr(ypos)+")")
            f.AddCmd(wgname+".shapes[1].width="+repr(self.width))
            f.AddCmd(wgname+".shapes[1].height="+repr(self.height))
        f.AddCmd(wgname+".shapes[1].xalign=1")
        f.AddCmd(wgname+".shapes[1].yalign=1")
        f.Exec("")

    def SetShapeProperties(self,appClass,wgname,ShapeInx,rix,alpha):
        f = appClass
        InxStr = repr(ShapeInx)
        f.AddCmd(wgname+".shapes["+InxStr+"].nr11="+repr(rix))
        f.AddCmd(wgname+".shapes["+InxStr+"].nr22="+repr(rix))
        f.AddCmd(wgname+".shapes["+InxStr+"].alpha11="+repr(alpha))
        f.AddCmd(wgname+".shapes["+InxStr+"].alpha22="+repr(alpha))
        f.Exec("")
               
    def CoordinateTransform(self,x,y,alpha):
        # USED BY RoundCorner()
        x_rot = math.cos(alpha)*x - math.sin(alpha)*y
        y_rot = math.sin(alpha)*x + math.cos(alpha)*y
        return x_rot,y_rot  
    def Getpolypoints(self):
        # This returns the co-ords of the vertices if the shape is a regular poly
        #
        xlist = []
        ylist = []
        # First work out the co-ords of the vertices for the regular polygon
        Omega = math.pi/2   #start at the top
        Delta_Omega = 2*math.pi/self.NumHexPoints
        r = self.width/2.0
        for i in range (0,self.NumHexPoints):
            x = r*math.cos(Omega)
            xlist.append(x)
            y = r*math.sin(Omega)
            ylist.append(y)
            Omega = Omega + Delta_Omega

        # if there is rounding then the points list will change......
        if self.IsRounding != 0:
            # We will create a new list of new points
            new_xlist = []
            new_ylist = []
            l = len(xlist)
            # This returns the co-ords of the vertices for the new polygon which has rounded edges.
            for i in range (0,l):
                # First get x and y coords of point
                Pt_x = xlist[i]
                Pt_y = ylist[i]
                
                # This is the theta between two original points
                theta = 2*math.pi/self.NumHexPoints
        
                # r = d/2 (d=is the distance to point of intersection of the circle along line)
                r = self.polyRoundingDiameter/2.0
        
                # We will take the point to be the origin
                CircleCentre_x = 0
                CircleCentre_y = 1.0*r/math.cos(theta/2.0)
                if self.polyroundingResolution != 1:
                    delta_theta = theta/(self.polyroundingResolution)
                    theta_temp = 1.0*theta/2
                    # Going anti-clockwise
                    while theta_temp >= -1.0*theta/2:
                        delta_x = r*math.sin(theta_temp)
                        delta_y = -1.0*CircleCentre_y + r*math.cos(theta_temp)
                        theta_temp = theta_temp - delta_theta
            
                        # Now transform coordinates
                        alpha = 2.0*math.pi*i/self.NumHexPoints
                        delta_x,delta_y = self.CoordinateTransform(delta_x,delta_y,alpha)
                        x = Pt_x + delta_x
                        y = Pt_y + delta_y

                        new_xlist.append(x)
                        new_ylist.append(y)
            
                else:
                    # If we only want 1 extra point we will use a different algorithm so that we conserve hole area
                    # A1 = area of circular segment
                    # Aseg = area of segment
                    # Atri = area of triangle
                    Aseg = 0.5*r*r*theta
                    Atri = r*r*math.cos(theta/2.0)*math.sin(theta/2.0)
                    A1 = Aseg - Atri
                    b1 = 2*r*math.sin(theta/2.0)
                    H = (r/math.cos(theta/2.0)) - (r*math.cos(theta/2.0))
                    AUpperTri = 0.5*b1*H
                    if A1 > AUpperTri:
                        print "A1 error"
                    # Atrap = area of trapezoid (1/2(b1+b2)h)
                    # b2 = (1-h/H)b1
                    # (b1/2H)h^2 - (b1)h + A1 = 0
                    a = b1/(2.0*H)
                    b = -1.0*b1
                    c = A1
                    h = (-1.0*b + math.sqrt((b*b)-(4*a*c)))/(2.0*a)
                    if h > H:
                        h = (-1.0*b - math.sqrt((b*b)-(4*a*c)))/(2.0*a)
                    if h > H:
                        print "H = "+repr(H)
                        print "h error"
            
                    # Now find the co-ordinates of the two points
                    b2 = (1-h/H)*b1
                    Atrap = 0.5*(b1+b2)*h
                    delta_x1 = b2/2.0
                    delta_x2 = -1.0*b2/2.0
                    theta_temp = theta/2
                    delta_y_orig = -1.0*H
                    delta_y = delta_y_orig + h

                    # Now transform coordinates
                    alpha = 2.0*math.pi*i/self.NumHexPoints
                    delta_x1,delta_y1 = self.CoordinateTransform(delta_x1,delta_y,alpha)
                    x1 = Pt_x + delta_x1
                    y1 = Pt_y + delta_y1
                    delta_x2,delta_y2 = self.CoordinateTransform(delta_x2,delta_y,alpha)
                    x2 = Pt_x + delta_x2
                    y2 = Pt_y + delta_y2
                
                    new_xlist.append(x1)
                    new_ylist.append(y1)
                    new_xlist.append(x2)
                    new_ylist.append(y2)

            xlist = new_xlist
            ylist = new_ylist

        return xlist,ylist      
class PCFibre:
    def __init__(self):
        self.nbulk = 2.5
        self.alphaBulk = 0
        self.nhole = 1
        self.alphaHole = 0
        self.latSpace = 0.5
        self.holeRad = 0.175
        self.NumHexPoints = 0                   # If 0 then circular holes
        self.nmodes = 1
        self.wavelength = 1.55
        self.mpr = 60
        self.nRings = 1             
        self.Dist2Edge = 0.0
        self.cladRad = 10
        self.bulk = PCFShape()
        self.latticeHole = PCFShape()
        self.IshollowCore = 0                  # If 1 then the user wants air core
        self.hollowCoreWidth = 5
        self.hollowCoreWallThickness = 0.2
        self.hollowCoreRoundingDiameter = 2.0
        self.hollowCoreHole = PCFShape()
        self.hollowCoreWall = PCFShape()
        self.AbsorbingJacket = PCFShape()
        self.QuarterWaveLayer = PCFShape()
        self.IsRounding = 0                     # If 1 then corners will be rounded
        self.polyroundingResolution = 2
        self.polyRoundingDiameter = 1
        self.jctSize = 5.0
        self.symmetry = 3                       # 1 = Vert, 2 = Hor, 3 = Both
        self.Wavelength = 1.55
        self.FillIn = 1
    def SetMOLABOptions(self,f,wgname,IshollowCore):
        # USED BY GENERATEFIBRE
        f.AddCmd(wgname+".evlist.svp.lambda=" + repr(self.Wavelength))
        f.AddCmd(wgname+".evlist.svp.hcurv=0")
        f.AddCmd(wgname+".evlist.svp.solvid=11")
        f.AddCmd(wgname+".evlist.svp.hsymmetry=0")
        f.AddCmd(wgname+".evlist.svp.vsymmetry=0")
        f.AddCmd(wgname+".evlist.svp.buff=V2 30 30 0 0.5")
        f.AddCmd(wgname+".evlist.mlp.autorun=1")
        f.AddCmd(wgname+".evlist.mlp.speed=0")
        f.AddCmd(wgname+".evlist.mlp.mintefrac=0")
        f.AddCmd(wgname+".evlist.mlp.maxtefrac=100")
        f.AddCmd(wgname+".evlist.mlp.maxnmodes=10")
        if IshollowCore ==1:
            f.AddCmd(wgname+".evlist.mlp.evstart=1.0")
            f.AddCmd(wgname+".evlist.mlp.evend=0.99")
        else:
            f.AddCmd(wgname+".evlist.mlp.evstart=1e+050")
            f.AddCmd(wgname+".evlist.mlp.evend=-1e+050")
        f.AddCmd(wgname+".evlist.mlp.nx=60")
        f.AddCmd(wgname+".evlist.mlp.ny=60")
        f.Exec("")

    def RingGenerator(self,f,wgname,NumringsCounter,CtrAtom_xpos,CtrAtom_ypos,delta_x,delta_y,WGWidth,WGHeight):
        # USED BY GENERATEFIBRE to determine the co-ords of ring in hex lattice
        #
        # First insert atom at right middle
        xpos = CtrAtom_xpos + NumringsCounter*delta_x
        ypos = CtrAtom_ypos + 0
        if xpos >= -1e-14 and xpos <= WGWidth and ypos >= -1e-14 and ypos <= WGHeight:
            self.latticeHole.InsertShapeInMWG(f,wgname,xpos,ypos)
            self.latticeHole.SetShapeProperties(f,wgname,1,self.nhole,self.alphaBulk)
        # Now place atoms up first side
        count = 0
        while count < NumringsCounter:
            count = count + 1
            xpos = xpos - 0.5*delta_x
            ypos = ypos + delta_y
            if xpos >= -1e-14 and xpos <= WGWidth and ypos >= -1e-14 and ypos <= WGHeight:
                self.latticeHole.InsertShapeInMWG(f,wgname,xpos,ypos) 
                self.latticeHole.SetShapeProperties(f,wgname,1,self.nhole,self.alphaBulk)      
        # Now place atoms along top
        count = 0
        while count < NumringsCounter:
            count = count + 1
            xpos = xpos - delta_x
            if xpos >= -1e-14 and xpos <= WGWidth and ypos >= -1e-14 and ypos <= WGHeight:
                self.latticeHole.InsertShapeInMWG(f,wgname,xpos,ypos)
                self.latticeHole.SetShapeProperties(f,wgname,1,self.nhole,self.alphaBulk)  
        # Now place atoms along top-left side
        count = 0
        while count < NumringsCounter:
            count = count + 1
            xpos = xpos - 0.5*delta_x
            ypos = ypos - delta_y
            if xpos >= -1e-14 and xpos <= WGWidth and ypos >= -1e-14 and ypos <= WGHeight:
                self.latticeHole.InsertShapeInMWG(f,wgname,xpos,ypos)
                self.latticeHole.SetShapeProperties(f,wgname,1,self.nhole,self.alphaBulk)  
        # Now place atoms along bottom-left side
        count = 0
        while count < NumringsCounter:
            count = count + 1
            xpos = xpos + 0.5*delta_x
            ypos = ypos - delta_y
            if xpos >= -1e-14 and xpos <= WGWidth and ypos >= -1e-14 and ypos <= WGHeight:
                self.latticeHole.InsertShapeInMWG(f,wgname,xpos,ypos)         
                self.latticeHole.SetShapeProperties(f,wgname,1,self.nhole,self.alphaBulk)    
        # Now place atoms along bottom side
        count = 0
        while count < NumringsCounter:
            count = count + 1
            xpos = xpos + delta_x
            if xpos >= -1e-14 and xpos <= WGWidth and ypos >= -1e-14 and ypos <= WGHeight:
                self.latticeHole.InsertShapeInMWG(f,wgname,xpos,ypos)
                self.latticeHole.SetShapeProperties(f,wgname,1,self.nhole,self.alphaBulk)  
        # Now place atoms along bottom-right side
        count = 1
        while count < NumringsCounter:
            count = count + 1
            xpos = xpos + 0.5*delta_x
            ypos = ypos + delta_y
            if xpos >= -1e-14 and xpos <= WGWidth and ypos >= -1e-14 and ypos <= WGHeight:
                self.latticeHole.InsertShapeInMWG(f,wgname,xpos,ypos)
                self.latticeHole.SetShapeProperties(f,wgname,1,self.nhole,self.alphaBulk)  
        
    def SetSymmetryPlane(self,f,wgname,symmetry,WGHeight,WGWidth):     
        # USED BY GENERATEFIBRE to set the symmetryplane
        #
        # New height and width:
        if symmetry == 1:    #Vert
            WGHeight = WGHeight/2.0
            WGWidth = WGWidth
        elif symmetry == 2:  #Hor
            WGHeight = WGHeight
            WGWidth = WGWidth/2.0
        elif symmetry == 3:  #Both
            WGHeight = WGHeight/2.0
            WGWidth = WGWidth/2.0
        else:
            WGHeight = WGHeight
            WGWidth = WGWidth

        # Re-size WG Window
        f.AddCmd(wgname+".width="+repr(WGWidth))
        f.AddCmd(wgname+".height="+repr(WGHeight))

        # Re-size and re-position bulk index if it is a rectangle
        n = f.Exec(wgname+".nshapes")
        n = int(n)
        if self.cladRad ==0:
            f.AddCmd(wgname+".shapes["+repr(n)+"].width="+repr(WGWidth))
            f.AddCmd(wgname+".shapes["+repr(n)+"].height="+repr(WGHeight))
            f.AddCmd(wgname+".shapes["+repr(n)+"].xposn="+repr(WGWidth/2.0))
            f.AddCmd(wgname+".shapes["+repr(n)+"].yposn="+repr(WGHeight/2.0))
            f.Exec("")
            
        # Now delete all shapes that are outside new boundary
        # We have to do this backwards as the number of shapes will change
        for i in range (n,0,-1):
            f.AddCmd(wgname+".shapes["+repr(i)+"].xposn")
            x,y = f.Exec(wgname+".shapes["+repr(i)+"].yposn")
            if x > (WGWidth+0.1) or y > (WGHeight+0.1):
                f.AddCmd(wgname+".deleteshape("+repr(i)+")") 
            
    def AddNotesSubnode(self,f,wgname):
        # This Notes subnode will store all specs of the waveguide.
        f.AddCmd(wgname+".addsubnode(notesTB,\"Specs\")")
        tempstring = "Bulk Index:   " + repr(self.nbulk)
        f.AddCmd(wgname+".subnodes[1].appendtext("+tempstring+")")
        f.AddCmd(wgname+".subnodes[1].appendtext(\"\\n\")")
        tempstring = "Hole Index:   " + repr(self.nhole)
        f.AddCmd(wgname+".subnodes[1].appendtext("+tempstring+")")
        f.AddCmd(wgname+".subnodes[1].appendtext(\"\\n\")")
        if self.NumHexPoints == 6:
            self.holeRad = self.holeRad / (2.0/math.sqrt(3.0))
        tempstring = "Hole radius:   " + repr(self.holeRad)
        f.AddCmd(wgname+".subnodes[1].appendtext("+tempstring+")")
        f.AddCmd(wgname+".subnodes[1].appendtext(\"\\n\")")
        tempstring = "Hole polygon resolution:   " + repr(self.NumHexPoints)
        f.AddCmd(wgname+".subnodes[1].appendtext("+tempstring+")")
        f.AddCmd(wgname+".subnodes[1].appendtext(\"\\n\")")
        tempstring = "Lattice Period:   " + repr(self.latSpace)
        f.AddCmd(wgname+".subnodes[1].appendtext("+tempstring+")")
        f.AddCmd(wgname+".subnodes[1].appendtext(\"\\n\")")
        if self.FillIn != 0:
            self.nRings = (self.nRings + 1)/2
        tempstring = "Number of rings:   " + repr(self.nRings)
        f.AddCmd(wgname+".subnodes[1].appendtext("+tempstring+")")
        f.AddCmd(wgname+".subnodes[1].appendtext(\"\\n\")")
        tempstring = "Cladding Radius:   " + repr(self.cladRad)
        f.AddCmd(wgname+".subnodes[1].appendtext("+tempstring+")")
        f.AddCmd(wgname+".subnodes[1].appendtext(\"\\n\")")
        tempstring = "Cladding Distance:   " + repr(self.Dist2Edge)
        f.AddCmd(wgname+".subnodes[1].appendtext("+tempstring+")")
        f.AddCmd(wgname+".subnodes[1].appendtext(\"\\n\")")
        if self.IshollowCore == 0:
            temp_string = "N"
        else:
            temp_string = "Y"
        tempstring = "Hollow Core:   " + temp_string
        f.AddCmd(wgname+".subnodes[1].appendtext("+tempstring+")")
        f.AddCmd(wgname+".subnodes[1].appendtext(\"\\n\")")
        tempstring = "Hollow Core Width:   " + repr(self.hollowCoreWidth)
        f.AddCmd(wgname+".subnodes[1].appendtext("+tempstring+")")
        f.AddCmd(wgname+".subnodes[1].appendtext(\"\\n\")")
        tempstring = "Hollow Core Wall Thickness:   " + repr(self.hollowCoreWallThickness)
        f.AddCmd(wgname+".subnodes[1].appendtext("+tempstring+")")
        f.AddCmd(wgname+".subnodes[1].appendtext(\"\\n\")")
        tempstring = "Hollow Core Rounding Diameter:   " + repr(self.hollowCoreRoundingDiameter)
        f.AddCmd(wgname+".subnodes[1].appendtext("+tempstring+")")
        f.AddCmd(wgname+".subnodes[1].appendtext(\"\\n\")")
        if self.IsRounding == 0:
            temp_string = "N"
        else:
            temp_string = "Y"
        tempstring = "Corner Rounding:   " + temp_string
        f.AddCmd(wgname+".subnodes[1].appendtext("+tempstring+")")
        f.AddCmd(wgname+".subnodes[1].appendtext(\"\\n\")")
        tempstring = "Rounding Resolution:   " + repr(self.polyroundingResolution)
        f.AddCmd(wgname+".subnodes[1].appendtext("+tempstring+")")
        f.AddCmd(wgname+".subnodes[1].appendtext(\"\\n\")")
        tempstring = "Rounding Diameter:   " + repr(self.polyRoundingDiameter)
        f.AddCmd(wgname+".subnodes[1].appendtext("+tempstring+")")
        f.AddCmd(wgname+".subnodes[1].appendtext(\"\\n\")")
        tempstring = "Absorbing Jacket Size:   " + repr(self.jctSize)
        f.AddCmd(wgname+".subnodes[1].appendtext("+tempstring+")")
        f.AddCmd(wgname+".subnodes[1].appendtext(\"\\n\")")
        tempstring = "Wavelength:   " + repr(self.Wavelength)
        f.AddCmd(wgname+".subnodes[1].appendtext("+tempstring+")")
        f.AddCmd(wgname+".subnodes[1].appendtext(\"\\n\")")
        f.Exec("")
    def GetAbsorptionParams(self,f,wgname):
        # Used by CreateAbsorbingLayer
        # returns the refractive indices of the absorbing layer and size of qw layer
        if self.IshollowCore ==1:
            ApproxModeIndex = 0.99*self.nhole
        else:
            ApproxModeIndex = 0.99*self.nbulk
            
        # This is the strength of the absorption - 0.01 = (1% reflection)
        reftol = 0.01     
                                                                     
        # Default: alpha is set so that power is reduced by 99% over 2*jctSize
        alpha = -1*math.log(math.sqrt(reftol))/self.jctSize
        jct_N_imag = (alpha*self.Wavelength)/(4*math.pi)

        # Need to determine refractive index of quarterwave plate
        n1 = self.nbulk + 0j
        n2 = complex(self.nbulk,jct_N_imag)
        n_m = cmath.sqrt(n1*n2)
        alpha_m = 4*math.pi*n_m.imag/self.Wavelength
                    
        # Need to determine the width needed for a quarter wave plate.
        # This depends on beta (i.e. angle of incidence)
        n_m_abs = abs(n_m)
        theta = math.acos(ApproxModeIndex/n_m_abs)
        d = (self.Wavelength/(4*n_m_abs))*math.sin(theta)

        return alpha,n_m,alpha_m,d
    
    def RealignShapes(self,f,wgname,start,end,delta_x,delta_y):
        # Used by CreateAbsorbingLayer
        for i in range (start,end):
            f.AddCmd(wgname + ".shapes[" + repr(i) + "].xposn")
            [xposn,yposn] = f.Exec(wgname + ".shapes[" + repr(i) + "].yposn")
            xposn = xposn + delta_x
            yposn = yposn + delta_y
            f.AddCmd(wgname + ".shapes[" + repr(i) + "].xposn=" + repr(xposn))
            f.AddCmd(wgname + ".shapes[" + repr(i) + "].yposn=" + repr(yposn))
        
    def CreateAbsorbingLayer(self,f,wgname):
        # Used by GenerateFibre()
        # This generates an absorbing jacket with a quarter wave layer to minimise reflections

        # First get parameters of absorbing layer
        alpha,n_m,alpha_m,d = self.GetAbsorptionParams(f,wgname)
        n_m_real = n_m.real
                
        # Resize Computational Window
        f.AddCmd(wgname+".width")
        [old_WGwidth,old_WGheight] = f.Exec(wgname+".height")
     
        # Note that self.Dist2Edge has a new use - when a jct is defined it describes the air gap
        if self.cladRad == 0:
            WGWidth = (old_WGwidth + 2.0*(d + self.jctSize + self.Dist2Edge))
            WGHeight = (old_WGheight + 2.0*(d + self.jctSize + self.Dist2Edge))
        else:
            WGWidth = 2.0*(self.cladRad + self.jctSize + d + self.Dist2Edge)
            WGHeight = 2.0*(self.cladRad + self.jctSize + d + self.Dist2Edge)

        # We want a square
        if WGWidth > WGHeight:
            WGHeight = WGWidth
        else:
            WGWidth = WGHeight
        f.AddCmd(wgname+".height="+repr(WGHeight))
        f.AddCmd(wgname+".width="+repr(WGWidth))

        # Insert quarterwave layer
        if self.cladRad !=0:
            self.QuarterWaveLayer.width = WGWidth - 2.0*self.Dist2Edge - 2.0*self.jctSize
            self.QuarterWaveLayer.height = WGHeight - 2.0*self.Dist2Edge - 2.0*self.jctSize
            self.QuarterWaveLayer.InsertShapeInMWG(f,wgname,WGWidth/2,WGHeight/2)
            self.QuarterWaveLayer.SetShapeProperties(f,wgname,1,n_m_real,alpha_m)
        else:
            self.QuarterWaveLayer.shapeType = 3     # Rectangle
            # Insert the two horizontal ones
            self.QuarterWaveLayer.width = d
            self.QuarterWaveLayer.height = WGHeight
            xposn = WGWidth - self.jctSize - self.QuarterWaveLayer.width/2.0
            yposn = WGHeight/2
            self.QuarterWaveLayer.InsertShapeInMWG(f,wgname,xposn,yposn)
            self.QuarterWaveLayer.SetShapeProperties(f,wgname,1,n_m_real,alpha_m)
            xposn = self.jctSize + self.QuarterWaveLayer.width/2.0
            self.QuarterWaveLayer.InsertShapeInMWG(f,wgname,xposn,yposn)
            self.QuarterWaveLayer.SetShapeProperties(f,wgname,1,n_m_real,alpha_m)
            # Insert the two vertical ones
            self.QuarterWaveLayer.width = WGWidth
            self.QuarterWaveLayer.height = d
            xposn = WGWidth/2
            yposn = WGHeight - self.jctSize - self.QuarterWaveLayer.height/2.0
            self.QuarterWaveLayer.InsertShapeInMWG(f,wgname,xposn,yposn)
            self.QuarterWaveLayer.SetShapeProperties(f,wgname,1,n_m_real,alpha_m)
            yposn = self.jctSize + self.QuarterWaveLayer.height/2.0
            self.QuarterWaveLayer.InsertShapeInMWG(f,wgname,xposn,yposn)
            self.QuarterWaveLayer.SetShapeProperties(f,wgname,1,n_m_real,alpha_m)
        
        # Insert absorbing layer
        if self.FillIn !=1 and self.cladRad !=0:
            self.AbsorbingJacket.shapeType = 1
            self.AbsorbingJacket.width = WGWidth-2.0*self.Dist2Edge
            self.AbsorbingJacket.height = WGHeight-2.0*self.Dist2Edge
            self.AbsorbingJacket.InsertShapeInMWG(f,wgname,WGWidth/2,WGHeight/2)
            self.AbsorbingJacket.SetShapeProperties(f,wgname,1,self.nbulk,alpha)
        else:
            # The absorbing layer will be 4 rectangles around the edges of the structure
            self.AbsorbingJacket.shapeType = 3      # Rectangle
            self.AbsorbingJacket.width = self.jctSize
            self.AbsorbingJacket.height = WGHeight
            # Insert the two horizontal ones
            xposn = WGWidth - self.AbsorbingJacket.width/2.0
            yposn = WGHeight/2
            self.AbsorbingJacket.InsertShapeInMWG(f,wgname,xposn,yposn)
            self.AbsorbingJacket.SetShapeProperties(f,wgname,1,n_m_real,alpha_m)
            xposn = self.AbsorbingJacket.width/2.0
            self.AbsorbingJacket.InsertShapeInMWG(f,wgname,xposn,yposn)
            self.AbsorbingJacket.SetShapeProperties(f,wgname,1,n_m_real,alpha_m)
            # Insert the two vertical ones
            self.AbsorbingJacket.width = WGWidth
            self.AbsorbingJacket.height = self.jctSize
            xposn = WGWidth/2
            yposn = WGHeight - self.AbsorbingJacket.height/2.0
            self.AbsorbingJacket.InsertShapeInMWG(f,wgname,xposn,yposn)
            self.AbsorbingJacket.SetShapeProperties(f,wgname,1,n_m_real,alpha_m)
            yposn = self.AbsorbingJacket.height/2.0
            self.AbsorbingJacket.InsertShapeInMWG(f,wgname,xposn,yposn)
            self.AbsorbingJacket.SetShapeProperties(f,wgname,1,n_m_real,alpha_m)

        # Re-Order shapes
        n = f.Exec(wgname + ".nshapes")
        n = int(n)
        if self.FillIn != 1 and self.cladRad !=0:
            f.AddCmd(wgname + ".ordershape(2," + repr(n) + ")")
            f.AddCmd(wgname + ".ordershape(1," + repr(n) + ")")
        
        # Now we need to re-align holes to jacket
        delta_x = (WGWidth - old_WGwidth)/2.0
        delta_y = (WGHeight - old_WGheight)/2.0
        if self.FillIn != 1 and self.cladRad !=0:
            start = 1
            end = n-1
        else:
            # Ignore the absorbing layers
            start = 9
            end = n        
        self.RealignShapes(f,wgname,start,end,delta_x,delta_y)
            
        return WGWidth,WGHeight       
        
    def CreateHollowCore(self,f,wgname,WGWidth,WGHeight):
        # Used by GenerateFibre() to insert the hollow core
        if self.IshollowCore == 1:
            # hollow Core is same type of shape as holes
            self.hollowCoreHole.shapeType = self.latticeHole.shapeType
            self.hollowCoreWall.shapeType = self.latticeHole.shapeType
            self.hollowCoreHole.NumHexPoints = self.NumHexPoints
            self.hollowCoreWall.NumHexPoints = self.NumHexPoints
            self.hollowCoreHole.IsRounding = self.IsRounding
            self.hollowCoreWall.IsRounding = self.IsRounding 
            self.hollowCoreHole.polyRoundingDiameter = self.hollowCoreRoundingDiameter
            self.hollowCoreWall.polyRoundingDiameter = self.hollowCoreRoundingDiameter
            self.hollowCoreHole.polyroundingResolution = self.polyroundingResolution
            self.hollowCoreWall.polyroundingResolution = self.polyroundingResolution
            
            # hollow Core will be placed in center
            xpos = WGWidth/2.0
            ypos = WGHeight/2.0

            # Insert the hollowCoreWall and set up properties.
            self.hollowCoreWall.width = self.hollowCoreWidth + 2.0*self.hollowCoreWallThickness
            self.hollowCoreWall.height = self.hollowCoreWidth + 2.0*self.hollowCoreWallThickness
            # If holes are polys then work out the points list for the vertices
            if self.hollowCoreWall.NumHexPoints != 0:
                self.hollowCoreWall.poly_xlist,self.hollowCoreWall.poly_ylist = self.hollowCoreWall.Getpolypoints()
                
            self.hollowCoreWall.InsertShapeInMWG(f,wgname,xpos,ypos)
            self.hollowCoreWall.SetShapeProperties(f,wgname,1,self.nbulk,self.alphaBulk)                


            # Insert the hollowCoreHole and set up properties.
            self.hollowCoreHole.width = self.hollowCoreWidth
            self.hollowCoreHole.height = self.hollowCoreWidth
            # If holes are polys then work out the points list for the vertices
            if self.hollowCoreHole.NumHexPoints != 0:
                self.hollowCoreHole.poly_xlist,self.hollowCoreHole.poly_ylist = self.hollowCoreHole.Getpolypoints()
                
            self.hollowCoreHole.InsertShapeInMWG(f,wgname,xpos,ypos)
            self.hollowCoreHole.SetShapeProperties(f,wgname,1,self.nhole,self.alphaHole)


            # By Default the Hole Core (and wall) should be rotated by 30 degrees
            # So that it is flat along the x direction
            # Only if the hollow core is hexagonal
            if self.NumHexPoints == 6:
                npoints = f.Exec(wgname+".shapes[1].npoints")
                npoints = int(npoints)
                alpha = math.pi/6.0
                for i in range (0,npoints):
                    j = i + 1
                    f.AddCmd(wgname+".shapes[1].plist["+repr(j)+"].x")
                    f.AddCmd(wgname+".shapes[2].plist["+repr(j)+"].x")
                    f.AddCmd(wgname+".shapes[1].plist["+repr(j)+"].y")
                    [x1,x2,y1,y2] = f.Exec(wgname+".shapes[2].plist["+repr(j)+"].y")
                    x1,y1 = self.hollowCoreWall.CoordinateTransform(x1,y1,alpha)
                    x2,y2 = self.hollowCoreHole.CoordinateTransform(x2,y2,alpha)
                    f.AddCmd(wgname+".shapes[1].plist["+repr(j)+"].x="+repr(x1))
                    f.AddCmd(wgname+".shapes[1].plist["+repr(j)+"].y="+repr(y1))
                    f.AddCmd(wgname+".shapes[2].plist["+repr(j)+"].x="+repr(x2))
                    f.AddCmd(wgname+".shapes[2].plist["+repr(j)+"].y="+repr(y2))
                
            # Delete holes underneath the hole core
            # Have only done this for the first ring
            d = self.hollowCoreWall.width
            nshapes = npoints = f.Exec(wgname+".nshapes")
            if d > 2.0*self.latSpace:
                # Delete holes in inner ring
                for i in range (0,6):
                    ShapeNo = nshapes - (i+1)
                    f.AddCmd(wgname+".deleteshape("+repr(ShapeNo)+")")
                f.Exec("")
                           
    def InsertLatticeHoles(self,f,wgname,WGWidth,WGHeight):
        # Used by GenerateFibre() to insert the lattice
        # Properties of all lattice atoms:
        if self.NumHexPoints == 0:
            self.latticeHole.shapeType = 1    # Circle
        else:
            self.latticeHole.shapeType = 2    # Regular Polygon
            self.latticeHole.NumHexPoints = self.NumHexPoints

        self.latticeHole.IsRounding = self.IsRounding      

        if self.latticeHole.IsRounding != 0:
            self.latticeHole.polyroundingResolution = self.polyroundingResolution
                    
        self.latticeHole.width = 2.0*self.holeRad
        self.latticeHole.height = 2.0*self.holeRad
        
        # Co-ordinates of the Centre Atom
        CtrAtom_xpos = WGWidth/2.0
        CtrAtom_ypos = WGHeight/2.0
        
        # BY DEFAULT THE CENTRE ATOM WILL NOT BE INSERTED
        # Now insert atoms for each ring.
        delta_x = self.latSpace
        delta_y = math.sqrt(3)/2*self.latSpace
        
        if self.FillIn != 0:
            # If we are filling in the whole of the region we simply insert more rings and then
            # delete the ones outside
            self.nRings = 2*self.nRings - 1

        # If holes are polys then work out the points list for the vertices
        if self.NumHexPoints != 0:
            self.latticeHole.poly_xlist,self.latticeHole.poly_ylist = self.latticeHole.Getpolypoints()
        # This will insert the holes one ring at a time
        NumringsCounter = 0
        while NumringsCounter < self.nRings:
            NumringsCounter = NumringsCounter + 1                       
            self.RingGenerator(f,wgname,NumringsCounter,CtrAtom_xpos,CtrAtom_ypos,delta_x,delta_y,WGWidth,WGHeight)
            
        f.Exec("")
            
    def InsertBulk(self,f,wgname):
        # Used by GenerateFibre() to insert the bulk material
        if self.cladRad == 0:
            # if no finte cladding then shape is a rectangle
            self.bulk.shapeType=3
            
            # Set width and height of the mwg
            WGWidth = 2.0*(self.nRings*self.latSpace + self.holeRad + self.Dist2Edge)
            if self.FillIn != 1:
                WGHeight = 2.0*(self.nRings*self.latSpace*math.cos(math.pi/6) + self.holeRad + self.Dist2Edge)
            else:
                # if we are "filling in" then we want to include half a row at top and bottom
                WGHeight = 2.0*((self.nRings + 0.5)*self.latSpace*math.cos(math.pi/6) + self.holeRad + self.Dist2Edge)
            f.AddCmd(wgname+".width="+repr(WGWidth)+"")
            f.AddCmd(wgname+".height="+repr(WGHeight)+"")
            
            # Bulk will be placed in center
            bulkxpos = WGWidth/2.0
            bulkypos = WGHeight/2.0 
            self.bulk.width=WGWidth
            self.bulk.height=WGHeight
            self.bulk.InsertShapeInMWG(f,wgname,bulkxpos,bulkypos)
            self.bulk.SetShapeProperties(f,wgname,1,self.nbulk,self.alphaBulk)
            
        else:
            # If a finite cladding is desired then bulk shape is a circle
            self.bulk.shapeType=1
            
            # Set width and height of the mwg
            WGWidth = 2.0*(self.cladRad)+2.0
            WGHeight = 2.0*(self.cladRad)+2.0
            f.AddCmd(wgname+".width="+repr(WGWidth)+"")
            f.AddCmd(wgname+".height="+repr(WGHeight)+"")
            
            # Bulk will be placed in center
            bulkxpos = WGWidth/2.0
            bulkypos = WGHeight/2.0

            # Insert the bulk material (cladding) and set up properties.
            self.bulk.width=2.0*self.cladRad
            self.bulk.height=2.0*self.cladRad
            self.bulk.InsertShapeInMWG(f,wgname,bulkxpos,bulkypos)
            self.bulk.SetShapeProperties(f,wgname,1,self.nbulk,self.alphaBulk)
            
        f.Exec("")
        return WGWidth,WGHeight
    
    def GenerateFibre(self,projectName,projectDir,pdAppClassName,wgname):
        # THIS IS THE MAIN FUNCTION OF PCFibre
        f = pdAppClassName
        f.AddCmd("app.setwdir("+projectDir+")")
        f.AddCmd("app.addsubnode(mwguideNode,"+wgname+")")
        f.AddCmd("Ref& "+wgname+"=app.findnode("+wgname+")")

        # To coincide with Saitoh's paper if the lattice holes are hexagonal then d is
        # the distance between the side walls, not point to point
        if self.NumHexPoints == 6:
            self.holeRad = self.holeRad * (2.0/math.sqrt(3.0))

        # To avoid conflict: if FillIn = 1, then we will override CladRad and Dist2Edge
        # Dist2Edge needs to be half a period
        if self.FillIn == 1:
            self.cladRad = 0.0
            self.Dist2Edge = 0.5*self.latSpace - self.holeRad  
        
        # Now Set size of MWG and bulk material          
        [WGWidth,WGHeight] = self.InsertBulk(f,wgname)    
        
        # Now Insert holes
        self.InsertLatticeHoles(f,wgname,WGWidth,WGHeight)
                
        # Now insert hollow core
        self.CreateHollowCore(f,wgname,WGWidth,WGHeight)
      
        # Now set default MOLAB Options
        self.SetMOLABOptions(f,wgname,self.IshollowCore)
        
        # Now insert absorbing Jacket
        if self.jctSize != 0:
            [WGHeight,WGWidth] = self.CreateAbsorbingLayer(f,wgname)

        # Now set symmetry plane if applicable
        self.SetSymmetryPlane(f,wgname,self.symmetry,WGHeight,WGWidth)            
            
        # Now create a Notes subnode and include all specs of the waveguide
        self.AddNotesSubnode(f,wgname)

        # prjectDir is currently "" and the file is saved in the working dir (with settings file)    
        f.Exec(wgname+".savetofile(" + projectDir + "\\" +projectName + ".mwg)") 
        
        

        
