epochpic · TomGoffrey · Feb 27, 2026 · Feb 27, 2026 · Feb 27, 2026 · Feb 27, 2026
diff --git a/epoch1d/Makefile b/epoch1d/Makefile
@@ -523,7 +523,7 @@ injectors.o: injectors.F90 evaluate.o file_injectors.o particle_temperature.o \
 ionise.o: ionise.F90 boundary.o calc_df.o numerics.o partlist.o \
   random_generator.o utilities.o
 iterators.o: iterators.F90 particle_id_hash.o particle_pointer_advance.o \
-  partlist.o
+  partlist.o photons.o
 laser.o: laser.f90 custom_laser.o evaluate.o
 mpi_routines.o: mpi_routines.F90 helper.o
 mpi_subtype_control.o: mpi_subtype_control.f90 shared_data.o

diff --git a/epoch1d/src/constants.F90 b/epoch1d/src/constants.F90
@@ -651,7 +651,8 @@ MODULE constants
   INTEGER, PARAMETER :: c_dump_part_rate_dr      = 76
   INTEGER, PARAMETER :: c_dump_part_rate_rr      = 77
   INTEGER, PARAMETER :: c_dump_part_rate_3br     = 78
-  INTEGER, PARAMETER :: num_vars_to_dump         = 78
+  INTEGER, PARAMETER :: c_dump_part_qed_chi      = 79
+  INTEGER, PARAMETER :: num_vars_to_dump         = 79
 
   INTEGER, PARAMETER :: c_subset_random     = 1
   INTEGER, PARAMETER :: c_subset_gamma_min  = 2

diff --git a/epoch1d/src/deck/deck_io_block.F90 b/epoch1d/src/deck/deck_io_block.F90
@@ -574,6 +574,9 @@ FUNCTION io_block_handle_element(element, value) RESULT(errcode)
 #ifdef PHOTONS
     ELSE IF (str_cmp(element, 'optical_depth')) THEN
       elementselected = c_dump_part_opdepth
+
+    ELSE IF (str_cmp(element, 'qed_chi')) THEN
+      elementselected = c_dump_part_qed_chi
 #endif
 
 #if defined(PHOTONS) || defined(BREMSSTRAHLUNG)

diff --git a/epoch1d/src/io/diagnostics.F90 b/epoch1d/src/io/diagnostics.F90
@@ -635,6 +635,9 @@ SUBROUTINE output_routines(step, force_write)   ! step = step index
 #ifdef PHOTONS
         CALL write_particle_variable(c_dump_part_opdepth, code, &
             'Optical depth', '', it_output_real)
+
+        CALL write_particle_variable(c_dump_part_qed_chi, code, &
+            'Chi', '', it_output_real)
 #endif
 #if defined(PHOTONS) || defined(BREMSSTRAHLUNG)
         CALL write_particle_variable(c_dump_part_qed_energy, code, &

diff --git a/epoch1d/src/io/iterators.F90 b/epoch1d/src/io/iterators.F90
@@ -18,6 +18,7 @@ MODULE iterators
   USE particle_pointer_advance
   USE partlist
   USE particle_id_hash_mod
+  USE photons
 
   IMPLICIT NONE
 
@@ -80,7 +81,9 @@ FUNCTION it_output_real(array, npoint_it, start, param)
     TYPE(particle_list), POINTER, SAVE :: current_list
     INTEGER :: part_count, ndim
     REAL(num) :: part_m, part_mc, part_mcc, part_mc2, gamma_mass, csqr, charge
-
+#ifdef PHOTONS
+    REAL(num) :: part_e, dir_x, dir_y, dir_z, norm, part_x
+#endif
     IF (start)  THEN
       CALL start_particle_list(current_species, current_list, cur)
     END IF
@@ -346,6 +349,27 @@ FUNCTION it_output_real(array, npoint_it, start, param)
           array(part_count) = cur%optical_depth
           cur => cur%next
         END DO
+      CASE (c_dump_part_qed_chi)
+        IF (current_species%species_type == c_species_id_photon) THEN
+          DO WHILE (ASSOCIATED(cur) .AND. (part_count < npoint_it))
+            part_x  = cur%part_pos - x_grid_min_local
+            norm  = c / cur%particle_energy
+            dir_x = cur%part_p(1) * norm
+            dir_y = cur%part_p(2) * norm
+            dir_z = cur%part_p(3) * norm
+            part_e  = cur%particle_energy / m0 / c**2
+            part_count = part_count + 1
+            array(part_count) = calculate_chi(part_x, dir_x, dir_y, &
+                dir_z, part_e)
+            cur => cur%next
+          END DO
+        ELSE
+          DO WHILE (ASSOCIATED(cur) .AND. (part_count < npoint_it))
+            part_count = part_count + 1
+            array(part_count) = 0.0_num
+            cur => cur%next
+          END DO
+        END IF
 #endif
 
 #if defined(PHOTONS) || defined(BREMSSTRAHLUNG)

diff --git a/epoch1d/src/physics_packages/photons.F90 b/epoch1d/src/physics_packages/photons.F90
@@ -606,7 +606,8 @@ SUBROUTINE qed_update_optical_depth
 
           current%optical_depth = current%optical_depth &
               - delta_optical_depth_photon(chi_val, part_e)
-          ! If optical depth dropped below zero generate pair...
+          ! If optical depth dropped below zero generate pair
+          ! Photon will be deleted
           IF (current%optical_depth < 0.0_num) THEN
             CALL generate_pair(current, chi_val, photon_species, &
                 breit_wheeler_electron_species, breit_wheeler_positron_species)
@@ -679,7 +680,10 @@ FUNCTION calculate_eta(part_x, part_ux, part_uy, part_uz, &
     REAL(num) :: e_at_part(3), b_at_part(3)
     REAL(num) :: beta_x, beta_y, beta_z
     REAL(num) :: flperp(3), i_e, tau0, roland_eta
-    REAL(num) :: lambdac, coeff_eta, moduclip, moduclip2, u_dot_e
+    REAL(num) :: moduclip, moduclip2, u_dot_e
+    REAL(num), PARAMETER :: lambdac = h_bar / mc0
+    REAL(num), PARAMETER :: coeff_eta = SQRT(3.0_num * lambdac &
+        / (2.0_num * alpha_f * m0 * c**3))
 
     CALL field_at_particle(part_x, e_at_part, b_at_part)
 
@@ -690,10 +694,6 @@ FUNCTION calculate_eta(part_x, part_ux, part_uy, part_uz, &
     beta_y = part_uy / gamma_rel
     beta_z = part_uz / gamma_rel
 
-    lambdac = h_bar / mc0
-
-    coeff_eta = SQRT(3.0_num * lambdac / (2.0_num * alpha_f * m0 * c**3))
-
     u_dot_e = (part_ux * e_at_part(1) + part_uy * e_at_part(2) &
         + part_uz * e_at_part(3)) / moduclip2
 

diff --git a/epoch2d/Makefile b/epoch2d/Makefile
@@ -523,7 +523,7 @@ injectors.o: injectors.F90 evaluate.o file_injectors.o particle_temperature.o \
 ionise.o: ionise.F90 boundary.o calc_df.o numerics.o partlist.o \
   random_generator.o utilities.o
 iterators.o: iterators.F90 particle_id_hash.o particle_pointer_advance.o \
-  partlist.o
+  partlist.o photons.o
 laser.o: laser.f90 custom_laser.o evaluate.o
 mpi_routines.o: mpi_routines.F90 helper.o
 mpi_subtype_control.o: mpi_subtype_control.f90 shared_data.o

diff --git a/epoch2d/src/constants.F90 b/epoch2d/src/constants.F90
@@ -653,7 +653,8 @@ MODULE constants
   INTEGER, PARAMETER :: c_dump_part_rate_dr      = 76
   INTEGER, PARAMETER :: c_dump_part_rate_rr      = 77
   INTEGER, PARAMETER :: c_dump_part_rate_3br     = 78
-  INTEGER, PARAMETER :: num_vars_to_dump         = 78
+  INTEGER, PARAMETER :: c_dump_part_qed_chi      = 79
+  INTEGER, PARAMETER :: num_vars_to_dump         = 79
 
   INTEGER, PARAMETER :: c_subset_random     = 1
   INTEGER, PARAMETER :: c_subset_gamma_min  = 2

diff --git a/epoch2d/src/deck/deck_io_block.F90 b/epoch2d/src/deck/deck_io_block.F90
@@ -575,6 +575,9 @@ FUNCTION io_block_handle_element(element, value) RESULT(errcode)
 #ifdef PHOTONS
     ELSE IF (str_cmp(element, 'optical_depth')) THEN
       elementselected = c_dump_part_opdepth
+
+    ELSE IF (str_cmp(element, 'qed_chi')) THEN
+      elementselected = c_dump_part_qed_chi
 #endif
 
 #if defined(PHOTONS) || defined(BREMSSTRAHLUNG)

diff --git a/epoch2d/src/io/diagnostics.F90 b/epoch2d/src/io/diagnostics.F90
@@ -644,6 +644,9 @@ SUBROUTINE output_routines(step, force_write)   ! step = step index
 #ifdef PHOTONS
         CALL write_particle_variable(c_dump_part_opdepth, code, &
             'Optical depth', '', it_output_real)
+
+        CALL write_particle_variable(c_dump_part_qed_chi, code, &
+            'Chi', '', it_output_real)
 #endif
 #if defined(PHOTONS) || defined(BREMSSTRAHLUNG)
         CALL write_particle_variable(c_dump_part_qed_energy, code, &

diff --git a/epoch2d/src/io/iterators.F90 b/epoch2d/src/io/iterators.F90
@@ -18,6 +18,7 @@ MODULE iterators
   USE particle_pointer_advance
   USE partlist
   USE particle_id_hash_mod
+  USE photons
 
   IMPLICIT NONE
 
@@ -80,7 +81,9 @@ FUNCTION it_output_real(array, npoint_it, start, param)
     TYPE(particle_list), POINTER, SAVE :: current_list
     INTEGER :: part_count, ndim
     REAL(num) :: part_m, part_mc, part_mcc, part_mc2, gamma_mass, csqr, charge
-
+#ifdef PHOTONS
+    REAL(num) :: part_e, dir_x, dir_y, dir_z, norm, part_x, part_y
+#endif
     IF (start)  THEN
       CALL start_particle_list(current_species, current_list, cur)
     END IF
@@ -346,6 +349,28 @@ FUNCTION it_output_real(array, npoint_it, start, param)
           array(part_count) = cur%optical_depth
           cur => cur%next
         END DO
+      CASE (c_dump_part_qed_chi)
+        IF (current_species%species_type == c_species_id_photon) THEN
+          DO WHILE (ASSOCIATED(cur) .AND. (part_count < npoint_it))
+            part_x = cur%part_pos(1) - x_grid_min_local
+            part_y = cur%part_pos(2) - y_grid_min_local
+            norm  = c / cur%particle_energy
+            dir_x = cur%part_p(1) * norm
+            dir_y = cur%part_p(2) * norm
+            dir_z = cur%part_p(3) * norm
+            part_e  = cur%particle_energy / m0 / c**2
+            part_count = part_count + 1
+            array(part_count) = calculate_chi(part_x, part_y, &
+                dir_x, dir_y, dir_z, part_e)
+            cur => cur%next
+          END DO
+        ELSE
+          DO WHILE (ASSOCIATED(cur) .AND. (part_count < npoint_it))
+            part_count = part_count + 1
+            array(part_count) = 0.0_num
+            cur => cur%next
+          END DO
+        END IF
 #endif
 
 #if defined(PHOTONS) || defined(BREMSSTRAHLUNG)

diff --git a/epoch2d/src/physics_packages/photons.F90 b/epoch2d/src/physics_packages/photons.F90
@@ -610,7 +610,8 @@ SUBROUTINE qed_update_optical_depth
 
           current%optical_depth = current%optical_depth &
               - delta_optical_depth_photon(chi_val, part_e)
-          ! If optical depth dropped below zero generate pair...
+          ! If optical depth dropped below zero generate pair
+          ! Photon will be deleted
           IF (current%optical_depth < 0.0_num) THEN
             CALL generate_pair(current, chi_val, photon_species, &
                 breit_wheeler_electron_species, breit_wheeler_positron_species)
@@ -683,7 +684,10 @@ FUNCTION calculate_eta(part_x, part_y, part_ux, part_uy, part_uz, &
     REAL(num) :: e_at_part(3), b_at_part(3)
     REAL(num) :: beta_x, beta_y, beta_z
     REAL(num) :: flperp(3), i_e, tau0, roland_eta
-    REAL(num) :: lambdac, coeff_eta, moduclip, moduclip2, u_dot_e
+    REAL(num) :: moduclip, moduclip2, u_dot_e
+    REAL(num), PARAMETER :: lambdac = h_bar / mc0
+    REAL(num), PARAMETER :: coeff_eta = SQRT(3.0_num * lambdac &
+        / (2.0_num * alpha_f * m0 * c**3))
 
     CALL field_at_particle(part_x, part_y, e_at_part, b_at_part)
 
@@ -694,10 +698,6 @@ FUNCTION calculate_eta(part_x, part_y, part_ux, part_uy, part_uz, &
     beta_y = part_uy / gamma_rel
     beta_z = part_uz / gamma_rel
 
-    lambdac = h_bar / mc0
-
-    coeff_eta = SQRT(3.0_num * lambdac / (2.0_num * alpha_f * m0 * c**3))
-
     u_dot_e = (part_ux * e_at_part(1) + part_uy * e_at_part(2) &
         + part_uz * e_at_part(3)) / moduclip2
 

diff --git a/epoch3d/Makefile b/epoch3d/Makefile
@@ -523,7 +523,7 @@ injectors.o: injectors.F90 evaluate.o file_injectors.o particle_temperature.o \
 ionise.o: ionise.F90 boundary.o calc_df.o numerics.o partlist.o \
   random_generator.o utilities.o
 iterators.o: iterators.F90 particle_id_hash.o particle_pointer_advance.o \
-  partlist.o
+  partlist.o photons.o
 laser.o: laser.f90 custom_laser.o evaluate.o
 mpi_routines.o: mpi_routines.F90 helper.o
 mpi_subtype_control.o: mpi_subtype_control.f90 shared_data.o

diff --git a/epoch3d/src/constants.F90 b/epoch3d/src/constants.F90
@@ -654,7 +654,8 @@ MODULE constants
   INTEGER, PARAMETER :: c_dump_part_rate_dr      = 76
   INTEGER, PARAMETER :: c_dump_part_rate_rr      = 77
   INTEGER, PARAMETER :: c_dump_part_rate_3br     = 78
-  INTEGER, PARAMETER :: num_vars_to_dump         = 78
+  INTEGER, PARAMETER :: c_dump_part_qed_chi      = 79
+  INTEGER, PARAMETER :: num_vars_to_dump         = 79
 
   INTEGER, PARAMETER :: c_subset_random     = 1
   INTEGER, PARAMETER :: c_subset_gamma_min  = 2

diff --git a/epoch3d/src/deck/deck_io_block.F90 b/epoch3d/src/deck/deck_io_block.F90
@@ -576,6 +576,9 @@ FUNCTION io_block_handle_element(element, value) RESULT(errcode)
 #ifdef PHOTONS
     ELSE IF (str_cmp(element, 'optical_depth')) THEN
       elementselected = c_dump_part_opdepth
+
+    ELSE IF (str_cmp(element, 'qed_chi')) THEN
+      elementselected = c_dump_part_qed_chi
 #endif
 
 #if defined(PHOTONS) || defined(BREMSSTRAHLUNG)

diff --git a/epoch3d/src/io/diagnostics.F90 b/epoch3d/src/io/diagnostics.F90
@@ -653,6 +653,9 @@ SUBROUTINE output_routines(step, force_write)   ! step = step index
 #ifdef PHOTONS
         CALL write_particle_variable(c_dump_part_opdepth, code, &
             'Optical depth', '', it_output_real)
+
+        CALL write_particle_variable(c_dump_part_qed_chi, code, &
+            'Chi', '', it_output_real)
 #endif
 #if defined(PHOTONS) || defined(BREMSSTRAHLUNG)
         CALL write_particle_variable(c_dump_part_qed_energy, code, &

diff --git a/epoch3d/src/io/iterators.F90 b/epoch3d/src/io/iterators.F90
@@ -18,6 +18,7 @@ MODULE iterators
   USE particle_pointer_advance
   USE partlist
   USE particle_id_hash_mod
+  USE photons
 
   IMPLICIT NONE
 
@@ -80,7 +81,9 @@ FUNCTION it_output_real(array, npoint_it, start, param)
     TYPE(particle_list), POINTER, SAVE :: current_list
     INTEGER :: part_count, ndim
     REAL(num) :: part_m, part_mc, part_mcc, part_mc2, gamma_mass, csqr, charge
-
+#ifdef PHOTONS
+    REAL(num) :: part_e, dir_x, dir_y, dir_z, norm, part_x, part_y, part_z
+#endif
     IF (start)  THEN
       CALL start_particle_list(current_species, current_list, cur)
     END IF
@@ -346,6 +349,29 @@ FUNCTION it_output_real(array, npoint_it, start, param)
           array(part_count) = cur%optical_depth
           cur => cur%next
         END DO
+      CASE (c_dump_part_qed_chi)
+        IF (current_species%species_type == c_species_id_photon) THEN
+          DO WHILE (ASSOCIATED(cur) .AND. (part_count < npoint_it))
+            part_x  = cur%part_pos(1) - x_grid_min_local
+            part_y  = cur%part_pos(2) - y_grid_min_local
+            part_z  = cur%part_pos(3) - z_grid_min_local
+            norm  = c / cur%particle_energy
+            dir_x = cur%part_p(1) * norm
+            dir_y = cur%part_p(2) * norm
+            dir_z = cur%part_p(3) * norm
+            part_e  = cur%particle_energy / m0 / c**2
+            part_count = part_count + 1
+            array(part_count) = calculate_chi(part_x, part_y, part_z, &
+                dir_x, dir_y, dir_z, part_e)
+            cur => cur%next
+          END DO
+        ELSE
+          DO WHILE (ASSOCIATED(cur) .AND. (part_count < npoint_it))
+            part_count = part_count + 1
+            array(part_count) = 0.0_num
+            cur => cur%next
+          END DO
+        END IF
 #endif
 
 #if defined(PHOTONS) || defined(BREMSSTRAHLUNG)

diff --git a/epoch3d/src/physics_packages/photons.F90 b/epoch3d/src/physics_packages/photons.F90
@@ -612,7 +612,8 @@ SUBROUTINE qed_update_optical_depth
 
           current%optical_depth = current%optical_depth &
               - delta_optical_depth_photon(chi_val, part_e)
-          ! If optical depth dropped below zero generate pair...
+          ! If optical depth dropped below zero generate pair
+          ! Photon will be deleted
           IF (current%optical_depth < 0.0_num) THEN
             CALL generate_pair(current, chi_val, photon_species, &
                 breit_wheeler_electron_species, breit_wheeler_positron_species)
@@ -685,7 +686,10 @@ FUNCTION calculate_eta(part_x, part_y, part_z, part_ux, part_uy, part_uz, &
     REAL(num) :: e_at_part(3), b_at_part(3)
     REAL(num) :: beta_x, beta_y, beta_z
     REAL(num) :: flperp(3), i_e, tau0, roland_eta
-    REAL(num) :: lambdac, coeff_eta, moduclip, moduclip2, u_dot_e
+    REAL(num) :: moduclip, moduclip2, u_dot_e
+    REAL(num), PARAMETER :: lambdac = h_bar / mc0
+    REAL(num), PARAMETER :: coeff_eta = SQRT(3.0_num * lambdac &
+        / (2.0_num * alpha_f * m0 * c**3))
 
     CALL field_at_particle(part_x, part_y, part_z, e_at_part, b_at_part)
 
@@ -696,10 +700,6 @@ FUNCTION calculate_eta(part_x, part_y, part_z, part_ux, part_uy, part_uz, &
     beta_y = part_uy / gamma_rel
     beta_z = part_uz / gamma_rel
 
-    lambdac = h_bar / mc0
-
-    coeff_eta = SQRT(3.0_num * lambdac / (2.0_num * alpha_f * m0 * c**3))
-
     u_dot_e = (part_ux * e_at_part(1) + part_uy * e_at_part(2) &
         + part_uz * e_at_part(3)) / moduclip2