diff --git a/.gitignore b/.gitignore index 70dd000b..570ba83b 100644 --- a/.gitignore +++ b/.gitignore @@ -5,6 +5,7 @@ docs/build/ **/__pycache__/** build/ docs/_build +**/.ipynb_checkpoints/** *.msh *.xdmf @@ -25,4 +26,4 @@ docs/_build Manifest.toml # bibtex-tidy -*.bib.original \ No newline at end of file +*.bib.original diff --git a/README.md b/README.md index 33931f9e..4edf733f 100644 --- a/README.md +++ b/README.md @@ -62,10 +62,6 @@ The documentation is lagging behind the state of code, so there's several featur - Calculation of overlap-integrals and confinement-factors - Bragg grating cells - Grating coupler cells -- Eigenmode of a coaxial cable and its specific impedance -- Eigenmodes of electric transmission lines - and determining their propagation constant (in work) -- Static electric fields - Overlap integrals between waveguide modes - Overlap integral between a waveguide mode and a fiber mode - Coupled mode theory - coupling between adjacent waveguides @@ -77,6 +73,14 @@ The documentation is lagging behind the state of code, so there's several featur - Static thermal profiles - Transient thermal behavior +### Electronics problems + +- Coplanar waveguide RF design +- Eigenmode of a coaxial cable and its specific impedance +- Eigenmodes of electric transmission lines + and determining their propagation constant (in work) +- Static electric fields + Something missing? Feel free to open an [issue](https://github.com/HelgeGehring/femwell/issues) :) ## Contributors @@ -93,6 +97,7 @@ Something missing? Feel free to open an [issue](https://github.com/HelgeGehring/ - Markus DeMartini (Google) - Lucas Grosjean (Google, Femto-ST Institute) - Eliza Leung (University of Adelaide) +- Duarte Silva (Eindhoven University of Technology) Happy about every form of contribution - pull requests, feature requests, issues, questions, ... :) diff --git a/docs/_toc.yml b/docs/_toc.yml index c3f5fc2d..00679db5 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -46,6 +46,7 @@ parts: - file: photonics/examples/refinement.py - file: electronics/examples/capacitor.py - file: electronics/examples/coax_cable.py + - file: electronics/examples/RF CPW transmission line tutorial/RF CPW transmission line tutorial.py - file: electronics/examples/coplanar_waveguide_vary_width.py - file: electronics/examples/coplanar_waveguide_vary_gap.py diff --git a/docs/electronics/examples/RF_CPW_transmission_line_tutorial/RF_CPW_transmission_line_tutorial.py b/docs/electronics/examples/RF_CPW_transmission_line_tutorial/RF_CPW_transmission_line_tutorial.py new file mode 100644 index 00000000..e9ddab7c --- /dev/null +++ b/docs/electronics/examples/RF_CPW_transmission_line_tutorial/RF_CPW_transmission_line_tutorial.py @@ -0,0 +1,1916 @@ +# --- +# jupyter: +# jupytext: +# formats: ipynb,py:percent +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.16.2 +# kernelspec: +# display_name: Python 3 (ipykernel) +# language: python +# name: python3 +# --- + +import time + +# %% +from collections import OrderedDict + +import enlighten +import numpy as np +from matplotlib import pyplot as plt +from matplotlib.gridspec import GridSpec +from pint import UnitRegistry +from shapely.geometry import LineString, MultiLineString, box +from shapely.ops import clip_by_rect, linemerge, unary_union +from skfem import Basis, ElementTriP0, Functional +from skfem.helpers import inner +from skfem.io.meshio import from_meshio + +from femwell.maxwell.waveguide import ( + calculate_overlap, + calculate_scalar_product, + compute_modes, +) +from femwell.mesh import mesh_from_OrderedDict +from femwell.visualization import plot_domains + +# %matplotlib inline +# To make things go smoother I advise to use %matplotlib widget so that you can inspect the mesh and other figures more clearly + +# %% [markdown] +# In this notebook we aim to study a simple CPW structure by repicating the measurements from [1] as in the image below. We wish to retrieve data from the microwave index and attenuation. This case is a good example of when the losses inside the metal will contribute greatly to the losses and, therefore, new conformal techniques were required to properly model the CPW. In our case, we will use FEMWELL to achieve the same results and confirm the theory and benchmark the software with the measurements. +# +#
+# +#
+# +# Furthermore, this notebook will also provide insight on how to use it for the design of RF waveguides, by exploring the fine line between waveguide and circuit theories [4, 5]. +# +# ## References +# +# [1] E. Tuncer, Beom-Taek Lee, M. S. Islam, and D. P. Neikirk, “Quasi-static conductor loss calculations in transmission lines using a new conformal mapping technique,” IEEE Trans. Microwave Theory Techn., vol. 42, no. 9, pp. 1807–1815, Sep. 1994, doi: 10.1109/22.310592. +# +# [2] G. W. Slade and K. J. Webb, “Computation of characteristic impedance for multiple microstrip transmission lines using a vector finite element method,” IEEE Trans. Microwave Theory Techn., vol. 40, no. 1, pp. 34–40, Jan. 1992, doi: 10.1109/22.108320. +# +# [3] - S. van Berkel, A. Garufo, A. Endo, N. LLombart, and A. Neto, �Characterization of printed transmission lines at high frequencies,� in The 9th European Conference on Antennas and Propagation (EuCAP 2015), (Lisbon, Portugal), Apr. 2015. (Software available at https://terahertz.tudelft.nl/Research/project.php?id=74&ti=27 ) +# +# [4] R. B. Marks and D. F. Williams, “A general waveguide circuit theory,” J. RES. NATL. INST. STAN., vol. 97, no. 5, p. 533, Sep. 1992, doi: 10.6028/jres.097.024. +# +# [5] D. F. Williams, L. A. Hayden, and R. B. Marks, “A complete multimode equivalent-circuit theory for electrical design,” J. Res. Natl. Inst. Stand. Technol., vol. 102, no. 4, p. 405, Jul. 1997, doi: 10.6028/jres.102.029. +# + +# %% [markdown] +# # Defining geometry and mesh + +# %% [markdown] +# The first step is to define our material parameters and geometry. To avoid any unit mistakes (AND THEY DO HAPPEN), we'll use `pint` to track every unit and conversions. Let us start by defining the frequency range we want to work with (in practice, in the following we will only use one frequency for now, but let's prepare everything to integrate seamlessly in a loop): + +# %% +reg = UnitRegistry() + +# Define frequency range +freq = np.linspace(0.2, 45, 10) * reg.GHz +omega = 2 * np.pi * freq + +# %% [markdown] +# Now some universal constants: + +# %% +## Define universal constants +mu0 = 4 * np.pi * 1e-7 * reg.henry / reg.meter # vacuum magnetic permeability +e0 = 8.854e-12 * reg.farad * reg.meter**-1 +c = 3e8 * reg.meter * reg.second**-1 # m s^-1 +e = 1.602176634e-19 * reg.coulomb # Coulombs +kb = 1.380649e-23 * reg.meter**2 * reg.kg * reg.second**-2 * reg.kelvin**-1 +T = 300 * reg.kelvin + +# %% [markdown] +# For the geometry we will follow the parametrization: +# +#
+# +#
+# +# Note that the port width can be smaller than the `w_sig+2*sep+2*w_gnd`. This, however, has to be chosen with a priori knowledge of the field profile. In this case, we know that the supported fields will be tightly confined in the gaps between the ground and signal pads. Therefore, we need to place the simulation boundaries far enough away from that region so as to guarantee that there is no coupling between the CPW mode and the boundary. + +# %% +w_sig = 7 * reg.micrometer +sep = 10 * reg.micrometer +w_gnd = 100 * reg.micrometer +t_metal = 0.8 * reg.micrometer + +port_width = (w_sig + 2 * sep + 2 * w_gnd) * 0.5 + 10 * reg.micrometer # um +bottom_height = 100 * reg.micrometer +top_height = 100 * reg.micrometer + + +eps_r_sub = 13 * reg.dimensionless # substrate relative permitivity +sig_metal = 6e5 * reg.siemens / reg.centimeter # Metal conductivity + +# %% [markdown] +# Next we will define the geometry that `skfem` needs to mesh. We also want to leave the option to use a symmetry plane or not. This will be useful to speed up computations and filter even and odd modes out. Therefore, we define 'use_symmetry_plane' and then tell which plane it is. The way we'll do it is to define the full geometry and then cut all the polygons and lines by that plane. +# +# Apart from the geometrical polygons that are included in the image above, we also want to include two additional things: +# +# The first one, is an integration path *just* outside the metal tracks. These will be used to integrate the magnetic field so as to retrieve the current flowing longitudinally in the transmission line. The second one will be two lines that will be used to calculate a voltage integral. Details will follow on the sections below. +# +#
+# +#
+ +# %% +######## Define FEM simulation ########### +use_symmetry_plane = True * 0 +symmetry_plane = box(0, -np.inf, np.inf, np.inf) + +near_field_width = 50 * reg.micrometer +near_field_height = 10 * reg.micrometer +########################################## + +metal_sig_box = box( + -w_sig.to(reg.micrometer).magnitude / 2, + 0, + w_sig.to(reg.micrometer).magnitude / 2, + t_metal.to(reg.micrometer).magnitude, +) + +metal_gnd_left_box = box( + max( + (-w_sig / 2 - sep - w_gnd).to(reg.micrometer).magnitude, + -(port_width / 2).to(reg.micrometer).magnitude, + ), + 0, + (-w_sig / 2 - sep).to(reg.micrometer).magnitude, + t_metal.to(reg.micrometer).magnitude, +) + +metal_gnd_right_box = box( + (w_sig / 2 + sep).to(reg.micrometer).magnitude, + 0, + min( + (w_sig / 2 + sep + w_gnd).to(reg.micrometer).magnitude, + (port_width / 2).to(reg.micrometer).magnitude, + ), + t_metal.to(reg.micrometer).magnitude, +) + + +air_box = box( + (-port_width / 2).to(reg.micrometer).magnitude, + 0, + (port_width / 2).to(reg.micrometer).magnitude, + top_height.to(reg.micrometer).magnitude, +) + + +substrate_box = box( + (-port_width / 2).to(reg.micrometer).magnitude, + -(bottom_height).to(reg.micrometer).magnitude, + (port_width / 2).to(reg.micrometer).magnitude, + 0, +) + +surface = unary_union([air_box, substrate_box]) + +##Make a line for a path integral +### Right conductor +xmin_1, ymin_1, xmax_1, ymax_1 = metal_sig_box.bounds +xmin_2, ymin_2, xmax_2, ymax_2 = metal_gnd_right_box.bounds + +points = [ + ((xmax_1 + xmin_1) / 2, (ymax_1 + ymin_1) / 2), + (xmax_1, (ymax_1 + ymin_1) / 2), + (xmin_2, (ymax_2 + ymin_2) / 2), + ((xmax_2 + xmin_2) / 2, (ymax_2 + ymin_2) / 2), +] + +path_integral_right = LineString(points) + +### left conductor +xmin_1, ymin_1, xmax_1, ymax_1 = metal_sig_box.bounds +xmin_2, ymin_2, xmax_2, ymax_2 = metal_gnd_left_box.bounds + +points = [ + ((xmax_1 + xmin_1) / 2, (ymax_1 + ymin_1) / 2), + (xmin_1, (ymax_1 + ymin_1) / 2), + (xmax_2, (ymax_2 + ymin_2) / 2), + ((xmax_2 + xmin_2) / 2, (ymax_2 + ymin_2) / 2), +] + +path_integral_left = LineString(points) + +polygons = OrderedDict( + surface=LineString(surface.exterior), + metal_sig_interface=LineString( + clip_by_rect( + metal_sig_box.buffer( + min(t_metal.to(reg.um).magnitude / 20, w_sig.to(reg.um).magnitude / 10), + join_style="bevel", + ), + *surface.bounds, + ).exterior + ), + metal_gnd_left_interface=LineString( + clip_by_rect( + metal_gnd_left_box.buffer( + min(t_metal.to(reg.um).magnitude / 20, w_gnd.to(reg.um).magnitude / 10), + join_style="bevel", + ), + *surface.bounds, + ).exterior + ), + metal_gnd_right_interface=LineString( + clip_by_rect( + metal_gnd_right_box.buffer( + min(t_metal.to(reg.um).magnitude / 20, w_gnd.to(reg.um).magnitude / 10), + join_style="bevel", + ), + *surface.bounds, + ).exterior + ), + path_integral_right=path_integral_right, + path_integral_left=path_integral_left, + metal_sig=metal_sig_box, + metal_gnd_left=metal_gnd_left_box, + metal_gnd_right=metal_gnd_right_box, + air=air_box, + substrate=substrate_box, +) + + +if use_symmetry_plane: + keys_to_pop = [] + for key, poly in polygons.items(): + if poly.intersects(symmetry_plane) and not symmetry_plane.contains(poly): + poly_tmp = clip_by_rect(poly, *symmetry_plane.bounds) + + if type(poly_tmp) == MultiLineString: + polygons[key] = linemerge(poly_tmp) + + elif poly_tmp.is_empty: + keys_to_pop.append(key) + else: + polygons[key] = poly_tmp + elif not poly.intersects(symmetry_plane): + keys_to_pop.append(key) + + for key in keys_to_pop: + polygons.pop(key) + +# Add the boundary polygons so that you can set custom boundary conditions +surf_bounds = polygons["surface"].bounds + +left = LineString([(surf_bounds[0], surf_bounds[1]), (surf_bounds[0], surf_bounds[3])]) + +bottom = LineString([(surf_bounds[0], surf_bounds[1]), (surf_bounds[2], surf_bounds[1])]) + +right = LineString([(surf_bounds[2], surf_bounds[1]), (surf_bounds[2], surf_bounds[3])]) + +top = LineString([(surf_bounds[0], surf_bounds[3]), (surf_bounds[2], surf_bounds[3])]) + +polygons["left"] = left +polygons["bottom"] = bottom +polygons["right"] = right +polygons["top"] = top + +polygons.move_to_end("top", last=False) +polygons.move_to_end("right", last=False) +polygons.move_to_end("bottom", last=False) +polygons.move_to_end("left", last=False) + +polygons.pop("surface") + +resolutions = dict( + surface={"resolution": 100, "distance": 1}, + metal_sig_interface={"resolution": 0.1, "distance": 10}, + metal_gnd_interface={"resolution": 0.5, "distance": 10}, + path_integral_right={"resolution": 0.2, "distance": 10}, + path_integral_left={"resolution": 0.2, "distance": 10}, + metal_sig={"resolution": 0.1, "distance": 0.1, "SizeMax": 0.1}, + metal_gnd_left={"resolution": 0.5, "distance": 0.2, "SizeMax": 0.5}, + metal_gnd_right={"resolution": 0.5, "distance": 0.2, "SizeMax": 0.5}, + air={"resolution": 10, "distance": 0.1}, + substrate={"resolution": 10, "distance": 1}, +) + + +# %% [markdown] +# We can now check our polygons: + +# %% +fig = plt.figure() +ax = fig.add_subplot(111) + +for key, polygon in polygons.items(): + if type(polygon) is not LineString: + x, y = polygon.exterior.xy + ax.plot(np.asarray(x), np.asarray(y), color="pink") + else: + ax.plot(*polygon.xy) + + +# %% [markdown] +# Notice how in the above we have used the argument `SizeMax` on the metal tracks. We do this so that at a distance of `distance` of the boundary of the polygons a resolution of `SizeMax` is assured. This is important as we expect the field to vary rapidly near the metal. +# +#
+# +#
+# +# Now we mesh it + +# %% +# Mesh it +mesh = from_meshio( + mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=100, verbose=True) +) + +print(mesh.nelements) + +# %% +fig = mesh.draw() +fig.axes.set_axis_on() + +# %% [markdown] +# Note that we are using a **very** fine mesh. Normally, if we were using thick metals then surface impedance would suffice and we wouldn't have to mesh the entire surface of the metal. However, in this case, the metal is far too thin and the skin depth of the frequencies we are using is far too big to neglect current flowing inside the metal. Sadly, the profile of $I(x,y)$ inside the metal is rapidly changing, so we need the very fine mesh. +# +# How do we know that it is fine enough? We won't do that here, but for a thorough study we advise to do a sweep on the resolution inside the metal and keep checking the absorption of the modes. Once it converges, you're good. +# +# For now, we'll move on to the definition of the materials. FEMWELL allows the definition of the $\epsilon(x,y)$ such that at each polygon you attribute a material constant as: +# +# $$ +# \epsilon = \epsilon^` +j\epsilon^{``} +# $$ +# +# where, for conductive materials: +# +# $$ +# \epsilon = -\frac{\sigma}{\omega \epsilon_0} +# $$ +# + +# %% +idx_freq = -1 +print(f"Frequency:{freq[idx_freq].magnitude:.2f} GHz") + +basis0 = Basis(mesh, ElementTriP0(), intorder=4) # Define the basis for the FEM +epsilon = basis0.ones(dtype=complex) + +epsilon[basis0.get_dofs(elements="air")] = 1 +epsilon[basis0.get_dofs(elements="substrate")] = eps_r_sub +epsilon[basis0.get_dofs(elements="metal_sig")] = ( + 1 - 1j * (sig_metal / omega / e0)[idx_freq].to(reg.dimensionless).magnitude +) +epsilon[basis0.get_dofs(elements="metal_gnd_left")] = ( + 1 - 1j * (sig_metal / omega / e0)[idx_freq].to(reg.dimensionless).magnitude +) +epsilon[basis0.get_dofs(elements="metal_gnd_right")] = ( + 1 - 1j * (sig_metal / omega / e0)[idx_freq].to(reg.dimensionless).magnitude +) + +# %% +fig = basis0.plot(epsilon.real, colorbar=True) +fig.axes.set_axis_on() + +fig = basis0.plot(epsilon.imag, colorbar=True) +fig.axes.set_axis_on() + +# %% [markdown] +# Let's compute the modes. Since we have 3 conductors, we expect only 2 modes to be present: an even and an odd mode. So let's just calculate 2 + +# %% +start = time.time() +modes = compute_modes( + basis0, + epsilon, + wavelength=(c / freq[idx_freq]).to(reg.micrometer).magnitude, + num_modes=2, + metallic_boundaries=False, +) +stop = time.time() + +modes = modes.sorted(key=lambda mode: mode.n_eff.real) + +print(stop - start) + +# %% +print(modes.n_effs.real) +print( + 20 + / np.log(10) + * (modes.n_effs.imag * 2 * np.pi * freq[idx_freq] / c).to(reg.centimeter**-1).magnitude, + "dB/cm", +) + +# %% [markdown] +# We should now inspect the modes. You can use `modes[i].plot(modes.E)` or even plot a single component with `modes[i].plot_component('E', 'x')`, but we're looking for something more custom, I will export the data onto a rectangular grid and plot a streamplot to visualize the field lines: + +# %% +from skfem import ElementDG, ElementTriP1, ElementVector + +# Trying interpolator +Nx = 50 +Ny = 50 +grid_data_E = np.zeros((2, Ny, Nx, 3), dtype=complex) +grid_data_H = np.zeros((2, Ny, Nx, 3), dtype=complex) + +xmin = -40 +xmax = 40 +ymin = -10 +ymax = 10 + +grid_x, grid_y = np.meshgrid(np.linspace(xmin, xmax, Nx), np.linspace(ymin, ymax, Ny)) + +for i, mode in enumerate(modes): + basis = mode.basis + basis_fix = basis.with_element(ElementVector(ElementDG(ElementTriP1()))) + + (et, et_basis), (ez, ez_basis) = basis.split(mode.E) + (et_x, et_x_basis), (et_y, et_y_basis) = basis_fix.split( + basis_fix.project(et_basis.interpolate(et)) + ) + + coordinates = np.array([grid_x.flatten(), grid_y.flatten()]) + + start = time.time() + grid_data = np.array( + ( + et_x_basis.interpolator(et_x)(coordinates), + et_y_basis.interpolator(et_y)(coordinates), + ez_basis.interpolator(ez)(coordinates), + ) + ).T + + grid_data_E[i] = grid_data.reshape((*grid_x.shape, -1)) + end = time.time() + print(end - start) + + (et, et_basis), (ez, ez_basis) = basis.split(mode.H) + (et_x, et_x_basis), (et_y, et_y_basis) = basis_fix.split( + basis_fix.project(et_basis.interpolate(et)) + ) + + coordinates = np.array([grid_x.flatten(), grid_y.flatten()]) + + grid_data = np.array( + ( + et_x_basis.interpolator(et_x)(coordinates), + et_y_basis.interpolator(et_y)(coordinates), + ez_basis.interpolator(ez)(coordinates), + ) + ).T + + grid_data_H[i] = grid_data.reshape((*grid_x.shape, -1)) + +# %% +fig = plt.figure(figsize=(7, 5)) +gs = GridSpec(2, 2) + +ax_even_E = fig.add_subplot(gs[0, 0]) +ax_odd_E = fig.add_subplot(gs[0, 1]) + +ax_even_H = fig.add_subplot(gs[1, 0]) +ax_odd_H = fig.add_subplot(gs[1, 1]) + +for ax, data, label in zip( + [ax_even_E, ax_odd_E, ax_even_H, ax_odd_H], + [grid_data_E[0], grid_data_E[1], grid_data_H[0], grid_data_H[1]], + [ + r"Mode 0 | $|E(x,y)|$", + r"Mode 1 | $|E(x,y)|$", + r"Mode 0 | $|H(x,y)|$", + r"Mode 1 | $|H(x,y)|$", + ], +): + + ax.imshow( + np.sum(np.abs(data), axis=2), + origin="lower", + extent=[xmin, xmax, ymin, ymax], + cmap="jet", + interpolation="bicubic", + aspect="auto", + ) + ax.streamplot( + grid_x, grid_y, data.real[:, :, 0], data.real[:, :, 1], color="black", linewidth=0.5 + ) + + for key, polygon in polygons.items(): + if type(polygon) is not LineString: + x, y = polygon.exterior.xy + ax.plot(np.asarray(x), np.asarray(y), color="pink") + + ax.set_xlim(xmin, xmax) + ax.set_ylim(ymin, ymax) + + ax.set_xlabel("x (um)") + ax.set_ylabel("y (um)") + + ax.set_title(label) + +fig.tight_layout() +# fig.savefig('modes.png', dpi = 400) + +# %% [markdown] +# # Transmission line characterization +# +# So far we have been foccused on studying the eigenmodes of maxwell's equations for a particular structure. This allows us to calculate $\mathbf{E(x,y)}_i$ and $\mathbf{H(x,y)}_i$ for each of the available eigenmodes. However, we are interested in studying the structure as a transmission line and for that we must start merging towards circuit theory. However, that is not as straightforward as it may appear. First, our modes are not TEM fields, and, therefore, cannot exist in a transmission line. Second, this is far from a lossless line, and it is also supports multiple modes, which not only leads to the requirement of the expressing the line as a multiconductor transmission line, but can also lead to a non-reciprocal line, depending on the formalism you decide to follow [4,5]. It is this ambiguity that can lead to a big confusion when attempting to translate a waveguiding system to a circuit. Here we will follow the formalism of [4] and eventually fully characterize this system as a transmission line. +# + +# %% [markdown] +# The first thing to check is if our modes are transverse enough. Fundamentally they are not, but if they're close enough, it's a good start. We can actually check this with FEMWELL by computing: +# +# $$ +# T = \int_S \frac{|E_x|^2 + |E_y|^2}{|E_z|^2} dS +# $$ + +# %% +print(modes[0].transversality) +print(modes[1].transversality) + + +# %% [markdown] +# We're good. Now we need to take care of the fact that we have a multimode system. More often than not, we will be interested in single mode operation. In the case of a CPW, we are interested in the ground-signal-ground excitation, which represents the fundamental mode. If we excite only this mode we can treat the system as a transmisison line, greatly simplifying our work. On that assumption let us now define the power integral as: +# +# $$ +# p(z) = \int_S \mathbf{E_t} \times \mathbf{H_t}^* dS +# $$ +# +# where $\mathbf{E_t}$, $\mathbf{H_t}$ denote the transverse components of the electric and magnetic field. These, in turn, can be defined as: +# +# $$ +# \mathbf{E_t} = (c_+ e^{-\gamma z} + c_- e^{\gamma z}) \mathbf{e}_t = \frac{v(z)}{v_0}\mathbf{e}_t +# $$ +# +# $$ +# \mathbf{H_t} = (c_+ e^{-\gamma z} - c_- e^{\gamma z}) \mathbf{h}_t = \frac{i(z)}{i_0}\mathbf{h}_t +# $$ +# +# where $i_0$, $v_0$, $i(z)$ and $v(z)$ all have units of current and voltage and the vectors have field units. With definition the power now becomes: +# +# $$ +# p(z) = \frac{v(z) i(z)^*}{v_0 i_0^*}p_0 +# $$ +# +# with: +# +# $$ +# p_0 = \int_S \mathbf{e_t} \times \mathbf{h_t}^* dS +# $$ +# +# Now we have quite some ambiguity left. What are the values of $v_0$, $i_0$ and $p_0$? There is no fundamental constraint telling us what the values should be and we can define any value (any 2 of the 3 variables) and still be correct in the sense that we will still end up with eigenmodes of the system. However, we have started this discussion wanting to transition to a circuit, therefore, we must start converging towards it. Since we can define some currents and voltages, let us identify these with the ones that make sense to be present in a circuit picture. Therefore, we can define: +# +# $$ +# p_0 = v_0 i_0^* +# $$ +# +# It is important to impose the constraint that $Re(p_0)>0$[4]. By fixing the power, we now have only one degree of freedom left to define. We can either define a current or a voltage. Since we want to have a clear and illustrative analogy with a circuit we must choose either $v_0$ or $i_0$ as circuit quantities. In the case of the CPW operating the GSG mode, we can make the analogy as in the picture below. +# +#
+# +#
+# +# With this model, it is now clear that we can define $i_0$ as: +# +# $$ +# i_0 = \int_C \mathbf{h_t} \cdot dl +# $$ +# +# where the integration is on the contour of the signal track. Having this quantity the voltage $v_0$ follows from the power definition. It is important that we define only two of the three unknown variables as they are not independent. This previous method is what is called in the literature as the PI model. Alternatively, we can define: +# +# $$ +# v_0 = -\int_{path} \mathbf{e_t} \cdot dl +# $$ +# +# where the integration path is defined so as to capture the voltage difference between the ground and signal tracks. This is called the PV model. Finally, it is also possible to define simultaneously the $i_0$ and the $v_0$ as above, and THEN define the power $p_0$. This is called the IV (or VI) model. +# +# +# With this in mind, it is clear that the characteristic impedance $Z_0$ follows from: +# +# $$ +# Z_0 = \frac{v_0}{i_0} = \frac{|v_0|^2}{p_0^*} = {p_0}{|i_0|^2} +# $$ +# +# Let us test this: + + +# %% +@Functional(dtype=np.complex64) +def current_form(w): + """ + What this does is it takes the normal vector to the boundary and rotates it 90deg + Then takes the inner product with the magnetic field in it's complex form + """ + return inner(np.array([w.n[1], -w.n[0]]), w.H) + + +@Functional(dtype=np.complex64) +def voltage_form(w): + """ + What this does is it takes the normal vector to the boundary and rotates it 90deg + Then takes the inner product with the electric field in it's complex form + """ + + return -inner(np.array([w.n[1], -w.n[0]]), w.E) + + +conductor = "metal_sig_interface" +line = "path_integral_right" +mode = modes[0] + +p0 = calculate_scalar_product( + mode.basis, np.conjugate(mode.E), mode.basis, np.conjugate(mode.H) +) # The conjugate is due to femwell internal definition as conj(E) x H + + +(ht, ht_basis), (hz, hz_basis) = mode.basis.split(mode.H) +facet_basis = ht_basis.boundary(facets=mesh.boundaries[conductor]) +i0 = current_form.assemble(facet_basis, H=facet_basis.interpolate(ht)) + +(et, et_basis), (ez, ez_basis) = mode.basis.split(mode.E) +facet_basis = et_basis.boundary(facets=mesh.boundaries[line]) +v0 = voltage_form.assemble(facet_basis, E=facet_basis.interpolate(et)) + + +print( + f"PI model : |Z0| = {np.abs(p0/np.abs(i0)**2):.2f} Ohm | angle(Z0) = {np.angle(p0/np.abs(i0)**2)/np.pi:.2f} pi radians" +) +print( + f"PV model : |Z0| = {np.abs(np.abs(v0)**2/p0.conj()):.2f} Ohm | angle(Z0) = {np.angle(np.abs(v0)**2/p0.conj())/np.pi:.2f} pi radians" +) +print( + f"VI model : |Z0| = {np.abs(v0/i0):.2f} Ohm | angle(Z0) = {np.angle(v0/i0)/np.pi:.2f} pi radians" +) + + +# %% [markdown] +# As you can see all the model give a very close result for the characteristic impedance. +# +# Another interesting aspect of waveguide circuit theory is the fact that it can be shown that, under the quasi-TEM approximation, it is possible to **analytically** extract the RLGC parameters of a circuit via the field solutions. This fits perfectly with the FEM solutions as the evaluation of the integrals is straightforward. The parameters are extracted as: +# +# $$ +# C = \frac{1}{|v_0|^2}\left[ \int_S \epsilon^\prime |\mathbf{e}_t|^2 dS - \int_S \mu^\prime |\mathbf{h}_z|^2 dS \right] +# $$ +# +# $$ +# L = \frac{1}{|i_0|^2}\left[ \int_S \mu^\prime |\mathbf{h}_t|^2 dS - \int_S \epsilon^\prime |\mathbf{e}_z|^2 dS \right] +# $$ +# +# $$ +# G = \frac{\omega}{|v_0|^2}\left[ \int_S \epsilon^{\prime\prime} |\mathbf{e}_t|^2 dS + \int_S \mu^{\prime\prime} |\mathbf{h}_z|^2 dS \right] +# $$ +# +# $$ +# R = \frac{\omega}{|i_0|^2}\left[ \int_S \mu^{\prime\prime} |\mathbf{h}_t|^2 dS + \int_S \epsilon^{\prime\prime} |\mathbf{e}_z|^2 dS \right] +# $$ + + +# %% +@Functional(dtype=np.complex64) +def C_form(w): + + return ( + 1 + / np.abs(w.v0) ** 2 + * ( + np.real(w.epsilon) * inner(w.et, np.conj(w.et)) + - np.real(w.mu) * inner(w.hz, np.conj(w.hz)) + ) + ) + + +@Functional(dtype=np.complex64) +def L_form(w): + + return ( + 1 + / np.abs(w.i0) ** 2 + * ( + np.real(w.mu) * inner(w.ht, np.conj(w.ht)) + - np.real(w.epsilon) * inner(w.ez, np.conj(w.ez)) + ) + ) + + +@Functional(dtype=np.complex64) +def G_form(w): + # The minus sign account for the fact that in the paper they define eps = eps_r-1j*eps_i + # whereas with python we have eps = eps_r+1j*eps_i. + return ( + -w.omega + / np.abs(w.v0) ** 2 + * ( + np.imag(w.epsilon) * inner(w.et, np.conj(w.et)) + + np.imag(w.mu) * inner(w.hz, np.conj(w.hz)) + ) + ) + + +@Functional(dtype=np.complex64) +def R_form(w): + # The minus sign account for the fact that in the paper they define eps = eps_r-1j*eps_i + # whereas with python we have eps = eps_r+1j*eps_i. + return ( + -w.omega + / np.abs(w.i0) ** 2 + * ( + np.imag(w.epsilon) * inner(w.ez, np.conj(w.ez)) + + np.imag(w.mu) * inner(w.ht, np.conj(w.ht)) + ) + ) + + +basis_t = et_basis +basis_z = ez_basis +basis_eps = basis0 + +# Be careful with the units!! +# In this case just make sure you adjust the length unit to micrometers +# everything else can stay as is. + +C = C_form.assemble( + mode.basis, + epsilon=basis_eps.interpolate(epsilon * e0.to(reg.farad / reg.micrometer).magnitude), + mu=mu0.to(reg.henry / reg.micrometer).magnitude, + i0=i0, + omega=omega[idx_freq].to(reg.second**-1).magnitude, + v0=v0, + et=basis_t.interpolate(et), + ez=basis_z.interpolate(ez), + ht=basis_t.interpolate(ht), + hz=basis_z.interpolate(hz), +) + +L = L_form.assemble( + mode.basis, + epsilon=basis_eps.interpolate(epsilon * e0.to(reg.farad / reg.micrometer).magnitude), + mu=mu0.to(reg.henry / reg.micrometer).magnitude, + i0=i0, + omega=omega[idx_freq].to(reg.second**-1).magnitude, + v0=v0, + et=basis_t.interpolate(et), + ez=basis_z.interpolate(ez), + ht=basis_t.interpolate(ht), + hz=basis_z.interpolate(hz), +) + +G = G_form.assemble( + mode.basis, + epsilon=basis_eps.interpolate(epsilon * e0.to(reg.farad / reg.micrometer).magnitude), + mu=mu0.to(reg.henry / reg.micrometer).magnitude, + i0=i0, + omega=omega[idx_freq].to(reg.second**-1).magnitude, + v0=v0, + et=basis_t.interpolate(et), + ez=basis_z.interpolate(ez), + ht=basis_t.interpolate(ht), + hz=basis_z.interpolate(hz), +) + +R = R_form.assemble( + mode.basis, + epsilon=basis_eps.interpolate(epsilon * e0.to(reg.farad / reg.micrometer).magnitude), + mu=mu0.to(reg.henry / reg.micrometer).magnitude, + i0=i0, + omega=omega[idx_freq].to(reg.second**-1).magnitude, + v0=v0, + et=basis_t.interpolate(et), + ez=basis_z.interpolate(ez), + ht=basis_t.interpolate(ht), + hz=basis_z.interpolate(hz), +) + +Z0 = np.sqrt( + (R + 1j * omega[idx_freq].to(reg.second**-1).magnitude * L) + / (G + 1j * omega[idx_freq].to(reg.second**-1).magnitude * C) +) + +gamma = np.sqrt( + (R + 1j * omega[idx_freq].to(reg.second**-1).magnitude * L) + * (G + 1j * omega[idx_freq].to(reg.second**-1).magnitude * C) +) + +print( + f"PI model : |Z0| = {np.abs(p0/np.abs(i0)**2):.2f} Ohm | angle(Z0) = {np.angle(p0/np.abs(i0)**2)/np.pi:.2f} pi radians" +) +print( + f"PV model : |Z0| = {np.abs(np.abs(v0)**2/p0.conj()):.2f} Ohm | angle(Z0) = {np.angle(np.abs(v0)**2/p0.conj())/np.pi:.2f} pi radians" +) +print( + f"VI model : |Z0| = {np.abs(v0/i0):.2f} Ohm | angle(Z0) = {np.angle(v0/i0)/np.pi:.2f} pi radians" +) +print(f"RLGC : |Z0| = {np.abs(Z0):.2f} Ohm | angle(Z0) = {np.angle(Z0)/np.pi:.2f} pi radians") + +print("=================== RLGC ==================") +print(f"R = {R.real/1e-3:.2f} mOhm/um") +print(f"L = {L.real/1e-12:.2f} picoHenry/um") +print(f"G = {G.real/1e-12:.2f} picoSiemens/um") +print(f"C = {C.real/1e-15:.2f} femtoFarad/um") + +print("================== propagation constant ===============") +print(f"beta : RLGC = {gamma.imag*1e6:.2f} 1/m | FEM = {mode.k.real*1e6:.2f} 1/m ") +print(f"alpha: RLGC = {gamma.real*1e6:.2f} 1/m | FEM = {-mode.k.imag*1e6:.2f} 1/m ") + +# %% [markdown] +# While there are some slight differences, the results match quite well!! Now, let us try to repeat the process but including a symmetry plane. + +# %% [markdown] +# Let's try to the same, but now with the symmetry plane: + +# %% +######## Define FEM simulation ########### +use_symmetry_plane = True +symmetry_plane = box(0, -np.inf, np.inf, np.inf) + +near_field_width = 50 * reg.micrometer +near_field_height = 10 * reg.micrometer +########################################## + + +metal_sig_box = box( + -w_sig.to(reg.micrometer).magnitude / 2, + 0, + w_sig.to(reg.micrometer).magnitude / 2, + t_metal.to(reg.micrometer).magnitude, +) + +metal_gnd_left_box = box( + max( + (-w_sig / 2 - sep - w_gnd).to(reg.micrometer).magnitude, + -(port_width / 2).to(reg.micrometer).magnitude, + ), + 0, + (-w_sig / 2 - sep).to(reg.micrometer).magnitude, + t_metal.to(reg.micrometer).magnitude, +) + +metal_gnd_right_box = box( + (w_sig / 2 + sep).to(reg.micrometer).magnitude, + 0, + min( + (w_sig / 2 + sep + w_gnd).to(reg.micrometer).magnitude, + (port_width / 2).to(reg.micrometer).magnitude, + ), + t_metal.to(reg.micrometer).magnitude, +) + + +air_box = box( + (-port_width / 2).to(reg.micrometer).magnitude, + 0, + (port_width / 2).to(reg.micrometer).magnitude, + top_height.to(reg.micrometer).magnitude, +) + + +substrate_box = box( + (-port_width / 2).to(reg.micrometer).magnitude, + -(bottom_height).to(reg.micrometer).magnitude, + (port_width / 2).to(reg.micrometer).magnitude, + 0, +) + +surface = unary_union([air_box, substrate_box]) + +##Make a line for a path integral +### Right conductor +xmin_1, ymin_1, xmax_1, ymax_1 = metal_sig_box.bounds +xmin_2, ymin_2, xmax_2, ymax_2 = metal_gnd_right_box.bounds + +points = [ + ((xmax_1 + xmin_1) / 2, (ymax_1 + ymin_1) / 2), + (xmax_1, (ymax_1 + ymin_1) / 2), + (xmin_2, (ymax_2 + ymin_2) / 2), + ((xmax_2 + xmin_2) / 2, (ymax_2 + ymin_2) / 2), +] + +path_integral_right = LineString(points) + +### left conductor +xmin_1, ymin_1, xmax_1, ymax_1 = metal_sig_box.bounds +xmin_2, ymin_2, xmax_2, ymax_2 = metal_gnd_left_box.bounds + +points = [ + ((xmax_1 + xmin_1) / 2, (ymax_1 + ymin_1) / 2), + (xmin_1, (ymax_1 + ymin_1) / 2), + (xmax_2, (ymax_2 + ymin_2) / 2), + ((xmax_2 + xmin_2) / 2, (ymax_2 + ymin_2) / 2), +] + +path_integral_left = LineString(points) + +polygons = OrderedDict( + surface=LineString(surface.exterior), + metal_sig_interface=LineString( + clip_by_rect( + metal_sig_box.buffer( + min(t_metal.to(reg.um).magnitude / 20, w_sig.to(reg.um).magnitude / 10), + join_style="bevel", + ), + *surface.bounds, + ).exterior + ), + metal_gnd_left_interface=LineString( + clip_by_rect( + metal_gnd_left_box.buffer( + min(t_metal.to(reg.um).magnitude / 20, w_gnd.to(reg.um).magnitude / 10), + join_style="bevel", + ), + *surface.bounds, + ).exterior + ), + metal_gnd_right_interface=LineString( + clip_by_rect( + metal_gnd_right_box.buffer( + min(t_metal.to(reg.um).magnitude / 20, w_gnd.to(reg.um).magnitude / 10), + join_style="bevel", + ), + *surface.bounds, + ).exterior + ), + path_integral_right=path_integral_right, + path_integral_left=path_integral_left, + metal_sig=metal_sig_box, + metal_gnd_left=metal_gnd_left_box, + metal_gnd_right=metal_gnd_right_box, + air=air_box, + substrate=substrate_box, +) + + +if use_symmetry_plane: + keys_to_pop = [] + for key, poly in polygons.items(): + if poly.intersects(symmetry_plane) and not symmetry_plane.contains(poly): + poly_tmp = clip_by_rect(poly, *symmetry_plane.bounds) + + if type(poly_tmp) == MultiLineString: + polygons[key] = linemerge(poly_tmp) + + elif poly_tmp.is_empty: + keys_to_pop.append(key) + else: + polygons[key] = poly_tmp + elif not poly.intersects(symmetry_plane): + keys_to_pop.append(key) + + for key in keys_to_pop: + polygons.pop(key) + + +# Add the boundary polygons so that you can set custom boundary conditions +surf_bounds = polygons["surface"].bounds + +left = LineString([(surf_bounds[0], surf_bounds[1]), (surf_bounds[0], surf_bounds[3])]) + +bottom = LineString([(surf_bounds[0], surf_bounds[1]), (surf_bounds[2], surf_bounds[1])]) + +right = LineString([(surf_bounds[2], surf_bounds[1]), (surf_bounds[2], surf_bounds[3])]) + +top = LineString([(surf_bounds[0], surf_bounds[3]), (surf_bounds[2], surf_bounds[3])]) + +polygons["left"] = left +polygons["bottom"] = bottom +polygons["right"] = right +polygons["top"] = top + +polygons.move_to_end("top", last=False) +polygons.move_to_end("right", last=False) +polygons.move_to_end("bottom", last=False) +polygons.move_to_end("left", last=False) + +polygons.pop("surface") + +resolutions = dict( + metal_sig_interface={"resolution": 0.1, "distance": 10}, + metal_gnd_interface={"resolution": 0.5, "distance": 10}, + path_integral_right={"resolution": 0.2, "distance": 10}, + path_integral_left={"resolution": 0.2, "distance": 10}, + metal_sig={"resolution": 0.1, "distance": 0.1, "SizeMax": 0.1}, + metal_gnd_left={"resolution": 0.5, "distance": 0.2, "SizeMax": 0.5}, + metal_gnd_right={"resolution": 0.5, "distance": 0.2, "SizeMax": 0.5}, + left={"resolution": 10, "distance": 0.1}, + right={"resolution": 10, "distance": 0.1}, + top={"resolution": 10, "distance": 0.1}, + bottom={"resolution": 10, "distance": 0.1}, + air={"resolution": 10, "distance": 0.1}, + substrate={"resolution": 10, "distance": 1}, +) + +# %% +fig = plt.figure() +ax = fig.add_subplot(111) + +for key, polygon in polygons.items(): + # print(polygon) + if type(polygon) is not LineString: + x, y = polygon.exterior.xy + ax.plot(np.asarray(x), np.asarray(y), color="pink") + else: + ax.plot(*polygon.xy) + + +# %% +# Mesh it +mesh = from_meshio( + mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=100, verbose=True) +) + +print(mesh.nelements) + +# %% +mesh.draw() + +# %% +idx_freq = -1 +print(f"Frequency:{freq[idx_freq].magnitude:.2f} GHz") + +basis0 = Basis(mesh, ElementTriP0(), intorder=4) # Define the basis for the FEM +epsilon = basis0.ones(dtype=complex) + +epsilon[basis0.get_dofs(elements="air")] = 1 +epsilon[basis0.get_dofs(elements="substrate")] = eps_r_sub +epsilon[basis0.get_dofs(elements="metal_sig")] = ( + 1 - 1j * (sig_metal / omega / e0)[idx_freq].to(reg.dimensionless).magnitude +) +# epsilon[basis0.get_dofs(elements = 'metal_gnd_left')] = (1 - 1j*(sig_metal/omega/e0)[idx_freq].to(reg.dimensionless).magnitude) +epsilon[basis0.get_dofs(elements="metal_gnd_right")] = ( + 1 - 1j * (sig_metal / omega / e0)[idx_freq].to(reg.dimensionless).magnitude +) + +# %% +start = time.time() +modes = compute_modes( + basis0, + epsilon, + wavelength=(c / freq[idx_freq]).to(reg.micrometer).magnitude, + num_modes=1, + metallic_boundaries=[], +) + +stop = time.time() + +print(stop - start) + +# %% +print(modes.n_effs.real) +print( + 20 + / np.log(10) + * (modes.n_effs.imag * 2 * np.pi * freq[idx_freq] / c).to(reg.centimeter**-1).magnitude, + "dB/cm", +) + +# %% [markdown] +# Much faster this time round + +# %% +from skfem import ElementDG, ElementTriP1, ElementVector + +# Trying interpolator +Nx = 25 +Ny = 50 +grid_data_E = np.zeros((len(modes), Ny, Nx, 3), dtype=complex) +grid_data_H = np.zeros((len(modes), Ny, Nx, 3), dtype=complex) + +xmin = 0 +xmax = 40 +ymin = -10 +ymax = 10 + +grid_x, grid_y = np.meshgrid(np.linspace(xmin, xmax, Nx), np.linspace(ymin, ymax, Ny)) + +for i, mode in enumerate(modes): + basis = mode.basis + basis_fix = basis.with_element(ElementVector(ElementDG(ElementTriP1()))) + + (et, et_basis), (ez, ez_basis) = basis.split(mode.E) + (et_x, et_x_basis), (et_y, et_y_basis) = basis_fix.split( + basis_fix.project(et_basis.interpolate(et)) + ) + + coordinates = np.array([grid_x.flatten(), grid_y.flatten()]) + + start = time.time() + grid_data = np.array( + ( + et_x_basis.interpolator(et_x)(coordinates), + et_y_basis.interpolator(et_y)(coordinates), + ez_basis.interpolator(ez)(coordinates), + ) + ).T + + grid_data_E[i] = grid_data.reshape((*grid_x.shape, -1)) + end = time.time() + print(end - start) + + (et, et_basis), (ez, ez_basis) = basis.split(mode.H) + (et_x, et_x_basis), (et_y, et_y_basis) = basis_fix.split( + basis_fix.project(et_basis.interpolate(et)) + ) + + coordinates = np.array([grid_x.flatten(), grid_y.flatten()]) + + grid_data = np.array( + ( + et_x_basis.interpolator(et_x)(coordinates), + et_y_basis.interpolator(et_y)(coordinates), + ez_basis.interpolator(ez)(coordinates), + ) + ).T + + grid_data_H[i] = grid_data.reshape((*grid_x.shape, -1)) + +# %% +fig = plt.figure(figsize=(4, 5)) +gs = GridSpec(2, 1) + +ax_even_E = fig.add_subplot(gs[0, 0]) +ax_even_H = fig.add_subplot(gs[1, 0]) + +mode_idx = 0 +for ax, data, label in zip( + [ax_even_E, ax_even_H], + [grid_data_E[mode_idx], grid_data_H[mode_idx]], + [r"Even mode | $|E(x,y)|$", r"Even mode | $|H(x,y)|$"], +): + + ax.imshow( + np.sum(np.abs(data), axis=2), + origin="lower", + extent=[xmin, xmax, ymin, ymax], + cmap="jet", + interpolation="bicubic", + aspect="auto", + ) + ax.streamplot( + grid_x, grid_y, data.real[:, :, 0], data.real[:, :, 1], color="black", linewidth=0.5 + ) + + for key, polygon in polygons.items(): + if type(polygon) is not LineString: + x, y = polygon.exterior.xy + ax.plot(np.asarray(x), np.asarray(y), color="pink") + + ax.set_xlim(xmin, xmax) + ax.set_ylim(ymin, ymax) + + ax.set_xlabel("x (um)") + ax.set_ylabel("y (um)") + + ax.set_title(label) + +fig.tight_layout() +# fig.savefig('modes_sym.png', dpi = 400) + +# %% [markdown] +# Now you have to be careful and account for a small artifact. Not only you are dealing with half the field, but the current integral that you do is half of the one calculated previously. Plus, the calculated fields also have different normalizations. Interestingly, if we apply the exact same algorithms without accounting for this fact, what we end up with the values corresponding to a circuit in parallel as below: +# +#
+# +#
+# +# + + +# %% +@Functional(dtype=np.complex64) +def current_form(w): + """ + What this does is it takes the normal vector to the boundary and rotates it 90deg + Then takes the inner product with the magnetic field in it's complex form + """ + return inner(np.array([w.n[1], -w.n[0]]), w.H) + + +@Functional(dtype=np.complex64) +def voltage_form(w): + """ + What this does is it takes the normal vector to the boundary and rotates it 90deg + Then takes the inner product with the electric field in it's complex form + """ + + return -inner(np.array([w.n[1], -w.n[0]]), w.E) + + +conductor = "metal_sig_interface" +line = "path_integral_right" +mode = modes[0] + +p0 = calculate_scalar_product(mode.basis, np.conjugate(mode.E), mode.basis, np.conjugate(mode.H)) +print(p0) +# The conjugate is due to femwell internal definition as conj(E) x H +# The factor of 2 is to account for the fact that we are working with half the field only + +(ht, ht_basis), (hz, hz_basis) = mode.basis.split(mode.H) +facet_basis = ht_basis.boundary(facets=mesh.boundaries[conductor]) +i0 = current_form.assemble(facet_basis, H=facet_basis.interpolate(ht)) + +(et, et_basis), (ez, ez_basis) = mode.basis.split(mode.E) +facet_basis = et_basis.boundary(facets=mesh.boundaries[line]) +v0 = voltage_form.assemble(facet_basis, E=facet_basis.interpolate(et)) + + +v0 = p0 / i0.conj() +print( + f"PI model : |Z0| = {np.abs(p0/np.abs(i0)**2):.2f} Ohm | angle(Z0) = {np.angle(p0/np.abs(i0)**2)/np.pi:.2f} pi radians" +) +print( + f"PV model : |Z0| = {np.abs(np.abs(v0)**2/p0.conj()):.2f} Ohm | angle(Z0) = {np.angle(np.abs(v0)**2/p0.conj())/np.pi:.2f} pi radians" +) +print( + f"VI model : |Z0| = {np.abs(v0/i0):.2f} Ohm | angle(Z0) = {np.angle(v0/i0)/np.pi:.2f} pi radians" +) + + +# %% [markdown] +# The fact that the Z0 is now double as before follows from the conclusion above. + + +# %% +@Functional(dtype=np.complex64) +def C_form(w): + + return ( + 1 + / np.abs(w.v0) ** 2 + * ( + np.real(w.epsilon) * inner(w.et, np.conj(w.et)) + - np.real(w.mu) * inner(w.hz, np.conj(w.hz)) + ) + ) + + +@Functional(dtype=np.complex64) +def L_form(w): + + return ( + 1 + / np.abs(w.i0) ** 2 + * ( + np.real(w.mu) * inner(w.ht, np.conj(w.ht)) + - np.real(w.epsilon) * inner(w.ez, np.conj(w.ez)) + ) + ) + + +@Functional(dtype=np.complex64) +def G_form(w): + # The minus sign account for the fact that in the paper they define eps = eps_r-1j*eps_i + # whereas with python we have eps = eps_r+1j*eps_i. + return ( + -w.omega + / np.abs(w.v0) ** 2 + * ( + np.imag(w.epsilon) * inner(w.et, np.conj(w.et)) + + np.imag(w.mu) * inner(w.hz, np.conj(w.hz)) + ) + ) + + +@Functional(dtype=np.complex64) +def R_form(w): + # The minus sign account for the fact that in the paper they define eps = eps_r-1j*eps_i + # whereas with python we have eps = eps_r+1j*eps_i. + return ( + -w.omega + / np.abs(w.i0) ** 2 + * ( + np.imag(w.epsilon) * inner(w.ez, np.conj(w.ez)) + + np.imag(w.mu) * inner(w.ht, np.conj(w.ht)) + ) + ) + + +basis_t = et_basis +basis_z = ez_basis +basis_eps = basis0 + +# Be careful with the units!! +# In this case just make sure you adjust the length unit to micrometers +# everything else can stay as is. + +C = C_form.assemble( + mode.basis, + epsilon=basis_eps.interpolate(epsilon * e0.to(reg.farad / reg.micrometer).magnitude), + mu=mu0.to(reg.henry / reg.micrometer).magnitude, + i0=i0, + omega=omega[idx_freq].to(reg.second**-1).magnitude, + v0=v0, + et=basis_t.interpolate(et), + ez=basis_z.interpolate(ez), + ht=basis_t.interpolate(ht), + hz=basis_z.interpolate(hz), +) + +L = L_form.assemble( + mode.basis, + epsilon=basis_eps.interpolate(epsilon * e0.to(reg.farad / reg.micrometer).magnitude), + mu=mu0.to(reg.henry / reg.micrometer).magnitude, + i0=i0, + omega=omega[idx_freq].to(reg.second**-1).magnitude, + v0=v0, + et=basis_t.interpolate(et), + ez=basis_z.interpolate(ez), + ht=basis_t.interpolate(ht), + hz=basis_z.interpolate(hz), +) + +G = G_form.assemble( + mode.basis, + epsilon=basis_eps.interpolate(epsilon * e0.to(reg.farad / reg.micrometer).magnitude), + mu=mu0.to(reg.henry / reg.micrometer).magnitude, + i0=i0, + omega=omega[idx_freq].to(reg.second**-1).magnitude, + v0=v0, + et=basis_t.interpolate(et), + ez=basis_z.interpolate(ez), + ht=basis_t.interpolate(ht), + hz=basis_z.interpolate(hz), +) + +R = R_form.assemble( + mode.basis, + epsilon=basis_eps.interpolate(epsilon * e0.to(reg.farad / reg.micrometer).magnitude), + mu=mu0.to(reg.henry / reg.micrometer).magnitude, + i0=i0, + omega=omega[idx_freq].to(reg.second**-1).magnitude, + v0=v0, + et=basis_t.interpolate(et), + ez=basis_z.interpolate(ez), + ht=basis_t.interpolate(ht), + hz=basis_z.interpolate(hz), +) + +Z0 = np.sqrt( + (R + 1j * omega[idx_freq].to(reg.second**-1).magnitude * L) + / (G + 1j * omega[idx_freq].to(reg.second**-1).magnitude * C) +) + +gamma = np.sqrt( + (R + 1j * omega[idx_freq].to(reg.second**-1).magnitude * L) + * (G + 1j * omega[idx_freq].to(reg.second**-1).magnitude * C) +) + +print( + f"PI model : |Z0| = {np.abs(p0/np.abs(i0)**2):.2f} Ohm | angle(Z0) = {np.angle(p0/np.abs(i0)**2)/np.pi:.2f} pi radians" +) +print( + f"PV model : |Z0| = {np.abs(np.abs(v0)**2/p0.conj()):.2f} Ohm | angle(Z0) = {np.angle(np.abs(v0)**2/p0.conj())/np.pi:.2f} pi radians" +) +print( + f"VI model : |Z0| = {np.abs(v0/i0):.2f} Ohm | angle(Z0) = {np.angle(v0/i0)/np.pi:.2f} pi radians" +) +print(f"RLGC : |Z0| = {np.abs(Z0):.2f} Ohm | angle(Z0) = {np.angle(Z0)/np.pi:.2f} pi radians") + +print("=================== RLGC ==================") +print(f"R = {R.real/1e-3:.2f} mOhm/um") +print(f"L = {L.real/1e-12:.2f} picoHenry/um") +print(f"G = {G.real/1e-12:.2f} picoSiemens/um") +print(f"C = {C.real/1e-15:.2f} femtoFarad/um") + +print("================== propagation constant ===============") +print(f"beta : RLGC = {gamma.imag*1e6:.2f} 1/m | FEM = {mode.k.real*1e6:.2f} 1/m ") +print(f"alpha: RLGC = {gamma.real*1e6:.2f} 1/m | FEM = {-mode.k.imag*1e6:.2f} 1/m ") + +# %% [markdown] +# To recover the RLGC values of the entire structure you simple half the Z and double the Y. + +# %% +R = R / 2 +L = L / 2 +C = 2 * C +G = 2 * G + +Z0 = np.sqrt( + (R + 1j * omega[idx_freq].to(reg.second**-1).magnitude * L) + / (G + 1j * omega[idx_freq].to(reg.second**-1).magnitude * C) +) + +gamma = np.sqrt( + (R + 1j * omega[idx_freq].to(reg.second**-1).magnitude * L) + * (G + 1j * omega[idx_freq].to(reg.second**-1).magnitude * C) +) + +print(f"RLGC : |Z0| = {np.abs(Z0):.2f} Ohm | angle(Z0) = {np.angle(Z0)/np.pi:.2f} pi radians") + +print("=================== RLGC ==================") +print(f"R = {R.real/1e-3:.2f} mOhm/um") +print(f"L = {L.real/1e-12:.2f} picoHenry/um") +print(f"G = {G.real/1e-12:.2f} picoSiemens/um") +print(f"C = {C.real/1e-15:.2f} femtoFarad/um") + +print("================== propagation constant ===============") +print(f"beta : RLGC = {gamma.imag*1e6:.2f} 1/m | FEM = {mode.k.real*1e6:.2f} 1/m ") +print(f"alpha: RLGC = {gamma.real*1e6:.2f} 1/m | FEM = {-mode.k.imag*1e6:.2f} 1/m ") + +# %% [markdown] +# Well, having all this set up nicely, all we are left to do is to make a frequency sweep and compare the frequency dependent values with other works and check the validity of it. Let's just define some helper functions + +# %% [markdown] +# # Frequency sweep + +# %% +reg = UnitRegistry() + +use_symmetry_plane = True +mesh_res = "fine" + +# Define frequency range +freq = np.logspace(np.log10(0.05), np.log10(45), 100) * reg.GHz +n_modes = 1 +omega = 2 * np.pi * freq +print("meshing") + + +################## MESHING ################################ +## Define universal constants +mu0 = 4 * np.pi * 1e-7 * reg.henry / reg.meter # vacuum magnetic permeability +e0 = 8.854e-12 * reg.farad * reg.meter**-1 +c = 3e8 * reg.meter * reg.second**-1 # m s^-1 +e = 1.602176634e-19 * reg.coulomb # Coulombs +kb = 1.380649e-23 * reg.meter**2 * reg.kg * reg.second**-2 * reg.kelvin**-1 +T = 300 * reg.kelvin + +w_sig = 7 * reg.micrometer +sep = 10 * reg.micrometer +w_gnd = 100 * reg.micrometer +t_metal = 0.8 * reg.micrometer + +h_sub = 500 * reg.micrometer +h_air = 500 * reg.micrometer + +eps_r_sub = 13 * reg.dimensionless + +sig_metal = 6e5 * reg.siemens / reg.centimeter + +######## Define FEM simulation ########### + +symmetry_plane = box(0, -np.inf, np.inf, np.inf) +port_width = (w_sig + 2 * sep) * 3 # um +bottom_height = (w_sig + 2 * sep) * 1 +top_height = (w_sig + 2 * sep) * 1 + +########################################## + +metal_sig_box = box( + -w_sig.to(reg.micrometer).magnitude / 2, + 0, + w_sig.to(reg.micrometer).magnitude / 2, + t_metal.to(reg.micrometer).magnitude, +) + +metal_gnd_left_box = box( + max( + (-w_sig / 2 - sep - w_gnd).to(reg.micrometer).magnitude, + -(port_width / 2).to(reg.micrometer).magnitude, + ), + 0, + (-w_sig / 2 - sep).to(reg.micrometer).magnitude, + t_metal.to(reg.micrometer).magnitude, +) + +metal_gnd_right_box = box( + (w_sig / 2 + sep).to(reg.micrometer).magnitude, + 0, + min( + (w_sig / 2 + sep + w_gnd).to(reg.micrometer).magnitude, + (port_width / 2).to(reg.micrometer).magnitude, + ), + t_metal.to(reg.micrometer).magnitude, +) + +air_box = box( + (-port_width / 2).to(reg.micrometer).magnitude, + 0, + (port_width / 2).to(reg.micrometer).magnitude, + top_height.to(reg.micrometer).magnitude, +) + + +substrate_box = box( + (-port_width / 2).to(reg.micrometer).magnitude, + -(bottom_height).to(reg.micrometer).magnitude, + (port_width / 2).to(reg.micrometer).magnitude, + 0, +) + +surface = unary_union([air_box, substrate_box]) + +##Make a line for a path integral +### Right conductor +xmin_1, ymin_1, xmax_1, ymax_1 = metal_sig_box.bounds +xmin_2, ymin_2, xmax_2, ymax_2 = metal_gnd_right_box.bounds + +points = [ + ((xmax_1 + xmin_1) / 2, (ymax_1 + ymin_1) / 2), + (xmax_1, (ymax_1 + ymin_1) / 2), + (xmin_2, (ymax_2 + ymin_2) / 2), + ((xmax_2 + xmin_2) / 2, (ymax_2 + ymin_2) / 2), +] + +path_integral_right = LineString(points) + +### left conductor +xmin_1, ymin_1, xmax_1, ymax_1 = metal_sig_box.bounds +xmin_2, ymin_2, xmax_2, ymax_2 = metal_gnd_left_box.bounds + +points = [ + ((xmax_1 + xmin_1) / 2, (ymax_1 + ymin_1) / 2), + (xmin_1, (ymax_1 + ymin_1) / 2), + (xmax_2, (ymax_2 + ymin_2) / 2), + ((xmax_2 + xmin_2) / 2, (ymax_2 + ymin_2) / 2), +] + +path_integral_left = LineString(points) + +polygons = OrderedDict( + surface=LineString(surface.exterior), + metal_sig_interface=LineString( + clip_by_rect( + metal_sig_box.buffer( + min(t_metal.to(reg.um).magnitude / 20, w_sig.to(reg.um).magnitude / 10), + join_style="bevel", + ), + *surface.bounds, + ).exterior + ), + metal_gnd_left_interface=LineString( + clip_by_rect( + metal_gnd_left_box.buffer( + min(t_metal.to(reg.um).magnitude / 20, w_gnd.to(reg.um).magnitude / 10), + join_style="bevel", + ), + *surface.bounds, + ).exterior + ), + metal_gnd_right_interface=LineString( + clip_by_rect( + metal_gnd_right_box.buffer( + min(t_metal.to(reg.um).magnitude / 20, w_gnd.to(reg.um).magnitude / 10), + join_style="bevel", + ), + *surface.bounds, + ).exterior + ), + path_integral_right=path_integral_right, + path_integral_left=path_integral_left, + metal_sig=metal_sig_box, + metal_gnd_left=metal_gnd_left_box, + metal_gnd_right=metal_gnd_right_box, + air=air_box, + substrate=substrate_box, +) + + +if use_symmetry_plane: + keys_to_pop = [] + for key, poly in polygons.items(): + if poly.intersects(symmetry_plane) and not symmetry_plane.contains(poly): + poly_tmp = clip_by_rect(poly, *symmetry_plane.bounds) + + if type(poly_tmp) == MultiLineString: + polygons[key] = linemerge(poly_tmp) + + elif poly_tmp.is_empty: + keys_to_pop.append(key) + else: + polygons[key] = poly_tmp + elif not poly.intersects(symmetry_plane): + keys_to_pop.append(key) + + for key in keys_to_pop: + polygons.pop(key) + +# Add the boundary polygons so that you can set custom boundary conditions +surf_bounds = polygons["surface"].bounds + +left = LineString([(surf_bounds[0], surf_bounds[1]), (surf_bounds[0], surf_bounds[3])]) + +bottom = LineString([(surf_bounds[0], surf_bounds[1]), (surf_bounds[2], surf_bounds[1])]) + +right = LineString([(surf_bounds[2], surf_bounds[1]), (surf_bounds[2], surf_bounds[3])]) + +top = LineString([(surf_bounds[0], surf_bounds[3]), (surf_bounds[2], surf_bounds[3])]) + +polygons["left"] = left +polygons["bottom"] = bottom +polygons["right"] = right +polygons["top"] = top + +polygons.move_to_end("top", last=False) +polygons.move_to_end("right", last=False) +polygons.move_to_end("bottom", last=False) +polygons.move_to_end("left", last=False) + +polygons.pop("surface") +# print(polygons.keys()) +if mesh_res == "fine": + resolutions = dict( + surface={"resolution": 100, "distance": 1}, + metal_sig_interface={"resolution": 0.1, "distance": 10}, + metal_gnd_interface={"resolution": 0.5, "distance": 10}, + path_integral_right={"resolution": 0.2, "distance": 10}, + path_integral_left={"resolution": 0.2, "distance": 10}, + metal_sig={"resolution": 0.1, "distance": 0.1, "SizeMax": 0.1}, + metal_gnd_left={"resolution": 0.5, "distance": 0.2, "SizeMax": 0.5}, + metal_gnd_right={"resolution": 0.5, "distance": 0.2, "SizeMax": 0.5}, + air={"resolution": 10, "distance": 0.1}, + substrate={"resolution": 10, "distance": 1}, + ) +else: + resolutions = dict( + surface={"resolution": 100, "distance": 1}, + metal_sig_interface={"resolution": 0.5, "distance": 10}, + metal_gnd_interface={"resolution": 0.5, "distance": 10}, + path_integral_right={"resolution": 0.5, "distance": 10}, + path_integral_left={"resolution": 0.5, "distance": 10}, + metal_sig={"resolution": 0.5, "distance": 0.1, "SizeMax": 0.1}, + metal_gnd_left={"resolution": 0.5, "distance": 0.2, "SizeMax": 0.5}, + metal_gnd_right={"resolution": 0.5, "distance": 0.2, "SizeMax": 0.5}, + air={"resolution": 20, "distance": 0.1}, + substrate={"resolution": 20, "distance": 1}, + ) + +# fig = plt.figure() +# ax = fig.add_subplot(111) + +# for key, polygon in polygons.items(): +# if type(polygon) is not LineString: +# x,y = polygon.exterior.xy +# ax.plot(np.asarray(x),np.asarray(y), color = 'pink') +# else: +# ax.plot(*polygon.xy) +# Mesh it +mesh = from_meshio( + mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=100, verbose=True) +) + + +print(mesh.nelements) +manager = enlighten.get_manager() +bar = manager.counter(total=len(freq), desc="Ticks", unit="ticks") + +neff_all = np.zeros(len(freq), dtype=complex) +atten = np.zeros(len(freq), dtype=float) +z0 = np.zeros(len(freq), dtype=complex) +R = np.zeros(len(freq), dtype=complex) +L = np.zeros(len(freq), dtype=complex) +G = np.zeros(len(freq), dtype=complex) +C = np.zeros(len(freq), dtype=complex) +z0_RLGC = np.zeros(len(freq), dtype=complex) + +for idx_freq in range(len(freq)): + # break + basis0 = Basis(mesh, ElementTriP0(), intorder=4) # Define the basis for the FEM + epsilon = basis0.ones(dtype=complex) + + epsilon[basis0.get_dofs(elements="air")] = 1 + epsilon[basis0.get_dofs(elements="substrate")] = eps_r_sub + epsilon[basis0.get_dofs(elements="metal_sig")] = ( + 1 - 1j * (sig_metal / omega / e0)[idx_freq].to(reg.dimensionless).magnitude + ) + # epsilon[basis0.get_dofs(elements = 'metal_gnd_left')] = (1 - 1j*(sig_metal/omega/e0)[idx_freq].to(reg.dimensionless).magnitude) + epsilon[basis0.get_dofs(elements="metal_gnd_right")] = ( + 1 - 1j * (sig_metal / omega / e0)[idx_freq].to(reg.dimensionless).magnitude + ) + # break + modes = compute_modes( + basis0, + epsilon, + wavelength=(c / freq[idx_freq]).to(reg.micrometer).magnitude, + num_modes=n_modes, + metallic_boundaries=False, + ) + + neff_all[idx_freq] = modes.n_effs + atten[idx_freq] = ( + 20 + / np.log(10) + * (modes.n_effs.imag * 2 * np.pi * freq[idx_freq] / c).to(reg.centimeter**-1).magnitude + ) # dB/cm + + ###### CALCULATE Z0 ############ + conductor = "metal_sig_interface" + line = "path_integral_right" + mode = modes[0] + + p0 = calculate_scalar_product( + mode.basis, np.conjugate(mode.E), mode.basis, np.conjugate(mode.H) + ) # The conjugate is due to femwell internal definition as conj(E) x H + + (ht, ht_basis), (hz, hz_basis) = mode.basis.split(mode.H) + (et, et_basis), (ez, ez_basis) = mode.basis.split(mode.E) + + facet_basis = ht_basis.boundary(facets=mesh.boundaries[conductor]) + i0 = current_form.assemble(facet_basis, H=facet_basis.interpolate(ht)) + + v0 = p0 / i0.conj() + + z0_tmp = v0 / i0 + + ########### CALCULATE RLGC PARAMETERS ############# + basis_t = et_basis + basis_z = ez_basis + basis_eps = basis0 + + C_tmp = C_form.assemble( + mode.basis, + epsilon=basis_eps.interpolate(epsilon * e0.to(reg.farad / reg.micrometer).magnitude), + mu=mu0.to(reg.henry / reg.micrometer).magnitude, + i0=i0, + omega=omega[idx_freq].to(reg.second**-1).magnitude, + v0=v0, + et=basis_t.interpolate(et), + ez=basis_z.interpolate(ez), + ht=basis_t.interpolate(ht), + hz=basis_z.interpolate(hz), + ) + + L_tmp = L_form.assemble( + mode.basis, + epsilon=basis_eps.interpolate(epsilon * e0.to(reg.farad / reg.micrometer).magnitude), + mu=mu0.to(reg.henry / reg.micrometer).magnitude, + i0=i0, + omega=omega[idx_freq].to(reg.second**-1).magnitude, + v0=v0, + et=basis_t.interpolate(et), + ez=basis_z.interpolate(ez), + ht=basis_t.interpolate(ht), + hz=basis_z.interpolate(hz), + ) + + G_tmp = G_form.assemble( + mode.basis, + epsilon=basis_eps.interpolate(epsilon * e0.to(reg.farad / reg.micrometer).magnitude), + mu=mu0.to(reg.henry / reg.micrometer).magnitude, + i0=i0, + omega=omega[idx_freq].to(reg.second**-1).magnitude, + v0=v0, + et=basis_t.interpolate(et), + ez=basis_z.interpolate(ez), + ht=basis_t.interpolate(ht), + hz=basis_z.interpolate(hz), + ) + + R_tmp = R_form.assemble( + mode.basis, + epsilon=basis_eps.interpolate(epsilon * e0.to(reg.farad / reg.micrometer).magnitude), + mu=mu0.to(reg.henry / reg.micrometer).magnitude, + i0=i0, + omega=omega[idx_freq].to(reg.second**-1).magnitude, + v0=v0, + et=basis_t.interpolate(et), + ez=basis_z.interpolate(ez), + ht=basis_t.interpolate(ht), + hz=basis_z.interpolate(hz), + ) + + z0_RLGC_tmp = np.sqrt( + (R_tmp + 1j * omega[idx_freq].to(reg.second**-1).magnitude * L_tmp) + / (G_tmp + 1j * omega[idx_freq].to(reg.second**-1).magnitude * C_tmp) + ) + + ## Save the parameters for the full device + + z0[idx_freq] = z0_tmp / 2 + z0_RLGC[idx_freq] = z0_RLGC_tmp / 2 + R[idx_freq] = R_tmp / 2 + L[idx_freq] = L_tmp / 2 + G[idx_freq] = G_tmp * 2 + C[idx_freq] = C_tmp * 2 + + bar.update() + +mesh.draw() + +# fig = basis0.plot(epsilon.real, colorbar = True) +# fig.axes.set_axis_on() + +# fig = basis0.plot(epsilon.imag, colorbar = True) +# fig.axes.set_axis_on() + +# %% +##DATA FROM PAPER +data_nu = np.asarray( + [ + [-1.317, 7.780], + [-1.218, 7.327], + [-1.106, 6.474], + [-0.983, 5.995], + [-0.800, 5.160], + [-0.647, 4.574], + [-0.499, 4.094], + [-0.270, 3.632], + [0.000, 3.188], + [0.311, 2.993], + [0.622, 2.940], + [0.958, 2.833], + [1.254, 2.833], + [1.559, 2.833], + ] +) + + +data_alpha = np.asarray( + [ + [-1.287, 0.728], + [-1.060, 0.950], + [-0.780, 1.181], + [-0.520, 1.465], + [-0.270, 1.785], + [-0.025, 2.087], + [0.204, 2.353], + [0.392, 2.567], + [0.632, 2.886], + [0.846, 3.117], + [1.070, 3.739], + [1.310, 4.680], + [1.478, 5.426], + [1.590, 5.977], + ] +) + +## DATA FROM TU DELFT SIM +freq_matlab = np.loadtxt("support/freq_delft.txt") # Frequency sampling from their simulator +ReZ0_delft = np.loadtxt("support/ReZ0_delft.txt") # Real part of the impedance +ImZ0_delft = np.loadtxt("support/ImZ0_delft.txt") # Imaginary part of the impedance +loss_delft = np.loadtxt("support/loss_delft.txt") # The loss of the mode +ReK_delft = np.loadtxt("support/ReK_delft.txt") # The real part of the wavevector + +# %% +fig = plt.figure() +ax = fig.add_subplot(111) +ax1 = ax.twinx() + +ax.plot(freq, neff_all, color="blue", label="FEMWELL") +ax.plot(freq_matlab, ReK_delft, color="blue", linestyle="dotted", label="TU Delft simulator") + +ax.plot( + 10 ** data_nu[:, 0], + data_nu[:, 1], + marker="^", + color="blue", + linestyle="dashed", + label="Tuncer et al. 1992", +) +ax.set_xscale("log") + +ax.legend(loc="upper center") +ax1.plot(freq, -atten, color="red") +ax1.plot(10 ** data_alpha[:, 0], data_alpha[:, 1], marker="^", color="red", linestyle="dashed") +ax1.plot(freq_matlab, loss_delft, color="red", linestyle="dotted") + +ax.set_xlabel("Frequency (GHz)") +ax.set_ylabel("Microwave index") +ax1.set_ylabel("Attenuation (dB/cm)") + +ax.arrow(0.1, 5, -0.05, 0, head_width=0.1, head_length=0.01, color="blue") +ax1.arrow(10, 3, 10, 0, head_width=0.1, head_length=5, color="red") + +# fig.savefig('Tuncer_et_al_1992_results.png', dpi = 400) + +# %% +fig = plt.figure() +ax = fig.add_subplot(111) + +ax.plot(freq, z0.real, label="FEMWELL - real", color="blue") +ax.plot(freq, z0.imag, label="FEMWELL - imag", color="red") +ax.plot(freq_matlab, ReZ0_delft, label="TU Delft sim - real", color="blue", linestyle="dashed") +ax.plot(freq_matlab, ImZ0_delft, label="TU Delft sim - imag", color="red", linestyle="dashed") + +ax.set_xlabel("Frequency (GHz)") +ax.set_ylabel("Z0 (Ohm)") + +ax.legend() + +ax.set_xscale("log") + +# fig.savefig('z0_delft.png', dpi = 400) + +# %% [markdown] +# The mismatch between the above curves can be explained due to TU Delft not being able to correctly estimate the propagation constant of the field, hence it will not correctly predict the characteristic impedance of the device. + +# %% +fig = plt.figure() +gs = GridSpec(2, 2) + +ax_C = fig.add_subplot(gs[0, 0]) +ax_G = fig.add_subplot(gs[0, 1]) +ax_L = fig.add_subplot(gs[1, 0]) +ax_R = fig.add_subplot(gs[1, 1]) + +for data, ax, label in zip( + [R / 1e-3, L / 1e-12, G / 1e-12, C / 1e-15], + [ax_R, ax_L, ax_G, ax_C], + ["kOhm/m", "uHenry/m", "uS/m", "nF/m"], +): + + ax.plot(freq, data, color="blue") + ax.set_xlabel("Frequency (GHz)") + ax.set_ylabel(label) + + ax.set_xscale("log") +fig.tight_layout() diff --git a/docs/electronics/examples/RF_CPW_transmission_line_tutorial/support/ImZ0_delft.txt b/docs/electronics/examples/RF_CPW_transmission_line_tutorial/support/ImZ0_delft.txt new file mode 100644 index 00000000..35c095d9 --- /dev/null +++ b/docs/electronics/examples/RF_CPW_transmission_line_tutorial/support/ImZ0_delft.txt @@ -0,0 +1,154 @@ +-1.157582306070029876e+00 +-1.182509787468440710e+00 +-1.208503087405807364e+00 +-1.235588034615594744e+00 +-1.263768931049162880e+00 +-1.293081539981793204e+00 +-1.323557946029030763e+00 +-1.355237558367177408e+00 +-1.388156740638433417e+00 +-1.422352034109460650e+00 +-1.457869823339989601e+00 +-1.494749245914654079e+00 +-1.533033300300297386e+00 +-1.572766357985726504e+00 +-1.613993723586580664e+00 +-1.656763166485924810e+00 +-1.701122794174423580e+00 +-1.747122599141999855e+00 +-1.794789653763162907e+00 +-1.844248824821063559e+00 +-1.895486117128052816e+00 +-1.948580531250056591e+00 +-2.003590981533915460e+00 +-2.060578434045229113e+00 +-2.119605996687136340e+00 +-2.180739009610220158e+00 +-2.244045135972282612e+00 +-2.309594452303071588e+00 +-2.377459540113373571e+00 +-2.447715577172612456e+00 +-2.520440428911534525e+00 +-2.595714741236326173e+00 +-2.673622032428730755e+00 +-2.754248786457532905e+00 +-2.837692988237921377e+00 +-2.924030895904208283e+00 +-3.013367263682980024e+00 +-3.105800566244479644e+00 +-3.201433503654893542e+00 +-3.300372292394565310e+00 +-3.402726754313660340e+00 +-3.508610402074802792e+00 +-3.618140522625397093e+00 +-3.731438255995030318e+00 +-3.848628670830473997e+00 +-3.969840833390715584e+00 +-4.095207873418067912e+00 +-4.224867035474326293e+00 +-4.358959731850598374e+00 +-4.497631578450308254e+00 +-4.641032423297423826e+00 +-4.789316363120608600e+00 +-4.942641745335079051e+00 +-5.101171159102533714e+00 +-5.265071405826591722e+00 +-5.434513460729482048e+00 +-5.609672409365172641e+00 +-5.790727371568686443e+00 +-5.977861403959793485e+00 +-6.171261384512746417e+00 +-6.371117874860692254e+00 +-6.577624967227535535e+00 +-6.790980107930280951e+00 +-7.011383902694460879e+00 +-7.239039907174851685e+00 +-7.474154394664134671e+00 +-7.716936112842202711e+00 +-7.967596024459555615e+00 +-8.226347038316509597e+00 +-8.493368069995694825e+00 +-8.769116357913990356e+00 +-9.053438798896502604e+00 +-9.346717875333606429e+00 +-9.649171930456445168e+00 +-9.961019532631356199e+00 +-1.028247920218426970e+01 +-1.061376914683452988e+01 +-1.095510702912872070e+01 +-1.130670973467646157e+01 +-1.166879318020991718e+01 +-1.204157213952655248e+01 +-1.242526010783475243e+01 +-1.282006918049309618e+01 +-1.322621001267457430e+01 +-1.364389172951264229e+01 +-1.407332198224052178e+01 +-1.451470696280678574e+01 +-1.496825148615143775e+01 +-1.543415912025571579e+01 +-1.591263233910674124e+01 +-1.640387271644796385e+01 +-1.690808117257032350e+01 +-1.742545822669816147e+01 +-1.795620431149833607e+01 +-1.850052007621407668e+01 +-1.905860677360249511e+01 +-1.963066660621158022e+01 +-2.021690314087108931e+01 +-2.081752169844617129e+01 +-2.143272980464999478e+01 +-2.206273759372864518e+01 +-2.271147997335348379e+01 +-2.337191642684674520e+01 +-2.404780852821738435e+01 +-2.473938119978217287e+01 +-2.544686414835368282e+01 +-2.617049227889619090e+01 +-2.691050610924286701e+01 +-2.766715216563436996e+01 +-2.844068336383317330e+01 +-2.923135937872348933e+01 +-3.003944700212987939e+01 +-3.086522046999905200e+01 +-3.170896180387996921e+01 +-3.257096109076084645e+01 +-3.345298297477481242e+01 +-3.435244289293947162e+01 +-3.527108288336980735e+01 +-3.620922845887920261e+01 +-3.716721422092346216e+01 +-3.814538437971605589e+01 +-3.914409291622719422e+01 +-4.016370393316071841e+01 +-4.120459168435123587e+01 +-4.226714085453937031e+01 +-4.335174669023051308e+01 +-4.445881518217638018e+01 +-4.558876328031279002e+01 +-4.674201895488760528e+01 +-4.791902143439809691e+01 +-4.912022133323333151e+01 +-5.034608082443121901e+01 +-5.159707373969254007e+01 +-5.287368580837868848e+01 +-5.417641474669940749e+01 +-5.550577042013694751e+01 +-5.686227506700453915e+01 +-5.824646334609037268e+01 +-5.965888261773633872e+01 +-6.110009302410490051e+01 +-6.257066775745295217e+01 +-6.407119312415588297e+01 +-6.560226883635576201e+01 +-6.716450814885490672e+01 +-6.875853803288995891e+01 +-7.038499945678564984e+01 +-7.204454750450928202e+01 +-7.373785166009164982e+01 +-7.546559597686869836e+01 +-7.722847942938810206e+01 +-7.902721592877759349e+01 +-8.086253483971837852e+01 +-8.273518104013440677e+01 +-8.464591516344316346e+01 diff --git a/docs/electronics/examples/RF_CPW_transmission_line_tutorial/support/ReK_delft.txt b/docs/electronics/examples/RF_CPW_transmission_line_tutorial/support/ReK_delft.txt new file mode 100644 index 00000000..8e328607 --- /dev/null +++ b/docs/electronics/examples/RF_CPW_transmission_line_tutorial/support/ReK_delft.txt @@ -0,0 +1,154 @@ +2.691485913925359075e+00 +2.692151926868994405e+00 +2.692826441486253053e+00 +2.693519551526634803e+00 +2.694226310540131841e+00 +2.694946955246150910e+00 +2.695675282891197622e+00 +2.696424752955969062e+00 +2.697189305790852032e+00 +2.697959524713574364e+00 +2.698755824544570636e+00 +2.699568884977982641e+00 +2.700399341492370819e+00 +2.701247861759427327e+00 +2.702114915749732837e+00 +2.703001684617745148e+00 +2.703908695441340271e+00 +2.704836729769424686e+00 +2.705772287716126545e+00 +2.706744804619915179e+00 +2.707740856974918842e+00 +2.708761340324274070e+00 +2.709807181408131260e+00 +2.710879339718949677e+00 +2.711978809320642725e+00 +2.713106620885914033e+00 +2.714263843944311194e+00 +2.715451589337029681e+00 +2.716671011877543673e+00 +2.717923313219855252e+00 +2.719209744938591022e+00 +2.720531611827278873e+00 +2.721890275423061389e+00 +2.723287157767730449e+00 +2.724698058981393167e+00 +2.726175871080242175e+00 +2.727696572665775587e+00 +2.729261870536163759e+00 +2.730873554773528333e+00 +2.732533504191862139e+00 +2.734243692115226665e+00 +2.736006192499397027e+00 +2.737823186408880449e+00 +2.739696968859475668e+00 +2.741629956034230542e+00 +2.743624692877720594e+00 +2.745683861070061571e+00 +2.747810287377837835e+00 +2.750006952374342006e+00 +2.752276999515936406e+00 +2.754623744555268594e+00 +2.757050685265315160e+00 +2.759561511440975057e+00 +2.762160115137277216e+00 +2.764850601095249161e+00 +2.767637297298348731e+00 +2.770524765594182259e+00 +2.773517812308225139e+00 +2.776621498768665486e+00 +2.779841151654383946e+00 +2.783182373071895643e+00 +2.786651050261772955e+00 +2.790253364831053684e+00 +2.793995801405521728e+00 +2.797885155594628248e+00 +2.801928541162580721e+00 +2.806133396301537353e+00 +2.810507488907405094e+00 +2.815058920765070916e+00 +2.819748064870253312e+00 +2.824679731147743844e+00 +2.829815066340673368e+00 +2.835163524812530866e+00 +2.840734896916673247e+00 +2.846539305188366065e+00 +2.852587199232630599e+00 +2.858889349231084065e+00 +2.865456838110005844e+00 +2.872301052431808888e+00 +2.879433672091868068e+00 +2.886866658921865447e+00 +2.894612244318917949e+00 +2.902682916045249417e+00 +2.911091404287765361e+00 +2.919850667325169624e+00 +2.928973876668762077e+00 +2.938474402164450883e+00 +2.948365797056039295e+00 +2.958661783256278266e+00 +2.969376236999468510e+00 +2.980523175052351537e+00 +2.992116741651275724e+00 +3.004171196321832049e+00 +3.016700902723012501e+00 +3.029720318641507593e+00 +3.043243987243854143e+00 +3.057286529674496656e+00 +3.071862639068193346e+00 +3.086987076024456211e+00 +3.102674665571999224e+00 +3.118940295631263293e+00 +3.135798921196399469e+00 +3.153265549201734075e+00 +3.171355265617371799e+00 +3.190083223819020652e+00 +3.209464653944571033e+00 +3.229514869597755578e+00 +3.250249275684259587e+00 +3.271683377290186101e+00 +3.293832789509170134e+00 +3.316713248122880930e+00 +3.340340621039453950e+00 +3.364730920395968550e+00 +3.389900315233658201e+00 +3.415865144658416952e+00 +3.442491762046922865e+00 +3.470096006880941886e+00 +3.498545824337362387e+00 +3.527858371546021665e+00 +3.558051046071218693e+00 +3.589141500381201055e+00 +3.621147656735406617e+00 +3.654087721093635999e+00 +3.687980198273299148e+00 +3.722843906066733766e+00 +3.758697989651845095e+00 +3.795561935826521349e+00 +3.833455587118438146e+00 +3.872399155761411471e+00 +3.912413237533613586e+00 +3.953518825456088681e+00 +3.995737323353187431e+00 +4.039090559279421733e+00 +4.083600798819269073e+00 +4.129290758268699513e+00 +4.176183617708841389e+00 +4.224303033983763278e+00 +4.273673153595286323e+00 +4.324318625529041782e+00 +4.376264614026507793e+00 +4.429536811318531164e+00 +4.484161450336041632e+00 +4.540165317414378343e+00 +4.597575765007191784e+00 +4.656420724426808100e+00 +4.716728718627032357e+00 +4.778528875045179092e+00 +4.841850938519123737e+00 +4.906725284295848866e+00 +4.973182931147079877e+00 +5.041255554607754874e+00 +5.110975500352574308e+00 +5.182375797725897293e+00 +5.255490173439739898e+00 diff --git a/docs/electronics/examples/RF_CPW_transmission_line_tutorial/support/ReZ0_delft.txt b/docs/electronics/examples/RF_CPW_transmission_line_tutorial/support/ReZ0_delft.txt new file mode 100644 index 00000000..c8d467b5 --- /dev/null +++ b/docs/electronics/examples/RF_CPW_transmission_line_tutorial/support/ReZ0_delft.txt @@ -0,0 +1,154 @@ +6.102385574926012879e+01 +6.104551690520960250e+01 +6.106660759490944201e+01 +6.108739654224218185e+01 +6.110778956018397423e+01 +6.112799942262114428e+01 +6.114826459881543030e+01 +6.116833493516440967e+01 +6.118842705717329267e+01 +6.120883484463112723e+01 +6.122911859359372500e+01 +6.124956324215520453e+01 +6.127020755411031416e+01 +6.129108784380473907e+01 +6.131224480027973556e+01 +6.133369835457791908e+01 +6.135548631509072237e+01 +6.137763908111366362e+01 +6.140021069839198731e+01 +6.142350545854514365e+01 +6.144692799167865616e+01 +6.147083094078191579e+01 +6.149524254550814106e+01 +6.152019106516704028e+01 +6.154570490011303718e+01 +6.157181271102260922e+01 +6.159854353110115710e+01 +6.162592687634925426e+01 +6.165399285405608509e+01 +6.168277227071557434e+01 +6.171229673690123718e+01 +6.174259877754125370e+01 +6.177371193951613293e+01 +6.180567090335308933e+01 +6.183912856860459328e+01 +6.187288786537260421e+01 +6.190760526989012646e+01 +6.194332090354865983e+01 +6.198007690987424922e+01 +6.201791731655830375e+01 +6.205688817727072859e+01 +6.209703772158894708e+01 +6.213841651481875061e+01 +6.218107761885220697e+01 +6.222507676685935252e+01 +6.227047254194565795e+01 +6.231732656393388226e+01 +6.236570368718663815e+01 +6.241567220340925815e+01 +6.246730405089373761e+01 +6.252067503472971310e+01 +6.257586505101465235e+01 +6.263295831715276307e+01 +6.269204360832888057e+01 +6.275321449515720929e+01 +6.281656959017098529e+01 +6.288221278669858094e+01 +6.295025349988044638e+01 +6.302080691851418237e+01 +6.309399422984090933e+01 +6.316994285179792001e+01 +6.324878666050017273e+01 +6.333066619439188116e+01 +6.341572885473902943e+01 +6.350412908553607849e+01 +6.359602853695201929e+01 +6.369159620807236166e+01 +6.379100856715874102e+01 +6.389444964897840862e+01 +6.400214211666279596e+01 +6.411531591714583556e+01 +6.423202558526564587e+01 +6.435357677811569488e+01 +6.448019178989409284e+01 +6.461210044097907712e+01 +6.474953996029825021e+01 +6.489275483799201538e+01 +6.504199664771238076e+01 +6.519752384104459964e+01 +6.535960150978695538e+01 +6.552850113473118654e+01 +6.570450029603257747e+01 +6.588788238115849083e+01 +6.607893626697959633e+01 +6.627795599223888701e+01 +6.648524042411435175e+01 +6.670109292007090573e+01 +6.692582099174367727e+01 +6.715973596512434085e+01 +6.740315266591022692e+01 +6.765638909882801499e+01 +6.791976615457770095e+01 +6.819360734285707792e+01 +6.847823852945849410e+01 +6.877398771243461795e+01 +6.908118483458150649e+01 +6.940016160583618898e+01 +6.973125137354213621e+01 +7.007478901776045177e+01 +7.043111089055379637e+01 +7.080055477563900013e+01 +7.118271584792721285e+01 +7.157945595311556986e+01 +7.199034528708668290e+01 +7.241572788516086234e+01 +7.285594952744975217e+01 +7.331135788694763278e+01 +7.378230272139265367e+01 +7.426913605136654439e+01 +7.477221240680763970e+01 +7.529188906327368613e+01 +7.582852628226599734e+01 +7.638248762111402357e+01 +7.695414018646303589e+01 +7.754385496229198793e+01 +7.815548911234613172e+01 +7.878248590484847114e+01 +7.942868489861962189e+01 +8.009447578825587755e+01 +8.078025360638723384e+01 +8.148641923749076454e+01 +8.221337959125223449e+01 +8.296154814674756039e+01 +8.373134509454395413e+01 +8.452319774390967666e+01 +8.533754085711237281e+01 +8.617481694247821622e+01 +8.703547657732279674e+01 +8.791997878040336900e+01 +8.882879123025921331e+01 +8.976239064182811944e+01 +9.072126308019812768e+01 +9.170590424863537748e+01 +9.271681974733117215e+01 +9.375452543905950620e+01 +9.481954772076542781e+01 +9.591242381350990343e+01 +9.703370205806540127e+01 +9.818394221262735755e+01 +9.936371571524450985e+01 +1.005736060313317495e+02 +1.018142088790417148e+02 +1.030861325579987806e+02 +1.043899981619856447e+02 +1.057264400949929239e+02 +1.070961059601505525e+02 +1.084996572906935484e+02 +1.099377695632773566e+02 +1.114111325680043336e+02 +1.129204506863547408e+02 +1.144664432939833745e+02 +1.160498448614285962e+02 +1.176714055121576905e+02 +1.193318911123788268e+02 diff --git a/docs/electronics/examples/RF_CPW_transmission_line_tutorial/support/freq_delft.txt b/docs/electronics/examples/RF_CPW_transmission_line_tutorial/support/freq_delft.txt new file mode 100644 index 00000000..1deac253 --- /dev/null +++ b/docs/electronics/examples/RF_CPW_transmission_line_tutorial/support/freq_delft.txt @@ -0,0 +1,154 @@ +5.000000000000000000e+01 +4.800000000000000000e+01 +4.607999999999999829e+01 +4.423679999999999524e+01 +4.246732799999999486e+01 +4.076863487999999336e+01 +3.913788948479999874e+01 +3.757237390540799282e+01 +3.606947894919166231e+01 +3.462669979122399866e+01 +3.324163179957503900e+01 +3.191196652759204255e+01 +3.063548786648835431e+01 +2.941006835182882284e+01 +2.823366561775566552e+01 +2.710431899304543890e+01 +2.602014623332361865e+01 +2.497934038399067447e+01 +2.398016676863104735e+01 +2.302096009788580488e+01 +2.210012169397037241e+01 +2.121611682621155381e+01 +2.036747215316308868e+01 +1.955277326703656726e+01 +1.877066233635510173e+01 +1.801983584290089624e+01 +1.729904240918486025e+01 +1.660708071281746712e+01 +1.594279748430476928e+01 +1.530508558493257709e+01 +1.469288216153527316e+01 +1.410516687507386280e+01 +1.354096020007090573e+01 +1.299932179206806993e+01 +1.247934892038534649e+01 +1.198017496356993306e+01 +1.150096796502713481e+01 +1.104092924642604778e+01 +1.059929207656900552e+01 +1.017532039350624729e+01 +9.768307577765995120e+00 +9.377575274655354676e+00 +9.002472263669142549e+00 +8.642373373122374858e+00 +8.296678438197480077e+00 +7.964811300669580874e+00 +7.646218848642796218e+00 +7.340370094697085790e+00 +7.046755290909201719e+00 +6.764885079272832336e+00 +6.494289676101919540e+00 +6.234518089057842438e+00 +5.985137365495527995e+00 +5.745731870875706804e+00 +5.515902596040677786e+00 +5.295266492199051100e+00 +5.083455832511088346e+00 +4.880117599210645807e+00 +4.684912895242218767e+00 +4.497516379432529021e+00 +4.317615724255229104e+00 +4.144911095285019442e+00 +3.979114651473618647e+00 +3.819950065414673190e+00 +3.667152062798086476e+00 +3.520465980286162644e+00 +3.379647341074716760e+00 +3.244461447431727930e+00 +3.114682989534458546e+00 +2.990095669953080293e+00 +2.870491843154956868e+00 +2.755672169428758256e+00 +2.645445282651607943e+00 +2.539627471345543341e+00 +2.438042372491721732e+00 +2.340520677592052845e+00 +2.246899850488370376e+00 +2.157023856468835987e+00 +2.070742902210081837e+00 +1.987913186121678555e+00 +1.908396658676811519e+00 +1.832060792329738952e+00 +1.758778360636549465e+00 +1.688427226211087318e+00 +1.620890137162643851e+00 +1.556054531676137831e+00 +1.493812350409092105e+00 +1.434059856392728483e+00 +1.376697462137019379e+00 +1.321629563651538719e+00 +1.268764381105477179e+00 +1.218013805861257737e+00 +1.169293253626807472e+00 +1.122521523481735173e+00 +1.077620662542465579e+00 +1.034515836040767134e+00 +9.931352025991362353e-01 +9.534097944951709369e-01 +9.152734027153639484e-01 +8.786624666067494038e-01 +8.435159679424794055e-01 +8.097753292247802381e-01 +7.773843160557889487e-01 +7.462889434135573374e-01 +7.164373856770150928e-01 +6.877798902499344669e-01 +6.602686946399370926e-01 +6.338579468543396045e-01 +6.085036289801658560e-01 +5.841634838209591729e-01 +5.607969444681207838e-01 +5.383650666893959880e-01 +5.168304640218202017e-01 +4.961572454609473604e-01 +4.763109556425094149e-01 +4.572585174168090072e-01 +4.389681767201366047e-01 +4.214094496513311472e-01 +4.045530716652778747e-01 +3.883709487986667930e-01 +3.728361108467200813e-01 +3.579226664128512536e-01 +3.436057597563372257e-01 +3.298615293660837056e-01 +3.166670681914403596e-01 +3.040003854637826630e-01 +2.918403700452313898e-01 +2.801667552434221609e-01 +2.689600850336852922e-01 +2.582016816323378094e-01 +2.478736143670443315e-01 +2.379586697923625205e-01 +2.284403230006679941e-01 +2.193027100806413110e-01 +2.105306016774156075e-01 +2.021093776103189965e-01 +1.940250025059062355e-01 +1.862640024056699539e-01 +1.788134423094431824e-01 +1.716609046170654296e-01 +1.647944684323827957e-01 +1.582026896950874950e-01 +1.518745821072839675e-01 +1.457995988229926210e-01 +1.399676148700728928e-01 +1.343689102752699815e-01 +1.289941538642591923e-01 +1.238343877096887979e-01 +1.188810122013012477e-01 +1.141257717132491967e-01 +1.095607408447192288e-01 +1.051783112109304641e-01 +1.009711787624932161e-01 +9.693233161199349412e-02 diff --git a/docs/electronics/examples/RF_CPW_transmission_line_tutorial/support/loss_delft.txt b/docs/electronics/examples/RF_CPW_transmission_line_tutorial/support/loss_delft.txt new file mode 100644 index 00000000..8b0efa22 --- /dev/null +++ b/docs/electronics/examples/RF_CPW_transmission_line_tutorial/support/loss_delft.txt @@ -0,0 +1,154 @@ +4.739199576599296471e+00 +4.647012687544633813e+00 +4.558218242030294576e+00 +4.472737258033426144e+00 +4.390419373715454654e+00 +4.311136387493473165e+00 +4.234734544214971308e+00 +4.161141192738623928e+00 +4.090205657526658101e+00 +4.021769103523255673e+00 +3.955792925471715460e+00 +3.892124139392544091e+00 +3.830653812404662073e+00 +3.771277277256577243e+00 +3.713892704675008538e+00 +3.658407302229667479e+00 +3.604727969944267141e+00 +3.552767548089233340e+00 +3.502400346173617507e+00 +3.453634805545053421e+00 +3.406351834173134030e+00 +3.360480679837588802e+00 +3.315954354723815722e+00 +3.272709479418081990e+00 +3.230686122419777373e+00 +3.189827637236431812e+00 +3.150080499055894379e+00 +3.111394142623802761e+00 +3.073720802634242855e+00 +3.037015357664051507e+00 +3.001235178442220874e+00 +2.966339981041535179e+00 +2.932291685406332782e+00 +2.899054279484842844e+00 +2.866546038454464984e+00 +2.834831497199962413e+00 +2.803830874363740566e+00 +2.773515189757886912e+00 +2.743856937864472822e+00 +2.714829985828722592e+00 +2.686409476587388045e+00 +2.658571736908973637e+00 +2.631294190115978537e+00 +2.604555273259819170e+00 +2.578334358526289449e+00 +2.552611678680725582e+00 +2.527368256232317645e+00 +2.502585836535490049e+00 +2.478246824038599971e+00 +2.454334222209606065e+00 +2.430831576665106741e+00 +2.407722921576634700e+00 +2.384992729311537385e+00 +2.362625863314483876e+00 +2.340607534264264711e+00 +2.318923259566478379e+00 +2.297558826265203802e+00 +2.276500257474750821e+00 +2.255733782445382030e+00 +2.235245810384023368e+00 +2.215022908151307313e+00 +2.195051781950370984e+00 +2.175319263109452006e+00 +2.155812298040462149e+00 +2.136517942428950878e+00 +2.117423359678053174e+00 +2.098515823590638707e+00 +2.079782725230922669e+00 +2.061211583859940077e+00 +2.042750476276841365e+00 +2.024467011801589944e+00 +2.006308978241504892e+00 +1.988264591073880316e+00 +1.970322304903539212e+00 +1.952470841442672578e+00 +1.934699221859771168e+00 +1.916996800625285324e+00 +1.899353300357262020e+00 +1.881758847150196923e+00 +1.864204005858222546e+00 +1.846679814799558450e+00 +1.829177819353787093e+00 +1.811690103957840137e+00 +1.794209321857834105e+00 +1.776728722600446453e+00 +1.759242176092288279e+00 +1.741744193592286383e+00 +1.724229944860523389e+00 +1.706695271340720188e+00 +1.689136695144437095e+00 +1.671551423691038929e+00 +1.653937349923479516e+00 +1.636293048086679081e+00 +1.618617765120738738e+00 +1.600911407783964124e+00 +1.583174525679961864e+00 +1.565408290417325876e+00 +1.547614471178377693e+00 +1.529795407015015707e+00 +1.511953976223662233e+00 +1.494093563177530504e+00 +1.476218023552321013e+00 +1.458331644867988164e+00 +1.440439112061591365e+00 +1.422545465328050485e+00 +1.404656061048842997e+00 +1.386776531806721113e+00 +1.368912746607069852e+00 +1.351070771637138179e+00 +1.333256831865709202e+00 +1.315477273754578347e+00 +1.297738529320441536e+00 +1.280047081752527127e+00 +1.262409432757672700e+00 +1.244832071771888815e+00 +1.227267033177987354e+00 +1.209830330469727988e+00 +1.192473026334426311e+00 +1.175201293632006649e+00 +1.158021169741950240e+00 +1.140938538244850031e+00 +1.123959112684016048e+00 +1.107088422285148477e+00 +1.090331799674984836e+00 +1.073694370408444732e+00 +1.057181044281076510e+00 +1.040796508308895030e+00 +1.024545221277803453e+00 +1.008431409758662189e+00 +9.924590654823298719e-01 +9.766319439688035953e-01 +9.609535643054692988e-01 +9.454272099713439959e-01 +9.300559306069960952e-01 +9.148425446331904798e-01 +8.997896426252972590e-01 +8.848995913548900472e-01 +8.701745384146549744e-01 +8.556164173476724644e-01 +8.412269532070437705e-01 +8.270076684771089370e-01 +8.129598892923184383e-01 +7.990847518949766304e-01 +7.853832092779218232e-01 +7.718560379627628087e-01 +7.585038448691042712e-01 +7.453270742342021915e-01 +7.323260145468937221e-01 +7.195008054632886507e-01 +7.068514446755712610e-01 +6.943777947084666113e-01 +6.820795896213168730e-01 +6.699564415964458508e-01 +6.580078473973850484e-01 diff --git a/femwell/mesh/mesh.py b/femwell/mesh/mesh.py index 5f880c05..4b056c87 100644 --- a/femwell/mesh/mesh.py +++ b/femwell/mesh/mesh.py @@ -14,14 +14,14 @@ Point, Polygon, ) -from shapely.ops import linemerge, polygonize, split, unary_union +from shapely.ops import linemerge, polygonize, snap, split, unary_union from femwell.mesh.meshtracker import MeshTracker initial_settings = np.seterr() # remove when shapely updated to more recent geos -def break_line_(line, other_line): +def break_line_(line, other_line, snap_tol): np.seterr(invalid="ignore") intersections = line.intersection(other_line) np.seterr(**initial_settings) @@ -31,7 +31,7 @@ def break_line_(line, other_line): ): # if type == "", intersection.type != 'Point': if intersection.geom_type == "Point": - line = linemerge(split(line, intersection)) + line = snap(line, intersection, snap_tol) else: new_coords_start, new_coords_end = intersection.boundary.geoms line = linemerge(split(line, new_coords_start)) @@ -73,7 +73,9 @@ def mesh_from_Dict( rings = [ LineString(list(object.exterior.coords)) for object in listpoly - if not (object.geom_type in ["Point", "LineString"] or object.is_empty) + if not ( + object.geom_type in ["Point", "LineString", "MultiLineString"] or object.is_empty + ) ] union = unary_union(rings) @@ -204,6 +206,7 @@ def mesh_from_OrderedDict( gmsh.option.setNumber("Mesh.Algorithm", gmsh_algorithm) model = geometry + meshtracker = MeshTracker(model=model) # Break up shapes in order so that plane is tiled with non-overlapping layers, overriding shapes according to an order shapes_tiled_dict = OrderedDict() @@ -250,7 +253,7 @@ def mesh_from_OrderedDict( else second_shape ) first_exterior_line = break_line_( - first_exterior_line, second_exterior_line + first_exterior_line, second_exterior_line, meshtracker.atol ) # Second line interiors for second_interior_line in ( @@ -260,7 +263,7 @@ def mesh_from_OrderedDict( ): second_interior_line = LineString(second_interior_line) first_exterior_line = break_line_( - first_exterior_line, second_interior_line + first_exterior_line, second_interior_line, meshtracker.atol ) # First line interiors if first_shape.geom_type in ["Polygon", "MultiPolygon"]: @@ -286,7 +289,7 @@ def mesh_from_OrderedDict( else second_shape ) first_interior_line = break_line_( - first_interior_line, second_exterior_line + first_interior_line, second_exterior_line, meshtracker.atol ) # Interiors for second_interior_line in ( @@ -301,7 +304,9 @@ def mesh_from_OrderedDict( ) np.seterr(**initial_settings) first_interior_line = break_line_( - first_interior_line, second_interior_line + first_interior_line, + second_interior_line, + meshtracker.atol, ) first_shape_interiors.append(first_interior_line) if first_shape.geom_type in ["Polygon", "MultiPolygon"]: @@ -318,7 +323,7 @@ def mesh_from_OrderedDict( ) # Add lines, reusing line segments - meshtracker = MeshTracker(model=model) + for line_name, line in lines_broken_dict.items(): meshtracker.add_get_xy_line(line, line_name) diff --git a/pyproject.toml b/pyproject.toml index 44bd60c3..39cb7621 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -41,6 +41,8 @@ flake8 = "*" [tool.poetry.group.docs.dependencies] python = ">=3.10" tqdm = "*" +enlighten = "*" +pint = ">0.20.1" # sphinx-book-theme = "*" jupytext = "*" # myst-parser = "*"