From 5beeb76d8759b2b3f64987158f74f10676c3bbe2 Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Tue, 15 Apr 2025 15:48:33 -0400 Subject: [PATCH 01/28] feat[FieldData]: export to ZBF --- CHANGELOG.md | 1 + tests/test_data/test_monitor_data.py | 226 ++++++++++++++++++++++++- tidy3d/components/data/monitor_data.py | 124 +++++++++++++- 3 files changed, 349 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 47db9ec7e..333788c47 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Methods `EMEExplicitGrid.from_structures` and `EMECompositeGrid.from_structure_groups` to place EME cell boundaries at structure bounds. - 'ModeSimulation' now supports 'PermittivityMonitor'. - Classmethod `from_frequency_range` in `GaussianPulse` for generating a pulse whose amplitude in the frequency_range [fmin, fmax] is maximized, which is particularly useful for running broadband simulations. +- `FieldData` supports exporting E fields to a Zemax Beam File (ZBF) with `FieldData.to_zbf()`. ### Changed - Performance enhancement for adjoint gradient calculations by optimizing field interpolation. diff --git a/tests/test_data/test_monitor_data.py b/tests/test_data/test_monitor_data.py index 9ca9a7051..f0974cddf 100644 --- a/tests/test_data/test_monitor_data.py +++ b/tests/test_data/test_monitor_data.py @@ -1,5 +1,7 @@ """Tests tidy3d/components/data/monitor_data.py""" +from struct import unpack + import matplotlib.pyplot as plt import numpy as np import pydantic.v1 as pydantic @@ -23,7 +25,7 @@ ) from tidy3d.exceptions import DataError -from ..utils import AssertLogLevel +from ..utils import AssertLogLevel, get_spatial_coords_dict from .test_data_arrays import ( AUX_FIELD_TIME_MONITOR, DIFFRACTION_MONITOR, @@ -865,3 +867,225 @@ def test_no_nans(): ) with pytest.raises(pydantic.ValidationError): td.CustomMedium(eps_dataset=eps_dataset_nan) + + +class TestZBF: + """Tests exporting field data to a zbf file""" + + simulation = SIM + monitor = FIELD_MONITOR_2D + + def read_beam_file(self, filename: str) -> tuple: + """Read in a Zemax Beam file + Adapted from https://community.zemax.com/code-exchange-10/python-class-to-open-zbf-file-4845 + + Parameters + ---------- + beamfilename : string + the filename of the beam file to read + + Returns + ------- + version : integer + the file format version number + n : 2-tuple, (nx, ny) + the number of samples in the x and y directions + ispol : boolean + is the beam polarized? + units : integer (0 or 1 or 2 or 3) + the units of the beam, 0 = mm, 1 = cm, 2 = in, 3 for m + d : 2-tuple, (dx, dy) + the x and y grid spacing + zposition : 2-tuple, (zpositionx, zpositiony) + the x and y z position of the beam + rayleigh : 2-tuple, (rayleighx, rayleighy) + the x and y rayleigh ranges of the beam + waist : 2-tuple, (waistx, waisty) + the x and y waists of the beam + lamda : double + the wavelength of the beam + index : double + the index of refraction in the current medium + receiver_eff : double + the receiver efficiency. Zero if fiber coupling is not computed + system_eff : double + the system efficiency. Zero if fiber coupling is not computed. + grid_pos : 2-tuple of lists, (x_matrix, y_matrix) + lists of x and y positions of the grid defining the beam + efield : 4-tuple of 2D lists, (Ex_real, Ex_imag, Ey_real, Ey_imag) + a tuple containing two dimensional lists with the real and + imaginary parts of the x and y polarizations of the beam + + """ + # will need to get rid of during automation. + # beamfilename = self.open_zbf() + f = open(filename, "rb") + # zemax version number + version = unpack("i", f.read(4))[0] + nx = unpack("i", f.read(4))[0] + ny = unpack("i", f.read(4))[0] + ispol = unpack("i", f.read(4))[0] + units = unpack("i", f.read(4))[0] + + if units == 0: + unit = "mm" + elif units == 1: + unit = "cm" + elif units == 2: + unit = "in" + elif units == 3: + unit = "m" + + f.read(16) + dx = unpack("d", f.read(8))[0] + dy = unpack("d", f.read(8))[0] + + # if version==1: + zposition_x = unpack("d", f.read(8))[0] + rayleigh_x = unpack("d", f.read(8))[0] + waist_x = unpack("d", f.read(8))[0] + # f.read(16) + zposition_y = unpack("d", f.read(8))[0] + rayleigh_y = unpack("d", f.read(8))[0] + waist_y = unpack("d", f.read(8))[0] + lamda = unpack("d", f.read(8))[0] + index = unpack("d", f.read(8))[0] + receiver_eff = unpack("d", f.read(8))[0] + system_eff = unpack("d", f.read(8))[0] + f.read(64) # 8 empty doubles + + rawx = [unpack("d", f.read(8))[0] for _ in range(2 * nx * ny)] + if ispol: + rawy = [unpack("d", f.read(8))[0] for _ in range(2 * nx * ny)] + f.close() + + # Get Image array + # Ex Real Value for point 1; + # Ex Imaginary Value for point 1; + # Ex Real Value for point 2; + # Ex Imaginary Value for point 2, etc + + Ex_real = np.asarray(rawx[0::2]).reshape(nx, ny) + Ex_imag = np.asarray(rawx[1::2]).reshape(nx, ny) + if ispol: + Ey_real = np.asarray(rawy[0::2]).reshape(nx, ny) + Ey_imag = np.asarray(rawy[1::2]).reshape(nx, ny) + else: + Ey_real = np.zeros((nx, ny)) + Ey_imag = np.zeros((nx, ny)) + + return ( + version, + nx, + ny, + ispol, + unit, + dx, + dy, + zposition_x, + zposition_y, + rayleigh_x, + rayleigh_y, + waist_x, + waist_y, + lamda, + index, + receiver_eff, + system_eff, + Ex_real, + Ex_imag, + Ey_real, + Ey_imag, + ) + + def make_data(self, coords: dict) -> td.components.data.data_array.DataArray: + """make a random DataArray out of supplied coordinates and data_type.""" + # get shape of the data array + data_shape = [len(coords[k]) for k in td.ScalarFieldDataArray._dims] + + # generate random data + np.random.seed(1) + data = np.random.random(data_shape) + 1j * np.random.random(data_shape) + + # return as a td.ScalarFieldDataArray + data_array = td.ScalarFieldDataArray(data, coords=coords) + return data_array + + @pytest.fixture(scope="class") + def field_data(self) -> td.FieldData: + """make a random FieldData from a FieldMonitor.""" + # get the grid for the FieldData + grid = self.simulation.discretize_monitor(self.monitor) + + # generate random complex field components + field_cmps = {} + for field_name in self.monitor.fields: + coords = get_spatial_coords_dict(self.simulation, self.monitor, field_name) + coords["f"] = list(self.monitor.freqs) + field_cmps[field_name] = self.make_data(coords=coords) + + return td.FieldData( + monitor=self.monitor, + symmetry=(0, 0, 0), + symmetry_center=self.simulation.center, + grid_expanded=grid, + **field_cmps, + ) + + @pytest.mark.parametrize("background_index", [1, 2, 3]) + @pytest.mark.parametrize("freq", list(monitor.freqs) + [None]) + @pytest.mark.parametrize("num_samples", [2, 4, 8, 16, 32, 64, None]) + def test_zbfread(self, tmp_path, field_data, background_index, freq, num_samples): + zbf_filename = tmp_path / "testzbf.zbf" + + # write to zbf and then load it back in + ex, ey = field_data.to_zbf( + fname=zbf_filename, + background_refractive_index=background_index, + freq=freq, + num_samples=num_samples, + ) + zbfdata = self.read_beam_file(zbf_filename) + ( + _version, + nx, + ny, + _ispol, + _unit, + _dx, + _dy, + _zposition_x, + _zposition_y, + _rayleigh_x, + _rayleigh_y, + _waist_x, + _waist_y, + wl, + index, + _receiver_eff, + _system_eff, + Ex_real, + Ex_imag, + Ey_real, + Ey_imag, + ) = zbfdata + + assert index == background_index + + if freq is not None: + assert np.isclose(wl * 1e3, td.C_0 / freq) + else: + assert np.isclose(wl * 1e3, td.C_0 / self.monitor.freqs[len(self.monitor.freqs) // 2]) + + if num_samples is not None: + assert nx == num_samples + assert ny == num_samples + else: + assert nx == 32 + assert ny == 64 + + # check that fields are close + assert np.allclose(ex.real.values, Ex_real) + assert np.allclose(ex.imag.values, Ex_imag) + assert np.allclose(ey.real.values, Ey_real) + assert np.allclose(ey.imag.values, Ey_imag) diff --git a/tidy3d/components/data/monitor_data.py b/tidy3d/components/data/monitor_data.py index 41e3bd2e2..ab2f6be9c 100644 --- a/tidy3d/components/data/monitor_data.py +++ b/tidy3d/components/data/monitor_data.py @@ -2,10 +2,11 @@ from __future__ import annotations +import struct import warnings from abc import ABC from math import isclose -from typing import Any, Callable, Dict, List, Tuple, Union +from typing import Any, Callable, Dict, List, Optional, Tuple, Union import autograd.numpy as np import pydantic.v1 as pd @@ -1224,6 +1225,127 @@ def make_adjoint_sources( return sources + def to_zbf( + self, + fname: str, + background_refractive_index: float, + freq: Optional[float] = None, + mode_index: Optional[int] = 0, + num_samples: Optional[int] = None, + ) -> Tuple[ScalarFieldDataArray, ScalarFieldDataArray]: + """For a 2D monitor, export the fields to a Zemax Beam File (zbf). + + The mode area is used to approximate the beam waist, which is only valid + if the beam profile approximates a Gaussian beam. + + Parameters + ---------- + fname : str + Full path to the .zbf file to be written. + background_refractive_index : float + Refractive index of the medium surrounding the monitor. + freq : float + Field frequency selection. If ``None``, the central one is used. + mode_index : int + Mode index selection. + num_samples: int + Number of samples to use. Must be a power of 2 + + Returns + ------- + Tuple[:class:`ScalarFieldDataArray`,:class:`ScalarFieldDataArray`] + The two E field components being exported to zbf. + """ + log.warning( + "Exporting to ZBF is currently an experimental feature. If any issues are encountered, please contact Flexcompute support." + ) + + # Mode area calculation ensures all E components are present + mode_area = self.mode_area + dim1, dim2 = self._tangential_dims + + # Using file-local coordinates x, y for the tangential components + e_x = self._tangential_fields["E" + dim1] + e_y = self._tangential_fields["E" + dim2] + x = e_x.coords[dim1].values + y = e_x.coords[dim2].values + + # Use the central frequency if freq is not specified + if freq is None: + freq = e_x.coords["f"].values[e_x.coords["f"].size // 2] + + mode_area = mode_area.sel(f=freq, drop=True) + e_x = e_x.sel(f=freq, drop=True) + e_y = e_y.sel(f=freq, drop=True) + + if "mode_index" in e_x.coords: + mode_area = mode_area.isel(mode_index=mode_index, drop=True) + e_x = e_x.isel(mode_index=mode_index, drop=True) + e_y = e_y.isel(mode_index=mode_index, drop=True) + + # Header info + version = 1 + polarized = 1 + units = 0 # "mm": 0, "cm": 1, "in": 2, "m": 3 + lda = C_0 / freq * 1e-3 + + # Pilot (reference) beam waist: use the mode area to approximate the expected value + w_x = (mode_area.item() / np.pi) ** 0.5 * 1e-3 + w_y = w_x + + # Pilot beam Rayleigh distance (ignored on input) + r_x = 0 + r_y = 0 + + # Pilot beam z position w.r.t. the waist + z_x = 0 * 1e-3 + z_y = 0 * 1e-3 + + # Number of samples: n_x = 2 ** a, a >= 5 and a <= 13 (32-bit), or a <= 14 (64-bit) + if num_samples is not None: + # custom user defined number of samples + # Raise an error if the number of samples is not a power of 2 + if (num_samples & (num_samples - 1)) != 0: + raise ValueError("num_samples must be a power of 2.") + n_x = num_samples + n_y = num_samples + else: + n_x = 2 ** min(13, max(5, int(np.log2(x.size) + 1))) + n_y = 2 ** min(13, max(5, int(np.log2(y.size) + 1))) + + # Interpolating coordinates + x = np.linspace(x.min(), x.max(), n_x) + y = np.linspace(y.min(), y.max(), n_y) + + # Interpolate fields + coords = {dim1: x, dim2: y} + e_x = e_x.interp(coords, assume_sorted=True) + e_y = e_y.interp(coords, assume_sorted=True) + + # Sampling distance + d_x = np.mean(np.diff(x)) * 1e-3 + d_y = np.mean(np.diff(y)) * 1e-3 + + # Receiver and system efficiencies + rec_efficiency = 0 # zero if fiber coupling is not computed + sys_efficiency = 0 # zero if fiber coupling is not computed + + with open(fname, "wb") as fout: + fout.write(struct.pack("=5I", version, n_x, n_y, polarized, units)) + fout.write(struct.pack("=4I", 0, 0, 0, 0)) # unused values + fout.write(struct.pack("=8d", d_x, d_y, z_x, r_x, w_x, z_y, r_y, w_y)) + fout.write( + struct.pack("=4d", lda, background_refractive_index, rec_efficiency, sys_efficiency) + ) + fout.write(struct.pack("=8d", 0, 0, 0, 0, 0, 0, 0, 0)) # unused values + for e in (e_x, e_y): + e_flat = e.values.flatten(order="C") + # Interweave real and imaginary parts + e_values = np.ravel(np.column_stack((e_flat.real, e_flat.imag))) + fout.write(struct.pack(f"={2 * n_x*n_y}d", *e_values)) + + return e_x, e_y + class FieldTimeData(FieldTimeDataset, ElectromagneticFieldData): """ From d9e2fd2290e0ae1ca50d7e1a858e36c9ef81e4fa Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Fri, 25 Apr 2025 13:49:26 -0400 Subject: [PATCH 02/28] use little endian when writing zbf --- tidy3d/components/data/monitor_data.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tidy3d/components/data/monitor_data.py b/tidy3d/components/data/monitor_data.py index ab2f6be9c..3c427c0b9 100644 --- a/tidy3d/components/data/monitor_data.py +++ b/tidy3d/components/data/monitor_data.py @@ -1331,18 +1331,18 @@ def to_zbf( sys_efficiency = 0 # zero if fiber coupling is not computed with open(fname, "wb") as fout: - fout.write(struct.pack("=5I", version, n_x, n_y, polarized, units)) - fout.write(struct.pack("=4I", 0, 0, 0, 0)) # unused values - fout.write(struct.pack("=8d", d_x, d_y, z_x, r_x, w_x, z_y, r_y, w_y)) + fout.write(struct.pack("<5I", version, n_x, n_y, polarized, units)) + fout.write(struct.pack("<4I", 0, 0, 0, 0)) # unused values + fout.write(struct.pack("<8d", d_x, d_y, z_x, r_x, w_x, z_y, r_y, w_y)) fout.write( - struct.pack("=4d", lda, background_refractive_index, rec_efficiency, sys_efficiency) + struct.pack("<4d", lda, background_refractive_index, rec_efficiency, sys_efficiency) ) - fout.write(struct.pack("=8d", 0, 0, 0, 0, 0, 0, 0, 0)) # unused values + fout.write(struct.pack("<8d", 0, 0, 0, 0, 0, 0, 0, 0)) # unused values for e in (e_x, e_y): e_flat = e.values.flatten(order="C") # Interweave real and imaginary parts e_values = np.ravel(np.column_stack((e_flat.real, e_flat.imag))) - fout.write(struct.pack(f"={2 * n_x*n_y}d", *e_values)) + fout.write(struct.pack(f"<{2 * n_x*n_y}d", *e_values)) return e_x, e_y From 8297bfb9cfa902bd3262545ad52a0fee3d6cb1c6 Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Fri, 25 Apr 2025 13:49:57 -0400 Subject: [PATCH 03/28] add read_zbf() method to td.components.data.utils --- tidy3d/components/data/utils.py | 115 ++++++++++++++++++++++++++++++++ 1 file changed, 115 insertions(+) diff --git a/tidy3d/components/data/utils.py b/tidy3d/components/data/utils.py index 9e64765fd..01e353945 100644 --- a/tidy3d/components/data/utils.py +++ b/tidy3d/components/data/utils.py @@ -2,6 +2,7 @@ from __future__ import annotations +from struct import unpack from typing import Union import numpy as np @@ -79,3 +80,117 @@ def _check_same_coordinates( return False return True + + +def read_zbf(filename: str): + """Reads a Zemax Beam File (zbf) + + Parameters + ---------- + filename : string + the filename of the beam file to read + + Returns + ------- + version : integer + the file format version number + n : 2-tuple, (nx, ny) + the number of samples in the x and y directions + ispol : boolean + is the beam polarized? + units : integer (0 or 1 or 2 or 3) + the units of the beam, 0 = mm, 1 = cm, 2 = in, 3 for m + d : 2-tuple, (dx, dy) + the x and y grid spacing + zposition : 2-tuple, (zpositionx, zpositiony) + the x and y z position of the beam + rayleigh : 2-tuple, (rayleighx, rayleighy) + the x and y rayleigh ranges of the beam + waist : 2-tuple, (waistx, waisty) + the x and y waists of the beam + lamda : double + the wavelength of the beam + index : double + the index of refraction in the current medium + receiver_eff : double + the receiver efficiency. Zero if fiber coupling is not computed + system_eff : double + the system efficiency. Zero if fiber coupling is not computed. + grid_pos : 2-tuple of lists, (x_matrix, y_matrix) + lists of x and y positions of the grid defining the beam + efield : 4-tuple of 2D lists, (Ex_real, Ex_imag, Ey_real, Ey_imag) + a tuple containing two dimensional lists with the real and + imaginary parts of the x and y polarizations of the beam + """ + + # Read the zbf file + with open(filename, "rb") as f: + # Load the header + version, nx, ny, ispol, units = unpack("<5I", f.read(20)) + f.read(16) # unused values + ( + dx, + dy, + zposition_x, + rayleigh_x, + waist_x, + zposition_y, + rayleigh_y, + waist_y, + wavelength, + background_refractive_index, + receiver_eff, + system_eff, + ) = unpack("<12d", f.read(96)) + f.read(64) # unused values + + # read E field + nsamps = 2 * nx * ny + rawx = [unpack(" Date: Fri, 25 Apr 2025 14:43:37 -0400 Subject: [PATCH 04/28] move read_zbf to its own file because i was running into circular dependencies --- tidy3d/components/data/utils.py | 115 ------------------------------- tidy3d/components/data/zbf.py | 118 ++++++++++++++++++++++++++++++++ 2 files changed, 118 insertions(+), 115 deletions(-) create mode 100644 tidy3d/components/data/zbf.py diff --git a/tidy3d/components/data/utils.py b/tidy3d/components/data/utils.py index 01e353945..9e64765fd 100644 --- a/tidy3d/components/data/utils.py +++ b/tidy3d/components/data/utils.py @@ -2,7 +2,6 @@ from __future__ import annotations -from struct import unpack from typing import Union import numpy as np @@ -80,117 +79,3 @@ def _check_same_coordinates( return False return True - - -def read_zbf(filename: str): - """Reads a Zemax Beam File (zbf) - - Parameters - ---------- - filename : string - the filename of the beam file to read - - Returns - ------- - version : integer - the file format version number - n : 2-tuple, (nx, ny) - the number of samples in the x and y directions - ispol : boolean - is the beam polarized? - units : integer (0 or 1 or 2 or 3) - the units of the beam, 0 = mm, 1 = cm, 2 = in, 3 for m - d : 2-tuple, (dx, dy) - the x and y grid spacing - zposition : 2-tuple, (zpositionx, zpositiony) - the x and y z position of the beam - rayleigh : 2-tuple, (rayleighx, rayleighy) - the x and y rayleigh ranges of the beam - waist : 2-tuple, (waistx, waisty) - the x and y waists of the beam - lamda : double - the wavelength of the beam - index : double - the index of refraction in the current medium - receiver_eff : double - the receiver efficiency. Zero if fiber coupling is not computed - system_eff : double - the system efficiency. Zero if fiber coupling is not computed. - grid_pos : 2-tuple of lists, (x_matrix, y_matrix) - lists of x and y positions of the grid defining the beam - efield : 4-tuple of 2D lists, (Ex_real, Ex_imag, Ey_real, Ey_imag) - a tuple containing two dimensional lists with the real and - imaginary parts of the x and y polarizations of the beam - """ - - # Read the zbf file - with open(filename, "rb") as f: - # Load the header - version, nx, ny, ispol, units = unpack("<5I", f.read(20)) - f.read(16) # unused values - ( - dx, - dy, - zposition_x, - rayleigh_x, - waist_x, - zposition_y, - rayleigh_y, - waist_y, - wavelength, - background_refractive_index, - receiver_eff, - system_eff, - ) = unpack("<12d", f.read(96)) - f.read(64) # unused values - - # read E field - nsamps = 2 * nx * ny - rawx = [unpack(" Date: Fri, 25 Apr 2025 16:54:47 -0400 Subject: [PATCH 05/28] fielddataset: remove kwargs, update docstring --- tidy3d/components/data/dataset.py | 87 +++++++++++++++++++++++++++++-- tidy3d/components/types.py | 4 ++ tidy3d/constants.py | 1 + 3 files changed, 89 insertions(+), 3 deletions(-) diff --git a/tidy3d/components/data/dataset.py b/tidy3d/components/data/dataset.py index 1e23ceda2..87f78f0ef 100644 --- a/tidy3d/components/data/dataset.py +++ b/tidy3d/components/data/dataset.py @@ -3,17 +3,17 @@ from __future__ import annotations from abc import ABC, abstractmethod -from typing import Any, Callable, Dict, Optional, Union +from typing import Any, Callable, Dict, Optional, Union, get_args import numpy as np import pydantic.v1 as pd import xarray as xr -from ...constants import PICOSECOND_PER_NANOMETER_PER_KILOMETER +from ...constants import C_0, PICOSECOND_PER_NANOMETER_PER_KILOMETER, UnitScaling from ...exceptions import DataError from ...log import log from ..base import Tidy3dBaseModel -from ..types import Axis +from ..types import Axis, xyz from .data_array import ( DataArray, EMEScalarFieldDataArray, @@ -28,6 +28,7 @@ TimeDataArray, TriangleMeshDataArray, ) +from .zbf import read_zbf DEFAULT_MAX_SAMPLES_PER_STEP = 10_000 DEFAULT_MAX_CELLS_PER_STEP = 10_000 @@ -257,6 +258,86 @@ class FieldDataset(ElectromagneticFieldDataset): description="Spatial distribution of the z-component of the magnetic field.", ) + def from_zbf(filename: str, dim1: xyz, dim2: xyz): + """Creates a :class:`.FieldDataset` from a Zemax Beam File (zbf) + + Parameters + ---------- + filename: str + Specification of the source time-dependence. + dim1: xyz + Tangential field component to map the x-axis of the zbf data to. + eg. `dim1 = "z"` sets `FieldDataset.Ez` to Ex of the zbf data. + dim2: xyz + Tangential field component to map the y-axis of the zbf data to. + + Returns + ------- + :class:`.FieldDataset` + A :class:`.FieldDataset` object with two tangential E field components populated + by zbf data. + """ + log.warning( + "`FieldDataset.from_zbf()` is currently an experimental feature." + "If any issues are encountered, please contact Flexcompute support https://www.flexcompute.com/tidy3d/technical-support/" + ) + + if dim1 not in get_args(xyz): + raise ValueError(f'`dim1` = {dim1} not allowed, must be one of `"x"`, `"y"`, or `"z"`.') + if dim2 not in get_args(xyz): + raise ValueError(f'`dim2` = {dim2} not allowed, must be one of `"x"`, `"y"`, or `"z"`.') + if dim1 == dim2: + raise ValueError("`dim1` and `dim2` must be different.") + + # get the third dimension + dim3 = list(set(get_args(xyz)) - {dim1, dim2})[0] + dims = {"x": 0, "y": 1, "z": 2} + dim2expand = dims[dim3] # this is for expanding E field arrays + + # load zbf data + zbfdata = read_zbf(filename) + + # Grab E fields, dimensions, wavelength + edim1 = zbfdata["Ex"] + edim2 = zbfdata["Ey"] + n1 = zbfdata["nx"] + n2 = zbfdata["ny"] + d1 = zbfdata["dx"] / UnitScaling[zbfdata["unit"]] + d2 = zbfdata["dy"] / UnitScaling[zbfdata["unit"]] + wavelength = zbfdata["wavelength"] / UnitScaling[zbfdata["unit"]] + + # make scalar field data arrays + len1 = d1 * (n1 - 1) + len2 = d2 * (n2 - 1) + coords1 = np.linspace(-len1 / 2, len1 / 2, n1) + coords2 = np.linspace(-len2 / 2, len2 / 2, n2) + f = [C_0 / wavelength] + Edim1 = ScalarFieldDataArray( + np.expand_dims(edim1, axis=(dim2expand, 3)), + coords={ + dim1: coords1, + dim2: coords2, + dim3: [0], + "f": f, + }, + ) + Edim2 = ScalarFieldDataArray( + np.expand_dims(edim2, axis=(dim2expand, 3)), + coords={ + dim1: coords1, + dim2: coords2, + dim3: [0], + "f": f, + }, + ) + + return FieldDataset( + **{ + f"E{dim1}": Edim1, + f"E{dim2}": Edim2, + } + ) + class FieldTimeDataset(ElectromagneticFieldDataset): """Dataset storing a collection of the scalar components of E and H fields in the time domain diff --git a/tidy3d/components/types.py b/tidy3d/components/types.py index 4f4c807be..9ff928d50 100644 --- a/tidy3d/components/types.py +++ b/tidy3d/components/types.py @@ -252,3 +252,7 @@ def __modify_schema__(cls, field_schema): """ lumped elements""" LumpDistType = Literal["off", "laterally_only", "on"] + +""" dataset """ + +xyz = Literal["x", "y", "z"] diff --git a/tidy3d/constants.py b/tidy3d/constants.py index b3439e83a..2610887e1 100644 --- a/tidy3d/constants.py +++ b/tidy3d/constants.py @@ -238,6 +238,7 @@ "mm": 1e-3, "cm": 1e-4, "m": 1e-6, + "in": 1.0 / 25400, } ) """Immutable dictionary for converting a unit specification to a scaling factor.""" From a0a5bba7bf25219271c3f10856e9be9b60a092b5 Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Fri, 25 Apr 2025 17:11:37 -0400 Subject: [PATCH 06/28] read_zbf: update docstring --- tidy3d/components/data/zbf.py | 62 ++++++++++++++++++++--------------- 1 file changed, 36 insertions(+), 26 deletions(-) diff --git a/tidy3d/components/data/zbf.py b/tidy3d/components/data/zbf.py index dbc807376..56707b5d5 100644 --- a/tidy3d/components/data/zbf.py +++ b/tidy3d/components/data/zbf.py @@ -15,35 +15,45 @@ def read_zbf(filename: str): Returns ------- - version : integer - the file format version number - n : 2-tuple, (nx, ny) - the number of samples in the x and y directions + A dictionary with the following fields: + version : int + File format version number. + nx : int + Number of samples in the x direction. + ny : int + Number of samples in the y direction. ispol : boolean - is the beam polarized? - units : integer (0 or 1 or 2 or 3) - the units of the beam, 0 = mm, 1 = cm, 2 = in, 3 for m - d : 2-tuple, (dx, dy) - the x and y grid spacing - zposition : 2-tuple, (zpositionx, zpositiony) - the x and y z position of the beam - rayleigh : 2-tuple, (rayleighx, rayleighy) - the x and y rayleigh ranges of the beam - waist : 2-tuple, (waistx, waisty) - the x and y waists of the beam - lamda : double - the wavelength of the beam - index : double - the index of refraction in the current medium + True if the beam is polarized. + unit : str + Spatial units, either "mm", "cm", "in", or "m". + dx : float + x grid spacing + dy : float + y grid spacing + zposition_x : float + Pilot beam z position with respect to the waist in x. + zposition_y : float + Pilot beam z position with respect to the waist in y. + rayleigh_x : float + Pilot beam Rayleigh distance in x. + rayleigh_y : float + Pilot beam Rayleigh distance in y. + waist_x : float + Pilot beam waist in x. + waist_y : float + Pilot beam waist in y. + wavelength : float + The wavelength of the beam. + background_refractive_index : float + The index of refraction of the material that the beam exists in. receiver_eff : double - the receiver efficiency. Zero if fiber coupling is not computed + The receiver efficiency. Zero if fiber coupling is not computed. system_eff : double - the system efficiency. Zero if fiber coupling is not computed. - grid_pos : 2-tuple of lists, (x_matrix, y_matrix) - lists of x and y positions of the grid defining the beam - efield : 4-tuple of 2D lists, (Ex_real, Ex_imag, Ey_real, Ey_imag) - a tuple containing two dimensional lists with the real and - imaginary parts of the x and y polarizations of the beam + The system efficiency. Zero if fiber coupling is not computed. + Ex : np.ndarray + Complex-valued electric field, x-component. + Ey : np.ndarray + Complex-valued electric field, x-component. """ # Read the zbf file From 6a4e1a0e9d877298cac88b132bea0af3c7564bf1 Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Mon, 28 Apr 2025 17:20:31 -0400 Subject: [PATCH 07/28] fixed up the tests to use emulated run and new read zbf function --- tests/test_data/test_monitor_data.py | 237 +++++---------------------- 1 file changed, 39 insertions(+), 198 deletions(-) diff --git a/tests/test_data/test_monitor_data.py b/tests/test_data/test_monitor_data.py index f0974cddf..4752fac85 100644 --- a/tests/test_data/test_monitor_data.py +++ b/tests/test_data/test_monitor_data.py @@ -1,7 +1,5 @@ """Tests tidy3d/components/data/monitor_data.py""" -from struct import unpack - import matplotlib.pyplot as plt import numpy as np import pydantic.v1 as pydantic @@ -23,9 +21,10 @@ ModeData, PermittivityData, ) +from tidy3d.components.data.zbf import read_zbf from tidy3d.exceptions import DataError -from ..utils import AssertLogLevel, get_spatial_coords_dict +from ..utils import AssertLogLevel, run_emulated from .test_data_arrays import ( AUX_FIELD_TIME_MONITOR, DIFFRACTION_MONITOR, @@ -872,169 +871,37 @@ def test_no_nans(): class TestZBF: """Tests exporting field data to a zbf file""" - simulation = SIM - monitor = FIELD_MONITOR_2D - - def read_beam_file(self, filename: str) -> tuple: - """Read in a Zemax Beam file - Adapted from https://community.zemax.com/code-exchange-10/python-class-to-open-zbf-file-4845 - - Parameters - ---------- - beamfilename : string - the filename of the beam file to read - - Returns - ------- - version : integer - the file format version number - n : 2-tuple, (nx, ny) - the number of samples in the x and y directions - ispol : boolean - is the beam polarized? - units : integer (0 or 1 or 2 or 3) - the units of the beam, 0 = mm, 1 = cm, 2 = in, 3 for m - d : 2-tuple, (dx, dy) - the x and y grid spacing - zposition : 2-tuple, (zpositionx, zpositiony) - the x and y z position of the beam - rayleigh : 2-tuple, (rayleighx, rayleighy) - the x and y rayleigh ranges of the beam - waist : 2-tuple, (waistx, waisty) - the x and y waists of the beam - lamda : double - the wavelength of the beam - index : double - the index of refraction in the current medium - receiver_eff : double - the receiver efficiency. Zero if fiber coupling is not computed - system_eff : double - the system efficiency. Zero if fiber coupling is not computed. - grid_pos : 2-tuple of lists, (x_matrix, y_matrix) - lists of x and y positions of the grid defining the beam - efield : 4-tuple of 2D lists, (Ex_real, Ex_imag, Ey_real, Ey_imag) - a tuple containing two dimensional lists with the real and - imaginary parts of the x and y polarizations of the beam - - """ - # will need to get rid of during automation. - # beamfilename = self.open_zbf() - f = open(filename, "rb") - # zemax version number - version = unpack("i", f.read(4))[0] - nx = unpack("i", f.read(4))[0] - ny = unpack("i", f.read(4))[0] - ispol = unpack("i", f.read(4))[0] - units = unpack("i", f.read(4))[0] - - if units == 0: - unit = "mm" - elif units == 1: - unit = "cm" - elif units == 2: - unit = "in" - elif units == 3: - unit = "m" - - f.read(16) - dx = unpack("d", f.read(8))[0] - dy = unpack("d", f.read(8))[0] - - # if version==1: - zposition_x = unpack("d", f.read(8))[0] - rayleigh_x = unpack("d", f.read(8))[0] - waist_x = unpack("d", f.read(8))[0] - # f.read(16) - zposition_y = unpack("d", f.read(8))[0] - rayleigh_y = unpack("d", f.read(8))[0] - waist_y = unpack("d", f.read(8))[0] - lamda = unpack("d", f.read(8))[0] - index = unpack("d", f.read(8))[0] - receiver_eff = unpack("d", f.read(8))[0] - system_eff = unpack("d", f.read(8))[0] - f.read(64) # 8 empty doubles - - rawx = [unpack("d", f.read(8))[0] for _ in range(2 * nx * ny)] - if ispol: - rawy = [unpack("d", f.read(8))[0] for _ in range(2 * nx * ny)] - f.close() - - # Get Image array - # Ex Real Value for point 1; - # Ex Imaginary Value for point 1; - # Ex Real Value for point 2; - # Ex Imaginary Value for point 2, etc - - Ex_real = np.asarray(rawx[0::2]).reshape(nx, ny) - Ex_imag = np.asarray(rawx[1::2]).reshape(nx, ny) - if ispol: - Ey_real = np.asarray(rawy[0::2]).reshape(nx, ny) - Ey_imag = np.asarray(rawy[1::2]).reshape(nx, ny) - else: - Ey_real = np.zeros((nx, ny)) - Ey_imag = np.zeros((nx, ny)) - - return ( - version, - nx, - ny, - ispol, - unit, - dx, - dy, - zposition_x, - zposition_y, - rayleigh_x, - rayleigh_y, - waist_x, - waist_y, - lamda, - index, - receiver_eff, - system_eff, - Ex_real, - Ex_imag, - Ey_real, - Ey_imag, - ) - - def make_data(self, coords: dict) -> td.components.data.data_array.DataArray: - """make a random DataArray out of supplied coordinates and data_type.""" - # get shape of the data array - data_shape = [len(coords[k]) for k in td.ScalarFieldDataArray._dims] - - # generate random data - np.random.seed(1) - data = np.random.random(data_shape) + 1j * np.random.random(data_shape) - - # return as a td.ScalarFieldDataArray - data_array = td.ScalarFieldDataArray(data, coords=coords) - return data_array + freq0 = td.C_0 / 0.75 + freqs = (freq0, freq0 * 1.01) @pytest.fixture(scope="class") def field_data(self) -> td.FieldData: - """make a random FieldData from a FieldMonitor.""" - # get the grid for the FieldData - grid = self.simulation.discretize_monitor(self.monitor) - - # generate random complex field components - field_cmps = {} - for field_name in self.monitor.fields: - coords = get_spatial_coords_dict(self.simulation, self.monitor, field_name) - coords["f"] = list(self.monitor.freqs) - field_cmps[field_name] = self.make_data(coords=coords) - - return td.FieldData( - monitor=self.monitor, - symmetry=(0, 0, 0), - symmetry_center=self.simulation.center, - grid_expanded=grid, - **field_cmps, + """Make random field data from an emulated simulation run.""" + source = td.PointDipole( + center=(-1.5, 0, 0), + source_time=td.GaussianPulse(freq0=self.freq0, fwidth=self.freq0 / 10.0), + polarization="Ey", + ) + monitor = td.FieldMonitor( + size=(td.inf, td.inf, 0), + freqs=self.freqs, + name="fields", + colocate=True, + ) + sim = td.Simulation( + size=(4, 3, 3), + grid_spec=td.GridSpec.auto(min_steps_per_wvl=10), + structures=[], + sources=[source], + monitors=[monitor], + run_time=120 / self.freq0, ) + simdata = run_emulated(sim) + return simdata["fields"] @pytest.mark.parametrize("background_index", [1, 2, 3]) - @pytest.mark.parametrize("freq", list(monitor.freqs) + [None]) - @pytest.mark.parametrize("num_samples", [2, 4, 8, 16, 32, 64, None]) + @pytest.mark.parametrize("freq", list(freqs) + [None]) + @pytest.mark.parametrize("num_samples", [2, 4, 8, 16, 32, 64]) def test_zbfread(self, tmp_path, field_data, background_index, freq, num_samples): zbf_filename = tmp_path / "testzbf.zbf" @@ -1045,47 +912,21 @@ def test_zbfread(self, tmp_path, field_data, background_index, freq, num_samples freq=freq, num_samples=num_samples, ) - zbfdata = self.read_beam_file(zbf_filename) - ( - _version, - nx, - ny, - _ispol, - _unit, - _dx, - _dy, - _zposition_x, - _zposition_y, - _rayleigh_x, - _rayleigh_y, - _waist_x, - _waist_y, - wl, - index, - _receiver_eff, - _system_eff, - Ex_real, - Ex_imag, - Ey_real, - Ey_imag, - ) = zbfdata - - assert index == background_index + zbfdata = read_zbf(zbf_filename) + + assert zbfdata["background_refractive_index"] == background_index if freq is not None: - assert np.isclose(wl * 1e3, td.C_0 / freq) + assert np.isclose(zbfdata["wavelength"] * 1e3, td.C_0 / freq) else: - assert np.isclose(wl * 1e3, td.C_0 / self.monitor.freqs[len(self.monitor.freqs) // 2]) + assert np.isclose( + zbfdata["wavelength"] * 1e3, + td.C_0 / field_data.monitor.freqs[len(field_data.monitor.freqs) // 2], + ) - if num_samples is not None: - assert nx == num_samples - assert ny == num_samples - else: - assert nx == 32 - assert ny == 64 + assert zbfdata["nx"] == num_samples + assert zbfdata["ny"] == num_samples # check that fields are close - assert np.allclose(ex.real.values, Ex_real) - assert np.allclose(ex.imag.values, Ex_imag) - assert np.allclose(ey.real.values, Ey_real) - assert np.allclose(ey.imag.values, Ey_imag) + assert np.allclose(ex.values, zbfdata["Ex"]) + assert np.allclose(ey.values, zbfdata["Ey"]) From a996edd3cf954a8af68360b45c5d864035210a18 Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Mon, 28 Apr 2025 18:29:23 -0400 Subject: [PATCH 08/28] testing out ModeData.to_zbf() turns out to_zbf should be moved to parent class ElectromagneticfieldData --- tests/test_data/test_monitor_data.py | 51 +++++++++++++++++++++------- 1 file changed, 39 insertions(+), 12 deletions(-) diff --git a/tests/test_data/test_monitor_data.py b/tests/test_data/test_monitor_data.py index 4752fac85..abe6cef56 100644 --- a/tests/test_data/test_monitor_data.py +++ b/tests/test_data/test_monitor_data.py @@ -874,20 +874,13 @@ class TestZBF: freq0 = td.C_0 / 0.75 freqs = (freq0, freq0 * 1.01) - @pytest.fixture(scope="class") - def field_data(self) -> td.FieldData: - """Make random field data from an emulated simulation run.""" + def simdata(self, monitor) -> td.SimulationData: + """Returns emulated simulation data""" source = td.PointDipole( center=(-1.5, 0, 0), source_time=td.GaussianPulse(freq0=self.freq0, fwidth=self.freq0 / 10.0), polarization="Ey", ) - monitor = td.FieldMonitor( - size=(td.inf, td.inf, 0), - freqs=self.freqs, - name="fields", - colocate=True, - ) sim = td.Simulation( size=(4, 3, 3), grid_spec=td.GridSpec.auto(min_steps_per_wvl=10), @@ -896,13 +889,35 @@ def field_data(self) -> td.FieldData: monitors=[monitor], run_time=120 / self.freq0, ) - simdata = run_emulated(sim) - return simdata["fields"] + return run_emulated(sim) + + @pytest.fixture(scope="class") + def field_data(self) -> td.FieldData: + """Make random field data from an emulated simulation run.""" + monitor = td.FieldMonitor( + size=(td.inf, td.inf, 0), + freqs=self.freqs, + name="fields", + colocate=True, + ) + return self.simdata(monitor)["fields"] + + @pytest.fixture(scope="class") + def mode_data(self) -> td.ModeData: + """Make random ModeData from an emulated simulation run.""" + monitor = td.ModeMonitor( + size=(td.inf, td.inf, 0), + freqs=self.freqs, + name="modes", + colocate=True, + mode_spec=td.ModeSpec(num_modes=2, target_neff=4.0), + ) + return self.simdata(monitor)["modes"] @pytest.mark.parametrize("background_index", [1, 2, 3]) @pytest.mark.parametrize("freq", list(freqs) + [None]) @pytest.mark.parametrize("num_samples", [2, 4, 8, 16, 32, 64]) - def test_zbfread(self, tmp_path, field_data, background_index, freq, num_samples): + def test_tozbf_readzbf(self, tmp_path, field_data, background_index, freq, num_samples): zbf_filename = tmp_path / "testzbf.zbf" # write to zbf and then load it back in @@ -930,3 +945,15 @@ def test_zbfread(self, tmp_path, field_data, background_index, freq, num_samples # check that fields are close assert np.allclose(ex.values, zbfdata["Ex"]) assert np.allclose(ey.values, zbfdata["Ey"]) + + def test_tozbf_modedata(self, tmp_path, mode_data): + zbf_filename = tmp_path / "testzbf_modedata.zbf" + + # write to zbf and then load it back in + ex, ey = mode_data.to_zbf( + fname=zbf_filename, + background_refractive_index=1, + freq=self.freq0, + num_samples=16, + ) + zbfdata = read_zbf(zbf_filename) From 840bcf832c2b5480eef7ac8d1334f74e5dee3640 Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Mon, 28 Apr 2025 19:01:31 -0400 Subject: [PATCH 09/28] ModeData tests --- tests/test_data/test_monitor_data.py | 26 ++++++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/tests/test_data/test_monitor_data.py b/tests/test_data/test_monitor_data.py index abe6cef56..0ef54f563 100644 --- a/tests/test_data/test_monitor_data.py +++ b/tests/test_data/test_monitor_data.py @@ -911,13 +911,17 @@ def mode_data(self) -> td.ModeData: name="modes", colocate=True, mode_spec=td.ModeSpec(num_modes=2, target_neff=4.0), + store_fields_direction="+", ) return self.simdata(monitor)["modes"] @pytest.mark.parametrize("background_index", [1, 2, 3]) @pytest.mark.parametrize("freq", list(freqs) + [None]) @pytest.mark.parametrize("num_samples", [2, 4, 8, 16, 32, 64]) - def test_tozbf_readzbf(self, tmp_path, field_data, background_index, freq, num_samples): + def test_fielddata_tozbf_readzbf( + self, tmp_path, field_data, background_index, freq, num_samples + ): + """Test that FieldData.to_zbf() works""" zbf_filename = tmp_path / "testzbf.zbf" # write to zbf and then load it back in @@ -946,7 +950,9 @@ def test_tozbf_readzbf(self, tmp_path, field_data, background_index, freq, num_s assert np.allclose(ex.values, zbfdata["Ex"]) assert np.allclose(ey.values, zbfdata["Ey"]) - def test_tozbf_modedata(self, tmp_path, mode_data): + @pytest.mark.parametrize("mode_index", [0, 1]) + def test_tozbf_modedata(self, tmp_path, mode_data, mode_index): + """Tests ModeData.to_zbf()""" zbf_filename = tmp_path / "testzbf_modedata.zbf" # write to zbf and then load it back in @@ -955,5 +961,21 @@ def test_tozbf_modedata(self, tmp_path, mode_data): background_refractive_index=1, freq=self.freq0, num_samples=16, + mode_index=mode_index, ) zbfdata = read_zbf(zbf_filename) + + # check that fields are close + assert np.allclose(ex.values, zbfdata["Ex"]) + assert np.allclose(ey.values, zbfdata["Ey"]) + + def test_tozbf_modedata_fails(self, tmp_path, mode_data): + """Asserts that Modedata.to_zbf() fails if mode_index is not specified""" + with pytest.raises(ValueError) as e: + _ = mode_data.to_zbf( + fname=tmp_path / "testzbf_modedata_fail.zbf", + background_refractive_index=1, + freq=self.freq0, + num_samples=16, + mode_index=None, + ) From ecd16640400a2f35caa7f239aace0efd9934f21c Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Mon, 28 Apr 2025 19:04:04 -0400 Subject: [PATCH 10/28] move to_zbf to electrmagneticfielddata (should have committed w previous --- tidy3d/components/data/monitor_data.py | 247 +++++++++++++------------ 1 file changed, 126 insertions(+), 121 deletions(-) diff --git a/tidy3d/components/data/monitor_data.py b/tidy3d/components/data/monitor_data.py index 3c427c0b9..ddb49d38f 100644 --- a/tidy3d/components/data/monitor_data.py +++ b/tidy3d/components/data/monitor_data.py @@ -1067,6 +1067,132 @@ def translated_copy(self, vector: Coordinate) -> ElectromagneticFieldData: **field_kwargs, ) + def to_zbf( + self, + fname: str, + background_refractive_index: float, + freq: Optional[float] = None, + mode_index: Optional[int] = None, + num_samples: Optional[int] = None, + ) -> Tuple[ScalarFieldDataArray, ScalarFieldDataArray]: + """For a 2D monitor, export the fields to a Zemax Beam File (zbf). + + The mode area is used to approximate the beam waist, which is only valid + if the beam profile approximates a Gaussian beam. + + TODO: add pilot beam parameter as user overrides + + Parameters + ---------- + fname : str + Full path to the .zbf file to be written. + background_refractive_index : float + Refractive index of the medium surrounding the monitor. + freq : float + Field frequency selection. If ``None``, the central one is used. + mode_index : Optional[int] + For :class:`ModeData`, choose which mode to save. + num_samples: int + Number of samples to use. Must be a power of 2 + + Returns + ------- + Tuple[:class:`ScalarFieldDataArray`,:class:`ScalarFieldDataArray`] + The two E field components being exported to zbf. + """ + log.warning( + "Exporting to ZBF is currently an experimental feature. If any issues are encountered, please contact Flexcompute support." + ) + + # Mode area calculation ensures all E components are present + mode_area = self.mode_area + dim1, dim2 = self._tangential_dims + + # Using file-local coordinates x, y for the tangential components + e_x = self._tangential_fields["E" + dim1] + e_y = self._tangential_fields["E" + dim2] + x = e_x.coords[dim1].values + y = e_x.coords[dim2].values + + # Use the central frequency if freq is not specified + if freq is None: + freq = e_x.coords["f"].values[e_x.coords["f"].size // 2] + + mode_area = mode_area.sel(f=freq, drop=True) + e_x = e_x.sel(f=freq, drop=True) + e_y = e_y.sel(f=freq, drop=True) + + # If the data is ModeData, choose one of the modes to save + if "mode_index" in e_x.coords: + if mode_index is None: + raise ValueError("`mode_index` is required for `ModeData.to_zbf()`") + mode_area = mode_area.isel(mode_index=mode_index, drop=True) + e_x = e_x.isel(mode_index=mode_index, drop=True) + e_y = e_y.isel(mode_index=mode_index, drop=True) + + # Header info + version = 1 + polarized = 1 + units = 0 # "mm": 0, "cm": 1, "in": 2, "m": 3 + lda = C_0 / freq * 1e-3 + + # Pilot (reference) beam waist: use the mode area to approximate the expected value + w_x = (mode_area.item() / np.pi) ** 0.5 * 1e-3 + w_y = w_x + + # Pilot beam Rayleigh distance (ignored on input) + r_x = 0 + r_y = 0 + + # Pilot beam z position w.r.t. the waist + z_x = 0 * 1e-3 + z_y = 0 * 1e-3 + + # Number of samples: n_x = 2 ** a, a >= 5 and a <= 13 (32-bit), or a <= 14 (64-bit) + if num_samples is not None: + # custom user defined number of samples + # Raise an error if the number of samples is not a power of 2 + if (num_samples & (num_samples - 1)) != 0: + raise ValueError("num_samples must be a power of 2.") + n_x = num_samples + n_y = num_samples + else: + n_x = 2 ** min(13, max(5, int(np.log2(x.size) + 1))) + n_y = 2 ** min(13, max(5, int(np.log2(y.size) + 1))) + + # Interpolating coordinates + x = np.linspace(x.min(), x.max(), n_x) + y = np.linspace(y.min(), y.max(), n_y) + + # Interpolate fields + coords = {dim1: x, dim2: y} + e_x = e_x.interp(coords, assume_sorted=True) + e_y = e_y.interp(coords, assume_sorted=True) + + # Sampling distance + d_x = np.mean(np.diff(x)) * 1e-3 + d_y = np.mean(np.diff(y)) * 1e-3 + + # Receiver and system efficiencies + rec_efficiency = 0 # zero if fiber coupling is not computed + sys_efficiency = 0 # zero if fiber coupling is not computed + + with open(fname, "wb") as fout: + fout.write(struct.pack("<5I", version, n_x, n_y, polarized, units)) + fout.write(struct.pack("<4I", 0, 0, 0, 0)) # unused values + fout.write(struct.pack("<8d", d_x, d_y, z_x, r_x, w_x, z_y, r_y, w_y)) + fout.write( + struct.pack("<4d", lda, background_refractive_index, rec_efficiency, sys_efficiency) + ) + fout.write(struct.pack("<8d", 0, 0, 0, 0, 0, 0, 0, 0)) # unused values + for e in (e_x, e_y): + e_flat = e.values.flatten(order="C") + # Interweave real and imaginary parts + e_values = np.ravel(np.column_stack((e_flat.real, e_flat.imag))) + fout.write(struct.pack(f"<{2 * n_x*n_y}d", *e_values)) + + return e_x, e_y + class FieldData(FieldDataset, ElectromagneticFieldData): """ @@ -1225,127 +1351,6 @@ def make_adjoint_sources( return sources - def to_zbf( - self, - fname: str, - background_refractive_index: float, - freq: Optional[float] = None, - mode_index: Optional[int] = 0, - num_samples: Optional[int] = None, - ) -> Tuple[ScalarFieldDataArray, ScalarFieldDataArray]: - """For a 2D monitor, export the fields to a Zemax Beam File (zbf). - - The mode area is used to approximate the beam waist, which is only valid - if the beam profile approximates a Gaussian beam. - - Parameters - ---------- - fname : str - Full path to the .zbf file to be written. - background_refractive_index : float - Refractive index of the medium surrounding the monitor. - freq : float - Field frequency selection. If ``None``, the central one is used. - mode_index : int - Mode index selection. - num_samples: int - Number of samples to use. Must be a power of 2 - - Returns - ------- - Tuple[:class:`ScalarFieldDataArray`,:class:`ScalarFieldDataArray`] - The two E field components being exported to zbf. - """ - log.warning( - "Exporting to ZBF is currently an experimental feature. If any issues are encountered, please contact Flexcompute support." - ) - - # Mode area calculation ensures all E components are present - mode_area = self.mode_area - dim1, dim2 = self._tangential_dims - - # Using file-local coordinates x, y for the tangential components - e_x = self._tangential_fields["E" + dim1] - e_y = self._tangential_fields["E" + dim2] - x = e_x.coords[dim1].values - y = e_x.coords[dim2].values - - # Use the central frequency if freq is not specified - if freq is None: - freq = e_x.coords["f"].values[e_x.coords["f"].size // 2] - - mode_area = mode_area.sel(f=freq, drop=True) - e_x = e_x.sel(f=freq, drop=True) - e_y = e_y.sel(f=freq, drop=True) - - if "mode_index" in e_x.coords: - mode_area = mode_area.isel(mode_index=mode_index, drop=True) - e_x = e_x.isel(mode_index=mode_index, drop=True) - e_y = e_y.isel(mode_index=mode_index, drop=True) - - # Header info - version = 1 - polarized = 1 - units = 0 # "mm": 0, "cm": 1, "in": 2, "m": 3 - lda = C_0 / freq * 1e-3 - - # Pilot (reference) beam waist: use the mode area to approximate the expected value - w_x = (mode_area.item() / np.pi) ** 0.5 * 1e-3 - w_y = w_x - - # Pilot beam Rayleigh distance (ignored on input) - r_x = 0 - r_y = 0 - - # Pilot beam z position w.r.t. the waist - z_x = 0 * 1e-3 - z_y = 0 * 1e-3 - - # Number of samples: n_x = 2 ** a, a >= 5 and a <= 13 (32-bit), or a <= 14 (64-bit) - if num_samples is not None: - # custom user defined number of samples - # Raise an error if the number of samples is not a power of 2 - if (num_samples & (num_samples - 1)) != 0: - raise ValueError("num_samples must be a power of 2.") - n_x = num_samples - n_y = num_samples - else: - n_x = 2 ** min(13, max(5, int(np.log2(x.size) + 1))) - n_y = 2 ** min(13, max(5, int(np.log2(y.size) + 1))) - - # Interpolating coordinates - x = np.linspace(x.min(), x.max(), n_x) - y = np.linspace(y.min(), y.max(), n_y) - - # Interpolate fields - coords = {dim1: x, dim2: y} - e_x = e_x.interp(coords, assume_sorted=True) - e_y = e_y.interp(coords, assume_sorted=True) - - # Sampling distance - d_x = np.mean(np.diff(x)) * 1e-3 - d_y = np.mean(np.diff(y)) * 1e-3 - - # Receiver and system efficiencies - rec_efficiency = 0 # zero if fiber coupling is not computed - sys_efficiency = 0 # zero if fiber coupling is not computed - - with open(fname, "wb") as fout: - fout.write(struct.pack("<5I", version, n_x, n_y, polarized, units)) - fout.write(struct.pack("<4I", 0, 0, 0, 0)) # unused values - fout.write(struct.pack("<8d", d_x, d_y, z_x, r_x, w_x, z_y, r_y, w_y)) - fout.write( - struct.pack("<4d", lda, background_refractive_index, rec_efficiency, sys_efficiency) - ) - fout.write(struct.pack("<8d", 0, 0, 0, 0, 0, 0, 0, 0)) # unused values - for e in (e_x, e_y): - e_flat = e.values.flatten(order="C") - # Interweave real and imaginary parts - e_values = np.ravel(np.column_stack((e_flat.real, e_flat.imag))) - fout.write(struct.pack(f"<{2 * n_x*n_y}d", *e_values)) - - return e_x, e_y - class FieldTimeData(FieldTimeDataset, ElectromagneticFieldData): """ From 77cf2ba965e6cce20bb8ec0f847f864ec9a15351 Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Tue, 29 Apr 2025 21:51:27 -0400 Subject: [PATCH 11/28] replace num_samples with n_x and n_y as required args --- tidy3d/components/data/monitor_data.py | 34 ++++++++++++++------------ 1 file changed, 19 insertions(+), 15 deletions(-) diff --git a/tidy3d/components/data/monitor_data.py b/tidy3d/components/data/monitor_data.py index ddb49d38f..f8ece1099 100644 --- a/tidy3d/components/data/monitor_data.py +++ b/tidy3d/components/data/monitor_data.py @@ -1071,9 +1071,10 @@ def to_zbf( self, fname: str, background_refractive_index: float, + n_x: int, + n_y: int, freq: Optional[float] = None, mode_index: Optional[int] = None, - num_samples: Optional[int] = None, ) -> Tuple[ScalarFieldDataArray, ScalarFieldDataArray]: """For a 2D monitor, export the fields to a Zemax Beam File (zbf). @@ -1090,10 +1091,12 @@ def to_zbf( Refractive index of the medium surrounding the monitor. freq : float Field frequency selection. If ``None``, the central one is used. + n_x : int + Number of field samples along x. Must be a power of 2, between 2^5 and 2^13 inclusive per Zemax's requirements. + n_y : int + Number of field samples along y. Must be a power of 2, between 2^5 and 2^13 inclusive per Zemax's requirements. mode_index : Optional[int] For :class:`ModeData`, choose which mode to save. - num_samples: int - Number of samples to use. Must be a power of 2 Returns ------- @@ -1104,6 +1107,19 @@ def to_zbf( "Exporting to ZBF is currently an experimental feature. If any issues are encountered, please contact Flexcompute support." ) + # Check that requirements are met for n_x and n_y + # TODO write tests for this + # n_x and n_y must be powers of 2 + if (n_x & (n_x - 1)) != 0: + raise ValueError("`n_x` must be a power of 2.") + if (n_y & (n_y - 1)) != 0: + raise ValueError("`n_y` must be a power of 2.") + # 32 <= n_x and n_y <= 2^13 + if n_x < 32 or n_x > 2**13: + raise ValueError("`n_x` must be between 2^5 and 2^13, inclusive.") + if n_y < 32 or n_y > 2**13: + raise ValueError("`n_y` must be between 2^5 and 2^13, inclusive.") + # Mode area calculation ensures all E components are present mode_area = self.mode_area dim1, dim2 = self._tangential_dims @@ -1148,18 +1164,6 @@ def to_zbf( z_x = 0 * 1e-3 z_y = 0 * 1e-3 - # Number of samples: n_x = 2 ** a, a >= 5 and a <= 13 (32-bit), or a <= 14 (64-bit) - if num_samples is not None: - # custom user defined number of samples - # Raise an error if the number of samples is not a power of 2 - if (num_samples & (num_samples - 1)) != 0: - raise ValueError("num_samples must be a power of 2.") - n_x = num_samples - n_y = num_samples - else: - n_x = 2 ** min(13, max(5, int(np.log2(x.size) + 1))) - n_y = 2 ** min(13, max(5, int(np.log2(y.size) + 1))) - # Interpolating coordinates x = np.linspace(x.min(), x.max(), n_x) y = np.linspace(y.min(), y.max(), n_y) From ff2ff1516a68f6141a6bbc37cb6ebe4cc4d59fb5 Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Tue, 29 Apr 2025 21:56:59 -0400 Subject: [PATCH 12/28] add tests for n_x and n_y --- tests/test_data/test_monitor_data.py | 27 ++++++++++++++++++++------- 1 file changed, 20 insertions(+), 7 deletions(-) diff --git a/tests/test_data/test_monitor_data.py b/tests/test_data/test_monitor_data.py index 0ef54f563..6c8fba3f5 100644 --- a/tests/test_data/test_monitor_data.py +++ b/tests/test_data/test_monitor_data.py @@ -917,10 +917,9 @@ def mode_data(self) -> td.ModeData: @pytest.mark.parametrize("background_index", [1, 2, 3]) @pytest.mark.parametrize("freq", list(freqs) + [None]) - @pytest.mark.parametrize("num_samples", [2, 4, 8, 16, 32, 64]) - def test_fielddata_tozbf_readzbf( - self, tmp_path, field_data, background_index, freq, num_samples - ): + @pytest.mark.parametrize("n_x", [2**5, 2**6]) + @pytest.mark.parametrize("n_y", [2**5, 2**6]) + def test_fielddata_tozbf_readzbf(self, tmp_path, field_data, background_index, freq, n_x, n_y): """Test that FieldData.to_zbf() works""" zbf_filename = tmp_path / "testzbf.zbf" @@ -929,7 +928,8 @@ def test_fielddata_tozbf_readzbf( fname=zbf_filename, background_refractive_index=background_index, freq=freq, - num_samples=num_samples, + n_x=n_x, + n_y=n_y, ) zbfdata = read_zbf(zbf_filename) @@ -943,8 +943,8 @@ def test_fielddata_tozbf_readzbf( td.C_0 / field_data.monitor.freqs[len(field_data.monitor.freqs) // 2], ) - assert zbfdata["nx"] == num_samples - assert zbfdata["ny"] == num_samples + assert zbfdata["nx"] == n_x + assert zbfdata["ny"] == n_y # check that fields are close assert np.allclose(ex.values, zbfdata["Ex"]) @@ -979,3 +979,16 @@ def test_tozbf_modedata_fails(self, tmp_path, mode_data): num_samples=16, mode_index=None, ) + + @pytest.mark.parametrize("n_x", [16, 2**14, 33]) + @pytest.mark.parametrize("n_y", [16, 2**14, 33]) + def test_tozbf_nxny_fails(self, tmp_path, field_data, n_x, n_y): + """Asserts that to_zbf() fails when n_x and n_y are invalid values.""" + with pytest.raises(ValueError) as e: + _ = field_data.to_zbf( + fname=tmp_path / "testzbf_nxny_fail.zbf", + background_refractive_index=1, + freq=self.freq0, + n_x=n_x, + n_y=n_y, + ) From 120f2433f22c5ef5a1f644f46444ea6df556044d Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Tue, 29 Apr 2025 21:57:12 -0400 Subject: [PATCH 13/28] remove todo msg --- tidy3d/components/data/monitor_data.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tidy3d/components/data/monitor_data.py b/tidy3d/components/data/monitor_data.py index f8ece1099..97e6550b8 100644 --- a/tidy3d/components/data/monitor_data.py +++ b/tidy3d/components/data/monitor_data.py @@ -1108,7 +1108,6 @@ def to_zbf( ) # Check that requirements are met for n_x and n_y - # TODO write tests for this # n_x and n_y must be powers of 2 if (n_x & (n_x - 1)) != 0: raise ValueError("`n_x` must be a power of 2.") From 16cd55a790be193e33289913c9349c9130ae2900 Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Tue, 29 Apr 2025 21:59:05 -0400 Subject: [PATCH 14/28] update warning msg with support url --- tidy3d/components/data/monitor_data.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tidy3d/components/data/monitor_data.py b/tidy3d/components/data/monitor_data.py index 97e6550b8..6d6507a83 100644 --- a/tidy3d/components/data/monitor_data.py +++ b/tidy3d/components/data/monitor_data.py @@ -1104,7 +1104,8 @@ def to_zbf( The two E field components being exported to zbf. """ log.warning( - "Exporting to ZBF is currently an experimental feature. If any issues are encountered, please contact Flexcompute support." + "`FieldData.to_zbf()` is currently an experimental feature." + "If any issues are encountered, please contact Flexcompute support https://www.flexcompute.com/tidy3d/technical-support/" ) # Check that requirements are met for n_x and n_y From 77f5e23e332303eb34888d00f8988be6d1ee4576 Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Tue, 29 Apr 2025 22:03:19 -0400 Subject: [PATCH 15/28] default frqeuency to the average of the recorded frequencies if not specified --- tests/test_data/test_monitor_data.py | 2 +- tidy3d/components/data/monitor_data.py | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/tests/test_data/test_monitor_data.py b/tests/test_data/test_monitor_data.py index 6c8fba3f5..276623ce2 100644 --- a/tests/test_data/test_monitor_data.py +++ b/tests/test_data/test_monitor_data.py @@ -940,7 +940,7 @@ def test_fielddata_tozbf_readzbf(self, tmp_path, field_data, background_index, f else: assert np.isclose( zbfdata["wavelength"] * 1e3, - td.C_0 / field_data.monitor.freqs[len(field_data.monitor.freqs) // 2], + td.C_0 / np.mean(field_data.monitor.freqs), ) assert zbfdata["nx"] == n_x diff --git a/tidy3d/components/data/monitor_data.py b/tidy3d/components/data/monitor_data.py index 6d6507a83..7318106e5 100644 --- a/tidy3d/components/data/monitor_data.py +++ b/tidy3d/components/data/monitor_data.py @@ -1090,7 +1090,7 @@ def to_zbf( background_refractive_index : float Refractive index of the medium surrounding the monitor. freq : float - Field frequency selection. If ``None``, the central one is used. + Field frequency selection. If ``None``, the average of the recorded frequencies is used. n_x : int Number of field samples along x. Must be a power of 2, between 2^5 and 2^13 inclusive per Zemax's requirements. n_y : int @@ -1130,13 +1130,13 @@ def to_zbf( x = e_x.coords[dim1].values y = e_x.coords[dim2].values - # Use the central frequency if freq is not specified + # Use the mean frequency if freq is not specified if freq is None: - freq = e_x.coords["f"].values[e_x.coords["f"].size // 2] + freq = np.mean(e_x.coords["f"].values) - mode_area = mode_area.sel(f=freq, drop=True) - e_x = e_x.sel(f=freq, drop=True) - e_y = e_y.sel(f=freq, drop=True) + mode_area = mode_area.interp(f=freq) + e_x = e_x.interp(f=freq) + e_y = e_y.interp(f=freq) # If the data is ModeData, choose one of the modes to save if "mode_index" in e_x.coords: From 277cd89dfb3eb9242d6ba681edd1f561f204623c Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Wed, 30 Apr 2025 17:25:15 -0400 Subject: [PATCH 16/28] let user pick units when writing zbf --- tests/test_data/test_monitor_data.py | 35 ++++++++++++++++++++++---- tidy3d/components/data/monitor_data.py | 32 +++++++++++++++-------- tidy3d/components/types.py | 1 + 3 files changed, 52 insertions(+), 16 deletions(-) diff --git a/tests/test_data/test_monitor_data.py b/tests/test_data/test_monitor_data.py index 276623ce2..29039715c 100644 --- a/tests/test_data/test_monitor_data.py +++ b/tests/test_data/test_monitor_data.py @@ -22,6 +22,7 @@ PermittivityData, ) from tidy3d.components.data.zbf import read_zbf +from tidy3d.constants import UnitScaling from tidy3d.exceptions import DataError from ..utils import AssertLogLevel, run_emulated @@ -919,7 +920,10 @@ def mode_data(self) -> td.ModeData: @pytest.mark.parametrize("freq", list(freqs) + [None]) @pytest.mark.parametrize("n_x", [2**5, 2**6]) @pytest.mark.parametrize("n_y", [2**5, 2**6]) - def test_fielddata_tozbf_readzbf(self, tmp_path, field_data, background_index, freq, n_x, n_y): + @pytest.mark.parametrize("units", ["mm", "cm", "in", "m"]) + def test_fielddata_tozbf_readzbf( + self, tmp_path, field_data, background_index, freq, n_x, n_y, units + ): """Test that FieldData.to_zbf() works""" zbf_filename = tmp_path / "testzbf.zbf" @@ -930,16 +934,19 @@ def test_fielddata_tozbf_readzbf(self, tmp_path, field_data, background_index, f freq=freq, n_x=n_x, n_y=n_y, + units=units, ) zbfdata = read_zbf(zbf_filename) assert zbfdata["background_refractive_index"] == background_index + unitscaling = UnitScaling[units] + if freq is not None: - assert np.isclose(zbfdata["wavelength"] * 1e3, td.C_0 / freq) + assert np.isclose(zbfdata["wavelength"] / unitscaling, td.C_0 / freq) else: assert np.isclose( - zbfdata["wavelength"] * 1e3, + zbfdata["wavelength"] / unitscaling, td.C_0 / np.mean(field_data.monitor.freqs), ) @@ -960,8 +967,10 @@ def test_tozbf_modedata(self, tmp_path, mode_data, mode_index): fname=zbf_filename, background_refractive_index=1, freq=self.freq0, - num_samples=16, mode_index=mode_index, + n_x=32, + n_y=32, + units="mm", ) zbfdata = read_zbf(zbf_filename) @@ -976,8 +985,10 @@ def test_tozbf_modedata_fails(self, tmp_path, mode_data): fname=tmp_path / "testzbf_modedata_fail.zbf", background_refractive_index=1, freq=self.freq0, - num_samples=16, mode_index=None, + n_x=32, + n_y=32, + units="mm", ) @pytest.mark.parametrize("n_x", [16, 2**14, 33]) @@ -991,4 +1002,18 @@ def test_tozbf_nxny_fails(self, tmp_path, field_data, n_x, n_y): freq=self.freq0, n_x=n_x, n_y=n_y, + units="mm", + ) + + @pytest.mark.parametrize("units", ["mmm", "123"]) + def test_tozbf_units_fails(self, tmp_path, field_data, units): + """Asserts that to_zbf() fails when units are invalid.""" + with pytest.raises(ValueError) as e: + _ = field_data.to_zbf( + fname=tmp_path / "testzbf_nxny_fail.zbf", + background_refractive_index=1, + freq=self.freq0, + n_x=32, + n_y=32, + units=units, ) diff --git a/tidy3d/components/data/monitor_data.py b/tidy3d/components/data/monitor_data.py index 7318106e5..e9b5e1787 100644 --- a/tidy3d/components/data/monitor_data.py +++ b/tidy3d/components/data/monitor_data.py @@ -6,7 +6,7 @@ import warnings from abc import ABC from math import isclose -from typing import Any, Callable, Dict, List, Optional, Tuple, Union +from typing import Any, Callable, Dict, List, Optional, Tuple, Union, get_args import autograd.numpy as np import pydantic.v1 as pd @@ -14,7 +14,7 @@ from pandas import DataFrame from xarray.core.types import Self -from ...constants import C_0, EPSILON_0, ETA_0, MICROMETER +from ...constants import C_0, EPSILON_0, ETA_0, MICROMETER, UnitScaling from ...exceptions import DataError, SetupError, Tidy3dNotImplementedError, ValidationError from ...log import log from ..base import TYPE_TAG_STR, cached_property, skip_if_fields_missing @@ -64,6 +64,7 @@ Size, Symmetry, TrackFreq, + UnitsZBF, ) from ..validators import enforce_monitor_fields_present, required_if_symmetry_present from .data_array import ( @@ -1073,6 +1074,7 @@ def to_zbf( background_refractive_index: float, n_x: int, n_y: int, + units: UnitsZBF, freq: Optional[float] = None, mode_index: Optional[int] = None, ) -> Tuple[ScalarFieldDataArray, ScalarFieldDataArray]: @@ -1089,12 +1091,14 @@ def to_zbf( Full path to the .zbf file to be written. background_refractive_index : float Refractive index of the medium surrounding the monitor. - freq : float - Field frequency selection. If ``None``, the average of the recorded frequencies is used. n_x : int Number of field samples along x. Must be a power of 2, between 2^5 and 2^13 inclusive per Zemax's requirements. n_y : int Number of field samples along y. Must be a power of 2, between 2^5 and 2^13 inclusive per Zemax's requirements. + units : `UnitsZBF` + Spatial units used for the zbf file. Options are `"mm"`, `"cm"`, `"in"`, or `"m"`. + freq : float + Field frequency selection. If ``None``, the average of the recorded frequencies is used. mode_index : Optional[int] For :class:`ModeData`, choose which mode to save. @@ -1108,6 +1112,10 @@ def to_zbf( "If any issues are encountered, please contact Flexcompute support https://www.flexcompute.com/tidy3d/technical-support/" ) + # Check that appropriate units are used + if units not in get_args(UnitsZBF): + raise ValueError("`units` must be either `'mm'`, `'cm'`, `'in'`, or `'m'`.") + # Check that requirements are met for n_x and n_y # n_x and n_y must be powers of 2 if (n_x & (n_x - 1)) != 0: @@ -1149,11 +1157,13 @@ def to_zbf( # Header info version = 1 polarized = 1 - units = 0 # "mm": 0, "cm": 1, "in": 2, "m": 3 - lda = C_0 / freq * 1e-3 + unit_mapping = {"mm": 0, "cm": 1, "in": 2, "m": 3} + unit_key = unit_mapping[units] + unit_scaling = UnitScaling[units] + lda = C_0 / freq * unit_scaling # Pilot (reference) beam waist: use the mode area to approximate the expected value - w_x = (mode_area.item() / np.pi) ** 0.5 * 1e-3 + w_x = (mode_area.item() / np.pi) ** 0.5 * unit_scaling w_y = w_x # Pilot beam Rayleigh distance (ignored on input) @@ -1164,7 +1174,7 @@ def to_zbf( z_x = 0 * 1e-3 z_y = 0 * 1e-3 - # Interpolating coordinates + # Interpolating coordinatesΠx = np.linspace(x.min(), x.max(), n_x) y = np.linspace(y.min(), y.max(), n_y) @@ -1174,15 +1184,15 @@ def to_zbf( e_y = e_y.interp(coords, assume_sorted=True) # Sampling distance - d_x = np.mean(np.diff(x)) * 1e-3 - d_y = np.mean(np.diff(y)) * 1e-3 + d_x = np.mean(np.diff(x)) * unit_scaling + d_y = np.mean(np.diff(y)) * unit_scaling # Receiver and system efficiencies rec_efficiency = 0 # zero if fiber coupling is not computed sys_efficiency = 0 # zero if fiber coupling is not computed with open(fname, "wb") as fout: - fout.write(struct.pack("<5I", version, n_x, n_y, polarized, units)) + fout.write(struct.pack("<5I", version, n_x, n_y, polarized, unit_key)) fout.write(struct.pack("<4I", 0, 0, 0, 0)) # unused values fout.write(struct.pack("<8d", d_x, d_y, z_x, r_x, w_x, z_y, r_y, w_y)) fout.write( diff --git a/tidy3d/components/types.py b/tidy3d/components/types.py index 9ff928d50..61697f87d 100644 --- a/tidy3d/components/types.py +++ b/tidy3d/components/types.py @@ -256,3 +256,4 @@ def __modify_schema__(cls, field_schema): """ dataset """ xyz = Literal["x", "y", "z"] +UnitsZBF = Literal["mm", "cm", "in", "m"] From 5d2fbbf3945b321897e7c5f167f70c7ea7844d98 Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Wed, 30 Apr 2025 17:32:42 -0400 Subject: [PATCH 17/28] add miscellaneous zemax options as kwargs if the user wants to define them --- tidy3d/components/data/monitor_data.py | 30 +++++++++++++++++++------- 1 file changed, 22 insertions(+), 8 deletions(-) diff --git a/tidy3d/components/data/monitor_data.py b/tidy3d/components/data/monitor_data.py index e9b5e1787..721227924 100644 --- a/tidy3d/components/data/monitor_data.py +++ b/tidy3d/components/data/monitor_data.py @@ -1077,6 +1077,12 @@ def to_zbf( units: UnitsZBF, freq: Optional[float] = None, mode_index: Optional[int] = None, + r_x: float = 0, + r_y: float = 0, + z_x: float = 0, + z_y: float = 0, + rec_efficiency: float = 0, + sys_efficiency: float = 0, ) -> Tuple[ScalarFieldDataArray, ScalarFieldDataArray]: """For a 2D monitor, export the fields to a Zemax Beam File (zbf). @@ -1101,6 +1107,18 @@ def to_zbf( Field frequency selection. If ``None``, the average of the recorded frequencies is used. mode_index : Optional[int] For :class:`ModeData`, choose which mode to save. + r_x : float + Pilot beam Rayleigh distance in x, um. Defaults to 0. + r_y : float + Pilot beam Rayleigh distance in y, um. Defaults to 0. + z_x : float + Pilot beam z position w.r.t. the waist in x, um. Defaults to 0. + z_y : float + Pilot beam z position w.r.t. the waist in y, um. Defaults to 0. + rec_efficiency : float + Receiver efficiency, zero if fiber coupling is not computed. Defaults to 0. + sys_efficiency : float + System efficiency, zero if fiber coupling is not computed. Defaults to 0. Returns ------- @@ -1167,12 +1185,12 @@ def to_zbf( w_y = w_x # Pilot beam Rayleigh distance (ignored on input) - r_x = 0 - r_y = 0 + r_x *= unit_scaling + r_y *= unit_scaling # Pilot beam z position w.r.t. the waist - z_x = 0 * 1e-3 - z_y = 0 * 1e-3 + z_x *= unit_scaling + z_y *= unit_scaling # Interpolating coordinatesΠx = np.linspace(x.min(), x.max(), n_x) @@ -1187,10 +1205,6 @@ def to_zbf( d_x = np.mean(np.diff(x)) * unit_scaling d_y = np.mean(np.diff(y)) * unit_scaling - # Receiver and system efficiencies - rec_efficiency = 0 # zero if fiber coupling is not computed - sys_efficiency = 0 # zero if fiber coupling is not computed - with open(fname, "wb") as fout: fout.write(struct.pack("<5I", version, n_x, n_y, polarized, unit_key)) fout.write(struct.pack("<4I", 0, 0, 0, 0)) # unused values From 7f27f05de67e7c5b218b1eb030679528b7e548fa Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Wed, 30 Apr 2025 17:35:06 -0400 Subject: [PATCH 18/28] add missing space in warning msg --- tidy3d/components/data/dataset.py | 2 +- tidy3d/components/data/monitor_data.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tidy3d/components/data/dataset.py b/tidy3d/components/data/dataset.py index 87f78f0ef..ca1e8ac57 100644 --- a/tidy3d/components/data/dataset.py +++ b/tidy3d/components/data/dataset.py @@ -279,7 +279,7 @@ def from_zbf(filename: str, dim1: xyz, dim2: xyz): """ log.warning( "`FieldDataset.from_zbf()` is currently an experimental feature." - "If any issues are encountered, please contact Flexcompute support https://www.flexcompute.com/tidy3d/technical-support/" + " If any issues are encountered, please contact Flexcompute support https://www.flexcompute.com/tidy3d/technical-support/" ) if dim1 not in get_args(xyz): diff --git a/tidy3d/components/data/monitor_data.py b/tidy3d/components/data/monitor_data.py index 721227924..74f6a0755 100644 --- a/tidy3d/components/data/monitor_data.py +++ b/tidy3d/components/data/monitor_data.py @@ -1127,7 +1127,7 @@ def to_zbf( """ log.warning( "`FieldData.to_zbf()` is currently an experimental feature." - "If any issues are encountered, please contact Flexcompute support https://www.flexcompute.com/tidy3d/technical-support/" + " If any issues are encountered, please contact Flexcompute support https://www.flexcompute.com/tidy3d/technical-support/" ) # Check that appropriate units are used From d40f1fdd9f8e311feab9a6a3f1cdd2b1f4c0f22c Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Wed, 30 Apr 2025 17:36:59 -0400 Subject: [PATCH 19/28] add double quotes around value of dim in error msg --- tidy3d/components/data/dataset.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/tidy3d/components/data/dataset.py b/tidy3d/components/data/dataset.py index ca1e8ac57..bb341d102 100644 --- a/tidy3d/components/data/dataset.py +++ b/tidy3d/components/data/dataset.py @@ -283,9 +283,13 @@ def from_zbf(filename: str, dim1: xyz, dim2: xyz): ) if dim1 not in get_args(xyz): - raise ValueError(f'`dim1` = {dim1} not allowed, must be one of `"x"`, `"y"`, or `"z"`.') + raise ValueError( + f'`dim1` = "{dim1}" not allowed, must be one of `"x"`, `"y"`, or `"z"`.' + ) if dim2 not in get_args(xyz): - raise ValueError(f'`dim2` = {dim2} not allowed, must be one of `"x"`, `"y"`, or `"z"`.') + raise ValueError( + f'`dim2` = "{dim2}" not allowed, must be one of `"x"`, `"y"`, or `"z"`.' + ) if dim1 == dim2: raise ValueError("`dim1` and `dim2` must be different.") From 0b174f1a85f30c87cd6486fa5a55b2c227e9ee65 Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Thu, 1 May 2025 14:09:16 -0400 Subject: [PATCH 20/28] add from_zbf tests --- tests/test_data/test_monitor_data.py | 38 ++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/tests/test_data/test_monitor_data.py b/tests/test_data/test_monitor_data.py index 29039715c..d14ab6507 100644 --- a/tests/test_data/test_monitor_data.py +++ b/tests/test_data/test_monitor_data.py @@ -1017,3 +1017,41 @@ def test_tozbf_units_fails(self, tmp_path, field_data, units): n_y=32, units=units, ) + + def test_from_zbf(self, tmp_path, field_data): + """Tests creating a field dataset from a zbf""" + zbf_filename = tmp_path / "testzbf.zbf" + # write to zbf and then load it back in + ex, ey = field_data.to_zbf( + fname=zbf_filename, + background_refractive_index=1, + n_x=32, + n_y=32, + units="mm", + ) + + # create a field dataset from the zbf file + fd = td.FieldDataset.from_zbf(filename=zbf_filename, dim1="x", dim2="y") + + # compare loaded field data to saved data + assert np.allclose(ex.values, fd.Ex.values.squeeze()) + assert np.allclose(ey.values, fd.Ey.values.squeeze()) + + @pytest.mark.parametrize( + "dim1,dim2", [("x", "x"), ("y", "y"), ("z", "z"), ("1", "2"), ("c", "z")] + ) + def test_from_zbf_dimsfail(self, tmp_path, field_data, dim1, dim2): + """Tests creating a field dataset from a zbf""" + zbf_filename = tmp_path / "testzbf.zbf" + # write to zbf and then load it back in + _, _ = field_data.to_zbf( + fname=zbf_filename, + background_refractive_index=1, + n_x=32, + n_y=32, + units="mm", + ) + + # create a field dataset from the zbf file + with pytest.raises(ValueError) as e: + _ = td.FieldDataset.from_zbf(filename=zbf_filename, dim1=dim1, dim2=dim2) From 88a7a20fb11e9a5f6cd84dd04815e8b1feccb747 Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Thu, 1 May 2025 14:15:29 -0400 Subject: [PATCH 21/28] remove todo comment --- tidy3d/components/data/monitor_data.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/tidy3d/components/data/monitor_data.py b/tidy3d/components/data/monitor_data.py index 74f6a0755..cc5e1c938 100644 --- a/tidy3d/components/data/monitor_data.py +++ b/tidy3d/components/data/monitor_data.py @@ -1089,8 +1089,6 @@ def to_zbf( The mode area is used to approximate the beam waist, which is only valid if the beam profile approximates a Gaussian beam. - TODO: add pilot beam parameter as user overrides - Parameters ---------- fname : str From fe6faac2f3d5b27c573866db970d2e5edd28fb83 Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Thu, 1 May 2025 14:21:24 -0400 Subject: [PATCH 22/28] update changelog --- CHANGELOG.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 333788c47..4db8c61f4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,7 +13,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Methods `EMEExplicitGrid.from_structures` and `EMECompositeGrid.from_structure_groups` to place EME cell boundaries at structure bounds. - 'ModeSimulation' now supports 'PermittivityMonitor'. - Classmethod `from_frequency_range` in `GaussianPulse` for generating a pulse whose amplitude in the frequency_range [fmin, fmax] is maximized, which is particularly useful for running broadband simulations. -- `FieldData` supports exporting E fields to a Zemax Beam File (ZBF) with `FieldData.to_zbf()`. +- `FieldData` and `ModeData` support exporting E fields to a Zemax Beam File (ZBF) with `.to_zbf()` (warning: experimental feature). +- `FieldDataset` supports reading E fields from a Zemax Beam File (ZBF) with `.from_zbf()` (warning: experimental feature). ### Changed - Performance enhancement for adjoint gradient calculations by optimizing field interpolation. From e1e1fbe1808b2f1f8c87ccb1c765c1af87d40ebd Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Mon, 5 May 2025 11:28:42 -0400 Subject: [PATCH 23/28] use pydantic model for zbf data --- tests/test_data/test_monitor_data.py | 31 ++- tidy3d/components/data/dataset.py | 50 ++--- tidy3d/components/data/zbf.py | 276 +++++++++++++++------------ 3 files changed, 200 insertions(+), 157 deletions(-) diff --git a/tests/test_data/test_monitor_data.py b/tests/test_data/test_monitor_data.py index d14ab6507..5ea9d6a6c 100644 --- a/tests/test_data/test_monitor_data.py +++ b/tests/test_data/test_monitor_data.py @@ -21,7 +21,7 @@ ModeData, PermittivityData, ) -from tidy3d.components.data.zbf import read_zbf +from tidy3d.components.data.zbf import ZBFData from tidy3d.constants import UnitScaling from tidy3d.exceptions import DataError @@ -924,7 +924,7 @@ def mode_data(self) -> td.ModeData: def test_fielddata_tozbf_readzbf( self, tmp_path, field_data, background_index, freq, n_x, n_y, units ): - """Test that FieldData.to_zbf() works""" + """Test that FieldData.to_zbf() -> ZBFData.read_zbf() works""" zbf_filename = tmp_path / "testzbf.zbf" # write to zbf and then load it back in @@ -936,26 +936,26 @@ def test_fielddata_tozbf_readzbf( n_y=n_y, units=units, ) - zbfdata = read_zbf(zbf_filename) + zbfdata = ZBFData.read_zbf(zbf_filename) - assert zbfdata["background_refractive_index"] == background_index + assert zbfdata.background_refractive_index == background_index unitscaling = UnitScaling[units] if freq is not None: - assert np.isclose(zbfdata["wavelength"] / unitscaling, td.C_0 / freq) + assert np.isclose(zbfdata.wavelength / unitscaling, td.C_0 / freq) else: assert np.isclose( - zbfdata["wavelength"] / unitscaling, + zbfdata.wavelength / unitscaling, td.C_0 / np.mean(field_data.monitor.freqs), ) - assert zbfdata["nx"] == n_x - assert zbfdata["ny"] == n_y + assert zbfdata.nx == n_x + assert zbfdata.ny == n_y # check that fields are close - assert np.allclose(ex.values, zbfdata["Ex"]) - assert np.allclose(ey.values, zbfdata["Ey"]) + assert np.allclose(ex.values, zbfdata.Ex) + assert np.allclose(ey.values, zbfdata.Ey) @pytest.mark.parametrize("mode_index", [0, 1]) def test_tozbf_modedata(self, tmp_path, mode_data, mode_index): @@ -972,11 +972,11 @@ def test_tozbf_modedata(self, tmp_path, mode_data, mode_index): n_y=32, units="mm", ) - zbfdata = read_zbf(zbf_filename) + zbfdata = ZBFData.read_zbf(zbf_filename) # check that fields are close - assert np.allclose(ex.values, zbfdata["Ex"]) - assert np.allclose(ey.values, zbfdata["Ey"]) + assert np.allclose(ex.values, zbfdata.Ex) + assert np.allclose(ey.values, zbfdata.Ey) def test_tozbf_modedata_fails(self, tmp_path, mode_data): """Asserts that Modedata.to_zbf() fails if mode_index is not specified""" @@ -1041,7 +1041,7 @@ def test_from_zbf(self, tmp_path, field_data): "dim1,dim2", [("x", "x"), ("y", "y"), ("z", "z"), ("1", "2"), ("c", "z")] ) def test_from_zbf_dimsfail(self, tmp_path, field_data, dim1, dim2): - """Tests creating a field dataset from a zbf""" + """Tests fail cases when the dimensions to populate are wrong.""" zbf_filename = tmp_path / "testzbf.zbf" # write to zbf and then load it back in _, _ = field_data.to_zbf( @@ -1051,7 +1051,6 @@ def test_from_zbf_dimsfail(self, tmp_path, field_data, dim1, dim2): n_y=32, units="mm", ) - - # create a field dataset from the zbf file + # this should fail with pytest.raises(ValueError) as e: _ = td.FieldDataset.from_zbf(filename=zbf_filename, dim1=dim1, dim2=dim2) diff --git a/tidy3d/components/data/dataset.py b/tidy3d/components/data/dataset.py index bb341d102..89dac6fd7 100644 --- a/tidy3d/components/data/dataset.py +++ b/tidy3d/components/data/dataset.py @@ -28,7 +28,7 @@ TimeDataArray, TriangleMeshDataArray, ) -from .zbf import read_zbf +from .zbf import ZBFData DEFAULT_MAX_SAMPLES_PER_STEP = 10_000 DEFAULT_MAX_CELLS_PER_STEP = 10_000 @@ -258,40 +258,42 @@ class FieldDataset(ElectromagneticFieldDataset): description="Spatial distribution of the z-component of the magnetic field.", ) - def from_zbf(filename: str, dim1: xyz, dim2: xyz): - """Creates a :class:`.FieldDataset` from a Zemax Beam File (zbf) + def from_zbf(filename: str, dim1: xyz, dim2: xyz) -> FieldDataset: + """Creates a :class:`.FieldDataset` from a Zemax Beam File (.zbf). Parameters ---------- filename: str - Specification of the source time-dependence. + The file name of the .zbf file to read. dim1: xyz - Tangential field component to map the x-axis of the zbf data to. - eg. `dim1 = "z"` sets `FieldDataset.Ez` to Ex of the zbf data. + Tangential field component to map the x-dimension of the zbf data to. + eg. ``dim1 = "z"`` sets ``FieldDataset.Ez`` to ``Ex`` of the zbf data. dim2: xyz - Tangential field component to map the y-axis of the zbf data to. + Tangential field component to map the y-dimension of the zbf data to. + eg. ``dim2 = "z"`` sets ``FieldDataset.Ez`` to ``Ey`` of the zbf data. Returns ------- :class:`.FieldDataset` A :class:`.FieldDataset` object with two tangential E field components populated by zbf data. + + See Also + -------- + :class:`.ZBFData': + A class containing data read in from a .zbf file. """ log.warning( - "`FieldDataset.from_zbf()` is currently an experimental feature." - " If any issues are encountered, please contact Flexcompute support https://www.flexcompute.com/tidy3d/technical-support/" + "'FieldDataset.from_zbf()' is currently an experimental feature." + " If any issues are encountered, please contact Flexcompute support 'https://www.flexcompute.com/tidy3d/technical-support/'" ) if dim1 not in get_args(xyz): - raise ValueError( - f'`dim1` = "{dim1}" not allowed, must be one of `"x"`, `"y"`, or `"z"`.' - ) + raise ValueError(f"'dim1' = '{dim1}' is not allowed, must be one of 'x', 'y', or 'z'.") if dim2 not in get_args(xyz): - raise ValueError( - f'`dim2` = "{dim2}" not allowed, must be one of `"x"`, `"y"`, or `"z"`.' - ) + raise ValueError(f"'dim2' = '{dim2}' is not allowed, must be one of 'x', 'y', or 'z'.") if dim1 == dim2: - raise ValueError("`dim1` and `dim2` must be different.") + raise ValueError("'dim1' and 'dim2' must be different.") # get the third dimension dim3 = list(set(get_args(xyz)) - {dim1, dim2})[0] @@ -299,16 +301,16 @@ def from_zbf(filename: str, dim1: xyz, dim2: xyz): dim2expand = dims[dim3] # this is for expanding E field arrays # load zbf data - zbfdata = read_zbf(filename) + zbfdata = ZBFData.read_zbf(filename) # Grab E fields, dimensions, wavelength - edim1 = zbfdata["Ex"] - edim2 = zbfdata["Ey"] - n1 = zbfdata["nx"] - n2 = zbfdata["ny"] - d1 = zbfdata["dx"] / UnitScaling[zbfdata["unit"]] - d2 = zbfdata["dy"] / UnitScaling[zbfdata["unit"]] - wavelength = zbfdata["wavelength"] / UnitScaling[zbfdata["unit"]] + edim1 = zbfdata.Ex + edim2 = zbfdata.Ey + n1 = zbfdata.nx + n2 = zbfdata.ny + d1 = zbfdata.dx / UnitScaling[zbfdata.unit] + d2 = zbfdata.dy / UnitScaling[zbfdata.unit] + wavelength = zbfdata.wavelength / UnitScaling[zbfdata.unit] # make scalar field data arrays len1 = d1 * (n1 - 1) diff --git a/tidy3d/components/data/zbf.py b/tidy3d/components/data/zbf.py index 56707b5d5..9dee4464b 100644 --- a/tidy3d/components/data/zbf.py +++ b/tidy3d/components/data/zbf.py @@ -1,128 +1,170 @@ """ZBF utilities""" +from __future__ import annotations + from struct import unpack import numpy as np +import pydantic.v1 as pd + +from ..base import Tidy3dBaseModel -def read_zbf(filename: str): - """Reads a Zemax Beam File (zbf) - - Parameters - ---------- - filename : string - the filename of the beam file to read - - Returns - ------- - A dictionary with the following fields: - version : int - File format version number. - nx : int - Number of samples in the x direction. - ny : int - Number of samples in the y direction. - ispol : boolean - True if the beam is polarized. - unit : str - Spatial units, either "mm", "cm", "in", or "m". - dx : float - x grid spacing - dy : float - y grid spacing - zposition_x : float - Pilot beam z position with respect to the waist in x. - zposition_y : float - Pilot beam z position with respect to the waist in y. - rayleigh_x : float - Pilot beam Rayleigh distance in x. - rayleigh_y : float - Pilot beam Rayleigh distance in y. - waist_x : float - Pilot beam waist in x. - waist_y : float - Pilot beam waist in y. - wavelength : float - The wavelength of the beam. - background_refractive_index : float - The index of refraction of the material that the beam exists in. - receiver_eff : double - The receiver efficiency. Zero if fiber coupling is not computed. - system_eff : double - The system efficiency. Zero if fiber coupling is not computed. - Ex : np.ndarray - Complex-valued electric field, x-component. - Ey : np.ndarray - Complex-valued electric field, x-component. +class ZBFData(Tidy3dBaseModel): """ + Contains data read in from a .zbf file + """ + + version: int = pd.Field(..., title="Version", description="File format version number.") + nx: int = pd.Field( + ..., title="Samples in X", description="Number of samples in the x direction." + ) + ny: int = pd.Field( + ..., title="Samples in Y", description="Number of samples in the y direction." + ) + ispol: bool = pd.Field( + ..., + title="Is Polarized", + description="``True`` if the beam is polarized, ``False`` otherwise.", + ) + unit: str = pd.Field( + ..., title="Spatial Units", description="Spatial units, either 'mm', 'cm', 'in', or 'm'." + ) + dx: float = pd.Field(..., title="Grid Spacing, X", description="Grid spacing in x.") + dy: float = pd.Field(..., title="Grid Spacing, Y", description="Grid spacing in y.") + zposition_x: float = pd.Field( + ..., + title="Z Position, X Direction", + description="The pilot beam z position with respect to the pilot beam waist, x direction.", + ) + zposition_y: float = pd.Field( + ..., + title="Z Position, Y Direction", + description="The pilot beam z position with respect to the pilot beam waist, y direction.", + ) + rayleigh_x: float = pd.Field( + ..., + title="Rayleigh Distance, X Direction", + description="The pilot beam Rayleigh distance in the x direction.", + ) + rayleigh_y: float = pd.Field( + ..., + title="Rayleigh Distance, Y Direction", + description="The pilot beam Rayleigh distance in the y direction.", + ) + waist_x: float = pd.Field( + ..., title="Beam Waist, X", description="The pilot beam waist in the x direction." + ) + waist_y: float = pd.Field( + ..., title="Beam Waist, Y", description="The pilot beam waist in the y direction." + ) + wavelength: float = pd.Field(..., title="Wavelength", description="The wavelength of the beam.") + background_refractive_index: float = pd.Field( + ..., + title="Background Refractive Index", + description="The index of refraction in the current medium.", + ) + receiver_eff: float = pd.Field( + ..., + title="Receiver Efficiency", + description="The receiver efficiency. Zero if fiber coupling is not computed.", + ) + system_eff: float = pd.Field( + ..., + title="System Efficiency", + description="The system efficiency. Zero if fiber coupling is not computed.", + ) + Ex: np.ndarray = pd.Field( + ..., + title="Electric Field, X Component", + description="Complex-valued electric field, x component.", + ) + Ey: np.ndarray = pd.Field( + ..., + title="Electric Field, Y Component", + description="Complex-valued electric field, y component.", + ) + + def read_zbf(filename: str) -> ZBFData: + """Reads a Zemax Beam File (.zbf) + + Parameters + ---------- + filename : str + The file name of the .zbf file to read. + + Returns + ------- + :class:`.ZBFData` + """ - # Read the zbf file - with open(filename, "rb") as f: - # Load the header - version, nx, ny, ispol, units = unpack("<5I", f.read(20)) - f.read(16) # unused values - ( - dx, - dy, - zposition_x, - rayleigh_x, - waist_x, - zposition_y, - rayleigh_y, - waist_y, - wavelength, - background_refractive_index, - receiver_eff, - system_eff, - ) = unpack("<12d", f.read(96)) - f.read(64) # unused values - - # read E field - nsamps = 2 * nx * ny - rawx = list(unpack(f"<{nsamps}d", f.read(8 * nsamps))) + # Read the zbf file + with open(filename, "rb") as f: + # Load the header + version, nx, ny, ispol, units = unpack("<5I", f.read(20)) + f.read(16) # unused values + ( + dx, + dy, + zposition_x, + rayleigh_x, + waist_x, + zposition_y, + rayleigh_y, + waist_y, + wavelength, + background_refractive_index, + receiver_eff, + system_eff, + ) = unpack("<12d", f.read(96)) + f.read(64) # unused values + + # read E field + nsamps = 2 * nx * ny + rawx = list(unpack(f"<{nsamps}d", f.read(8 * nsamps))) + if ispol: + rawy = list(unpack(f"<{nsamps}d", f.read(8 * nsamps))) + + # convert unit key to unit string + map_units = {0: "mm", 1: "cm", 2: "in", 3: "m"} + try: + unit = map_units[units] + except KeyError: + raise KeyError( + f"Invalid units specified in the zbf file (expected `0`, `1`, `2`, or `3`, got `{units}`)." + ) + + # load E field + Ex_real = np.asarray(rawx[0::2]).reshape(nx, ny) + Ex_imag = np.asarray(rawx[1::2]).reshape(nx, ny) if ispol: - rawy = list(unpack(f"<{nsamps}d", f.read(8 * nsamps))) - - # convert unit key to unit string - map_units = {0: "mm", 1: "cm", 2: "in", 3: "m"} - try: - unit = map_units[units] - except KeyError: - raise KeyError( - f"Invalid units specified in the zbf file (expected `0`, `1`, `2`, or `3`, got `{units}`)." - ) + Ey_real = np.asarray(rawy[0::2]).reshape(nx, ny) + Ey_imag = np.asarray(rawy[1::2]).reshape(nx, ny) + else: + Ey_real = np.zeros((nx, ny)) + Ey_imag = np.zeros((nx, ny)) + + Ex = Ex_real + 1j * Ex_imag + Ey = Ey_real + 1j * Ey_imag - # load E field - Ex_real = np.asarray(rawx[0::2]).reshape(nx, ny) - Ex_imag = np.asarray(rawx[1::2]).reshape(nx, ny) - if ispol: - Ey_real = np.asarray(rawy[0::2]).reshape(nx, ny) - Ey_imag = np.asarray(rawy[1::2]).reshape(nx, ny) - else: - Ey_real = np.zeros((nx, ny)) - Ey_imag = np.zeros((nx, ny)) - - Ex = Ex_real + 1j * Ex_imag - Ey = Ey_real + 1j * Ey_imag - - return { - "version": version, - "nx": nx, - "ny": ny, - "ispol": ispol, - "unit": unit, - "dx": dx, - "dy": dy, - "zposition_x": zposition_x, - "zposition_y": zposition_y, - "rayleigh_x": rayleigh_x, - "rayleigh_y": rayleigh_y, - "waist_x": waist_x, - "waist_y": waist_y, - "wavelength": wavelength, - "background_refractive_index": background_refractive_index, - "receiver_eff": receiver_eff, - "system_eff": system_eff, - "Ex": Ex, - "Ey": Ey, - } + return ZBFData( + version=version, + nx=nx, + ny=ny, + ispol=ispol, + unit=unit, + dx=dx, + dy=dy, + zposition_x=zposition_x, + zposition_y=zposition_y, + rayleigh_x=rayleigh_x, + rayleigh_y=rayleigh_y, + waist_x=waist_x, + waist_y=waist_y, + wavelength=wavelength, + background_refractive_index=background_refractive_index, + receiver_eff=receiver_eff, + system_eff=system_eff, + Ex=Ex, + Ey=Ey, + ) From 6f4455b9ea9cb9acf7056e22c3a09b71055b854c Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Mon, 5 May 2025 11:35:42 -0400 Subject: [PATCH 24/28] update docstring and console print formatting --- tidy3d/components/data/dataset.py | 4 ++-- tidy3d/components/data/monitor_data.py | 32 +++++++++++++------------- tidy3d/components/data/zbf.py | 6 ++--- 3 files changed, 21 insertions(+), 21 deletions(-) diff --git a/tidy3d/components/data/dataset.py b/tidy3d/components/data/dataset.py index 89dac6fd7..df6ed974d 100644 --- a/tidy3d/components/data/dataset.py +++ b/tidy3d/components/data/dataset.py @@ -259,7 +259,7 @@ class FieldDataset(ElectromagneticFieldDataset): ) def from_zbf(filename: str, dim1: xyz, dim2: xyz) -> FieldDataset: - """Creates a :class:`.FieldDataset` from a Zemax Beam File (.zbf). + """Creates a :class:`.FieldDataset` from a Zemax Beam File (`.zbf`). Parameters ---------- @@ -281,7 +281,7 @@ def from_zbf(filename: str, dim1: xyz, dim2: xyz) -> FieldDataset: See Also -------- :class:`.ZBFData': - A class containing data read in from a .zbf file. + A class containing data read in from a `.zbf` file. """ log.warning( "'FieldDataset.from_zbf()' is currently an experimental feature." diff --git a/tidy3d/components/data/monitor_data.py b/tidy3d/components/data/monitor_data.py index cc5e1c938..635bac473 100644 --- a/tidy3d/components/data/monitor_data.py +++ b/tidy3d/components/data/monitor_data.py @@ -1084,7 +1084,7 @@ def to_zbf( rec_efficiency: float = 0, sys_efficiency: float = 0, ) -> Tuple[ScalarFieldDataArray, ScalarFieldDataArray]: - """For a 2D monitor, export the fields to a Zemax Beam File (zbf). + """For a 2D monitor, export the fields to a Zemax Beam File (`.zbf`). The mode area is used to approximate the beam waist, which is only valid if the beam profile approximates a Gaussian beam. @@ -1099,20 +1099,20 @@ def to_zbf( Number of field samples along x. Must be a power of 2, between 2^5 and 2^13 inclusive per Zemax's requirements. n_y : int Number of field samples along y. Must be a power of 2, between 2^5 and 2^13 inclusive per Zemax's requirements. - units : `UnitsZBF` - Spatial units used for the zbf file. Options are `"mm"`, `"cm"`, `"in"`, or `"m"`. + units : ``UnitsZBF`` + Spatial units used for the zbf file. Options are ``"mm"``, ``"cm"``, ``"in"``, or ``"m"``. freq : float Field frequency selection. If ``None``, the average of the recorded frequencies is used. mode_index : Optional[int] - For :class:`ModeData`, choose which mode to save. + For :class:`.ModeData`, choose which mode to save. r_x : float Pilot beam Rayleigh distance in x, um. Defaults to 0. r_y : float Pilot beam Rayleigh distance in y, um. Defaults to 0. z_x : float - Pilot beam z position w.r.t. the waist in x, um. Defaults to 0. + Pilot beam z position with respect to the waist in x, um. Defaults to 0. z_y : float - Pilot beam z position w.r.t. the waist in y, um. Defaults to 0. + Pilot beam z position with respect to the waist in y, um. Defaults to 0. rec_efficiency : float Receiver efficiency, zero if fiber coupling is not computed. Defaults to 0. sys_efficiency : float @@ -1120,29 +1120,29 @@ def to_zbf( Returns ------- - Tuple[:class:`ScalarFieldDataArray`,:class:`ScalarFieldDataArray`] - The two E field components being exported to zbf. + Tuple[:class:`.ScalarFieldDataArray`,:class:`.ScalarFieldDataArray`] + The two E field components being exported to ``.zbf``. """ log.warning( - "`FieldData.to_zbf()` is currently an experimental feature." - " If any issues are encountered, please contact Flexcompute support https://www.flexcompute.com/tidy3d/technical-support/" + "'FieldData.to_zbf()' is currently an experimental feature." + " If any issues are encountered, please contact Flexcompute support 'https://www.flexcompute.com/tidy3d/technical-support/'" ) # Check that appropriate units are used if units not in get_args(UnitsZBF): - raise ValueError("`units` must be either `'mm'`, `'cm'`, `'in'`, or `'m'`.") + raise ValueError("'units' must be either 'mm', 'cm', 'in', or 'm'.") # Check that requirements are met for n_x and n_y # n_x and n_y must be powers of 2 if (n_x & (n_x - 1)) != 0: - raise ValueError("`n_x` must be a power of 2.") + raise ValueError("'n_x' must be a power of 2.") if (n_y & (n_y - 1)) != 0: - raise ValueError("`n_y` must be a power of 2.") + raise ValueError("'n_y' must be a power of 2.") # 32 <= n_x and n_y <= 2^13 if n_x < 32 or n_x > 2**13: - raise ValueError("`n_x` must be between 2^5 and 2^13, inclusive.") + raise ValueError("'n_x' must be between 2^5 and 2^13, inclusive.") if n_y < 32 or n_y > 2**13: - raise ValueError("`n_y` must be between 2^5 and 2^13, inclusive.") + raise ValueError("'n_y' must be between 2^5 and 2^13, inclusive.") # Mode area calculation ensures all E components are present mode_area = self.mode_area @@ -1165,7 +1165,7 @@ def to_zbf( # If the data is ModeData, choose one of the modes to save if "mode_index" in e_x.coords: if mode_index is None: - raise ValueError("`mode_index` is required for `ModeData.to_zbf()`") + raise ValueError("'mode_index' is required for 'ModeData.to_zbf()'") mode_area = mode_area.isel(mode_index=mode_index, drop=True) e_x = e_x.isel(mode_index=mode_index, drop=True) e_y = e_y.isel(mode_index=mode_index, drop=True) diff --git a/tidy3d/components/data/zbf.py b/tidy3d/components/data/zbf.py index 9dee4464b..5ec823cfa 100644 --- a/tidy3d/components/data/zbf.py +++ b/tidy3d/components/data/zbf.py @@ -86,12 +86,12 @@ class ZBFData(Tidy3dBaseModel): ) def read_zbf(filename: str) -> ZBFData: - """Reads a Zemax Beam File (.zbf) + """Reads a Zemax Beam File (`.zbf`) Parameters ---------- filename : str - The file name of the .zbf file to read. + The file name of the `.zbf` file to read. Returns ------- @@ -131,7 +131,7 @@ def read_zbf(filename: str) -> ZBFData: unit = map_units[units] except KeyError: raise KeyError( - f"Invalid units specified in the zbf file (expected `0`, `1`, `2`, or `3`, got `{units}`)." + f"Invalid units specified in the zbf file (expected '0', '1', '2', or '3', got '{units}')." ) # load E field From 3ea8adb91f9e1df1111272d3ae7c591bf234c1a3 Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Mon, 5 May 2025 12:08:56 -0400 Subject: [PATCH 25/28] add default values for background index and ny, ny --- tidy3d/components/data/monitor_data.py | 85 +++++++++++++++----------- 1 file changed, 49 insertions(+), 36 deletions(-) diff --git a/tidy3d/components/data/monitor_data.py b/tidy3d/components/data/monitor_data.py index 635bac473..50fb6d14d 100644 --- a/tidy3d/components/data/monitor_data.py +++ b/tidy3d/components/data/monitor_data.py @@ -1071,10 +1071,10 @@ def translated_copy(self, vector: Coordinate) -> ElectromagneticFieldData: def to_zbf( self, fname: str, - background_refractive_index: float, - n_x: int, - n_y: int, units: UnitsZBF, + background_refractive_index: float = 1, + n_x: Optional[int] = None, + n_y: Optional[int] = None, freq: Optional[float] = None, mode_index: Optional[int] = None, r_x: float = 0, @@ -1093,30 +1093,34 @@ def to_zbf( ---------- fname : str Full path to the .zbf file to be written. - background_refractive_index : float - Refractive index of the medium surrounding the monitor. - n_x : int - Number of field samples along x. Must be a power of 2, between 2^5 and 2^13 inclusive per Zemax's requirements. - n_y : int - Number of field samples along y. Must be a power of 2, between 2^5 and 2^13 inclusive per Zemax's requirements. units : ``UnitsZBF`` Spatial units used for the zbf file. Options are ``"mm"``, ``"cm"``, ``"in"``, or ``"m"``. - freq : float + background_refractive_index : float = 1 + Refractive index of the medium surrounding the monitor. Defaults to ``1``. + n_x : Optional[int] = None + Number of field samples along x. + Must be a power of 2, between 2^5 and 2^13 inclusive per Zemax's requirements. + Defaults to ``None``, in which case a value is chosen for the user depending on the coordinates in the field data. + n_y : Optional[int] = None + Number of field samples along y. + Must be a power of 2, between 2^5 and 2^13 inclusive per Zemax's requirements. + Defaults to ``None``, in which case a value is chosen for the user depending on the coordinates in the field data. + freq : Optional[float] = None Field frequency selection. If ``None``, the average of the recorded frequencies is used. - mode_index : Optional[int] + mode_index : Optional[int] = None For :class:`.ModeData`, choose which mode to save. - r_x : float - Pilot beam Rayleigh distance in x, um. Defaults to 0. - r_y : float - Pilot beam Rayleigh distance in y, um. Defaults to 0. - z_x : float - Pilot beam z position with respect to the waist in x, um. Defaults to 0. - z_y : float - Pilot beam z position with respect to the waist in y, um. Defaults to 0. - rec_efficiency : float - Receiver efficiency, zero if fiber coupling is not computed. Defaults to 0. - sys_efficiency : float - System efficiency, zero if fiber coupling is not computed. Defaults to 0. + r_x : float = 0 + Pilot beam Rayleigh distance in x, um. Defaults to ``0``. + r_y : float = 0 + Pilot beam Rayleigh distance in y, um. Defaults to ``0``. + z_x : float = 0 + Pilot beam z position with respect to the waist in x, um. Defaults to ``0``. + z_y : float = 0 + Pilot beam z position with respect to the waist in y, um. Defaults to ``0``. + rec_efficiency : float = 0 + Receiver efficiency, zero if fiber coupling is not computed. Defaults to ``0``. + sys_efficiency : float = 0 + System efficiency, zero if fiber coupling is not computed. Defaults to ``0``. Returns ------- @@ -1132,18 +1136,6 @@ def to_zbf( if units not in get_args(UnitsZBF): raise ValueError("'units' must be either 'mm', 'cm', 'in', or 'm'.") - # Check that requirements are met for n_x and n_y - # n_x and n_y must be powers of 2 - if (n_x & (n_x - 1)) != 0: - raise ValueError("'n_x' must be a power of 2.") - if (n_y & (n_y - 1)) != 0: - raise ValueError("'n_y' must be a power of 2.") - # 32 <= n_x and n_y <= 2^13 - if n_x < 32 or n_x > 2**13: - raise ValueError("'n_x' must be between 2^5 and 2^13, inclusive.") - if n_y < 32 or n_y > 2**13: - raise ValueError("'n_y' must be between 2^5 and 2^13, inclusive.") - # Mode area calculation ensures all E components are present mode_area = self.mode_area dim1, dim2 = self._tangential_dims @@ -1156,6 +1148,9 @@ def to_zbf( # Use the mean frequency if freq is not specified if freq is None: + log.warning( + "'freq' was not specified for 'FieldData.to_zbf()'. Defaulting to the mean frequency of the dataset." + ) freq = np.mean(e_x.coords["f"].values) mode_area = mode_area.interp(f=freq) @@ -1190,7 +1185,25 @@ def to_zbf( z_x *= unit_scaling z_y *= unit_scaling - # Interpolating coordinatesΠ+ # defaults for n_x and n_y + if n_x is None: + n_x = 2 ** min(13, max(5, int(np.log2(x.size) + 1))) + if n_y is None: + n_y = 2 ** min(13, max(5, int(np.log2(y.size) + 1))) + + # Check that requirements are met for n_x and n_y + # n_x and n_y must be powers of 2 + if (n_x & (n_x - 1)) != 0: + raise ValueError("'n_x' must be a power of 2.") + if (n_y & (n_y - 1)) != 0: + raise ValueError("'n_y' must be a power of 2.") + # 32 <= n_x and n_y <= 2^13 + if n_x < 32 or n_x > 2**13: + raise ValueError("'n_x' must be between 2^5 and 2^13, inclusive.") + if n_y < 32 or n_y > 2**13: + raise ValueError("'n_y' must be between 2^5 and 2^13, inclusive.") + + # Interpolating coordinates x = np.linspace(x.min(), x.max(), n_x) y = np.linspace(y.min(), y.max(), n_y) From 71b61f13e3cab9760b444375652e139fe9b17087 Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Mon, 5 May 2025 12:20:30 -0400 Subject: [PATCH 26/28] add warnings that default index, nx, ny are being used --- tidy3d/components/data/monitor_data.py | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/tidy3d/components/data/monitor_data.py b/tidy3d/components/data/monitor_data.py index 50fb6d14d..70a2757bd 100644 --- a/tidy3d/components/data/monitor_data.py +++ b/tidy3d/components/data/monitor_data.py @@ -1072,7 +1072,7 @@ def to_zbf( self, fname: str, units: UnitsZBF, - background_refractive_index: float = 1, + background_refractive_index: Optional[float] = None, n_x: Optional[int] = None, n_y: Optional[int] = None, freq: Optional[float] = None, @@ -1095,7 +1095,7 @@ def to_zbf( Full path to the .zbf file to be written. units : ``UnitsZBF`` Spatial units used for the zbf file. Options are ``"mm"``, ``"cm"``, ``"in"``, or ``"m"``. - background_refractive_index : float = 1 + background_refractive_index : Optional[float] = None Refractive index of the medium surrounding the monitor. Defaults to ``1``. n_x : Optional[int] = None Number of field samples along x. @@ -1132,6 +1132,13 @@ def to_zbf( " If any issues are encountered, please contact Flexcompute support 'https://www.flexcompute.com/tidy3d/technical-support/'" ) + # Sets default background refractive index, issues warning if it was not set. + if background_refractive_index is None: + background_refractive_index = 1 + log.warning( + "'background_refractive_index' was not set for 'FieldData.to_zbf()'. Defaulting to 1." + ) + # Check that appropriate units are used if units not in get_args(UnitsZBF): raise ValueError("'units' must be either 'mm', 'cm', 'in', or 'm'.") @@ -1188,8 +1195,14 @@ def to_zbf( # defaults for n_x and n_y if n_x is None: n_x = 2 ** min(13, max(5, int(np.log2(x.size) + 1))) + log.warning( + f"'n_x' was not specified for 'FieldData.to_zbf()'. Defaulting to 'n_x' = {n_x}." + ) if n_y is None: n_y = 2 ** min(13, max(5, int(np.log2(y.size) + 1))) + log.warning( + f"'n_y' was not specified for 'FieldData.to_zbf()'. Defaulting to 'n_y' = {n_y}." + ) # Check that requirements are met for n_x and n_y # n_x and n_y must be powers of 2 From 38a4ce545a2d91fbc126731c7d29508f4d80ff63 Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Mon, 5 May 2025 12:26:33 -0400 Subject: [PATCH 27/28] update docstring for UnitScaling so its more clear --- tidy3d/constants.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tidy3d/constants.py b/tidy3d/constants.py index 2610887e1..73c93b6cd 100644 --- a/tidy3d/constants.py +++ b/tidy3d/constants.py @@ -241,4 +241,4 @@ "in": 1.0 / 25400, } ) -"""Immutable dictionary for converting a unit specification to a scaling factor.""" +"""Immutable dictionary for converting microns to another spatial unit, eg. nm = um * UnitScaling["nm"].""" From 7533e425c22dfcd540a9d4514c5ca1cc3a2f7cd8 Mon Sep 17 00:00:00 2001 From: Bohan Zhang Date: Mon, 5 May 2025 17:55:03 -0400 Subject: [PATCH 28/28] set defaults for units, default background index to 1 --- tidy3d/components/data/monitor_data.py | 16 +++++----------- 1 file changed, 5 insertions(+), 11 deletions(-) diff --git a/tidy3d/components/data/monitor_data.py b/tidy3d/components/data/monitor_data.py index 70a2757bd..43b52bb42 100644 --- a/tidy3d/components/data/monitor_data.py +++ b/tidy3d/components/data/monitor_data.py @@ -1071,8 +1071,8 @@ def translated_copy(self, vector: Coordinate) -> ElectromagneticFieldData: def to_zbf( self, fname: str, - units: UnitsZBF, - background_refractive_index: Optional[float] = None, + units: UnitsZBF = "mm", + background_refractive_index: float = 1, n_x: Optional[int] = None, n_y: Optional[int] = None, freq: Optional[float] = None, @@ -1093,9 +1093,10 @@ def to_zbf( ---------- fname : str Full path to the .zbf file to be written. - units : ``UnitsZBF`` + units : UnitsZBF = "mm" Spatial units used for the zbf file. Options are ``"mm"``, ``"cm"``, ``"in"``, or ``"m"``. - background_refractive_index : Optional[float] = None + Defaults to ``"mm"``. + background_refractive_index : float = 1 Refractive index of the medium surrounding the monitor. Defaults to ``1``. n_x : Optional[int] = None Number of field samples along x. @@ -1132,13 +1133,6 @@ def to_zbf( " If any issues are encountered, please contact Flexcompute support 'https://www.flexcompute.com/tidy3d/technical-support/'" ) - # Sets default background refractive index, issues warning if it was not set. - if background_refractive_index is None: - background_refractive_index = 1 - log.warning( - "'background_refractive_index' was not set for 'FieldData.to_zbf()'. Defaulting to 1." - ) - # Check that appropriate units are used if units not in get_args(UnitsZBF): raise ValueError("'units' must be either 'mm', 'cm', 'in', or 'm'.")