# -*- coding: utf-8 -*-
from __future__ import division, print_function, absolute_import
import numpy as np
import itertools
from numpy import sqrt, pi, cos, sin, log, exp, sinh
# from scipy.interpolate import interp1d

import gdspy
from gdspy import clipper
from phidl.device_layout import Device, Port, Polygon, CellArray
from phidl.device_layout import _parse_layer, DeviceReference
import copy as python_copy
from collections import OrderedDict
import pickle
import json
import warnings
from functools import update_wrapper
from phidl.constants import _glyph,_width,_indent


##### Categories:
# Polygons / shapes
# Boolean functions
# Lithography test structures
# Utility functions (copying, importing, extracting)
# Pads
# Taper
# Text
# Wafer / Die
# Waveguide
# Packer tool / Fill tool
# Photonics
# Optimal (current-crowding) curves
# Superconducting devices


#==============================================================================
#
# Polygons / Shapes
#
#==============================================================================


def rectangle(size = (4,2), layer = 0):
    """Generate rectangle geometry.

    Parameters
    ----------
    size : tuple
        Width and height of rectangle.
    layer : int, array-like[2], or set
        Specific layer(s) to put polygon geometry on.

    Returns
    -------
    A Device with a single rectangle in it
    """

    D = Device(name = 'rectangle')
    points = [[size[0], size[1]], [size[0], 0], [0, 0], [0, size[1]]]
    D.add_polygon(points, layer = layer)
    return D


def bbox(bbox = [(-1,-1),(3,4)], layer = 0):
    """ Creates a bounding box rectangle from coordinates, to allow
    creation of a rectangle bounding box directly form another shape.

    Parameters
    ----------
    bbox : list of tuples
        Coordinates of the box [(x1,y1),(x2,y2)].
    layer : int, array-like[2], or set
        Specific layer(s) to put polygon geometry on.

    Returns
    -------
    A Device with a single rectangle in it

    Examples
    --------
    >>> D = pg.bbox(anothershape.bbox)
    """

    D = Device(name = 'bbox')
    (a,b),(c,d)  = bbox
    points = ((a,b), (c,b), (c,d), (a,d))
    D.add_polygon(points, layer = layer)
    return D


def cross(length = 10, width = 3, layer = 0):
    """Generates a right-angle cross (+ shape, symmetric) from two
    rectangles of specified length and width.

    Parameters
    ----------
    length : float
        Length of the cross from one end to the other.
    width : float
        Width of the arms of the cross.
    layer : int, array-like[2], or set
        Specific layer(s) to put polygon geometry on.

    Returns
    -------
    A Device with a cross in it
    """

    D = Device(name = 'cross')
    R = rectangle(size = (width, length), layer = layer)
    r1 = D.add_ref(R).rotate(90)
    r2 = D.add_ref(R)
    r1.center = (0,0)
    r2.center = (0,0)
    return D


def ellipse(radii = (10,5), angle_resolution = 2.5, layer = 0):
    """Generate an ellipse geometry.

    Parameters
    ----------
    radii : tuple
        Semimajor (x) and semiminor (y) axis lengths of the ellipse.
    angle_resolution : float
        Resolution of the curve of the ring (# of degrees per point).
    layer : int, array-like[2], or set
        Specific layer(s) to put polygon geometry on.

    Returns
    -------
    A Device with an ellipse polygon in it
    """

    D = Device(name = 'ellipse')
    a = radii[0]
    b = radii[1]
    t = np.linspace(0, 360, int(np.ceil(360/angle_resolution) + 1))*pi/180
    r = a*b/(sqrt((b*cos(t))**2 + (a*sin(t))**2))
    xpts = r*cos(t)
    ypts = r*sin(t)
    D.add_polygon(points = (xpts,ypts), layer = layer)
    return D


def circle(radius = 10, angle_resolution = 2.5, layer = 0):
    """Generate a circle geometry.

    Parameters
    ----------
    radius : float
        Radius of the circle.
    angle_resolution : float
        Resolution of the curve of the ring (# of degrees per point).
    layer : int, array-like[2], or set
        Specific layer(s) to put polygon geometry on.

    Returns
    -------
    A Device with an circle polygon in it

    """

    D = Device(name = 'circle')
    t = np.linspace(0, 360, int(np.ceil(360/angle_resolution) + 1))*pi/180
    xpts = (radius*cos(t)).tolist()
    ypts = (radius*sin(t)).tolist()
    D.add_polygon(points = (xpts,ypts), layer = layer)
    return D


def ring(radius = 10, width = 0.5, angle_resolution = 2.5, layer = 0):
    """Generate a ring geometry.

    Parameters
    ----------
    radius : float
        Middle radius of the ring.
    width : float
        Width of the ring.
    angle_resolution : float
        Resolution of the curve of the ring (# of degrees per point).
    layer : int, array-like[2], or set
        Specific layer(s) to put polygon geometry on.

    Returns
    -------
    A Device with an ring polygon in it

    Notes
    -----
    The ring is formed by taking the radius out to the specified value, and then constructing the thickness by dividing the width in half and adding that value to either side of the radius.

    The angle_resolution alters the precision of the curve of the ring. Larger values yield lower resolution.
    """

    D = Device(name = 'ring')
    inner_radius = radius - width/2
    outer_radius = radius + width/2
    n = int(np.round(360/angle_resolution))
    t = np.linspace(0, 360, n+1)*pi/180
    inner_points_x = (inner_radius*cos(t)).tolist()
    inner_points_y = (inner_radius*sin(t)).tolist()
    outer_points_x = (outer_radius*cos(t)).tolist()
    outer_points_y = (outer_radius*sin(t)).tolist()
    xpts = inner_points_x + outer_points_x[::-1]
    ypts = inner_points_y + outer_points_y[::-1]
    D.add_polygon(points = (xpts,ypts), layer = layer)
    return D


def arc(radius = 10, width = 0.5, theta = 45, start_angle = 0, angle_resolution = 2.5, layer = 0):
    """ Creates an arc of arclength ``theta`` starting at angle ``start_angle``

    Parameters
    ----------
    radius : float
        Radius of the arc centerline.
    width : float
        Width of the arc.
    theta : float
        Total angle coverage of the arc.
    start_angle : float
        Starting angle.
    angle_resolution : float
        Resolution of the curve of the arc.
    layer : int, array-like[2], or set
        Specific layer(s) to put polygon geometry on.

    Returns
    -------
    A Device with an arc polygon in it, and two ports (`1` and `2`) on either end


    Notes
    -----
    Theta = 0 is located along the positive x-axis relative to the centre of the arc.
    Ports are added to each end of the arc to facilitate connecting those ends to other geometries.
    """

    inner_radius = radius-width/2
    outer_radius = radius+width/2
    angle1 = (start_angle)*pi/180
    angle2 = (start_angle + theta)*pi/180
    t = np.linspace(angle1, angle2, int(np.ceil(abs(theta)/angle_resolution)))
    inner_points_x = (inner_radius*cos(t)).tolist()
    inner_points_y = (inner_radius*sin(t)).tolist()
    outer_points_x = (outer_radius*cos(t)).tolist()
    outer_points_y = (outer_radius*sin(t)).tolist()
    xpts = inner_points_x + outer_points_x[::-1]
    ypts = inner_points_y + outer_points_y[::-1]

    D = Device('arc')
    D.add_polygon(points = (xpts,ypts), layer = layer)
    D.add_port(name = 1, midpoint = (radius*cos(angle1), radius*sin(angle1)),  width = width, orientation = start_angle - 90 + 180*(theta<0))
    D.add_port(name = 2, midpoint = (radius*cos(angle2), radius*sin(angle2)),  width = width, orientation = start_angle + theta + 90 - 180*(theta<0))
    D.info['length'] = (abs(theta)*pi/180)*radius
    return D


def turn(port, radius = 10, angle = 270, angle_resolution = 2.5, layer = 0):
    """ Starting form a port, create a arc which connects to the port
    """

    D = arc(radius = radius, width = port.width, theta = angle, start_angle = 0,
            angle_resolution = angle_resolution, layer = layer)
    D.rotate(angle =  180 + port.orientation - D.ports[1].orientation, center = D.ports[1].midpoint)
    D.move(origin = D.ports[1], destination = port)
    return D


def straight(size = (4,2), layer = 0):
    """Generates a rectangular wire geometry with ports on the length edges.

    Parameters
    ----------
    size : tuple
        The length and width of the rectangle.
    layer : int, array-like[2], or set
        Specific layer(s) to put polygon geometry on.

    Returns
    -------
    A Device with an rectangle polygon in it, and two ports (`1` and `2`) on either end

    Notes
    -----
    Ports are included on both sides of the length edge (i.e. size[0]) of the geometry.
    """

    D = Device(name = 'wire')
    points = [[size[0], size[1]], [size[0], 0], [0, 0], [0, size[1]]]
    D.add_polygon(points, layer = layer)
    D.add_port(name = 1, midpoint = (size[0]/2, size[1]),  width = size[0], orientation = 90)
    D.add_port(name = 2, midpoint = (size[0]/2, 0),  width = size[0], orientation = -90)
    return D


def L(width = 1, size = (10,20) , layer = 0):
    """Generates an "L" geometry with ports on both ends.

    Parameters
    ----------
    width : float
        Width of the L.
    size : tuple
        Lengths of the base and height of the L, respectively.
    layer : int, array-like[2], or set
        Specific layer(s) to put polygon geometry on.

    Returns
    -------
    A Device with an L-shaped polygon in it, and two ports (`1` and `2`) on
    either end of the L
    """

    D = Device(name = 'L')
    w = width/2
    s1, s2 = size
    points = [(-w,-w), (s1,-w), (s1,w), (w,w), (w,s2), (-w,s2), (-w,-w)]
    D.add_polygon(points, layer = layer)
    D.add_port(name = 1, midpoint = (0,s2),  width = width, orientation = 90)
    D.add_port(name = 2, midpoint = (s1, 0),  width = width, orientation = 0)
    return D


def C(width = 1, size = (10,20) , layer = 0):
    """Generates a "C" geometry with ports on both ends.

    Parameters
    ----------
    width : float
        Width of the C.
    size : tuple
        Lengths of the base + top edges and the height of the C, respectively.
    layer : int, array-like[2], or set
        Specific layer(s) to put polygon geometry on.

    Returns
    -------
    A Device with an [-bracket-shaped polygon in it, and two ports (`1` and `2`) on
    either end of the [ shape
    """

    D = Device(name = 'C')
    w = width/2
    s1, s2 = size
    points = [(-w,-w), (s1,-w), (s1,w), (w,w), (w,s2-w), (s1,s2-w), (s1,s2+w), (-w, s2+w), (-w,-w)]
    D.add_polygon(points, layer = layer)
    D.add_port(name = 1, midpoint = (s1,s2),  width = width, orientation = 0)
    D.add_port(name = 2, midpoint = (s1, 0),  width = width, orientation = 0)
    return D



#==============================================================================
#
# Boolean functions
#
#==============================================================================

def offset(elements, distance = 0.1, join_first = True, precision = 1e-4, 
        num_divisions = [1,1],  join='miter', tolerance=2,
        max_points = 4000, layer = 0):
    if type(elements) is not list: elements = [elements]
    polygons_to_offset = []
    for e in elements:
        if isinstance(e, (Device, DeviceReference)): polygons_to_offset += e.get_polygons(by_spec = False)
        elif isinstance(e, (Polygon, gdspy.Polygon)): polygons_to_offset.append(e)
    if len(polygons_to_offset) == 0:
        return Device('offset')
    polygons_to_offset = _merge_floating_point_errors(polygons_to_offset, tol = precision/1000)
    gds_layer, gds_datatype = _parse_layer(layer)
    if all(np.array(num_divisions) == np.array([1,1])):
        p = gdspy.offset(polygons_to_offset, distance = distance, join=join, tolerance=tolerance,
                         precision=precision, join_first=join_first,
                         max_points=max_points, layer=gds_layer, datatype = gds_datatype)
    else:
        p = _offset_polygons_parallel(
            polygons_to_offset,
            distance = distance,
            num_divisions = num_divisions,
            join_first = join_first,
            precision = precision,
            join = join,
            tolerance = tolerance,
            )

    D = Device('offset')
    polygons = D.add_polygon(p, layer=layer)
    [polygon.fracture(max_points = max_points, precision = precision) for polygon in polygons]
    return D


def boolean(A, B, operation, precision = 1e-4, num_divisions = [1,1],
            max_points=4000, layer = 0):
    """
    Performs boolean operations between 2 Device/DeviceReference objects,
    or lists of Devices/DeviceReferences.

    ``operation`` should be one of {'not', 'and', 'or', 'xor', 'A-B', 'B-A', 'A+B'}.  
    Note that 'A+B' is equivalent to 'or', 'A-B' is equivalent to 'not', and
    'B-A' is equivalent to 'not' with the operands switched
    """
    D = Device('boolean')

    A_polys = []
    B_polys = []
    if type(A) is not list: A = [A]
    if type(B) is not list: B = [B]
    for e in A:
        if isinstance(e, Device): A_polys += e.get_polygons()
        elif isinstance(e, DeviceReference): A_polys += e.get_polygons()
    for e in B:
        if isinstance(e, Device): B_polys += e.get_polygons()
        elif isinstance(e, DeviceReference): B_polys += e.get_polygons()

    gds_layer, gds_datatype = _parse_layer(layer)

    operation = operation.lower().replace(' ','')
    if operation == 'a-b':
        operation = 'not'
    elif operation == 'b-a':
        operation = 'not'
        A_polys, B_polys = B_polys, A_polys
    elif operation == 'a+b':
        operation = 'or'
    elif operation not in ['not', 'and', 'or', 'xor', 'a-b', 'b-a', 'a+b']:
        raise ValueError("[PHIDL] phidl.geometry.boolean() `operation` parameter not recognized, must be one of the following:  'not', 'and', 'or', 'xor', 'A-B', 'B-A', 'A+B'")

    # Check for trivial solutions
    if (len(A_polys) == 0) or (len(B_polys) == 0):
        if (operation == 'not'):
            if len(A_polys) == 0:   p = None
            elif len(B_polys) == 0: p = A_polys
        elif (operation == 'and'):
            p = None
        elif (operation == 'or') or (operation == 'xor'):
            if (len(A_polys) == 0) and (len(B_polys) == 0): p = None
            elif len(A_polys) == 0: p = B_polys
            elif len(B_polys) == 0: p = A_polys
    else:
        # If no trivial solutions, run boolean operation either in parallel or straight
        if all(np.array(num_divisions) == np.array([1,1])):
            p = gdspy.boolean(operand1 = A_polys, operand2 = B_polys, operation = operation, precision=precision,
                         max_points=max_points, layer=gds_layer, datatype=gds_datatype)
        else:
            p = _boolean_polygons_parallel(polygons_A = A_polys, polygons_B = B_polys,
                       num_divisions = num_divisions, operation = operation,
                       precision = precision)

    if p is not None:
        polygons = D.add_polygon(p, layer = layer)
        [polygon.fracture(max_points = max_points, precision = precision) for polygon in polygons]
    return D


def outline(elements, distance = 1, precision = 1e-4, num_divisions = [1,1],
    join = 'miter', tolerance = 2, join_first = True, max_points = 4000, layer = 0):
    """ Creates an outline around all the polygons passed in the `elements`
    argument.  `elements` may be a Device, Polygon, or list of Devices
    """
    D = Device('outline')
    if type(elements) is not list: elements = [elements]
    for e in elements:
        if isinstance(e, Device): D.add_ref(e)
        else: D.add(e)
    gds_layer, gds_datatype = _parse_layer(layer)

    D_bloated = offset(D, distance = distance, join_first = join_first,
        num_divisions = num_divisions, precision = precision, max_points = max_points,
        join = join, tolerance = tolerance, layer = layer)
    Outline = boolean(A = D_bloated, B = D, operation = 'A-B', num_divisions = num_divisions,
         max_points = max_points, precision = precision, layer = layer)
    return Outline


def inset(elements, distance = 0.1, join_first = True, precision = 1e-4, layer = 0):
    raise ValueError('[PHIDL] pg.inset() is deprecated, please use pg.offset()')


def invert(elements, border = 10, precision = 1e-4, num_divisions = [1,1], max_points = 4000, layer = 0):
    """ Creates an inverted version of the input shapes with an additional
    border around the edges """
    Temp = Device()
    if type(elements) is not list: elements = [elements]
    for e in elements:
        if isinstance(e, Device): Temp.add_ref(e)
        else: Temp.add(e)
    gds_layer, gds_datatype = _parse_layer(layer)

    # Build the rectangle around the device D
    R = rectangle(size = (Temp.xsize + 2*border, Temp.ysize + 2*border))
    R.center = Temp.center

    D = boolean(A = R, B = Temp, operation = 'A-B', precision = precision,
                num_divisions = num_divisions, max_points = max_points, layer = layer)
    return D

def xor_diff(A,B, precision = 1e-4):
    """ Given two Devices A and B, performs the layer-by-layer XOR
    difference between A and B, and returns polygons representing
    the differences between A and B.
    """
    D = Device()
    A_polys = A.get_polygons(by_spec = True)
    B_polys = B.get_polygons(by_spec = True)
    A_layers = A_polys.keys()
    B_layers = B_polys.keys()
    all_layers = set()
    all_layers.update(A_layers )
    all_layers.update(B_layers)
    for layer in all_layers:
        if (layer in A_layers) and (layer in B_layers):
            p = gdspy.boolean(operand1 = A_polys[layer], operand2 = B_polys[layer],
                                   operation = 'xor', precision=precision,
                                   max_points=4000, layer=layer[0], datatype=layer[1])
        elif (layer in A_layers):
            p = A_polys[layer]
        elif (layer in B_layers):
            p = B_polys[layer]
        if p is not None:
            D.add_polygon(p, layer = layer)
    return D


def union(D, by_layer = False, precision = 1e-4, join_first = True, max_points = 4000, layer = 0):
    U = Device()

    if by_layer == True:
        all_polygons = D.get_polygons(by_spec = True)
        for layer, polygons in all_polygons.items():
            unioned_polygons = _union_polygons(polygons, precision = precision, max_points=max_points)
            U.add_polygon(unioned_polygons, layer = layer)
    else:
        all_polygons = D.get_polygons(by_spec = False)
        unioned_polygons = _union_polygons(all_polygons, precision = precision, max_points=max_points)
        U.add_polygon(unioned_polygons, layer = layer)
    return U

def _union_polygons(polygons, precision = 1e-4, max_points = 4000):
    polygons = _merge_floating_point_errors(polygons, tol = precision/1000)
    unioned = gdspy.boolean(polygons, [], operation = 'or',
                                 precision=precision, max_points=max_points)
    return unioned


def _merge_floating_point_errors(polygons, tol = 1e-10):
    stacked_polygons = np.vstack(polygons)
    x = stacked_polygons[:,0]
    y = stacked_polygons[:,1]
    polygon_indices = np.cumsum([len(p) for p in polygons])

    xfixed = _merge_nearby_floating_points(x, tol = tol)
    yfixed = _merge_nearby_floating_points(y, tol = tol)
    stacked_polygons_fixed = np.vstack([xfixed, yfixed]).T
    polygons_fixed = np.vsplit(stacked_polygons_fixed, polygon_indices[:-1])
    return polygons_fixed


def _merge_nearby_floating_points(x, tol = 1e-10):
    """ Takes an array `x` and merges any values within the tolerance `tol`
    So if given
    >>> x = [-2, -1, 0, 1.0001, 1.0002, 1.0003, 4, 5, 5.003, 6, 7, 8]
    >>> _merge_nearby_floating_points(x, tol = 1e-3)
    will then return:
    >>> [-2, -1, 0, 1.0001, 1.0001, 1.0001, 4, 5, 5.003, 6, 7, 8] """
    xargsort = np.argsort(x)
    xargunsort = np.argsort(xargsort)
    xsort = x[xargsort]
    xsortthreshold = (np.diff(xsort) < tol)
    xsortthresholdind = np.argwhere(xsortthreshold)

    # Merge nearby floating point values
    for xi in xsortthresholdind:
         xsort[xi+1] = xsort[xi]
    return xsort[xargunsort]


def _crop_region(polygons, left, bottom, right, top, precision):
    """ Given a rectangular boundary defined by left/bottom/right/top, this
    takes a list of polygons and cuts them at the boundary, discarding parts
    of the polygons outside the rectangle. """
    cropped_polygons = []
    for p in polygons:
        clipped_polys = clipper._chop(p, [top, bottom], 1, 1 / precision) # polygon, [cuts], axis, scale
        for cp in clipped_polys[1]:
            result = clipper._chop(cp, [left, right], 0, 1 / precision)
            cropped_polygons += list(result[1])
    return cropped_polygons



def _crop_edge_polygons(all_polygons, bboxes, left, bottom, right, top, precision):
    """ Parses out which polygons are along the edge of the rectangle and need 
     to be  cropped and which are deep inside the rectangle region and can be
      left alone, then crops only those polygons along the edge """
    polygons_in_rect_i = _find_bboxes_in_rect(bboxes, left, bottom, right, top)
    polygons_edge_i = _find_bboxes_on_rect_edge(bboxes, left, bottom, right, top)
    polygons_in_rect_no_edge_i = polygons_in_rect_i & (~polygons_edge_i)
    
    # Crop polygons along the edge and recombine them with polygons inside the rectangle
    polygons_edge = all_polygons[polygons_edge_i]
    polygons_in_rect_no_edge = all_polygons[polygons_in_rect_no_edge_i].tolist()
    polygons_edge_cropped = _crop_region(polygons_edge, left, bottom, right, top, precision = precision)
    polygons_to_process = polygons_in_rect_no_edge + polygons_edge_cropped
    
    return polygons_to_process

def _find_bboxes_in_rect(bboxes, left, bottom, right, top):
    """ Given a list of polygon bounding boxes and a rectangle defined by
    left/bottom/right/top, this function returns those polygons which overlap
    the rectangle. """
    result = (bboxes[:,0] <= right) & (bboxes[:,2] >= left) & \
             (bboxes[:,1] <= top)   & (bboxes[:,3] >= bottom)
    return result

# _find_bboxes_on_rect_edge
def _find_bboxes_on_rect_edge(bboxes, left, bottom, right, top):
    """ Given a list of polygon bounding boxes and a rectangular boundary defined by
    left/bottom/right/top, this function returns those polygons which intersect
    the rectangular boundary. """
    bboxes_left   = _find_bboxes_in_rect(bboxes, left, bottom, left, top)
    bboxes_right  = _find_bboxes_in_rect(bboxes, right, bottom, right, top)
    bboxes_top    = _find_bboxes_in_rect(bboxes, left, top, right, top)
    bboxes_bottom = _find_bboxes_in_rect(bboxes, left, bottom, right, bottom)
    result = bboxes_left | bboxes_right | bboxes_top | bboxes_bottom
    return result


def _offset_region(all_polygons, bboxes, left, bottom, right, top,
                distance = 5,
                join_first = True,
                precision = 1e-4,
                join = 'miter',
                tolerance = 2,
                ):
    """ Taking a region of e.g. size (x,y) which needs to be offset by distance d,
    this function crops out a region (x+2*d, y+2*d) large, offsets that region,
    then crops it back to size (x,y) to create a valid result """
        
    # Mark out a region slightly larger than the final desired region
    # FIXME: Necessary?
    d = distance*1.01
    
    polygons_to_offset = _crop_edge_polygons(all_polygons, bboxes, left-d, bottom-d, right+d, top+d, precision = precision)
    
    # Offset the resulting cropped polygons and recrop to final desired size
    polygons_offset = clipper.offset(polygons_to_offset, distance, join, tolerance, 1/precision, int(join_first))
    polygons_offset_cropped = _crop_region(polygons_offset, left, bottom, right, top, precision = precision)
    
    return polygons_offset_cropped


def _polygons_to_bboxes(polygons):
#    Build bounding boxes
    bboxes = np.empty([len(polygons),4])
    for n,p in enumerate(polygons):
        try:
            left,bottom = np.min(p,axis = 0)
        except:
            import pdb
            pdb.set_trace()
        right,top = np.max(p,axis = 0)
        bboxes[n] = [left, bottom, right, top]
    return bboxes

def _offset_polygons_parallel(
    polygons,
    distance = 5,
    num_divisions = [10,10],
    join_first = True,
    precision = 1e-4,
    join = 'miter',
    tolerance = 2,
    ):
    
#    Build bounding boxes
    polygons = np.asarray(polygons)
    bboxes = _polygons_to_bboxes(polygons)
    
    xmin,ymin = np.min(bboxes[:,0:2], axis = 0) - distance
    xmax,ymax = np.max(bboxes[:,2:4], axis = 0) + distance
    
    xsize = xmax - xmin
    ysize = ymax - ymin
    xdelta = xsize/num_divisions[0]
    ydelta = ysize/num_divisions[1]
    xcorners = xmin + np.arange(num_divisions[0])*xdelta
    ycorners = ymin + np.arange(num_divisions[1])*ydelta

    offset_polygons = []
    for n,xc in enumerate(xcorners):
        for m,yc in enumerate(ycorners):
            left = xc
            right = xc+xdelta
            bottom = yc
            top = yc+ydelta
            _offset_region_polygons = _offset_region(polygons, bboxes,
                                            left, bottom, right, top,
                                            distance = distance,
                                            join_first = join_first,
                                            precision = precision,
                                            join = join,
                                            tolerance = tolerance,
                                            )
            offset_polygons += _offset_region_polygons
            
    return offset_polygons


def _boolean_region(all_polygons_A, all_polygons_B,
                    bboxes_A, bboxes_B,
                    left, bottom, right, top,
                operation = 'and',
                precision = 1e-4,
                ):
    """ Taking a region of e.g. size (x,y) which needs to be booleaned,
    this function crops out a region (x, y) large from each set of polygons
    (A and B), booleans that cropped region and returns the result"""
        
    polygons_to_boolean_A = _crop_edge_polygons(all_polygons_A, bboxes_A, left, bottom, right, top, precision)
    polygons_to_boolean_B = _crop_edge_polygons(all_polygons_B, bboxes_B, left, bottom, right, top, precision)
    polygons_boolean = clipper.clip(polygons_to_boolean_A, polygons_to_boolean_B,
                                         operation, 1/precision)
    return polygons_boolean



def _boolean_polygons_parallel(
        polygons_A, polygons_B,
        num_divisions = [10,10],
        operation = 'and',
        precision = 1e-4,
        ):
    
    #    Build bounding boxes
    polygons_A = np.asarray(polygons_A)
    polygons_B = np.asarray(polygons_B)
    bboxes_A = _polygons_to_bboxes(polygons_A)
    bboxes_B = _polygons_to_bboxes(polygons_B)
    
    xmin,ymin = np.min([np.min(bboxes_A[:,0:2], axis = 0), np.min(bboxes_B[:,0:2], axis = 0)], axis = 0)
    xmax,ymax = np.max([np.max(bboxes_A[:,2:4], axis = 0), np.max(bboxes_B[:,2:4], axis = 0)], axis = 0)
    
    xsize = xmax - xmin
    ysize = ymax - ymin
    xdelta = xsize/num_divisions[0]
    ydelta = ysize/num_divisions[1]
    xcorners = xmin + np.arange(num_divisions[0])*xdelta
    ycorners = ymin + np.arange(num_divisions[1])*ydelta
    
    boolean_polygons = []
    for n,xc in enumerate(xcorners):
        for m,yc in enumerate(ycorners):
            left = xc
            right = xc+xdelta
            bottom = yc
            top = yc+ydelta
            _boolean_region_polygons = _boolean_region(polygons_A, polygons_B, bboxes_A, bboxes_B,
                                            left, bottom, right, top,
                                            operation = operation,
                                            precision = precision,
                                            )
            boolean_polygons += _boolean_region_polygons
        
    return boolean_polygons

#==============================================================================
#
# Lithography test structures
#
#==============================================================================



def litho_steps(
        line_widths = [1,2,4,8,16],
        line_spacing = 10,
        height = 100,
        layer = 0
        ):
    """ Produces a positive + negative tone linewidth test, used for
    lithography resolution test patterning """
    D = Device('litho_steps')

    height = height / 2
    T1 = text(text = '%s' % str(line_widths[-1]),
        size = height, justify = 'center', layer = layer)
    t1 = D.add_ref(T1).rotate(90).movex(-height/10)
    R1 = rectangle(size = (line_spacing, height), layer = layer)
    r1 = D.add_ref(R1).movey(-height)
    count = 0
    for i in reversed(line_widths):
        count += line_spacing + i
        R2 = rectangle(size = (i, height), layer = layer)
        r1 = D.add_ref(R1).movex(count).movey(-height)
        r2 = D.add_ref(R2).movex(count - i)

    return(D)


def litho_star(
        num_lines = 20,
        line_width = 2,
        diameter = 200,
        layer = 0
        ):
    """ Creates a circular-star shape from lines, used as a lithographic
    resolution test pattern """
    D = Device('litho_star')

    degree = 180 / num_lines
    R1 = rectangle(size = (line_width, diameter), layer = layer)
    for i in range(num_lines):
        r1 = D.add_ref(R1).rotate(degree * i)
        r1.center = (0,0)

    return(D)


def litho_calipers(
        notch_size = [2,5],
        notch_spacing = 2,
        num_notches = 11,
        offset_per_notch = 0.1,
        row_spacing = 0,
        layer1 = 1,
        layer2 = 2):
    """ Creates a vernier caliper structure for lithography alignment
    tests.  Vernier structure is made horizontally. """

    D = Device('litho_calipers')
    num_notches_total = num_notches*2+1
    centre_notch = num_notches
    R1 = rectangle(size = (notch_size), layer = layer1)
    R2 = rectangle(size = (notch_size), layer = layer2)
    for i in range(num_notches_total):
        if i == centre_notch:
            r1 = D.add_ref(R1).movex(i * (notch_size[0] + notch_spacing)).movey(notch_size[1])
            r2 = D.add_ref(R2).movex(i * (notch_size[0] + notch_spacing) + offset_per_notch * (centre_notch - i)).movey(-2 * notch_size[1] - row_spacing)
        r1 = D.add_ref(R1).movex(i * (notch_size[0] + notch_spacing))
        r2 = D.add_ref(R2).movex(i * (notch_size[0] + notch_spacing) + offset_per_notch * (centre_notch - i)).movey(-notch_size[1] - row_spacing)

    return(D)



#==============================================================================
#
# Utility functions
#
#==============================================================================



def extract(D, layers = [0,1]):
    D_extracted = Device('extract')
    if type(layers) not in (list, tuple):
        raise ValueError('[PHIDL] pg.extract() Argument `layers` needs to be passed a list or tuple')
    poly_dict = D.get_polygons(by_spec = True)
    parsed_layer_list = [_parse_layer(layer) for layer in layers]
    for layer, polys in poly_dict.items():
        if _parse_layer(layer) in parsed_layer_list:
            D_extracted.add_polygon(polys, layer = layer)
    return D_extracted


def copy(D):
    D_copy = Device(name = D._internal_name)
    D_copy.info = python_copy.deepcopy(D.info)
    for ref in D.references:
        new_ref = DeviceReference(device = ref.parent,
                                origin = ref.origin,
                                rotation = ref.rotation,
                                magnification = ref.magnification,
                                x_reflection = ref.x_reflection)
        new_ref.owner = D_copy
        D_copy.add(new_ref)
        for alias_name, alias_ref in D.aliases.items():
            if alias_ref == ref: D_copy.aliases[alias_name] = new_ref

    for port in D.ports.values():      D_copy.add_port(port = port)
    for poly in D.polygons:   D_copy.add_polygon(poly)
    for label in D.labels:    D_copy.add_label(text = label.text,
                                           position = label.position,
                                           layer = (label.layer, label.texttype))
    return D_copy


def deepcopy(D):
    D_copy = python_copy.deepcopy(D)
    D_copy.uid = Device._next_uid
    Device._next_uid += 1
    D_copy._internal_name = D._internal_name
    D_copy.name = '%s%06d' % (D_copy._internal_name[:20], D_copy.uid) # Write name e.g. 'Unnamed000005'
    # Make sure _bb_valid is set to false for these new objects so new
    # bounding boxes are created in the cache
    for D in D_copy.get_dependencies(True):
        D._bb_valid = False
    D_copy._bb_valid = False
    
    return D_copy


def copy_layer(D, layer = 1, new_layer = 2):
    D_copied_layer = extract(D, layers = [layer])
    D_copied_layer.flatten(single_layer = new_layer)
    return D_copied_layer


def import_gds(filename, cellname = None, flatten = False):
    gdsii_lib = gdspy.GdsLibrary()
    gdsii_lib.read_gds(filename)
    top_level_cells = gdsii_lib.top_level()
    if cellname is not None:
        if cellname not in gdsii_lib.cells:
            raise ValueError('[PHIDL] import_gds() The requested cell (named %s) \
                        is not present in file %s' % (cellname,filename))
        topcell = gdsii_lib.cells[cellname]
    elif cellname is None and len(top_level_cells) == 1:
        topcell = top_level_cells[0]
    elif cellname is None and len(top_level_cells) > 1:
        raise ValueError('[PHIDL] import_gds() There are multiple top-level cells, \
                        you must specify `cellname` to select of one of them')

    if flatten == False:
        D_list = []
        c2dmap = {}
        for cell in gdsii_lib.cells.values():
            D = Device(name = cell.name)
            D.polygons = cell.polygons
            D.references = cell.references
            D.name = cell.name
            D.labels = cell.labels
            c2dmap.update({cell:D})
            D_list += [D]

        for D in D_list:
            # First convert each reference so it points to the right Device
            converted_references = []
            for e in D.references:
                ref_device = c2dmap[e.ref_cell]
                if isinstance(e, gdspy.CellReference):
                    dr = DeviceReference(
                        device = ref_device,
                        origin = e.origin,
                        rotation = e.rotation,
                        magnification = e.magnification,
                        x_reflection = e.x_reflection,
                        )
                    dr.owner = D
                    converted_references.append(dr)
                elif isinstance(e, gdspy.CellArray):
                    dr = CellArray(
                        device = ref_device,
                        columns = e.columns,
                        rows = e.rows,
                        spacing = e.spacing,
                        origin = e.origin,
                        rotation = e.rotation,
                        magnification = e.magnification,
                        x_reflection = e.x_reflection,
                        )
                    dr.owner = D
                    converted_references.append(dr)
            D.references = converted_references
            # Next convert each Polygon
            temp_polygons = list(D.polygons)
            D.polygons = []
            for p in temp_polygons:
                D.add_polygon(p)

        topdevice = c2dmap[topcell]
        return topdevice

    elif flatten == True:
        D = Device('import_gds')
        polygons = topcell.get_polygons(by_spec = True)

        for layer_in_gds, polys in polygons.items():
            D.add_polygon(polys, layer = layer_in_gds)
        return D


def _translate_cell(c):
    D = Device(name = c.name)
    for e in c.elements:
        if isinstance(e, gdspy.PolygonSet):
            for n, points in enumerate(e.polygons):
                polygon_layer = _parse_layer((e.layers[n], e.datatypes[n]))
                D.add_polygon(points = points, layer = polygon_layer)
        elif isinstance(e, gdspy.CellReference):
            dr = DeviceReference(device = _translate_cell(e.ref_cell),
                            origin = e.origin,
                            rotation = e.rotation, magnification = None,
                            x_reflection = e.x_reflection)
            D.elements.append(dr)
    D.labels = c.labels
    return D


def preview_layerset(ls, size = 100, spacing = 100):
    """ Generates a preview Device with representations of all the layers,
    used for previewing LayerSet color schemes in quickplot or saved .gds
    files """
    D = Device()
    scale = size/100
    num_layers = len(ls._layers)
    matrix_size = int(np.ceil(np.sqrt(num_layers)))
    sorted_layers = sorted(ls._layers.values(), key = lambda x: (x.gds_layer, x.gds_datatype))
    for n, layer in enumerate(sorted_layers):
        R = rectangle(size = (100*scale, 100*scale), layer = layer)
        T = text(
                text = '%s\n%s / %s' % (layer.name, layer.gds_layer, layer.gds_datatype),
                size = 20*scale,
                justify = 'center',
                layer = layer)

        T.move((50*scale,-20*scale))
        xloc = n % matrix_size
        yloc = int(n // matrix_size)
        D.add_ref(R).movex((100+spacing) * xloc *scale).movey(-(100+spacing) * yloc*scale)
        D.add_ref(T).movex((100+spacing) * xloc *scale).movey(-(100+spacing) * yloc*scale)
    return D

class device_lru_cache:
    def __init__(self, fn):
        self.maxsize = 32
        self.fn = fn
        self.memo = OrderedDict()
        update_wrapper(self, fn)

    def __call__(self, *args, **kwargs):
        pickle_str = pickle.dumps(args, 1) + pickle.dumps(kwargs, 1)
        if pickle_str not in self.memo.keys():
            new_cache_item = self.fn(*args, **kwargs)
            if not isinstance(new_cache_item, Device):
                raise ValueError('[PHIDL] @device_lru_cache can only be used on functions which return a Device')
            if len(self.memo) > self.maxsize:
                self.memo.popitem(last = False) # Remove oldest item from cache
            # Add a deepcopy of new item to cache so that if we change the
            # returned device, our stored cache item is not changed
            self.memo[pickle_str] = python_copy.deepcopy(new_cache_item)
            return new_cache_item
        else: # if found in cache
            # Pop cache item out and put it back on the top of the cache
            cached_output = self.memo.pop(pickle_str)
            self.memo[pickle_str] = cached_output
            # Then return a copy of the cached Device
            return deepcopy(cached_output)


def _convert_port_to_geometry(port, layer = 0):
    ''' Converts a Port to a label and a triangle Device that are then added to the parent.
        The Port must start with a parent.
    '''
    if port.parent is None:
        raise ValueError('Port {}: Port needs a parent in which to draw'.format(port.name))
    if isinstance(port.parent, DeviceReference):
        device = port.parent.parent
    else:
        device = port.parent

    # A visual marker
    triangle_points = [[0, 0]] * 3
    triangle_points[0] = port.endpoints[0]
    triangle_points[1] = port.endpoints[1]
    triangle_points[2] = (port.midpoint + (port.normal - port.midpoint) * port.width / 10)[1]
    device.add_polygon(triangle_points, layer)

    # Label carrying actual information that will be recovered
    label_contents = (str(port.name),
                      # port.midpoint,  # rather than put this in the text, use the label position
                      float(np.round(port.width, decimals=3)),  # this can have rounding errors that are less than a nanometer
                      float(port.orientation),
                      # device,  # this is definitely not serializable
                      # port.info,  # would like to include, but it might go longer than 1024 characters
                      # port.uid,  # not including because it is part of the build process, not the port state
                     )
    label_text = json.dumps(label_contents)
    device.add_label(text = label_text, position = port.midpoint + _calculate_label_offset(port),
                      magnification = .04 * port.width, rotation = (90 + port.orientation) % 360,
                      layer = layer)


def _calculate_label_offset(port):
    ''' Used to put the label in a pretty position.
        It is added when drawing and substracted when extracting.
    '''
    offset_position = np.array((-np.cos(np.pi / 180 * port.orientation),
                                -np.sin(np.pi / 180 * port.orientation)))
    offset_position *= port.width * .05
    return offset_position


def _convert_geometry_to_port(label, layer = 0):
    ''' Converts a label into a Port in the parent Device.
        The label contains name, width, orientation.
        Does not remove that label from the parent.
        Returns the new port.
    '''
    name, width, orientation = json.loads(label.text)
    new_port = Port(name=name, width=width, orientation=orientation)
    new_port.midpoint = label.position - _calculate_label_offset(new_port)
    return new_port


def ports_to_geometry(device, layer = 0):
    ''' Converts Port objects over the whole Device hierarchy to geometry and labels.
        layer: the special port record layer
        Does not change the device used as argument. Returns a new one lacking all Ports.
    '''
    temp_device = deepcopy(device)
    all_cells = list(temp_device.get_dependencies(recursive=True))
    all_cells.append(temp_device)
    for subcell in all_cells:
        for port in subcell.ports.values():
            _convert_port_to_geometry(port, layer=layer)
            subcell.remove(port)
    return temp_device


def geometry_to_ports(device, layer = 0):
    ''' Converts geometry representing ports over the whole Device hierarchy into Port objects.
        layer: the special port record layer
        Does not mutate the device in the argument. Returns a new one lacking all port geometry (incl. labels)
    '''
    temp_device = deepcopy(device)
    all_cells = list(temp_device.get_dependencies(recursive=True))
    all_cells.append(temp_device)
    for subcell in all_cells: # Walk through cells
        for lab in subcell.labels:
            if lab.layer == layer:
                the_port = _convert_geometry_to_port(lab)
                subcell.add_port(name=the_port.name, port=the_port)
    temp_device.remove_layers(layers=[layer], include_labels=True)
    return temp_device


#==============================================================================
#
# Connectors
#
#==============================================================================

def connector(midpoint = (0,0), width = 1, orientation = 0):
    """ Creates a Device which has back-to-back ports """
    D = Device(name = 'connector')
    D.add_port(name = 1, midpoint = [midpoint[0], midpoint[1]],  width = width, orientation = orientation)
    D.add_port(name = 2, midpoint = [midpoint[0], midpoint[1]],  width = width, orientation = orientation-180)
    return D



#==============================================================================
#
# Contact pads
#
#==============================================================================


def compass(size = (4,2), layer = 0):
    """ Creates a rectangular contact pad with centered ports on edges of the
    rectangle (north, south, east, and west)
    """

    D = Device(name = 'compass')
    r = D.add_ref( rectangle(size, layer = layer) )
    r.center = (0,0)

    dx = size[0]
    dy = size[1]
    D.add_port(name = 'N', midpoint = [0, dy/2],  width = dx, orientation = 90)
    D.add_port(name = 'S', midpoint = [0, -dy/2], width = dx, orientation = -90)
    D.add_port(name = 'E', midpoint = [dx/2, 0],  width = dy, orientation = 0)
    D.add_port(name = 'W', midpoint = [-dx/2, 0], width = dy, orientation = 180)

    return D


def compass_multi(size = (4,2), ports = {'N':3,'S':4}, layer = 0):
    """ Creates a rectangular contact pad with multiple ports along the edges
    rectangle (north, south, east, and west).
    """

    D = Device(name = 'compass_multi')
    r = D.add_ref( rectangle(size, layer = layer) )
    r.center = (0,0)

    dx = size[0]/2
    dy = size[1]/2

    if 'N' in ports:
        num_ports = ports['N']
        m = dx-dx/num_ports
        p_list = np.linspace(-m, m, num_ports)
        [D.add_port(name = ('N%s' % (n+1)), midpoint = [p, dy],  width = dx/num_ports*2, orientation = 90) for n,p in enumerate(p_list)]
    if 'S' in ports:
        num_ports = ports['S']
        m = dx-dx/num_ports
        p_list = np.linspace(-m, m, num_ports)
        [D.add_port(name = ('S%s' % (n+1)), midpoint = [p, -dy],  width = dx/num_ports*2, orientation = -90) for n,p in enumerate(p_list)]
    if 'E' in ports:
        num_ports = ports['E']
        m = dy-dy/num_ports
        p_list = np.linspace(-m, m, num_ports)
        [D.add_port(name = ('E%s' % (n+1)), midpoint = [dx, p],  width = dy/num_ports*2, orientation = 0) for n,p in enumerate(p_list)]
    if 'W' in ports:
        num_ports = ports['W']
        m = dy-dy/num_ports
        p_list = np.linspace(-m, m, num_ports)
        [D.add_port(name = ('W%s' % (n+1)), midpoint = [-dx, p],  width = dy/num_ports*2, orientation = 180) for n,p in enumerate(p_list)]

    return D



# TODO: Fix the fillet here, right now only goes halfway down
def flagpole(size = (4,2), stub_size = (2,1), shape = 'p', taper_type = 'straight', layer = 0):
    f = np.array(size)
    p = np.array(stub_size)
    shape = shape.lower()

    assert shape in 'pqbd', '[DEVICE]  flagpole() shape must be p, q, b, or d'
    assert taper_type in ['straight','fillet',None], '[DEVICE]  flagpole() taper_type must "straight" or "fillet" or None'

    if shape ==   'p':
        orientation = -90
    elif shape == 'q':
        f[0], p[0] = -size[0], -stub_size[0]
        orientation = -90
    elif shape == 'b':
        f[1], p[1] = -size[1], -stub_size[1]
        orientation = 90
    elif shape == 'd':
        f[1], p[1] = -size[1], -stub_size[1]
        f[0], p[0] = -size[0], -stub_size[0]
        orientation = 90
    xpts = [0, 0, f[0], f[0], p[0], p[0], 0]
    ypts = [0, f[1], f[1], 0, 0, -p[1], -p[1]]

    D = Device(name = 'flagpole')
    pad_poly = D.add_polygon([xpts,ypts], layer = layer)
    if taper_type == 'fillet':
        taper_amount = min([abs(f[0]-p[0]), abs(p[1])])
        pad_poly.fillet([0,0,0,0,taper_amount,0,0])
    elif taper_type == 'straight':
        D.add_polygon([xpts[3:6],ypts[3:6]], layer = layer)

    D.add_port(name = 1, midpoint = [p[0]/2, -p[1]],  width = abs(p[0]), orientation = orientation)
    D.add_port(name = 2, midpoint = [f[0]/2, f[1]],  width = abs(f[0]), orientation = orientation-180)
    return D


def tee(size = (4,2), stub_size = (2,1), taper_type = None, layer = 0):
    f = np.array(size)
    p = np.array(stub_size)

    xpts = np.array([f[0], f[0], p[0], p[0], -p[0], -p[0], -f[0], -f[0]])/2
    ypts = [f[1], 0, 0, -p[1], -p[1], 0, 0, f[1]]

    D = Device(name = 'tee')
    pad_poly = D.add_polygon([xpts,ypts], layer = layer)
    if taper_type == 'fillet':
        taper_amount = min([abs(f[0]-p[0]), abs(p[1])])
        pad_poly.fillet([0,0,taper_amount,0,0,taper_amount,0,0])
    elif taper_type == 'straight':
        taper_poly1 = D.add_polygon([xpts[1:4],ypts[1:4]], layer = layer)
        taper_poly2 = D.add_polygon([xpts[4:7],ypts[4:7]], layer = layer)

    D.add_port(name = 1, midpoint = [f[0]/2, f[1]/2],  width = f[1], orientation = 0)
    D.add_port(name = 2, midpoint = [-f[0]/2, f[1]/2],  width = f[1], orientation = 180)
    D.add_port(name = 3, midpoint = [0, -p[1]],  width = p[0], orientation = -90)
    return D


#==============================================================================
# Example code
#==============================================================================

#cp = compass(size = [4,2])
#quickplot(cp)


#cpm = compass_multi(size = [40,20], ports = {'N':3,'S':4, 'E':1, 'W':8}, layer = 0)
#quickplot(cpm)


#cpm = compass_multi(size = [40,20], ports = {'N':3,'S':4, 'E':1, 'W':8}, layer = 0)
#inset_polygon = offset(cpm, distance = -2, layer = 1)
#cpm.add_polygon(inset_polygon)
#quickplot(cpm)

#fp = flagpole(size = [4,2], stub_size = [2,1], shape = 'p', taper_type = 'straight', layer = 0)
#quickplot(fp)


#tp = tee(size = [4,2], stub_size = [2,1], taper_type = 'fillet', layer = 0)
#quickplot(tp)







#==============================================================================
#
# Tapers
#
#==============================================================================

# TODO change this so "width1" and "width2" arguments can accept Port directly
def taper(length = 10, width1 = 5, width2 = None, port = None, layer = 0):
    if type(port) is Port and width1 is None: width1 = port.width
    if width2 is None: width2 = width1
    xpts = [0, length, length, 0]
    ypts = [width1/2, width2/2, -width2/2, -width1/2]

    D = Device('taper')
    D.add_polygon([xpts,ypts], layer = layer)
    D.add_port(name = 1, midpoint = [0, 0], width = width1, orientation = 180)
    D.add_port(name = 2, midpoint = [length, 0], width = width2, orientation = 0)
    if type(port) is Port:
        D.rotate(angle = port.orientation, center = [0,0])
        D.move(origin = [0,0], destination = port.midpoint)
    return D


def ramp(length = 10, width1 = 5, width2 = 8, layer = 0):
    if width2 is None: width2 = width1
    xpts = [0, length, length, 0]
    ypts = [width1, width2, 0, 0]
    D = Device('ramp')
    D.add_polygon([xpts,ypts], layer = layer)
    D.add_port(name = 1, midpoint = [0, width1/2], width = width1, orientation = 180)
    D.add_port(name = 2, midpoint = [length, width2/2], width = width2, orientation = 0)
    return D


# Equations taken from
# Hammerstad, E., & Jensen, O. (1980). Accurate Models for Microstrip
# Computer-Aided Design.  http://doi.org/10.1109/MWSYM.1980.1124303
def _microstrip_Z(wire_width, dielectric_thickness, eps_r):
    # Note these equations can be further corrected for thick films (Hammersted Eqs 6-9)
    # and also for frequency since microstrips are dispersive  (Hammersted Eqs 10-12)

    u = wire_width/dielectric_thickness
    eta = 376.73 # Vacuum impedance

    a = 1 + log((u**4 + (u/52)**2)/(u**4 + 0.432))/49 + log(1 + (u/18.1)**3)/18.7;
    b = 0.564*((eps_r-0.9)/(eps_r+3))**0.053;
    F = 6 + (2*pi-6)*exp(-(30.666/u)**0.7528);
    eps_eff = 0.5*(eps_r+1) + 0.5*(eps_r-1)*(1 + 10/u)**(-a*b);
    Z = eta/(2*pi) * log(F/u + sqrt(1+(2/u)**2)) /sqrt(eps_eff);
    return Z,eps_eff


def _microstrip_LC_per_meter(wire_width, dielectric_thickness, eps_r):
    # Use the fact that v = 1/sqrt(L_m*C_m) = 1/sqrt(eps*mu) and
    # Z = sqrt(L_m/C_m)   [Where L_m is inductance per meter]

    Z, eps_eff =  _microstrip_Z(wire_width, dielectric_thickness, eps_r)
    eps0 =  8.854e-12
    mu0 = 4*pi*1e-7

    eps = eps_eff*eps0
    mu = mu0
    L_m = sqrt(eps*mu)*Z
    C_m = sqrt(eps*mu)/Z
    return L_m, C_m


def _microstrip_Z_with_Lk(wire_width, dielectric_thickness, eps_r, Lk_per_sq):
    # Add a kinetic inductance and recalculate the impedance, be careful
    # to input Lk as a per-meter inductance

    L_m, C_m = _microstrip_LC_per_meter(wire_width, dielectric_thickness, eps_r)
    Lk_m = Lk_per_sq*(1.0/wire_width)
    Z = sqrt((L_m+Lk_m)/C_m)
    return Z

def _microstrip_v_with_Lk(wire_width, dielectric_thickness, eps_r, Lk_per_sq):
    L_m, C_m = _microstrip_LC_per_meter(wire_width, dielectric_thickness, eps_r)
    Lk_m = Lk_per_sq*(1.0/wire_width)
    v = 1/sqrt((L_m+Lk_m)*C_m)
    return v

def _find_microstrip_wire_width(Z_target, dielectric_thickness, eps_r, Lk_per_sq):

    def error_fun(wire_width):
        Z_guessed = _microstrip_Z_with_Lk(wire_width, dielectric_thickness, eps_r, Lk_per_sq)
        return (Z_guessed-Z_target)**2 # The error

    x0 = dielectric_thickness
    try:
        from scipy.optimize import fmin
    except:
        raise ImportError(""" [PHIDL] To run the microsctrip functions you need scipy,
          please install it with `pip install scipy` """)
    w = fmin(error_fun, x0, args=(), disp=False)
    return w[0]

def _G_integrand(xip, B):
    try:
        from scipy.special import iv as besseli
    except:
        """ [PHIDL] To run this function you need scipy, please install it with
        pip install scipy """
    return besseli(0, B*sqrt(1-xip**2))


def _G(xi, B):
    try:
        import scipy.integrate
    except:
        raise ImportError(""" [PHIDL] To run the microsctrip functions you need scipy,
          please install it with `pip install scipy` """)
    return B/sinh(B)*scipy.integrate.quad(_G_integrand, 0, xi, args = (B))[0]

@device_lru_cache
def hecken_taper(length = 200, B = 4.0091, dielectric_thickness = 0.25, eps_r = 2,
                 Lk_per_sq = 250e-12, Z1 = None, Z2 = None, width1 = None, width2 = None,
                 num_pts = 100, layer = 0):
    if width1 is not None:  Z1 = _microstrip_Z_with_Lk(width1*1e-6, dielectric_thickness*1e-6, eps_r, Lk_per_sq)
    if width2 is not None:  Z2 = _microstrip_Z_with_Lk(width2*1e-6, dielectric_thickness*1e-6, eps_r, Lk_per_sq)
    xi_list = np.linspace(-1,1, num_pts) # Normalized length of the wire [-1 to +1]
    Z = [np.exp( 0.5*log(Z1*Z2) + 0.5*log(Z2/Z1)*_G(xi, B) ) for xi in xi_list]
    widths = np.array([_find_microstrip_wire_width(z, dielectric_thickness*1e-6, eps_r, Lk_per_sq)*1e6 for z in Z])
    x = ((xi_list/2)*length)

    # Compensate for varying speed of light in the microstrip by shortening
    # and lengthening sections according to the speed of light in that section
    v = np.array([_microstrip_v_with_Lk(w*1e-6, dielectric_thickness*1e-6, eps_r, Lk_per_sq) for w in widths])
    dx = np.diff(x)
    dx_compensated = dx*v[:-1]
    x_compensated = np.cumsum(dx_compensated)
    x = np.hstack([0,x_compensated])/max(x_compensated)*length

    # Create blank device and add taper polygon
    D = Device('hecken')
    xpts = np.concatenate([x, x[::-1]])
    ypts = np.concatenate([widths/2, -widths[::-1]/2])
    D.add_polygon((xpts,ypts), layer = layer)
    D.add_port(name = 1, midpoint = (0,0), width = widths[0], orientation = 180)
    D.add_port(name = 2, midpoint = (length,0), width = widths[-1], orientation = 0)

    # Add meta information about the taper
    D.info['num_squares'] = np.sum(np.diff(x)/widths[:-1])
    D.info['width1'] = widths[0]
    D.info['width2'] = widths[-1]
    D.info['Z1'] = Z[0]
    D.info['Z2'] = Z[-1]
    # Note there are two values for v/c (and f_cutoff) because the speed of
    # light is different at the beginning and end of the taper
    D.info['w'] = widths
    D.info['x'] = x
    D.info['Z'] = Z
    D.info['v/c'] = v/3e8
    D.info['time_length'] = np.sum(np.diff(D.info['x']*1e-6)/(D.info['v/c'][:-1]*3e8))
    BetaLmin = np.sqrt(B**2 + 6.523)
    D.info['f_cutoff'] = 1/(2*D.info['time_length'])
    D.info['length'] = length

    return D


@device_lru_cache
def meander_taper(x_taper, w_taper, meander_length = 1000, spacing_factor = 3,
                  min_spacing = 0.5, layer = 0):

    def taper_width(x):
        return np.interp(x, x_taper, w_taper)


    def taper_section(x_start, x_end, num_pts = 30, layer = 0):
        D = Device('tapersec')
        length =  x_end - x_start
        x = np.linspace(0, length, num_pts)
        widths = np.linspace(taper_width(x_start), taper_width(x_end), num_pts)
        xpts = np.concatenate([x, x[::-1]])
        ypts = np.concatenate([widths/2, -widths[::-1]/2])
        D.add_polygon((xpts,ypts), layer = layer)
        D.add_port(name = 1, midpoint = (0,0), width = widths[0], orientation = 180)
        D.add_port(name = 2, midpoint = (length,0), width = widths[-1], orientation = 0)
        return D

    def arc_tapered(radius = 10, width1 = 1, width2 = 2, theta = 45, angle_resolution = 2.5, layer = 0):
        D = Device('arctaper')
        path1 = gdspy.Path(width = width1, initial_point = (0, 0))
        path1.turn(radius = radius, angle = theta*np.pi/180, number_of_points=int(abs(2*theta/angle_resolution)), final_width = width2)
        [D.add_polygon(p, layer = layer) for p in path1.polygons]
        D.add_port(name = 1, midpoint = (0, 0), width = width1, orientation = 180)
        D.add_port(name = 2, midpoint = (path1.x, path1.y), width = width2, orientation = path1.direction*180/np.pi)
        return D

    D = Device('meander-taper')
    xpos1 = min(x_taper)
    xpos2 = min(x_taper) + meander_length
    t = D.add_ref( taper_section(x_start = xpos1, x_end = xpos2, num_pts = 50, layer = layer) )
    D.add_port(t.ports[1])
    dir_toggle = -1
    while xpos2 < max(x_taper):
        arc_width1 = taper_width(xpos2)
        arc_radius = max(spacing_factor*arc_width1, min_spacing)
        arc_length = np.pi*arc_radius
        arc_width2 = taper_width(xpos2 + arc_length)
        A = arc_tapered(radius = arc_radius, width1 = arc_width1,
                        width2 = arc_width2, theta = 180*dir_toggle, layer = layer)
        a = D.add_ref(A)
        a.connect(port = 1, destination = t.ports[2])
        dir_toggle = -dir_toggle
        xpos1 = xpos2 + arc_length
        xpos2 = xpos1 + meander_length
        t = D.add_ref( taper_section(x_start = xpos1, x_end = xpos2, num_pts = 30, layer = layer) )
        t.connect(port = 1, destination = a.ports[2])
    D.add_port(t.ports[2])

    return D

#==============================================================================
# Example code
#==============================================================================

#D = racetrack_gradual(width, R = 5, N=3)
#quickplot(D)

# D = hecken_taper(length = 200, B = 4.0091, dielectric_thickness = 0.25, eps_r = 2,
#                  Lk_per_sq = 250e-12, Z1 = 50, width2 = 0.3,
#                  num_pts = 100, layer = 0)
# quickplot(D)

#t = np.linspace(0,1)
#x,y = _racetrack_gradual_parametric(t, R = 5, N = 3)
#plt.plot(x,y)
#plt.axis('equal')





#==============================================================================
#
# Text
#
#==============================================================================


def text(text = 'abcd', size = 10, justify = 'left', layer = 0):
    scaling = size/1000
    position = (0,0)
    xoffset = 0
    yoffset = 0
    t = Device('text')
    for line in text.split('\n'):
        l = Device(name = 'textline')
        for c in line:
            ascii_val = ord(c)
            if c == ' ':
                xoffset += 500*scaling
            elif (33 <= ascii_val <= 126) or (ascii_val == 181):
                for poly in _glyph[ascii_val]:
                    xpts = np.array(poly)[:,0]*scaling
                    ypts = np.array(poly)[:,1]*scaling
                    l.add_polygon([xpts + xoffset,ypts + yoffset], layer=layer)
                xoffset += (_width[ascii_val] + _indent[ascii_val])*scaling
            else:
                valid_chars = '!"#$%&\'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\\]^_`abcdefghijklmnopqrstuvwxyz{|}~µ'
                warnings.warn('[PHIDL] text(): Warning, some characters ignored, no geometry for character "%s" with ascii value %s. ' \
                'Valid characters: %s'  % (chr(ascii_val), ascii_val,valid_chars))
        t.add_ref(l)
        yoffset -= 1500*scaling
        xoffset = position[0]
    justify = justify.lower()
    for l in t.references:
        if justify == 'left':   pass
        if justify == 'right':  l.xmax = position[0]
        if justify == 'center': l.move(origin = l.center, destination = position, axis = 'x')
    t.flatten()
    return t


#==============================================================================
# Example code
#==============================================================================

#D = text('the quick brown\n fox jumped over\nthe lazy dog', justify = 'center', size = 800)
#quickplot(D)




#==============================================================================
#
# Wafer and die
#
#==============================================================================



def basic_die(
              size = (10000, 10000),
              street_width = 100,
              street_length = 1000,
              die_name = 'chip99',
              text_size = 100,
              text_location = 'SW',
              layer = 0,
              draw_bbox = True,
              bbox_layer = 99,
              ):
    D = Device(name = 'die')
    sx, sy = size[0]/2, size[1]/2
    xpts = np.array([sx, sx, sx-street_width, sx-street_width,
        sx-street_length, sx-street_length])
    ypts = np.array([sy, sy-street_length, sy-street_length,
        sy-street_width, sy-street_width, sy])
    D.add_polygon([xpts,ypts], layer = layer)
    D.add_polygon([-xpts,ypts], layer = layer)
    D.add_polygon([xpts,-ypts], layer = layer)
    D.add_polygon([-xpts,-ypts], layer = layer)

    if draw_bbox is True:
        D.add_polygon([[sx,sy], [sx,-sy],[-sx,-sy],[-sx,sy]], layer = bbox_layer)
    D.center = (0,0)
    t = D.add_ref(text(text = die_name, size = text_size, layer=layer))

    d = street_width + 20
    if type(text_location) is str:
        text_location = text_location.upper()
        if text_location == 'NW':
            t.xmin, t.ymax = [-sx + d, sy - d]
        elif text_location == 'N':
            t.x, t.ymax = [0, sy - d]
        elif text_location == 'NE':
            t.xmax, t.ymax = [sx - d, sy - d]
        if text_location == 'SW':
            t.xmin, t.ymin = [-sx + d, -sy + d]
        elif text_location == 'S':
            t.x, t.ymin = [0, -sy + d]
        elif text_location == 'SE':
            t.xmax, t.ymin = [sx - d, -sy + d]
    else:
        t.x, t.y = text_location

    return D



#==============================================================================
# Example code
#==============================================================================

# D = basic_die(size = (10000, 10000), street_width = 100, street_length = 1000,
#               die_name = 'chip99', text_size = 300, text_location = 'SW',  layer = 0,
#               draw_bbox = True,  bbox_layer = 99)
# quickplot(D)





#==============================================================================
#
# Waveguide curves
#
#==============================================================================


def racetrack_gradual(width = 0.3, R = 5, N = 3, layer = 0):
    curve_fun = lambda t: _racetrack_gradual_parametric(t, R = 5, N = 3)
    route_path = gdspy.Path(width = width, initial_point = [0,0])
    route_path.parametric(curve_fun, number_of_evaluations=99,\
            max_points=4000,  final_distance=None, layer=layer)
    D = Device('racetrack')
    D.add(route_path)
    return D


def _racetrack_gradual_parametric(t, R, N):
    """ Takes in a parametric value ``t`` on (0,1), returns the x,y coordinates
    of a racetrack bent according to 20090810_EOS4_modulator_designs_excerptForJasonGradualBends.ppt """
    x0 = R/2**(1/N)
    Rmin = 2**(0.5-1/N)/(N-1)*R
    R0 = R-(x0-Rmin/sqrt(2))
    t = np.array(t)
    x,y = np.zeros(t.shape), np.zeros(t.shape)

    # Doing the math
    x = cos(t*pi/2)*R0 # t (0-1) while x (0 to R0)
    ii =  (Rmin/sqrt(2) < x) & (x <= R0)
    jj =  (0 < x) & (x <= Rmin/sqrt(2))
    y[ii] = (R**N - (x[ii]+(x0-Rmin/sqrt(2)))**N)**(1/N)
    y[jj] = (x0-Rmin/sqrt(2))+sqrt(Rmin**2-x[jj]**2)
    return x,y



#==============================================================================
# Example code
#==============================================================================

# D = racetrack_gradual(width = 0.3, R = 5, N = 3)
# quickplot(D)



# t = np.linspace(0,1)
# x,y = _racetrack_gradual_parametric(t, R = 5, N = 3)
# plt.plot(x,y)
# plt.axis('equal')







#==============================================================================
#
# Packing and fill
#
#==============================================================================


def _pack_single_bin(
    rect_dict,
    aspect_ratio,
    max_size,
    sort_by_area,
    density,
    precision,
    verbose,
    ):
    """ Takes a `rect_dict` argument of the form {id:(w,h)} and tries to 
    pack it into a bin as small as possible with aspect ratio `aspect_ratio`
    Will iteratively grow the bin size until everything fits or the bin size 
    reaches `max_size`.  Returns a dictionary of of the packed rectangles
    in the form {id:(x,y,w,h)}, and a dictionary of remaining unpacked rects """

    try:
        import rectpack
    except:
        raise ImportError('[PHIDL] The packer() function requires the module "rectpack"' + 
              "to operate.  Please retry after installing rectpack:" +
              "\n\n$ pip install rectpack")

    # Compute total area and use it for an initial estimate of the bin size
    total_area = 0
    for r in rect_dict.values():
        total_area += r[0]*r[1]
    aspect_ratio = np.asarray(aspect_ratio)/np.linalg.norm(aspect_ratio) # Normalize
    
    # Setup variables
    box_size = np.asarray(aspect_ratio*np.sqrt(total_area), dtype = np.float64)
    box_size = np.clip(box_size, None, max_size)
    if sort_by_area == True: rp_sort = rectpack.SORT_AREA
    else:                    rp_sort = rectpack.SORT_NONE
    
    
    # Repeatedly run the rectangle-packing algorithm with increasingly larger
    # areas until everything fits or we've reached the maximum size
    while True:
        # Create the packer object
        rect_packer = rectpack.newPacker(
                                    mode = rectpack.PackingMode.Offline,
                                    pack_algo = rectpack.MaxRectsBlsf,
                                    sort_algo = rp_sort,
                                    bin_algo= rectpack.PackingBin.BBF,
                                    rotation=False,)
        
        # Add each rectangle to the packer, create a single bin, and pack
        for rid, r in rect_dict.items():
            rect_packer.add_rect(width = r[0], height = r[1], rid = rid)
        rect_packer.add_bin(width = box_size[0], height = box_size[1])
        rect_packer.pack()
        
        # Adjust the box size for next time
        box_size *= density  # Increase area to try to fit
        box_size = np.clip(box_size, None, max_size)
        if verbose == True:  print('Trying to pack in bin size (%0.2f, %0.2f)' % tuple(box_size*precision))
        
        # Quit the loop if we've packed all the rectangles or reached the max size
        if (len(rect_packer.rect_list()) == len(rect_dict)):
            if verbose == True: print('Success!')
            break
        elif all(box_size >= max_size):
            if verbose == True: print('Reached max_size, creating an additional bin')
            break

    
    # Separate packed from unpacked rectangles, make dicts of form {id:(x,y,w,h)}
    packed_rect_dict = {r[-1]:r[:-1] for r in rect_packer[0].rect_list()}
    unpacked_rect_dict = {}
    for k,v in rect_dict.items():
        if k not in packed_rect_dict:
            unpacked_rect_dict[k] = v
            
    return (packed_rect_dict, unpacked_rect_dict)



def packer(
        D_list,
        spacing = 10,
        aspect_ratio = (1,1),
        max_size = (None,None),
        sort_by_area = True,
        density = 1.1,
        precision = 1e-2,
        verbose = False,
        ):
    if density < 1.01:
        raise ValueError("[PHIDL] packer() was given a `density` argument that is" + 
              " too small.  The density argument must be >= 1.01")

    # Santize max_size variable
    max_size = [np.inf if v is None else v for v in max_size]
    max_size = np.asarray(max_size, dtype = np.float64) # In case it's integers
    max_size = max_size/precision
    
    # Convert Devices to rectangles
    rect_dict = {}
    for n, D in enumerate(D_list):
        w,h = (D.size + spacing)/precision
        w,h = int(w), int(h)
        if (w > max_size[0]) or (h > max_size[1]):
            raise ValueError("[PHIDL] packer() failed because one of the objects " + 
                  "in `D_list` is has an x or y dimension larger than `max_size` and " +
                  "so cannot be packed")
        rect_dict[n] = (w,h)
        
    packed_list = []
    while len(rect_dict) > 0:
        (packed_rect_dict, rect_dict) = _pack_single_bin(rect_dict,
                                                         aspect_ratio = aspect_ratio,
                                                         max_size = max_size,
                                                         sort_by_area = sort_by_area,
                                                         density = density,
                                                         precision = precision,
                                                         verbose = verbose,)
        packed_list.append(packed_rect_dict)
    
    D_packed_list = []
    for rect_dict in packed_list:
        D_packed = Device()
        for n, rect in rect_dict.items():
            x, y, w, h = rect
            xcenter = x + w/2 + spacing/2
            ycenter = y + h/2 + spacing/2
            d = D_packed.add_ref(D_list[n])
            d.center = (xcenter*precision,ycenter*precision)
        D_packed_list.append(D_packed)

    return D_packed_list



def _rasterize_polygons(polygons, bounds = [[-100, -100], [100, 100]], dx = 1, dy = 1):
    try:
        from skimage import draw
    except:
        raise ImportError("""
            The fill function requires the module "scikit-image"
            to operate.  Please retry after installing scikit-image:

            $ pip install --upgrade scikit-image """)


    # Prepare polygon array by shifting all points into the first quadrant and
    # separating points into x and y lists
    xpts = []
    ypts = []
    for p in polygons:
        p_array = np.asarray(p)
        x = p_array[:,0]
        y = p_array[:,1]
        xpts.append((x-bounds[0][0])/dx-0.5)
        ypts.append((y-bounds[0][1])/dy-0.5)

    # Initialize the raster matrix we'll be writing to
    xsize = int(np.ceil((bounds[1][0]-bounds[0][0]))/dx)
    ysize = int(np.ceil((bounds[1][1]-bounds[0][1]))/dy)
    raster = np.zeros((ysize, xsize), dtype=np.bool)

    # TODO: Replace polygon_perimeter with the supercover version
    for n in range(len(xpts)):
        rr, cc = draw.polygon(ypts[n], xpts[n], shape=raster.shape)
        rrp, ccp = draw.polygon_perimeter(ypts[n], xpts[n], shape=raster.shape, clip=False)
        raster[rr, cc] = 1
        raster[rrp, ccp] = 1

    return raster

def _raster_index_to_coords(i, j, bounds = [[-100, -100], [100, 100]], dx = 1, dy = 1):
    x = (j+0.5)*dx + bounds[0][0]
    y = (i+0.5)*dy + bounds[0][1]
    return x,y


def _expand_raster(raster, distance = (4,2)):
    try:
        from skimage import draw, morphology
    except:
        raise ImportError("""
            The fill function requires the module "scikit-image"
            to operate.  Please retry after installing scikit-image:

            $ pip install --upgrade scikit-image """)
    if distance[0] <= 0.5 and distance[1] <= 0.5: return raster

    num_pixels = np.array(np.ceil(distance), dtype = int)
    neighborhood = np.zeros((num_pixels[1]*2+1, num_pixels[0]*2+1), dtype=np.bool)
    rr, cc = draw.ellipse(num_pixels[1], num_pixels[0], distance[1]+0.5, distance[0]+0.5)
    neighborhood[rr, cc] = 1

    return morphology.binary_dilation(image = raster, selem=neighborhood)



def _fill_cell_rectangle(size = (20,20), layers = (0,1,3),
                         densities = (0.5, 0.25, 0.7), inverted = (False, False, False)):
    D = Device('fillcell')
    for layer, density, inv in zip(layers, densities, inverted):
        rectangle_size = np.array(size)*np.sqrt(density)
#        r = D.add_ref(rectangle(size = rectangle_size, layer = layer))
        R = rectangle(size = rectangle_size, layer = layer)
        R.center = (0,0)
        if inv is True:
            A = rectangle(size = size)
            A.center = (0,0)
            A = A.get_polygons()
            B = R.get_polygons()
            p = gdspy.boolean(A, B, operation = 'not')
            D.add_polygon(p, layer = layer)
        else:
            D.add_ref(R)
    return D

def _loop_over(var):
    # Checks if a variable is in the form of an iterable (list/tuple)
    # and if not, returns it as a list.  Useful for allowing argument
    # inputs to be either lists (e.g. [1,3,4]) or single-valued (e.g. 3)
    if hasattr(var,"__iter__"):
        return var
    else:
        return [var]

def fill_rectangle(D, fill_size = (40,10), avoid_layers = 'all', include_layers = None,
                    margin = 100, fill_layers = (0,1,3),
                   fill_densities = (0.5, 0.25, 0.7), fill_inverted = None, bbox = None):

    # Create the fill cell.  If fill_inverted is not specified, assume all False
    fill_layers = _loop_over(fill_layers)
    fill_densities = _loop_over(fill_densities)
    if fill_inverted is None: fill_inverted = [False]*len(fill_layers)
    fill_inverted = _loop_over(fill_inverted)
    if len(fill_layers) != len(fill_densities):
        raise ValueError("[PHIDL] phidl.geometry.fill_rectangle() `fill_layers` and" +
        " `fill_densities` parameters must be lists of the same length")
    if len(fill_layers) != len(fill_inverted):
        raise ValueError("[PHIDL] phidl.geometry.fill_rectangle() `fill_layers` and" +
        " `fill_inverted` parameters must be lists of the same length")

    fill_cell = _fill_cell_rectangle(size = fill_size, layers = fill_layers,
                                     densities = fill_densities, inverted = fill_inverted)
    F = Device(name = 'fill_pattern')

    if avoid_layers == 'all':
        exclude_polys = D.get_polygons(by_spec=False, depth=None)
    else:
        avoid_layers = [_parse_layer(l) for l in _loop_over(avoid_layers)]
        exclude_polys = D.get_polygons(by_spec=True, depth=None)
        exclude_polys = {key:exclude_polys[key] for key in exclude_polys if key in avoid_layers}
        exclude_polys = itertools.chain.from_iterable(exclude_polys.values())

    if include_layers is None:
        include_polys = []
    else:
        include_layers = [_parse_layer(l) for l in _loop_over(include_layers)]
        include_polys = D.get_polygons(by_spec=True, depth=None)
        include_polys = {key:include_polys[key] for key in include_polys if key in include_layers}
        include_polys = itertools.chain.from_iterable(include_polys.values())



    if bbox is None:  bbox = D.bbox

    raster = _rasterize_polygons(polygons = exclude_polys, bounds = bbox, dx = fill_size[0], dy = fill_size[1])
    raster = raster & ~_rasterize_polygons(polygons = include_polys, bounds = bbox, dx = fill_size[0], dy = fill_size[1])
    raster = _expand_raster(raster, distance = margin/np.array(fill_size))

    for i in range(np.size(raster,0)):
        sub_rasters = [list(g) for k, g in itertools.groupby(raster[i])]
        j = 0
        for s in sub_rasters:
            if s[0] == 0:
                x,y = _raster_index_to_coords(i, j, bbox, fill_size[0], fill_size[1])
                # F.add(gdspy.CellArray(ref_cell = fill_cell, columns = len(s), rows = 1, spacing = fill_size, ))
                a = F.add_array(fill_cell, columns = len(s), rows = 1, spacing = fill_size)
                a.move((x, y))
            j += len(s)

    return F





#==============================================================================
#
# Photonics
#
#==============================================================================


def polygon_ports(xpts=[-1,-1, 0, 0],
            ypts = [0, 1, 1, 0],
            layer = 0):
    # returns a polygon with ports on all edges
    P = Device('polygon')
    P.add_polygon([xpts, ypts], layer = layer)
    n = len(xpts)
    xpts.append(xpts[0])
    ypts.append(ypts[0])
    #determine if clockwise or counterclockwise
    cc = 0
    for i in range(0,n):
        cc += ((xpts[i+1]-xpts[i])*(ypts[i+1]+ypts[i]))

    for i in range(0,n):
        midpoint_n = [(xpts[i+1]+xpts[i])/2, (ypts[i+1]+ypts[i])/2]
        orientation_n = np.arctan2(np.sign(cc)*(xpts[i+1]-xpts[i]),np.sign(cc)*(ypts[i]-ypts[i+1]))*180/np.pi
        width_n = np.abs(np.sqrt((xpts[i+1]-xpts[i])**2+(ypts[i+1]-ypts[i])**2))
        P.add_port(name = str(i+1), midpoint = midpoint_n, width = width_n, orientation = orientation_n)

    return P

#==============================================================================
# Example code
#==============================================================================

# P = polygon(xpts=[-1,-3, 0, 0], ypts = [0, 1, 2, 0], layer = 3)
# quickplot(P)


@device_lru_cache
def grating(num_periods = 20, period = 0.75, fill_factor = 0.5, width_grating = 5, length_taper = 10, width = 0.4, partial_etch = False):
    #returns a fiber grating
    G = Device('grating')

# make the deep etched grating
    if partial_etch is False:
        # make the grating teeth
        for i in range(num_periods):
            cgrating = G.add_ref(compass(size=[period*fill_factor,width_grating], layer = 0))
            cgrating.x+=i*period

        # make the taper
        tgrating = G.add_ref(taper(length = length_taper, width1 = width_grating, width2 = width, port = None, layer = 0))
        tgrating.xmin = cgrating.xmax
        # define the port of the grating
        p = G.add_port(port = tgrating.ports[2], name = 1)
# make a partially etched grating
    if partial_etch is True:
        # hard coded overlap
            partetch_overhang = 5
            # make the etched areas (opposite to teeth)
            for i in range(num_periods):
                cgrating = G.add_ref(compass(size=[period*(1-fill_factor),width_grating+partetch_overhang*2]), layer = 1)
                cgrating.x+=i*period
                        # define the port of the grating
            p = G.add_port(port = cgrating.ports['E'], name = 1)
            p.midpoint=p.midpoint+np.array([(1-fill_factor)*period,0])

        #draw the deep etched square around the grating
            deepbox = G.add_ref(compass(size=[num_periods*period, width_grating]), layer=0)
    return G

#==============================================================================
# Example code
#==============================================================================

# G = grating(num_periods = 20, period = 0.75, fill_factor = 0.5, width_grating = 20, length_taper = 10, width = 0.4, partial_etch = False)
# quickplot(G)




#==============================================================================
#
# Test Structures
#
#==============================================================================


# Via Route ----------------------------------------
def _via_iterable(via_spacing, wire_width, wiring1_layer, wiring2_layer, via_layer, via_width):
    VI = Device('test_via_iter')
    wire1 = VI.add_ref(compass(size=(via_spacing, wire_width), layer=wiring1_layer))
    wire2 = VI.add_ref(compass(size=(via_spacing, wire_width), layer=wiring2_layer))
    via1 = VI.add_ref(compass(size=(via_width, via_width), layer=via_layer))
    via2 = VI.add_ref(compass(size=(via_width, via_width), layer=via_layer))
    wire1.connect(port='E', destination = wire2.ports['W'], overlap=wire_width)
    via1.connect(port='W', destination = wire1.ports['E'], overlap = (wire_width + via_width)/2)
    via2.connect(port='W', destination = wire2.ports['E'], overlap = (wire_width + via_width)/2)
    VI.add_port(name='W', port = wire1.ports['W'])
    VI.add_port(name='E', port = wire2.ports['E'])
    VI.add_port(name='S', midpoint = [(1*wire_width)+ wire_width/2,-wire_width/2], width = wire_width, orientation = -90)
    VI.add_port(name='N', midpoint = [(1*wire_width)+ wire_width/2,wire_width/2], width = wire_width, orientation = 90)

    return VI

def test_via(num_vias = 100, wire_width = 10, via_width = 15, via_spacing = 40, pad_size = (300,300), min_pad_spacing = 0,
                             pad_layer = 0, wiring1_layer = 1, wiring2_layer = 2, via_layer = 3):
    """
    Usage:
        Call via_route_test_structure() by indicating the number of vias you want drawn. You can also change the other parameters however
        if you do not specifiy a value for a parameter it will just use the default value
        Ex::

            via_route_test_structure(num_vias=54)

        - or -::

            via_route_test_structure(num_vias=12, pad_size=(100,100),wire_width=8)

        total requested vias (num_vias) -> this needs to be even
        pad size (pad_size) -> given in a pair (width, height)
        wire_width -> how wide each wire should be
        pad_layer -> GDS layer number of the pads
        wiring1_layer -> GDS layer number of the top wiring
        wiring2_layer -> GDS layer number of the bottom wiring
        via_layer -> GDS layer number of the vias
        ex: via_route(54, min_pad_spacing=300)
    """

    VR = Device('test_via')
    pad1 = VR.add_ref(rectangle(size=pad_size, layer=pad_layer))
    pad1_overlay = VR.add_ref(rectangle(size=pad_size, layer=wiring1_layer))
    pad2 = VR.add_ref(rectangle(size=pad_size, layer=pad_layer))
    pad2_overlay = VR.add_ref(rectangle(size=pad_size, layer=wiring1_layer))
    nub = VR.add_ref(compass(size=(3*wire_width,wire_width),layer=pad_layer))
    nub_overlay = VR.add_ref(compass(size=(3*wire_width,wire_width),layer=wiring1_layer))
    head = VR.add_ref(compass(size=(wire_width,wire_width),layer=pad_layer))
    head_overlay = VR.add_ref(compass(size=(wire_width,wire_width),layer=wiring1_layer))
    nub.ymax = pad1.ymax-5
    nub.xmin = pad1.xmax
    nub_overlay.ymax = pad1.ymax-5
    nub_overlay.xmin = pad1.xmax
    head.connect(port = "W", destination = nub.ports["E"])
    head_overlay.connect(port = "W", destination = nub_overlay.ports["E"])
    pad1_overlay.xmin = pad1.xmin
    pad1_overlay.ymin = pad1.ymin

    old_port = head.ports['S']
    count = 0
    width_via_iter = 2*via_spacing - 2*wire_width

    pad2.xmin = pad1.xmax + min_pad_spacing
    up = False
    down = True
    edge = True
    current_width = 3*wire_width + wire_width #width of nub and 1 overlap
    obj_old = head
    obj = head
    via_iterable = _via_iterable(via_spacing, wire_width, wiring1_layer, wiring2_layer, via_layer, via_width)
    while( (count+2) <= num_vias):
        obj = VR.add_ref(via_iterable)
        obj.connect(port = 'W', destination = old_port, overlap = wire_width)
        old_port = obj.ports['E']
        edge = False
        if(obj.ymax > pad1.ymax):
            obj.connect(port = 'W', destination = obj_old.ports['S'], overlap = wire_width)
            old_port = obj.ports['S']
            current_width += width_via_iter
            down = True
            up = False
            edge = True

        elif(obj.ymin < pad1.ymin):
            obj.connect(port = 'W', destination = obj_old.ports['N'], overlap = wire_width)
            old_port = obj.ports['N']
            current_width += width_via_iter
            up = True
            down = False
            edge = True
        count = count + 2
        obj_old = obj

    if(current_width < min_pad_spacing and (min_pad_spacing - current_width) > 3*wire_width):
        tail = VR.add_ref(compass(size=(min_pad_spacing-current_width+wire_width,wire_width),layer=wiring1_layer))
        tail_overlay = VR.add_ref(compass(size=(min_pad_spacing-current_width+wire_width,wire_width),layer=pad_layer))
    else:
        tail = VR.add_ref(compass(size=(3*wire_width,wire_width),layer=wiring1_layer))
        tail_overlay = VR.add_ref(compass(size=(3*wire_width,wire_width),layer=wiring1_layer))

    if(up == True and edge != True):
        tail.connect(port = 'W', destination = obj.ports['S'], overlap = wire_width)
        tail_overlay.connect(port = 'W', destination = obj.ports['S'], overlap = wire_width)
    elif(down == True and edge != True):
        tail.connect(port = 'W', destination = obj.ports['N'], overlap = wire_width)
        tail_overlay.connect(port = 'W', destination = obj.ports['N'], overlap = wire_width)
    else:
        tail.connect(port = 'W', destination = obj.ports['E'], overlap = wire_width)
        tail_overlay.connect(port = 'W', destination = obj.ports['E'], overlap = wire_width)


    pad2.xmin = tail.xmax
    pad2_overlay.xmin = pad2.xmin
    pad2_overlay.ymin = pad2.ymin

    return VR





def test_comb(pad_size = (200,200), wire_width = 1, wire_gap = 3,
              comb_layer = 0, overlap_zigzag_layer = 1,
              comb_pad_layer = None, comb_gnd_layer = None, overlap_pad_layer = None):
    """
    Usage:

    Call comb_insulation_test_structure() with any of the
    parameters shown below which you'd like to change. You
    only need to supply the parameters which you intend on
    changing You can alternatively call it with no parameters
    and it will take all the default alues shown below.
    Ex::

        comb_insulation_test_structure(pad_size=(175,175), wire_width=2, wire_gap=5)

    - or -::

        comb_insulation_test_structure()
    """
    CI = Device("test_comb")

    if comb_pad_layer is None:  comb_pad_layer = comb_layer
    if comb_gnd_layer is None:  comb_gnd_layer = comb_layer
    if overlap_pad_layer is None:  overlap_pad_layer = overlap_zigzag_layer
    wire_spacing = wire_width + wire_gap*2



    #%% pad overlays
    overlay_padb = CI.add_ref(rectangle(size=(pad_size[0]*9/10,pad_size[1]*9/10), layer=overlap_pad_layer))
    overlay_padl = CI.add_ref(rectangle(size=(pad_size[0]*9/10,pad_size[1]*9/10), layer=comb_pad_layer ) )
    overlay_padt = CI.add_ref(rectangle(size=(pad_size[0]*9/10,pad_size[1]*9/10), layer=comb_pad_layer ) )
    overlay_padr = CI.add_ref(rectangle(size=(pad_size[0]*9/10,pad_size[1]*9/10), layer=comb_gnd_layer))

    overlay_padl.xmin = 0
    overlay_padl.ymin = 0
    overlay_padb.ymax = 0
    overlay_padb.xmin = overlay_padl.xmax + pad_size[1]/5
    overlay_padr.ymin = overlay_padl.ymin
    overlay_padr.xmin = overlay_padb.xmax + pad_size[1]/5
    overlay_padt.xmin = overlay_padl.xmax + pad_size[1]/5
    overlay_padt.ymin = overlay_padl.ymax

    #%% pads
    padl = CI.add_ref(rectangle(size=pad_size, layer=comb_layer))
    padt = CI.add_ref(rectangle(size=pad_size, layer=comb_layer))
    padr = CI.add_ref(rectangle(size=pad_size, layer=comb_layer))
    padb = CI.add_ref(rectangle(size=pad_size, layer=overlap_zigzag_layer))
    padl_nub = CI.add_ref(rectangle(size=(pad_size[0]/4,pad_size[1]/2), layer=comb_layer))
    padr_nub = CI.add_ref(rectangle(size=(pad_size[0]/4,pad_size[1]/2), layer=comb_layer))

    padl.xmin = overlay_padl.xmin
    padl.center = [padl.center[0],overlay_padl.center[1]]
    padt.ymax = overlay_padt.ymax
    padt.center = [overlay_padt.center[0],padt.center[1]]
    padr.xmax = overlay_padr.xmax
    padr.center = [padr.center[0],overlay_padr.center[1]]
    padb.ymin = overlay_padb.ymin
    padb.center = [overlay_padb.center[0],padb.center[1]]
    padl_nub.xmin = padl.xmax
    padl_nub.center = [padl_nub.center[0],padl.center[1]]
    padr_nub.xmax = padr.xmin
    padr_nub.center = [padr_nub.center[0],padr.center[1]]

    #%% connected zig

    head = CI.add_ref(compass(size=(pad_size[0]/12, wire_width), layer=comb_layer))
    head.xmin = padl_nub.xmax
    head.ymax = padl_nub.ymax
    connector = CI.add_ref(compass(size=(wire_width, wire_width), layer=comb_layer))
    connector.connect(port = 'W', destination=head.ports['E'])
    old_port = connector.ports['S']
    top = True
    obj = connector
    while(obj.xmax + pad_size[0]/12 < padr_nub.xmin):
        #long zig zag rectangle
        obj = CI.add_ref(compass(size=(pad_size[1]/2 - 2*wire_width, wire_width), layer=comb_layer))
        obj.connect(port = 'W', destination=old_port)
        old_port = obj.ports['E']
        if(top):
            #zig zag edge rectangle
            obj = CI.add_ref(compass(size=(wire_width, wire_width), layer=comb_layer))
            obj.connect(port = 'N', destination=old_port)
            top = False
        else:
            #zig zag edge rectangle
            obj = CI.add_ref(compass(size=(wire_width, wire_width), layer=comb_layer))
            obj.connect(port = 'S', destination=old_port)
            top = True
            # comb rectange
            comb = CI.add_ref(rectangle(size=((padt.ymin-head.ymax)+pad_size[1]/2  - (wire_spacing + wire_width)/2, wire_width), layer=comb_layer))
            comb.rotate(90)
            comb.ymax = padt.ymin
            comb.xmax = obj.xmax - (wire_spacing+wire_width)/2
        old_port = obj.ports['E']
        obj = CI.add_ref(compass(size=(wire_spacing, wire_width), layer=comb_layer))
        obj.connect(port = 'W', destination=old_port)
        old_port = obj.ports['E']
        obj = CI.add_ref(compass(size=(wire_width, wire_width), layer=comb_layer))
        obj.connect(port = 'W', destination=old_port)
        if(top):
            old_port = obj.ports['S']
        else:
            old_port = obj.ports['N']
    old_port = obj.ports['E']
    if(padr_nub.xmin-obj.xmax > 0):
        tail = CI.add_ref(compass(size=(padr_nub.xmin-obj.xmax, wire_width), layer=comb_layer))
    else:
        tail = CI.add_ref(compass(size=(wire_width, wire_width), layer=comb_layer))
    tail.connect(port = 'W', destination=old_port)

    #%% disconnected zig

    dhead = CI.add_ref(compass(size=(padr_nub.ymin -padb.ymax - wire_width, wire_width), layer=overlap_zigzag_layer))
    dhead.rotate(90)
    dhead.ymin = padb.ymax
    dhead.xmax = tail.xmin-(wire_spacing+wire_width)/2
    connector = CI.add_ref(compass(size=(wire_width, wire_width), layer=overlap_zigzag_layer))
    connector.connect(port = 'S', destination=dhead.ports['E'])
    old_port = connector.ports['N']
    right = True
    obj = connector
    while(obj.ymax + wire_spacing + wire_width < head.ymax):
        obj = CI.add_ref(compass(size=(wire_spacing, wire_width), layer=overlap_zigzag_layer))
        obj.connect(port = 'W', destination=old_port)
        old_port = obj.ports['E']
        if(right):
            obj = CI.add_ref(compass(size=(wire_width, wire_width), layer=overlap_zigzag_layer))
            obj.connect(port = 'W', destination=old_port)
            right = False
        else:
            obj = CI.add_ref(compass(size=(wire_width, wire_width), layer=overlap_zigzag_layer))
            obj.connect(port = 'E', destination=old_port)
            right = True
        old_port = obj.ports['N']
        obj = CI.add_ref(compass(size=(dhead.xmin-(head.xmax+head.xmin+wire_width)/2, wire_width), layer=overlap_zigzag_layer))
        obj.connect(port = 'E', destination=old_port)
        old_port = obj.ports['W']
        obj = CI.add_ref(compass(size=(wire_width, wire_width), layer=overlap_zigzag_layer))
        obj.connect(port = 'S', destination=old_port)
        if(right):
            old_port = obj.ports['W']
        else:
            old_port = obj.ports['E']

    return CI




#This is a helper function to make the Ic step wire structure
def _test_ic_wire_step(thick_width = 10, thin_width = 1, wire_layer = 2):
    WS4 = Device('test_ic_step')
    wire_stepa = WS4.add_ref(optimal_step(thick_width/2, thin_width/2, layer=wire_layer))
    wire_stepb = WS4.add_ref(optimal_step(thin_width/2, thick_width/2, layer=wire_layer))
    wire_stepc = WS4.add_ref(optimal_step(thick_width/2, thin_width/2, layer=wire_layer))
    wire_stepd = WS4.add_ref(optimal_step(thin_width/2, thick_width/2, layer=wire_layer))
    wire_stepb.rotate(180)
    wire_stepb.xmin = wire_stepa.xmin
    wire_stepc.rotate(180)
    wire_stepc.xmin = wire_stepa.xmax
    wire_stepd.xmin = wire_stepc.xmin
    return WS4


def test_ic(wire_widths = [0.25, 0.5,1,2,4], wire_widths_wide = [0.75, 1.5, 3, 4, 6], pad_size = (200,200), pad_gap=75,
            wire_layer = 0, pad_layer = 1, gnd_layer = None):
    """
    Usage:

    Call ic_test_structure() with either a list of widths for the thickest part of each wire to test and a list for the
    thinnest parts of each wire. Alternatively, specify a list of widths for the thinnest part of each wire and ignore the
    wire_widths parameter. Instead you should specify the width_growth_factor which indicates by what factor the thick
    part of the wire will be larger than the thin part.
    Ex::

        ic_test_structure(wire_widths = [5,10,10,10,10], thin_width=[0.5,1,2,3,4])

    - or -::

        ic_test_structure(width_growth_factor = 5, thin_width=[0.5,1,2,3,4])
    """
    ICS = Device('test_ic')

    if gnd_layer is None: gnd_layer = pad_layer
    translation = 0
    padb = ICS.add_ref(rectangle(size=(np.size(wire_widths) * (pad_size[0]*6/5), pad_size[1]), layer=wire_layer))
    padb_overlay = ICS.add_ref(rectangle(size=((np.size(wire_widths) * (pad_size[0]*6/5))*9/10, pad_size[1]*9/10), layer=gnd_layer))
    padb_overlay.center = padb.center
    padb_overlay.ymin = padb.ymin
    for i, x in enumerate(wire_widths_wide):
        padt = ICS.add_ref(rectangle(pad_size, wire_layer))
        padt.xmin = padb.xmin + translation
        padt.ymin = padb.ymax + pad_gap
        padt_overlay = ICS.add_ref(rectangle(size=(pad_size[0]*9/10, pad_size[1]*9/10), layer=pad_layer))
        padt_overlay.center = padt.center
        padt_overlay.ymax = padt.ymax
        difference = padt.ymin-padb.ymax
        wire_step = ICS.add_ref(_test_ic_wire_step(wire_widths_wide[i], wire_widths[i], wire_layer=wire_layer))
        wire_step.rotate(90)
        wire_step.center = (padt.center[0], padb.ymax + difference/2)
        translation = translation + pad_size[0]*12/10
        conn_wire_top = ICS.add_ref(rectangle(size=(wire_widths_wide[i], padt.ymin-wire_step.ymax), layer=wire_layer))
        conn_wire_bottom = ICS.add_ref(rectangle(size=(wire_widths_wide[i], wire_step.ymin-padb.ymax), layer=wire_layer))
        conn_wire_top.ymax = padt.ymin
        conn_wire_top.xmin = wire_step.xmin
        conn_wire_bottom.ymin = padb.ymax
        conn_wire_bottom.xmin = wire_step.xmin
    return ICS

def test_res(pad_size = [50,50],
                     num_squares = 1000,
                     width = 1,
                     res_layer = 0,
                     pad_layer = None,
                     gnd_layer = None):

    """ Creates an efficient resonator structure for a wafer layout.

    Keyword arguments:
    pad_size    -- Size of the two matched impedance pads (microns)
    num_squares -- Number of squares comprising the resonator wire
    width       -- The width of the squares (microns)
    """

    x = pad_size[0]
    z = pad_size[1]

    # Checking validity of input
    if x <= 0 or z <= 0:
        raise ValueError('Pad must have positive, real dimensions')
    elif width > z:
        raise ValueError('Width of cell cannot be greater than height of pad')
    elif num_squares <= 0:
        raise ValueError('Number of squares must be a positive real number')
    elif width <= 0:
        raise ValueError('Width of cell must be a positive real number')

    # Performing preliminary calculations
    num_rows = int(np.floor(z / (2 * width)))
    if num_rows % 2 == 0:
        num_rows -= 1
    num_columns = num_rows - 1
    squares_in_row = (num_squares - num_columns - 2) / num_rows

    # Compensating for weird edge cases
    if squares_in_row < 1:
        num_rows = int(round(num_rows / 2) - 2)
        squares_in_row = 1
    if width * 2 > z:
        num_rows = 1
        squares_in_row = num_squares - 2

    length_row = squares_in_row * width

    # Creating row/column corner combination structure
    T = Device()
    Row = rectangle(size = (length_row, width), layer = res_layer)
    Col = rectangle(size = (width, width), layer = res_layer)

    row = T.add_ref(Row)
    col = T.add_ref(Col)
    col.move([length_row - width, -width])

    # Creating entire waveguide net
    N = Device('net')
    n = 1
    for i in range(num_rows):
        if i != num_rows - 1:
            d = N.add_ref(T)
        else:
            d = N.add_ref(Row)
        if n % 2 == 0:
            d.mirror(p1 = (d.x, d.ymax), p2 = (d.x, d.ymin))
        d.movey(-(n - 1) * T.ysize)
        n += 1
    d = N.add_ref(Col).movex(-width)
    d = N.add_ref(Col).move([length_row, -(n - 2) * T.ysize])

    # Creating pads
    P = Device('pads')
    Pad1 = rectangle(size = (x,z), layer = pad_layer)
    Pad2 = rectangle(size = (x + 5, z), layer = pad_layer)
    Gnd1 = offset(Pad1, distance = -5, layer = gnd_layer)
    Gnd2 = offset(Pad2, distance = -5, layer = gnd_layer)
    pad1 = P.add_ref(Pad1).movex(-x - width)
    pad2 = P.add_ref(Pad1).movex(length_row + width)
    gnd1 = P.add_ref(Gnd1).center = pad1.center
    gnd2 = P.add_ref(Gnd2)
    nets = P.add_ref(N).y = pad1.y
    gnd2.center = pad2.center
    gnd2.movex(2.5)

    return P


#==============================================================================
#
# Optimal current-crowding superconducting structures
#
#==============================================================================


@device_lru_cache
def optimal_hairpin(width = 0.2, pitch = 0.6, length = 10,
    turn_ratio = 4, num_pts = 50, layer = 0):
    # Optimal structure from https://doi.org/10.1103/PhysRevB.84.174510
    # Clem, J., & Berggren, K. (2011). Geometry-dependent critical currents in 
    # superconducting nanocircuits. Physical Review B, 84(17), 1–27. 
    #==========================================================================
    #  Create the basic geometry
    #==========================================================================
    a = (pitch + width)/2
    y = -(pitch - width)/2
    x = -pitch
    dl = width/(num_pts*2)
    n = 0

    # Get points of ideal curve from conformal mapping
    # TODO This is an inefficient way of finding points that you need
    xpts = [x]; ypts = [y]
    while (y < 0) & (n<1e6):
        s = x + 1j*y
        w = np.sqrt(1 - np.exp(pi*s/a))
        wx = np.real(w); wy = np.imag(w)
        wx = wx/np.sqrt(wx**2+wy**2); wy = wy/np.sqrt(wx**2+wy**2)
        x = x + wx*dl; y = y + wy*dl
        xpts.append(x); ypts.append(y)
        n = n+1
    ypts[-1] = 0 # Set last point be on the x=0 axis for sake of cleanliness
    ds_factor = int(len(xpts)/num_pts) # Downsample the total number of points
    xpts = xpts[::-ds_factor]; xpts = xpts[::-1]    # This looks confusing, but it's just flipping the arrays around
    ypts = ypts[::-ds_factor]; ypts = ypts[::-1]    # so the last point is guaranteed to be included when downsampled

    # Add points for the rest of meander
    xpts.append(xpts[-1] + turn_ratio*width); ypts.append(0)
    xpts.append(xpts[-1]); ypts.append(-a)
    xpts.append(xpts[0]); ypts.append(-a)
    xpts.append(max(xpts)-length); ypts.append(-a)
    xpts.append(xpts[-1]); ypts.append(-a + width)
    xpts.append(xpts[0]); ypts.append(ypts[0])

    xpts = np.array(xpts)
    ypts = np.array(ypts)

    #==========================================================================
    #  Create a blank device, add the geometry, and define the ports
    #==========================================================================
    D = Device(name = 'hairpin')
    D.add_polygon([xpts,ypts], layer = layer)
    D.add_polygon([xpts,-ypts], layer = layer)

    xports = min(xpts)
    yports = -a + width/2
    D.add_port(name = 1, midpoint = [xports,-yports], width = width, orientation = 180)
    D.add_port(name = 2, midpoint = [xports,yports], width = width, orientation = 180)

    return D


# TODO Include parameter which specifies "half" (one edge flat) vs "full" (both edges curved)
@device_lru_cache
def optimal_step(start_width = 10, end_width = 22, num_pts = 50, width_tol = 1e-3,
                 anticrowding_factor = 1.2, symmetric = False, layer = 0):
    # Optimal structure from https://doi.org/10.1103/PhysRevB.84.174510
    # Clem, J., & Berggren, K. (2011). Geometry-dependent critical currents in 
    # superconducting nanocircuits. Physical Review B, 84(17), 1–27. 
    #==========================================================================
    #  Create the basic geometry
    #==========================================================================
    def step_points(eta, W, a):
        # Returns points from a unit semicircle in the w (= u + iv) plane to
        # the optimal curve in the zeta (= x + iy) plane which transitions
        # a wire from a width of 'W' to a width of 'a'
        # eta takes value 0 to pi

        W = np.complex(W)
        a = np.complex(a)

        gamma = (a*a + W*W)/(a*a - W*W)

        w = np.exp(1j*eta)

        zeta = 4*1j/pi*(W*np.arctan(np.sqrt((w-gamma)/(gamma+1))) \
                           + a*np.arctan(np.sqrt((gamma-1)/(w-gamma))))

        x = np.real(zeta)
        y = np.imag(zeta)
        return x,y


    def invert_step_point(x_desired = -10, y_desired = None, W = 1, a = 2):
        # Finds the eta associated with the value x_desired along the optimal curve
        def fh(eta):
            guessed_x, guessed_y = step_points(eta, W = W, a = a)
            if y_desired is None:   return (guessed_x-x_desired)**2 # The error
            else:                   return (guessed_y-y_desired)**2

        try:
            from scipy.optimize import fminbound
        except:
            raise ImportError(""" [PHIDL] To run the optimal-curve geometry functions you need scipy,
              please install it with `pip install scipy` """)
        found_eta = fminbound(fh, x1 = 0, x2 = pi, args=())
        return step_points(found_eta, W = W, a = a)

    if start_width > end_width:
        reverse = True
        start_width, end_width = end_width, start_width
    else:
        reverse = False

    if start_width == end_width: # Just return a square
        ypts = [0, start_width, start_width,           0]
        xpts = [0,           0, start_width, start_width]
    else:
        xmin,ymin = invert_step_point(y_desired = start_width*(1+width_tol), W = start_width, a = end_width)
        xmax,ymax = invert_step_point(y_desired = end_width*(1-width_tol), W = start_width, a = end_width)

        xpts = np.linspace(xmin, xmax, num_pts).tolist()
        ypts = []
        for x in xpts:
            x,y = invert_step_point(x_desired = x, W = start_width, a = end_width)
            ypts.append(y)

        ypts[-1] = end_width
        ypts[0] =  start_width
        if symmetric == False:
            xpts.append(xpts[-1])
            ypts.append(0)
            xpts.append(xpts[0])
            ypts.append(0)
        else:
            xpts += [x for x in xpts[::-1]]
            ypts += [-y for y in ypts[::-1]]
            xpts = [x/2 for x in xpts]
            ypts = [y/2 for y in ypts]

        # anticrowding_factor stretches the wire out; a stretched wire is a gentler
        # transition, so there's less chance of current crowding if the fabrication
        # isn't perfect but as a result, the wire isn't as short as it could be
        xpts = (np.array(xpts)*anticrowding_factor).tolist()

        if reverse is True:
            xpts = (-np.array(xpts)).tolist()
            start_width, end_width = end_width, start_width

    #==========================================================================
    #  Create a blank device, add the geometry, and define the ports
    #==========================================================================
    D = Device(name = 'step')
    D.add_polygon([xpts,ypts], layer = layer)

    if symmetric == False:
        D.add_port(name = 1, midpoint = [min(xpts),start_width/2], width = start_width, orientation = 180)
        D.add_port(name = 2, midpoint = [max(xpts),end_width/2], width = end_width, orientation = 0)
    if symmetric == True:
        D.add_port(name = 1, midpoint = [min(xpts),0], width = start_width, orientation = 180)
        D.add_port(name = 2, midpoint = [max(xpts),0], width = end_width, orientation = 0)

    return D


def optimal_90deg(width = 100.0, num_pts = 15, length_adjust = 1, layer = 0):
    # Optimal structure from https://doi.org/10.1103/PhysRevB.84.174510
    # Clem, J., & Berggren, K. (2011). Geometry-dependent critical currents in 
    # superconducting nanocircuits. Physical Review B, 84(17), 1–27. 
    D = Device()

    # Get points of ideal curve
    a = 2*width
    v = np.logspace(-length_adjust,length_adjust,num_pts)
    xi = a/2.0*((1+2/np.pi*np.arcsinh(1/v)) + 1j*(1+2/np.pi*np.arcsinh(v)))
    xpts = list(np.real(xi)); ypts = list(np.imag(xi))

    # Add points for the rest of curve
    d = 2*xpts[0] # Farthest point out * 2, rounded to nearest 100
    xpts.append(width); ypts.append(d)
    xpts.append(0); ypts.append(d)
    xpts.append(0); ypts.append(0)
    xpts.append(d); ypts.append(0)
    xpts.append(d); ypts.append(width)
    xpts.append(xpts[0]); ypts.append(ypts[0])

    D.add_polygon([xpts, ypts], layer = layer)

    D.add_port(name = 1, midpoint = [a/4,d], width = a/2, orientation = 90)
    D.add_port(name = 2, midpoint = [d,a/4], width = a/2, orientation = 0)
    return D



#==============================================================================
# Example code
#==============================================================================

#hairpin = optimal_hairpin(width = 1, pitch = 3, length = 30, num_pts = 20)
#quickplot(hairpin)


#step = optimal_step(start_width = 5, end_width = 1, num_pts = 80, width_tol = 1e-3)
#quickplot(step)


#turn = optimal_90deg(width = 90, length_adjust = 1)
#quickplot(turn)





#==============================================================================
#
# Superconducting devices
#
#==============================================================================


@device_lru_cache
def snspd(wire_width = 0.2, wire_pitch = 0.6, size = (10,8),
        num_squares = None, turn_ratio = 4,
        terminals_same_side = False, layer = 0):
    # Convenience tests to auto-shape the size based
    # on the number of squares
    if num_squares is not None and ((size is None) or ((size[0] is None) and (size[1]) is None)):
        xy = np.sqrt(num_squares*wire_pitch*wire_width)
        size = [xy,xy]
        num_squares = None
    if ([size[0], size[1], num_squares].count(None) != 1):
        raise ValueError('[PHIDL] snspd() requires that exactly ONE value of' +
                         ' the arguments ``num_squares`` and ``size`` be None'+
                         ' to prevent overconstraining, for example:\n' +
                         '>>> snspd(size = (3, None), num_squares = 2000)')
    if size[0] is None:
        ysize = size[1]
        xsize = num_squares*wire_pitch*wire_width/ysize
    elif size[1] is None:
        xsize = size[0]
        ysize = num_squares*wire_pitch*wire_width/xsize
    else:
        xsize = size[0]
        ysize = size[1]

    num_meanders = int(np.ceil(ysize/wire_pitch))

    D = Device(name = 'snspd')
    hairpin = optimal_hairpin(width = wire_width, pitch = wire_pitch,
        turn_ratio = turn_ratio, length = xsize/2, num_pts = 20, layer = layer)


    if (terminals_same_side is False) and ((num_meanders % 2) == 0):
        num_meanders += 1
    elif (terminals_same_side is True) and ((num_meanders % 2) == 1):
        num_meanders += 1

    start_nw = D.add_ref(compass(size = [xsize/2 ,wire_width], layer = layer))

    hp_prev = D.add_ref(hairpin)
    hp_prev.connect(1, start_nw.ports['E'])
    alternate = True
    for n in range(2,num_meanders):
        hp = D.add_ref(hairpin)
        if alternate:
            hp.connect(2, hp_prev.ports[2])
            last_port = hp.ports[1]
        else:
            hp.connect(1, hp_prev.ports[1])
            last_port = hp.ports[2]
        hp_prev = hp
        alternate = not alternate

    finish_se = D.add_ref(compass(size = [xsize/2 ,wire_width], layer = layer))
    finish_se.connect('E', last_port)

    D.add_port(port = start_nw.ports['W'], name = 1)
    D.add_port(port = finish_se.ports['W'], name = 2)

    D.info['num_squares'] = num_meanders*(xsize/wire_width)
    D.info['area'] = xsize*ysize
    D.info['size'] = (xsize, ysize)

    return D


def snspd_expanded(wire_width = 0.2, wire_pitch = 0.6, size = (10,8),
           num_squares = None, connector_width = 1, connector_symmetric = False,
            turn_ratio = 4, terminals_same_side = False, layer = 0):
    """ Creates an optimally-rounded SNSPD with wires coming out of it that expand"""
    D = Device('snspd_expanded')
    S = snspd(wire_width = wire_width, wire_pitch = wire_pitch,
                        size = size, num_squares = num_squares, turn_ratio = turn_ratio,
                        terminals_same_side = terminals_same_side, layer = layer)
    s = D.add_ref(S)
    step_device = optimal_step(start_width = wire_width, end_width = connector_width,
                            num_pts = 100, anticrowding_factor = 2, width_tol = 1e-3,
                            symmetric = connector_symmetric, layer = layer)
    step1 = D.add_ref(step_device)
    step2 = D.add_ref(step_device)
    step1.connect(port = 1, destination = s.ports[1])
    step2.connect(port = 1, destination = s.ports[2])
    D.add_port(name = 1, port = step1.ports[2])
    D.add_port(name = 2, port = step2.ports[2])

    D.info = S.info
    S.info = {}

    return D


#==============================================================================
# Example code
#==============================================================================

#s = snspd(wire_width = 0.2, wire_pitch = 0.6, size = [10,3], terminals_same_side = True)
#quickplot(s)

#s = snspd(wire_width = 0.2, wire_pitch = 0.6, size = [10, None],
#          num_squares = 1000, terminals_same_side = True)
#quickplot(s)

def ytron_round(rho = 1, arm_lengths = (500,300),  source_length = 500,
                arm_widths = (200, 200), theta = 2.5, theta_resolution = 10,
                layer = 0):

    #==========================================================================
    #  Create the basic geometry
    #==========================================================================
    theta = theta*pi/180
    theta_resolution = theta_resolution*pi/180
    thetalist = np.linspace(-(pi-theta),-theta, int((pi-2*theta)/theta_resolution) + 2)
    semicircle_x = rho*cos(thetalist)
    semicircle_y = rho*sin(thetalist)+rho

    # Rest of yTron
    xc = rho*cos(theta)
    yc = rho*sin(theta)
    arm_x_left  = arm_lengths[0]*sin(theta)
    arm_y_left  = arm_lengths[0]*cos(theta)
    arm_x_right = arm_lengths[1]*sin(theta)
    arm_y_right = arm_lengths[1]*cos(theta)

    # Write out x and y coords for yTron polygon
    xpts = semicircle_x.tolist() + [xc+arm_x_right, xc+arm_x_right+arm_widths[1], xc+arm_widths[1], \
           xc+arm_widths[1], 0, -(xc+arm_widths[0]), -(xc+arm_widths[0]), -(xc+arm_x_left+arm_widths[0]), -(xc+arm_x_left)]
    ypts = semicircle_y.tolist() + [yc+arm_y_right,      yc+arm_y_right,      yc,   yc-source_length, yc-source_length,  \
            yc-source_length,        yc,        yc+arm_y_left,    yc+arm_y_left]

    #==========================================================================
    #  Create a blank device, add the geometry, and define the ports
    #==========================================================================
    D = Device(name = 'ytron')
    D.add_polygon([xpts,ypts], layer = layer)
    D.add_port(name = 'left', midpoint = [-(xc+arm_x_left+arm_widths[0]/2), yc+arm_y_left],  width = arm_widths[0], orientation = 90)
    D.add_port(name = 'right', midpoint = [xc+arm_x_right+arm_widths[1]/2, yc+arm_y_right],  width = arm_widths[1], orientation = 90)
    D.add_port(name = 'source', midpoint = [0+(arm_widths[1]-arm_widths[0])/2, -source_length+yc],  width = arm_widths[0] + arm_widths[1] + 2*xc, orientation = -90)

    #==========================================================================
    #  Record any parameters you may want to access later
    #==========================================================================
    D.info['rho'] = rho
    D.info['left_width'] =   arm_widths[0]
    D.info['right_width'] =  arm_widths[1]
    D.info['source_width'] = arm_widths[0] + arm_widths[1] + 2*xc

    return D


#==============================================================================
# Example code
#==============================================================================

#y = ytron_round(rho = 1, arm_lengths = (500,300),  source_length = 500,
                # arm_widths = (200, 200), theta = 2.5, theta_resolution = 10,
                # layer = 0)
#quickplot(y)
