# -*- 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)
