Files

455 lines
19 KiB
Python

"""Curve and spiral geometry primitives."""
from typing import Any, Optional
from cmath import pi
import nazca as nd
import numpy as np
from .polygons import _my_polygon
class Conchoid:
def __init__(self,R0: Any,kR: Any,T: Any,w: float,layer: str,w_end: Optional[float]=None,res: float=0.1,final_flat: Any=None,begin_flat: Any=None,xs: Optional[str]=None) -> None:
## with half circle to be one cycle
if (w_end==None):
w_end = w
with nd.Cell(instantiate=False)as C:
n_sects = int(np.floor(T/np.pi)) ## intersecting into different semi-circle
L = R0*T + 1/2*kR*T**2 ## The total length of the Conchoid center line
n_points = int(np.floor(L/res))+1
## calculating sections
if (np.abs(T-n_sects*np.pi) <0.0001):
n_sects = n_sects
else:
n_sects = n_sects+1
# res_sect = int(np.floor(res/n_sects))+1
# res = 0
if (layer!=None):
nd.add_xsection(name='temp')
nd.add_layer2xsection(xsection='temp',layer=layer,growx=0,growy=0)
xs = 'temp'
for _n_ in range(0,n_sects):
# phi_start = _n_*pi
# phi_end = min(T,_n_*pi+pi) ## forward placement
phi_end = T - _n_*pi
phi_start = max(0,T - _n_*pi - pi)
L_sect = R0*(phi_end-phi_start) + 1/2*kR*(phi_end**2 - phi_start**2)
n_points = int(L_sect/res)+1
if (layer!=None):
Theta = np.linspace(phi_start, phi_end ,n_points)
R = (Theta*kR+R0) ## conchoid function
# res = kR/2*T*T + R0*T
# res = res + np.sum(R[0:-1]*np.diff(Theta))
w_cur = w
if (_n_==0):
w_cur = np.linspace(w,w_end,n_points)
# print("Loading Taper area")
vtx_cx = R*np.cos(Theta)
vtx_cy = R*np.sin(Theta)
vtx_center = np.c_[vtx_cx,vtx_cy]
e_theta = -1/((R0/kR)+Theta) ## actuall norm towards spiral
# e_rou = np.ones(len(e_theta))
ey = np.sin(Theta) - np.cos(Theta)*kR/R
ex = np.cos(Theta) + np.sin(Theta)*kR/R
if (final_flat!=None and _n_==0):
ey[-1] = np.sin(final_flat/180*np.pi)
ex[-1] = np.cos(final_flat/180*np.pi)
if (begin_flat!=None and _n_==n_sects-1):
ey[0] = np.sin(begin_flat/180*np.pi)
ex[0] = np.cos(begin_flat/180*np.pi)
# if (final_flat!=None and _n_==0):
# e_theta[-1] = final_flat
# if (begin_flat!=None and _n_==n_sects-1):
# e_theta[0] = begin_flat
if (_n_==0):
self.Atilt = np.arctan(ey[0]/ex[0])/np.pi*180
# print("Atilt_conchoid = %.5f" % self.Atilt)
Lnorm = np.sqrt(np.power(ex,2)+np.power(ey,2))
vtx_x = R*np.cos(Theta)
vtx_y = R*np.sin(Theta)
vtx_out_x = vtx_x + w_cur/2*ex/Lnorm
vtx_out_y = vtx_y + w_cur/2*ey/Lnorm
vtx_in_x = vtx_x - w_cur/2*ex/Lnorm
vtx_in_y = vtx_y - w_cur/2*ey/Lnorm
vtx_in = np.c_[np.flip(vtx_in_x),np.flip(vtx_in_y)]
vtx_out = np.c_[vtx_out_x,vtx_out_y]
vtx = np.r_[vtx_out,vtx_in]
_my_polygon(layer_wg=layer,vtx=vtx).put(0,0,0)
elif(layer==None and xs!=None):
for layers,growx,growy,acc in nd.layeriter(xs=xs):
(a1,b1), (a2,b2),c1,c2 = growx
Theta = np.linspace(phi_start, phi_end ,n_points)
R = (Theta*kR+R0)
# res = kR/2*T*T + R0*T
# res = res + np.sum(R[0:-1]*np.diff(Theta)
w_cur = w*(a1-a2) + (b1-b2)
if (_n_==0):
w_cur = np.linspace(w,w_end,n_points)
w_cur = w_cur*(a1-a2) + (b1-b2)
vtx_cx = R*np.cos(Theta)
vtx_cy = R*np.sin(Theta)
vtx_center = np.c_[vtx_cx,vtx_cy]
# e_theta = -1/((R0/kR)+Theta)
# e_rou = np.ones(len(e_theta))
# if (final_flat!=None and _n_==0):
# e_theta[-1] = final_flat
# if (begin_flat!=None and _n_==n_sects-1):
# e_theta[0] = begin_flat
ey = np.sin(Theta)*R - np.sin(Theta)*kR
ex = np.cos(Theta)*R + np.cos(Theta)*kR
if (final_flat!=None and _n_==0):
ey[-1] = np.sin(final_flat/180*np.pi)
ex[-1] = np.cos(final_flat/180*np.pi)
if (begin_flat!=None and _n_==n_sects-1):
ey[0] = np.sin(begin_flat/180*np.pi)
ex[0] = np.cos(begin_flat/180*np.pi)
Lnorm = np.sqrt(np.power(ex,2)+np.power(ey,2))
vtx_x = R*np.cos(Theta)
vtx_y = R*np.sin(Theta)
vtx_out_x = vtx_x + w_cur/2*ex/Lnorm
vtx_out_y = vtx_y + w_cur/2*ey/Lnorm
vtx_in_x = vtx_x - w_cur/2*ex/Lnorm
vtx_in_y = vtx_y - w_cur/2*ey/Lnorm
vtx_in = np.c_[np.flip(vtx_in_x),np.flip(vtx_in_y)]
vtx_out = np.c_[vtx_out_x,vtx_out_y]
vtx = np.r_[vtx_out,vtx_in]
_my_polygon(layer_wg=layers,vtx=vtx).put(0,0,0)
Rmax = T*kR+R0
## revised in 2026.06.07 by Qin Yue
# legacy: nd.Pin(name="a1").put(R0,0,-90)
nd.Pin(name="opt_a1",type="optical:").put(R0,0,-90)
## revised in 2026.06.07 by Qin Yue
# legacy: nd.Pin(name="b1").put(Rmax*np.cos(T),Rmax*np.sin(T),(T/np.pi*180+90))
nd.Pin(name="opt_b1",type="optical:").put(Rmax*np.cos(T),Rmax*np.sin(T),(T/np.pi*180+90))
self.L = L
self.cell =C
self.vtx_center = vtx_center
self.vtx = vtx
self.K_end = (np.power(np.max(R),2) + 2*np.power(kR,2)) / np.power((np.power(np.max(R),2) + np.power(kR,2)),1.5)
self.R_end = 1/self.K_end
def _line2wg_(x,y,wu,wd,theta,n_points):
""" building waveguide with center line and side expansion
Args:
vtx_line (list[float]): the location of points, [x,y]
width (list[float]): the expansion width of points vertical to the pointing vector, [wu,wd]
theta (list[float]): the pointing angle, [theta], 0 represent right, 180 represent left
"""
theta = theta*np.pi/180
x_u = x+wu*np.cos(theta+pi/2)
x_d = x+wd*np.cos(theta-pi/2)
y_u = y+wu*np.sin(theta+pi/2)
y_d = y+wd*np.sin(theta-pi/2)
### polygon section, reducing resolution
sect = np.linspace(start= 0,stop= len(x_u)-1,num= n_points)
sect = np.asarray(sect, dtype = int)
x_u = x_u[sect]
x_d = x_d[sect]
y_u = y_u[sect]
y_d = y_d[sect]
vtx_u = np.c_[x_u,y_u]
vtx_d = np.c_[x_d,y_d]
vtx = np.r_[vtx_u,np.flip(vtx_d,0)]
return vtx
def _my_poly_spiral(r,theta,order,res,R_max,sz_restrict=None):
''' generating a poly spiral curve
Args
r (2*1 list) :r[0] is the begining
theta (2*1 list) :theta[0] is the begining [in degree]
Return
frame (nazca.cell):
'''
theta[0] = theta[0]/180*np.pi ## angle format changing
theta[1] = theta[1]/180*np.pi ## angle format changing
K_ends = np.array([1/r[0],1/r[1]]) ## definition of the curvature, r[0] is the beginnin and r[1] is the ending
L0 = np.abs(theta[0]-theta[1])/(K_ends[0] + (K_ends[1]-K_ends[0])*order/(order+1))
L = np.linspace(0,L0,int(np.floor(L0/res)+1)) ## L = [0:res:L0];
K = K_ends[0] + (K_ends[1] - K_ends[0])/np.power(L0,order)*(np.power(L0,order) - np.power(np.abs(L-L0),order))
R = 1/K
dir = np.sign(theta[1] - theta[0])
dt = dir*res/R
theta_temp = np.cumsum(dt) + theta[0]
""" 2023.08.01 updated, using array calculation instead of for loop"""
dx = dir*R[1:]*( np.sin(theta_temp[1:]) - np.sin(theta_temp[0:-1]))
dy = -dir*R[1:]*( np.cos(theta_temp[1:]) - np.cos(theta_temp[0:-1]))
x = np.r_[0,np.cumsum(dx)]
y = np.r_[0,np.cumsum(dy)]
# x = np.zeros(len(L))
# y = np.zeros(len(L))
# idx = np.linspace(1,len(L)-1,len(L)-1)
# for _idx_ in idx :
# _idx_ = int(_idx_)
# x[_idx_] = x[_idx_-1] + dir*R[_idx_]*( np.sin(theta_temp[_idx_]) - np.sin(theta_temp[_idx_-1]))
# y[_idx_] = y[_idx_-1] - dir*R[_idx_]*( np.cos(theta_temp[_idx_]) - np.cos(theta_temp[_idx_-1]))
vector = np.c_[x,y,theta_temp,L]
return (vector,L0)
class Clothoid:
def __init__(self,
name:str=None,
R: 'list|np.ndarray'=[10,20],
w: 'list|np.ndarray|float'=[0.4,0.5], ## w either has the length as R, or 1 element, or 2 element
A: 'list|np.ndarray'=[0,45],
width_type: str='sine',
spiral_order: float=1,
Rmax:float=10000,
dL_cal:float=0.001,
dL_wg: float = 0.1,
# n_points:int=1024,
xs:str='strip',
layer:str=None,
sharp_patch:bool=True,
end_patch : bool=True,
show_pins:bool=False) -> None:
"""_summary_
Args:
R (list|np.ndarray, optional): Curvature radius in each attaching point. Defaults to [10,20].
w (list|np.ndarray|float, optional): Width at each attaching point corresponding to R, or it can be set to one or two element. Defaults to [0.4,0.5].
A (list|np.ndarray, optional): Angle at each attaching point. Defaults to [0,45].
width_type (str, optional): The width function with length or angle 'linear' 'linear2' 'sine' 'sine2'. Defaults to 'sine'.
spiral_order (float, optional): The curvature order of spiral. Defaults to 1.
Rmax (float, optional): Maxmum radius. Defaults to 10000.
res (float, optional): Resolution in calculation. Defaults to 0.001.
n_points (int, optional): Resolution in GDS. Defaults to 1024.
xs (str, optional): XSection of the devices. Defaults to 'strip'.
layer (str, optional): Layer of the devices. Defaults to None.
sharp_patch (bool, optional): Either to patch. Defaults to True.
show_pins (bool, optional): Either to show pins. Defaults to False.
Raises:
Exception: _description_
Exception: _description_
"""
if (isinstance(w,int) or isinstance(w,float)):
w= np.array([w,w])
self.name = name
self.R = R
self.A = A
self.width_type = width_type
self.spiral_order = spiral_order
self.dL_cal = dL_cal
# self.n_points = n_points
self.xs =xs
self.layer = layer
self.dL_wg = dL_wg
if (len(R) != len(A)):
raise Exception("ERROR: In <mxpic::structures::Colthoid>, <A> and <R> are not matched in length, please keep len(A) = len(R)")
if (isinstance(spiral_order,int) or isinstance(spiral_order,float)):
spiral_order = spiral_order*np.ones(len(R)-1)
elif(isinstance(spiral_order,list)):
spiral_order = np.array(spiral_order)
## center curve routing
_idx_act_=0
for _idx_ in range(0,len(R)-1):
if ( abs( A[_idx_] - A[_idx_+1] )<0.001 ):
continue
vec_cur,L0_cur = _my_poly_spiral([R[_idx_],R[_idx_+1]],[A[_idx_],A[_idx_+1]],spiral_order[_idx_],dL_cal,Rmax)
_idx_act_ = _idx_act_+1
x_cur = vec_cur[:,0]
y_cur = vec_cur[:,1]
theta_cur = vec_cur[:,2]/np.pi*180 ## pointing vector
L_cur = vec_cur[:,3]
if (_idx_act_==1):
L = L_cur
x = x_cur
y = y_cur
theta = theta_cur
L0 = L0_cur
else :
L = np.r_[L,L_cur+L[-1]]
x = np.r_[x,x_cur+x[-1]]
y = np.r_[y,y_cur+y[-1]]
theta = np.r_[theta,theta_cur]
L0 = L0 + L0_cur
if (len(w)>2 and len(w)==len(R)):
w_cur = (w[_idx_+1]-w[_idx_])/L0_cur*L_cur + (w[_idx_])
if (_idx_act_==1):
w_fianl = w_cur
else :
w_fianl = np.r_[w_fianl,w_cur]
self.x = x
self.y = y
self.L = L
self.L0 = L0
self.theta = theta
self.vtx_center = np.c_[x,y]
self.end_patch = end_patch
self.sz = [np.abs(max(self.x) - min(self.x)),np.abs(max(self.y) - min(self.y))]
if (dL_wg!=None):
self.n_points = int(np.floor(self.L0/self.dL_wg)+1) ## overwrite n_points
# print("n points",self.n_points)
if (len(w)==2):
## width winding
if (width_type=='linear'):
w = (w[1]-w[0])/L0*L + w[0]
elif (width_type=='dual_linear'):
w = (w[1]-w[0])/L0/2*np.abs(L-L0/2) + w[0]
elif (width_type=='sine'):
w = (w[0]-w[1])*np.cos(theta/180*pi)*np.cos(theta/180*pi) + w[1]
elif (width_type=='dual_sine'):
w = (w[0]-w[1])*np.cos(theta/2/180*pi)*np.cos(theta/2/180*pi) + w[1]
elif (width_type=='crow_customize' or width_type=='pumpkin'):
dw = (w[1]-w[0])
z = theta/180*np.pi
z = np.sqrt(z)*np.sqrt(np.pi/2)
z = np.sin(z)**2*np.pi/2
w = dw*np.sin(z)**2 + w[0]
else :
w = (w[1]-w[0])/L0*L + w[0]
self.w = np.array(w)
elif (len(w)==len(R)):
self.w = w_fianl
else:
raise Exception("ERROR, In <mxpic::structures::Clothoid>, <w> is not matched with <R>, please keep len(w)=2 or len(w)=R or w=int")
self.cell = self.generate_gds(sharp_patch=sharp_patch,show_pins=show_pins)
def generate_gds(self,sharp_patch,show_pins):
if (self.name is None):
self.instantiate = False
else:
self.instantiate = True
with nd.Cell(name=self.name,instantiate=self.instantiate) as C:
if (self.layer==None and self.xs!=None): ## if definition is in layers
for layers,growx,growy,acc in nd.layeriter(xs=self.xs):
(a1,b1), (a2,b2),c1,c2 = growx
if (b1!=0 and b2!=0):
vtx_wg = _line2wg_(x=self.x,y=self.y,wu=self.w*a1+b1,wd= -self.w*a2-b2,theta=self.theta,n_points=self.n_points)
dX = np.max(vtx_wg[:,0]) - np.min(vtx_wg[:,0])
dY = np.max(vtx_wg[:,1]) - np.min(vtx_wg[:,1])
cX = np.max(vtx_wg[:,0])/2 + np.min(vtx_wg[:,0])/2
cY = np.max(vtx_wg[:,1])/2 + np.min(vtx_wg[:,1])/2
if (sharp_patch):
if (self.end_patch):
nd.strt(length = dX+(b1-b2),width = dY,layer=layers).put(cX-dX/2-b1,cY,0)
else:
nd.strt(length = dX,width = dY,layer=layers).put(cX-dX/2,cY,0)
else :
_my_polygon(layers,vtx_wg).put(0,0,0)
else :
vtx_wg = _line2wg_(x=self.x,y=self.y,wu=self.w*a1+b1,wd= -self.w*a2-b2,theta=self.theta,n_points=self.n_points)
self.vtx =vtx_wg
_my_polygon(layers,vtx_wg).put(0,0,0)
nd.Pin(name='a0',width=self.w[0],type="Optical:").put(self.x[0],self.y[0],self.A[0]+180)
nd.Pin(name='b0',width=self.w[-1],type="Optical:").put(self.x[-1],self.y[-1],self.A[-1])
## revised in 2026.06.07 by Qin Yue
# legacy: nd.Pin(name='a1',width=self.w[0],type="Optical:").put(self.x[0],self.y[0],self.A[0]+180)
nd.Pin(name='opt_a1',width=self.w[0],type="optical:").put(self.x[0],self.y[0],self.A[0]+180)
## revised in 2026.06.07 by Qin Yue
# legacy: nd.Pin(name='b1',width=self.w[-1],type="Optical:").put(self.x[-1],self.y[-1],self.A[-1])
nd.Pin(name='opt_b1',width=self.w[-1],type="optical:").put(self.x[-1],self.y[-1],self.A[-1])
elif(self.layer!=None) : ## if definition is in xsections
vtx_wg = _line2wg_(x=self.x,y=self.y,wu=self.w/2,wd= self.w/2,theta=self.theta,n_points=self.n_points)
_my_polygon(self.layer,vtx_wg).put(0,0,0)
nd.Pin(name='a0',width=self.w[0],type="Optical:").put(self.x[0],self.y[0],self.A[0]+180)
nd.Pin(name='b0',width=self.w[-1],type="Optical:").put(self.x[-1],self.y[-1],self.A[-1])
## revised in 2026.06.07 by Qin Yue
# legacy: nd.Pin(name='a1',width=self.w[0],type="Optical:").put(self.x[0],self.y[0],self.A[0]+180)
nd.Pin(name='opt_a1',width=self.w[0],type="optical:").put(self.x[0],self.y[0],self.A[0]+180)
## revised in 2026.06.07 by Qin Yue
# legacy: nd.Pin(name='b1',width=self.w[-1],type="Optical:").put(self.x[-1],self.y[-1],self.A[-1])
nd.Pin(name='opt_b1',width=self.w[-1],type="optical:").put(self.x[-1],self.y[-1],self.A[-1])
else:
raise Exception("ERROR: In <mxpic::structures::Colthoid>, <layer | xs> not defined")
if (show_pins):
nd.put_stub()
self.sz_p2p = [np.abs(self.x[-1] - self.x[0]),np.abs(self.y[-1] - self.y[0])]
return C