Add custom laser profile injection to epoch2d/3d

epoch2d/example_decks/custom_laser_profile.deck

@@ -0,0 +1,66 @@
+begin:control
+  nx = 500
+  ny = 500
+
+  # Long enough for the 30fs custom pulse to enter and cross the domain
+  t_end = 60 * femto
+
+  # Size of domain
+  x_min = 0
+  x_max = 20 * micron
+  y_min = -10 * micron
+  y_max =  10 * micron
+
+  stdout_frequency = 10
+end:control
+
+
+begin:boundaries
+  bc_x_min = simple_laser
+  bc_x_max = open
+  bc_y_min = periodic
+  bc_y_max = periodic
+end:boundaries
+
+
+# Custom laser profile injection: the amplitude and (optionally) phase are
+# read from raw binary files instead of an analytic deck expression. Run
+# generate_custom_laser_profile.py first to create amplitude.dat and
+# phase.dat; both must sit alongside this deck in the data directory.
+begin:laser
+  boundary = x_min
+  intensity_w_cm2 = 1e22
+  lambda = 1.0 * micron
+  pol = 0
+
+  use_custom_profile = T
+  use_spatiotemporal_profile = T
+  profile_data_file = amplitude.dat
+
+  n_t = 50
+  n_y = 30
+  y_min = -10.0 * micron
+  y_max =  10.0 * micron
+  t_start = 0.0
+  t_end = 30e-15
+
+  use_phase_from_file = T
+  phase_data_file = phase.dat
+end:laser
+
+
+begin:output
+  name = normal
+
+  # Simulated time between output dumps
+  dt_snapshot = 2 * femto
+
+  # Properties on grid
+  grid = always
+  ex = always
+  ey = always
+  ez = always
+  bx = always
+  by = always
+  bz = always
+end:output

epoch2d/example_decks/generate_custom_laser_profile.py

@@ -0,0 +1,30 @@
+# Generates amplitude.dat and phase.dat for custom_laser_profile.deck.
+# Run this once, with both output files placed in the same data directory
+# as the deck, before starting epoch2d.
+import numpy as np
+
+N_T = 50
+N_Y = 30
+Y_MIN, Y_MAX = -10e-6, 10e-6
+T_START, T_END = 0.0, 30e-15
+
+y = np.linspace(Y_MIN, Y_MAX, N_Y)
+t = np.linspace(T_START, T_END, N_T)
+
+# meshgrid with indexing="ij" gives arr[t_index, y_index]
+Tg, Yg = np.meshgrid(t, y, indexing="ij")
+
+# --- amplitude envelope (values in [0, 1]) ---
+t0, tau = 15e-15, 5e-15
+w0 = 3e-6
+amplitude = np.exp(-((Tg - t0) / tau) ** 2) * np.exp(-(Yg / w0) ** 2)
+
+assert amplitude.shape == (N_T, N_Y)
+amplitude.astype(np.float64).tofile("amplitude.dat")
+
+# --- phase envelope (radians) ---
+phase = np.zeros((N_T, N_Y))  # flat phase for a plane wave
+phase.astype(np.float64).tofile("phase.dat")
+
+print("Wrote amplitude.dat and phase.dat "
+      f"({N_T} x {N_Y} points, {amplitude.nbytes} bytes each)")

epoch2d/src/deck/deck_laser_block.f90

@@ -58,6 +58,18 @@ SUBROUTINE laser_block_start
     working_laser%use_profile_function = .TRUE.
     working_laser%use_omega_function = .FALSE.
 
+    ! Explicitly initialise custom-profile flags and filename for each new
+    ! laser block. Blank profile_data_file triggers the default filename
+    ! ('temporal_spatial_profile.dat' or 'spatial_profile.dat') at load time.
+    working_laser%use_custom_profile = .FALSE.
+    working_laser%use_spatiotemporal = .TRUE.
+    working_laser%profile_data_file = ' '
+
+    ! Phase-from-file defaults: disabled, with a blank filename that triggers
+    ! the default 'phase_profile.dat' at load time.
+    working_laser%use_phase_from_file = .FALSE.
+    working_laser%phase_data_file = ' '
+
   END SUBROUTINE laser_block_start
 
 
@@ -251,6 +263,80 @@ FUNCTION laser_block_handle_element(element, value) RESULT(errcode)
       RETURN
     END IF
 
+    ! Custom laser profile: read the amplitude (and optionally phase) from
+    ! an external binary file instead of a deck expression. See
+    ! src/shared_data.F90 for the laser_block fields these set.
+    IF (str_cmp(element, 'use_custom_profile')) THEN
+      working_laser%use_custom_profile = as_logical_print(value, element, &
+          errcode)
+      RETURN
+    END IF
+
+    IF (str_cmp(element, 'use_spatiotemporal_profile')) THEN
+      working_laser%use_spatiotemporal = as_logical_print(value, element, &
+          errcode)
+      RETURN
+    END IF
+
+    ! Parse the custom profile data filename. The value can be:
+    !   - A plain filename (e.g. 'TS01.dat'), resolved relative to data_dir
+    !   - An absolute path (e.g. '/home/user/profiles/TS01.dat'), used as-is
+    ! If omitted, the default filenames are used for backward compatibility.
+    IF (str_cmp(element, 'profile_data_file')) THEN
+      working_laser%profile_data_file = TRIM(ADJUSTL(value))
+      RETURN
+    END IF
+
+    ! Enable reading the spatiotemporal phase from file. When true, EPOCH
+    ! ignores any 'phase = ...' deck expression and instead interpolates the
+    ! phase from phase_data_file at every time step (handled in
+    ! laser_update_phase / custom_laser_phase).
+    IF (str_cmp(element, 'use_phase_from_file')) THEN
+      working_laser%use_phase_from_file = as_logical_print(value, &
+          element, errcode)
+      RETURN
+    END IF
+
+    ! Parse the phase data filename (same resolution rules as profile_data_file:
+    ! plain filenames resolve relative to data_dir, absolute paths used as-is).
+    ! If omitted, the default 'phase_profile.dat' is used.
+    IF (str_cmp(element, 'phase_data_file')) THEN
+      working_laser%phase_data_file = TRIM(ADJUSTL(value))
+      RETURN
+    END IF
+
+    ! Shape and bounds of the spatiotemporal profile/phase binary files
+    ! (use_spatiotemporal_profile = T only). These files carry no embedded
+    ! header (per EPOCH's documented binary-file convention), so the grid
+    ! must be declared here. The temporal extent reuses t_start/t_end rather
+    ! than a separate pair of elements, since the laser is only ever active
+    ! within that window anyway.
+    IF (str_cmp(element, 'n_t_points') .OR. str_cmp(element, 'n_t')) THEN
+      working_laser%n_t_points = as_integer_print(value, element, errcode)
+      RETURN
+    END IF
+
+    IF (str_cmp(element, 'n_transverse_points') &
+        .OR. str_cmp(element, 'n_y')) THEN
+      working_laser%n_transverse_points = &
+          as_integer_print(value, element, errcode)
+      RETURN
+    END IF
+
+    IF (str_cmp(element, 'profile_transverse_min') &
+        .OR. str_cmp(element, 'y_min')) THEN
+      working_laser%profile_transverse_min = &
+          as_real_print(value, element, errcode)
+      RETURN
+    END IF
+
+    IF (str_cmp(element, 'profile_transverse_max') &
+        .OR. str_cmp(element, 'y_max')) THEN
+      working_laser%profile_transverse_max = &
+          as_real_print(value, element, errcode)
+      RETURN
+    END IF
+
     errcode = c_err_unknown_element
 
   END FUNCTION laser_block_handle_element