From 2772501795a912fe27ce1ae604d49128be85f231 Mon Sep 17 00:00:00 2001 From: Tom Goffrey Date: Fri, 27 Feb 2026 13:03:03 +0000 Subject: [PATCH 1/6] Use constants where possible in QED routines --- epoch1d/src/physics_packages/photons.F90 | 9 ++++----- epoch2d/src/physics_packages/photons.F90 | 9 ++++----- epoch3d/src/physics_packages/photons.F90 | 9 ++++----- 3 files changed, 12 insertions(+), 15 deletions(-) diff --git a/epoch1d/src/physics_packages/photons.F90 b/epoch1d/src/physics_packages/photons.F90 index e4cfbf5ca..b20ffa0a5 100644 --- a/epoch1d/src/physics_packages/photons.F90 +++ b/epoch1d/src/physics_packages/photons.F90 @@ -679,7 +679,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 +693,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/src/physics_packages/photons.F90 b/epoch2d/src/physics_packages/photons.F90 index 6e7f40d6b..9b9298e41 100644 --- a/epoch2d/src/physics_packages/photons.F90 +++ b/epoch2d/src/physics_packages/photons.F90 @@ -683,7 +683,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 +697,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/src/physics_packages/photons.F90 b/epoch3d/src/physics_packages/photons.F90 index d45420cf3..ec3597c75 100644 --- a/epoch3d/src/physics_packages/photons.F90 +++ b/epoch3d/src/physics_packages/photons.F90 @@ -685,7 +685,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 +699,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 From 73d5ee06671434e5a01e047125e2cb18da966bba Mon Sep 17 00:00:00 2001 From: Tom Goffrey Date: Fri, 27 Feb 2026 13:09:46 +0000 Subject: [PATCH 2/6] Improve comment in QED routines --- epoch1d/src/physics_packages/photons.F90 | 3 ++- epoch2d/src/physics_packages/photons.F90 | 3 ++- epoch3d/src/physics_packages/photons.F90 | 3 ++- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/epoch1d/src/physics_packages/photons.F90 b/epoch1d/src/physics_packages/photons.F90 index b20ffa0a5..bb8e83680 100644 --- 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) diff --git a/epoch2d/src/physics_packages/photons.F90 b/epoch2d/src/physics_packages/photons.F90 index 9b9298e41..3ff3cc741 100644 --- 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) diff --git a/epoch3d/src/physics_packages/photons.F90 b/epoch3d/src/physics_packages/photons.F90 index ec3597c75..ae0ab0d62 100644 --- 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) From 8d77fcb99cc34cac2c3463fba1ab6e74fe5184d5 Mon Sep 17 00:00:00 2001 From: Tom Goffrey Date: Fri, 27 Feb 2026 14:28:26 +0000 Subject: [PATCH 3/6] Add QED Chi diagnostic --- epoch1d/Makefile | 2 +- epoch1d/src/constants.F90 | 3 ++- epoch1d/src/deck/deck_io_block.F90 | 3 +++ epoch1d/src/io/diagnostics.F90 | 3 +++ epoch1d/src/io/iterators.F90 | 23 +++++++++++++++++++++++ epoch2d/Makefile | 2 +- epoch2d/src/constants.F90 | 3 ++- epoch2d/src/deck/deck_io_block.F90 | 3 +++ epoch2d/src/io/diagnostics.F90 | 3 +++ epoch2d/src/io/iterators.F90 | 24 ++++++++++++++++++++++++ epoch3d/Makefile | 2 +- epoch3d/src/constants.F90 | 3 ++- epoch3d/src/deck/deck_io_block.F90 | 3 +++ epoch3d/src/io/diagnostics.F90 | 3 +++ epoch3d/src/io/iterators.F90 | 25 +++++++++++++++++++++++++ 15 files changed, 99 insertions(+), 6 deletions(-) diff --git a/epoch1d/Makefile b/epoch1d/Makefile index dc633fe65..790450480 100644 --- 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 index 46bfba58a..56ed134e4 100644 --- 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 index d0952b3dd..6d6a43f0c 100644 --- 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 index 85a3e1d94..53577bedb 100644 --- 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 index 52654ad7a..41b6b83b0 100644 --- 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,6 +81,7 @@ 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 + REAL(num) :: part_e, dir_x, dir_y, dir_z, norm, part_x IF (start) THEN CALL start_particle_list(current_species, current_list, cur) @@ -346,6 +348,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/epoch2d/Makefile b/epoch2d/Makefile index 84212ce4d..e84ad7dff 100644 --- 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 index 810d9b48b..1bcdcf0af 100644 --- 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 index c9c392d1b..0ba447bd3 100644 --- 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 index 0637579f2..1c367e6ca 100644 --- 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 index e9d6941ac..9c213196f 100644 --- 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,6 +81,7 @@ 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 + REAL(num) :: part_e, dir_x, dir_y, dir_z, norm, part_x, part_y IF (start) THEN CALL start_particle_list(current_species, current_list, cur) @@ -346,6 +348,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/epoch3d/Makefile b/epoch3d/Makefile index 29f71489f..70cdca79b 100644 --- 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 index 8fe305fad..d28c46d4e 100644 --- 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 index 0b8e85069..d545aca7a 100644 --- 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 index 640ca14fc..0fe270567 100644 --- 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 index e9d6941ac..1e1f6730a 100644 --- 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,6 +81,7 @@ 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 + REAL(num) :: part_e, dir_x, dir_y, dir_z, norm, part_x, part_y, part_z IF (start) THEN CALL start_particle_list(current_species, current_list, cur) @@ -346,6 +348,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 - x_grid_min_local + part_y = current%part_pos(2) - y_grid_min_local + part_z = current%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) From 7f1a1035945e9e0912347757eccea012f72197ee Mon Sep 17 00:00:00 2001 From: Tom Goffrey Date: Fri, 27 Feb 2026 15:47:06 +0000 Subject: [PATCH 4/6] IFDEF some variable declarations --- epoch1d/src/io/iterators.F90 | 3 ++- epoch2d/src/io/iterators.F90 | 3 ++- epoch3d/src/io/iterators.F90 | 3 ++- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/epoch1d/src/io/iterators.F90 b/epoch1d/src/io/iterators.F90 index 41b6b83b0..91cd776c9 100644 --- a/epoch1d/src/io/iterators.F90 +++ b/epoch1d/src/io/iterators.F90 @@ -81,8 +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 diff --git a/epoch2d/src/io/iterators.F90 b/epoch2d/src/io/iterators.F90 index 9c213196f..d9c290b2a 100644 --- a/epoch2d/src/io/iterators.F90 +++ b/epoch2d/src/io/iterators.F90 @@ -81,8 +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 diff --git a/epoch3d/src/io/iterators.F90 b/epoch3d/src/io/iterators.F90 index 1e1f6730a..ce99d49a7 100644 --- a/epoch3d/src/io/iterators.F90 +++ b/epoch3d/src/io/iterators.F90 @@ -81,8 +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 From 4f0028382e3476880208292cc0ff2230c09d7dee Mon Sep 17 00:00:00 2001 From: Tom Goffrey Date: Fri, 27 Feb 2026 16:06:10 +0000 Subject: [PATCH 5/6] Fix typos --- epoch3d/src/io/iterators.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/epoch3d/src/io/iterators.F90 b/epoch3d/src/io/iterators.F90 index ce99d49a7..56b4ccdd0 100644 --- a/epoch3d/src/io/iterators.F90 +++ b/epoch3d/src/io/iterators.F90 @@ -352,9 +352,9 @@ FUNCTION it_output_real(array, npoint_it, start, param) 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 - part_y = current%part_pos(2) - y_grid_min_local - part_z = current%part_pos(3) - z_grid_min_local + 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 From 9db034ac0f3af244825e770507386faa5a70756b Mon Sep 17 00:00:00 2001 From: Tom Goffrey Date: Wed, 6 May 2026 19:43:52 +0100 Subject: [PATCH 6/6] Add chi diagnostic for electrons and positrons --- epoch1d/src/io/iterators.F90 | 16 +++++++++++++++- epoch2d/src/io/iterators.F90 | 17 ++++++++++++++++- epoch3d/src/io/iterators.F90 | 16 ++++++++++++++++ 3 files changed, 47 insertions(+), 2 deletions(-) diff --git a/epoch1d/src/io/iterators.F90 b/epoch1d/src/io/iterators.F90 index 91cd776c9..f1a724b94 100644 --- a/epoch1d/src/io/iterators.F90 +++ b/epoch1d/src/io/iterators.F90 @@ -82,7 +82,8 @@ FUNCTION it_output_real(array, npoint_it, start, param) 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 + REAL(num) :: part_e, dir_x, dir_y, dir_z, norm, part_x, gamma_rel + REAL(num) :: part_ux, part_uy, part_uz #endif IF (start) THEN CALL start_particle_list(current_species, current_list, cur) @@ -363,6 +364,19 @@ FUNCTION it_output_real(array, npoint_it, start, param) dir_z, part_e) cur => cur%next END DO + ELSE IF (current_species%species_type == c_species_id_electron .OR. & + current_species%species_type == c_species_id_positron) THEN + DO WHILE (ASSOCIATED(cur) .AND. (part_count < npoint_it)) + part_x = cur%part_pos - x_grid_min_local + part_ux = cur%part_p(1) / mc0 + part_uy = cur%part_p(2) / mc0 + part_uz = cur%part_p(3) / mc0 + gamma_rel = SQRT(part_ux**2 + part_uy**2 + part_uz**2 + 1.0_num) + part_count = part_count + 1 + array(part_count) = calculate_eta(part_x, part_ux, part_uy, & + part_uz, gamma_rel) + cur => cur%next + END DO ELSE DO WHILE (ASSOCIATED(cur) .AND. (part_count < npoint_it)) part_count = part_count + 1 diff --git a/epoch2d/src/io/iterators.F90 b/epoch2d/src/io/iterators.F90 index d9c290b2a..92ac9e091 100644 --- a/epoch2d/src/io/iterators.F90 +++ b/epoch2d/src/io/iterators.F90 @@ -82,7 +82,8 @@ FUNCTION it_output_real(array, npoint_it, start, param) 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 + REAL(num) :: part_e, dir_x, dir_y, dir_z, norm, part_x, part_y, gamma_rel + REAL(num) :: part_ux, part_uy, part_uz #endif IF (start) THEN CALL start_particle_list(current_species, current_list, cur) @@ -364,6 +365,20 @@ FUNCTION it_output_real(array, npoint_it, start, param) dir_x, dir_y, dir_z, part_e) cur => cur%next END DO + ELSE IF (current_species%species_type == c_species_id_electron .OR. & + current_species%species_type == c_species_id_positron) 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_ux = cur%part_p(1) / mc0 + part_uy = cur%part_p(2) / mc0 + part_uz = cur%part_p(3) / mc0 + gamma_rel = SQRT(part_ux**2 + part_uy**2 + part_uz**2 + 1.0_num) + part_count = part_count + 1 + array(part_count) = calculate_eta(part_x, part_y,part_ux, part_uy, & + part_uz, gamma_rel) + cur => cur%next + END DO ELSE DO WHILE (ASSOCIATED(cur) .AND. (part_count < npoint_it)) part_count = part_count + 1 diff --git a/epoch3d/src/io/iterators.F90 b/epoch3d/src/io/iterators.F90 index 56b4ccdd0..8a464cec7 100644 --- a/epoch3d/src/io/iterators.F90 +++ b/epoch3d/src/io/iterators.F90 @@ -83,6 +83,7 @@ FUNCTION it_output_real(array, npoint_it, start, param) 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 + REAL(num) :: part_ux, part_uy, part_uz, gamma_rel #endif IF (start) THEN CALL start_particle_list(current_species, current_list, cur) @@ -365,6 +366,21 @@ FUNCTION it_output_real(array, npoint_it, start, param) dir_x, dir_y, dir_z, part_e) cur => cur%next END DO + ELSE IF (current_species%species_type == c_species_id_electron .OR. & + current_species%species_type == c_species_id_positron) 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 + part_ux = cur%part_p(1) / mc0 + part_uy = cur%part_p(2) / mc0 + part_uz = cur%part_p(3) / mc0 + gamma_rel = SQRT(part_ux**2 + part_uy**2 + part_uz**2 + 1.0_num) + part_count = part_count + 1 + array(part_count) = calculate_eta(part_x, part_y, part_z, & + part_ux, part_uy, part_uz, gamma_rel) + cur => cur%next + END DO ELSE DO WHILE (ASSOCIATED(cur) .AND. (part_count < npoint_it)) part_count = part_count + 1