WI add tests

nidn/fdtd/detectors.py

@@ -14,7 +14,7 @@
 # relative
 from .grid import Grid
 from .backend import backend as bd
-from .constants import X, Y, Z
+from .constants import FDTD_X, FDTD_Y, FDTD_Z
 
 ## Detector
 class LineDetector:
@@ -441,20 +441,22 @@ def single_point_current(self, px, py, pz):
         # option for unit factor here? Units will get very complicated otherwise
 
         current_vector_1 = (
-            (self.grid.H[px, py - 1, pz, X]) - (self.grid.H[px, py, pz, X])
+            (self.grid.H[px, py - 1, pz, FDTD_X]) - (self.grid.H[px, py, pz, FDTD_X])
         ) * self.grid.grid_spacing
         current_vector_2 = (
-            (self.grid.H[px, py, pz, Y]) - (self.grid.H[px - 1, py, pz, Y])
+            (self.grid.H[px, py, pz, FDTD_Y]) - (self.grid.H[px - 1, py, pz, FDTD_Y])
         ) * self.grid.grid_spacing
 
         current_1 = current_vector_1 + current_vector_2
         # current_1 = float(current.cpu())
 
         current_vector_1 = (
-            (self.grid.H[px, py - 1, pz - 1, X]) - (self.grid.H[px, py, pz - 1, X])
+            (self.grid.H[px, py - 1, pz - 1, FDTD_X])
+            - (self.grid.H[px, py, pz - 1, FDTD_X])
         ) * self.grid.grid_spacing
         current_vector_2 += (
-            (self.grid.H[px, py, pz - 1, Y]) - (self.grid.H[px - 1, py, pz - 1, Y])
+            (self.grid.H[px, py, pz - 1, FDTD_Y])
+            - (self.grid.H[px - 1, py, pz - 1, FDTD_Y])
         ) * self.grid.grid_spacing
         # current_2 = float(current_2.cpu())
         current_2 = current_vector_1 + current_vector_2

nidn/fdtd_integration/calculate_transmission_reflection_coefficients.py

@@ -1,6 +1,8 @@
 from dotmap import DotMap
 from torch.fft import rfft, rfftfreq
-from torch import sqrt, tensor
+import torch
+from loguru import logger
+
 
 from nidn.fdtd_integration.constants import FDTD_GRID_SCALE
 
@@ -49,6 +51,20 @@ def calculate_transmission_reflection_coefficients(
         reflection_coefficient = _fft(true_reflection, cfg) / _fft(
             reflection_signals[0], cfg
         )
+
+    if not (transmission_coefficient >= 0 and transmission_coefficient <= 1):
+        raise ValueError(
+            f"The transmission coefficient is outsied of the possible range between 0 and 1"
+        )
+
+    if not (reflection_coefficient >= 0 and reflection_coefficient <= 1):
+        raise ValueError(
+            f"The reflection coefficient is outsied of the possible range between 0 and 1"
+        )
+    if transmission_coefficient + reflection_coefficient > 1:
+        logger.warning(
+            f"The sum of the transmission and reflection coefficient is greater than 1, which is physically impossible"
+        )
     return transmission_coefficient, reflection_coefficient
 
 

nidn/fdtd_integration/compute_spectrum_fdtd.py

@@ -1,6 +1,6 @@
 from dotmap import DotMap
 from tqdm import tqdm
-from torch import sqrt, tensor
+import torch
 from loguru import logger
 
 from ..trcwa.get_frequency_points import get_frequency_points
@@ -23,9 +23,10 @@ def compute_spectrum_fdtd(permittivity, cfg: DotMap):
     transmission_spectrum = []
     reflection_spectrum = []
     physical_wavelengths, norm_freq = get_frequency_points(cfg)
-    logger.debug("Wavelenghts in spectrum")
+    logger.debug("Wavelenghts in spectrum : ")
     logger.debug(physical_wavelengths)
-
+    logger.debug("Number of layers: ")
+    logger.debug(len(permittivity[0, 0, :, 0]))
     # For each wavelength, calculate transmission and reflection coefficents
     disable_progress_bar = logger._core.min_level >= 20
     for i in tqdm(range(len(physical_wavelengths)), disable=disable_progress_bar):
@@ -59,8 +60,6 @@ def compute_spectrum_fdtd(permittivity, cfg: DotMap):
         )
         transmission_signal.append(transmission_material)
         reflection_signal.append(reflection_material)
-        time = [i for i in range(len(transmission_signal[0]))]
-
         # Calculate transmission and reflection coefficients,
         # by using the signals from the free space simulation and the material simulation
         (
@@ -77,7 +76,7 @@ def compute_spectrum_fdtd(permittivity, cfg: DotMap):
     logger.debug("Reflection spectrum")
     logger.debug(reflection_spectrum)
 
-    return transmission_spectrum, reflection_spectrum
+    return torch.tensor(transmission_spectrum), torch.tensor(reflection_spectrum)
 
 
 def _get_detector_values(transmission_detector, reflection_detector):
@@ -120,7 +119,7 @@ def _get_abs_value_from_3D_signal(signal):
     abs_value = []
     for i in range(len(signal)):
         abs_value.append(
-            sqrt(tensor(signal[i][0] ** 2 + signal[i][1] ** 2 + signal[i][2] ** 2))
+            torch.sqrt(signal[i][0] ** 2 + signal[i][1] ** 2 + signal[i][2] ** 2)
         )
     return abs_value
 
@@ -136,7 +135,7 @@ def _average_along_detector(signal):
     """
     avg = []
     for e in signal:
-        s = [tensor(0.0), tensor(0.0), tensor(0.0)]
+        s = [torch.tensor(0.0), torch.tensor(0.0), torch.tensor(0.0)]
         for p in e:
             s[0] += p[0] / len(e)
             s[1] += p[1] / len(e)

nidn/fdtd_integration/init_fdtd.py

@@ -1,6 +1,6 @@
 from dotmap import DotMap
-import fdtd
 
+from ..fdtd import *
 from ..utils.global_constants import EPS_0, SPEED_OF_LIGHT, UNIT_MAGNITUDE
 from .constants import FDTD_GRID_SCALE
 
@@ -17,13 +17,13 @@ def init_fdtd(cfg: DotMap, include_object, wavelength, permittivity):
     Returns:
         fdtd:Grid: Grid with all the added object, ready to be run
     """
-    fdtd.set_backend("torch")
+    set_backend("torch")
     scaling = UNIT_MAGNITUDE / (cfg.physical_wavelength_range[0] * FDTD_GRID_SCALE)
     x_grid_size = int(
         scaling * (cfg.FDTD_grid[0] + cfg.N_layers * cfg.PER_LAYER_THICKNESS[0])
     )
     y_grid_size = int(cfg.FDTD_grid[1] * scaling)
-    grid = fdtd.Grid(
+    grid = Grid(
         (x_grid_size, y_grid_size, 1),
         grid_spacing=cfg.physical_wavelength_range[0] * FDTD_GRID_SCALE,
         permittivity=1.0,
@@ -38,16 +38,24 @@ def init_fdtd(cfg: DotMap, include_object, wavelength, permittivity):
                 cfg.FDTD_reflection_detector_x
                 + cfg.N_layers * cfg.PER_LAYER_THICKNESS[0]
             )
+            + cfg.FDTD_min_gridpoints_between_detectors
         ),
         int(cfg.FDTD_reflection_detector_x * scaling),
     )
+    use_pulse = True
+    if cfg.FDTD_pulse_type == "pulse":
+        pass
+    elif cfg.FDTD_pulse_type == "continuous":
+        use_pulse = False
+    else:
+        raise ValueError(f' FDTD_pulse_type must either be "pulse" or "continuous"')
     grid = _add_source(
         grid,
-        int(cfg.FDTD_source[0] * scaling),
-        int(cfg.FDTD_source[1] * scaling),
+        int(cfg.FDTD_source_position[0] * scaling),
+        int(cfg.FDTD_source_position[1] * scaling),
         wavelength / SPEED_OF_LIGHT,
-        cfg.FDTD_use_pulsesource,
-        cfg.FDTD_use_pointsource,
+        use_pulse,
+        cfg.FDTD_source_type,
     )
     if include_object:
         for i in range(cfg.N_layers):
@@ -88,17 +96,17 @@ def _add_boundaries(grid, pml_thickness):
         fdtd.Grid: The grid with the added boundaries
     """
     # Add PML boundary to the left side of the grid
-    grid[0:pml_thickness, :, :] = fdtd.PML(name="pml_xlow")
+    grid[0:pml_thickness, :, :] = PML(name="pml_xlow")
     # Add PML boundary to the right side of the grid
-    grid[-pml_thickness:, :, :] = fdtd.PML(name="pml_xhigh")
+    grid[-pml_thickness:, :, :] = PML(name="pml_xhigh")
     # Add periodic boundaries at both sides in y-direction
-    grid[:, 0, :] = fdtd.PeriodicBoundary(name="ybounds")
+    grid[:, 0, :] = PeriodicBoundary(name="ybounds")
     # Add periodic boundaries on both sides in z-direction. Only applicable for 3D grids
     # grid[:, :, 0] = fdtd.PeriodicBoundary(name="zbounds")
     return grid
 
 
-def _add_source(grid, source_x, source_y, period, use_pulse_source, use_point_source):
+def _add_source(grid, source_x, source_y, period, use_pulse_source, source_type):
     """Add a specified source to the fdtd grid
 
     Args:
@@ -113,17 +121,18 @@ def _add_source(grid, source_x, source_y, period, use_pulse_source, use_point_so
         fdtd.Grid: The grid with the added source
     """
 
-    if use_point_source:
-        grid[source_x, source_y, 0] = fdtd.PointSource(
+    if source_type == "point":
+        grid[source_x, source_y, 0] = PointSource(
             period=period,
-            name="linesource",
+            name="pointsource",
             pulse=use_pulse_source,
             cycle=1,
             hanning_dt=1e-15,
         )
+    elif source_type == "line":
+        grid[source_x, :, 0] = LineSource(period=period, name="linesource")
     else:
-        grid[source_x, :, 0] = fdtd.LineSource(period=period, name="linesource")
-
+        raise ValueError(f'FDTD_source_type must be either "line" or "point"')
     return grid
 
 
@@ -139,8 +148,8 @@ def _add_detectors(grid, transmission_detector_x, reflection_detector_x):
         fdtd.Grid: The grid with the added detectors
     """
 
-    transmission_detector = fdtd.LineDetector(name="t_detector")
-    reflection_detector = fdtd.LineDetector(name="r_detector")
+    transmission_detector = LineDetector(name="t_detector")
+    reflection_detector = LineDetector(name="r_detector")
     grid[transmission_detector_x, :, 0] = transmission_detector
     grid[reflection_detector_x, :, 0] = reflection_detector
     return grid, transmission_detector, reflection_detector
@@ -158,9 +167,9 @@ def _add_object(grid, object_start_x, object_end_x, permittivity, frequency):
     Returns:
         fdtd.Grid: The grid with the added object
     """
-    # Not sure whether the conductivity should be relative or absolute, i.e. if it should be multiplied with EPS_0.
+    # Not sure whether the conductivity should be relative or absolute, i.e. if it should be multiplied with EPS_0. Multiplied with 2pi to get w(angular frequency)?
     # Since the permittivity is set to 1 for the free space grid, I'll leave it at an relative value for now. Also, permittivity for object is relative.
-    grid[object_start_x:object_end_x, :, :] = fdtd.AbsorbingObject(
+    grid[object_start_x:object_end_x, :, :] = AbsorbingObject(
         permittivity=permittivity.real,
         conductivity=permittivity.imag * frequency,
     )

nidn/plots/plot_spectra.py

@@ -3,7 +3,7 @@
 from matplotlib import pyplot as plt
 from matplotlib.ticker import FormatStrFormatter
 from ..utils.convert_units import freq_to_wl, wl_to_phys_wl
-from ..trcwa.compute_spectrum import compute_spectrum
+from ..utils.compute_spectrum import compute_spectrum
 from ..training.model.model_to_eps_grid import model_to_eps_grid
 
 

nidn/tests/calculate_coefficient_test.py

@@ -0,0 +1,26 @@
+import numpy as np
+from ..fdtd_integration.calculate_transmission_reflection_coefficients import (
+    calculate_transmission_reflection_coefficients,
+)
+from ..utils.load_default_cfg import load_default_cfg
+
+
+def test_calculate_coefficient():
+    cfg = load_default_cfg()
+    time_points = [0.002 * np.pi * i for i in range(1000)]
+    signal_a = 2 * np.sin(time_points)
+    signal_b = np.sin(time_points)
+    signal_array = [signal_a, signal_b]
+    (
+        transmission_coefficient_ms,
+        reflection_coefficient_ms,
+    ) = calculate_transmission_reflection_coefficients(
+        signal_array, signal_array, "mean square", cfg
+    )
+    # TODO: Add test for fft, when the method is complete
+    assert transmission_coefficient_ms - 0.25 == 0
+    assert reflection_coefficient_ms - 0.25 == 0
+
+
+if __name__ == "__main__":
+    test_calculate_coefficient()