Add Linear Compton Scattering

epoch1d/Makefile

@@ -536,7 +536,7 @@ particle_temperature.o: particle_temperature.F90 constants.o evaluate.o \
   random_generator.o
 particles.o: particles.F90 boundary.o partlist.o prefetch.o
 partlist.o: partlist.F90 particle_id_hash.o random_generator.o shared_data.o
-photons.o: photons.F90 partlist.o
+photons.o: photons.F90 collisions.o partlist.o
 prefetch.o: prefetch.F90 shared_data.o
 probes.o: probes.F90 partlist.o $(SDFMOD)
 random_generator.o: random_generator.f90

epoch1d/src/constants.F90

@@ -229,6 +229,15 @@ MODULE constants
   REAL(num), PARAMETER :: alpha_f = 7.297352575523020256850802729527158e-3_num
   ! tau_c = h_bar / (m0 * c**2)
   REAL(num), PARAMETER :: tau_c = 1.288088667367242662108649212042082e-21_num
+
+  REAL(num), PARAMETER :: classical_re = 0.25_num / pi / epsilon0 / m0 &
+      * (q0 / c)**2
+  REAL(num), PARAMETER :: sigma_thomson = 8.0_num * pi / 3.0_num &
+      * classical_re**2
+  REAL(num), PARAMETER :: inv_mc0 = 1.0_num / mc0
+  REAL(num), PARAMETER :: inv_m0c2 = 1.0_num / m0c2
+  REAL(num), PARAMETER :: pire2 = pi * classical_re**2
+  REAL(num), PARAMETER :: quarter_pire2 = 0.25_num * pire2
 #endif
 
   ! Constants used for bremsstrahlung with plasma screening

epoch1d/src/deck/deck_qed_block.F90

@@ -45,6 +45,9 @@ SUBROUTINE qed_deck_initialise
       use_radiation_reaction = .TRUE.
       produce_photons = .FALSE.
       photon_dynamics = .FALSE.
+      use_binary_collisions = .FALSE.
+      use_LCS = .FALSE.
+      use_LCS_diff = .TRUE.
     END IF
 #endif
 
@@ -57,6 +60,7 @@ SUBROUTINE qed_deck_finalise
     INTEGER :: io, iu
 #ifdef PHOTONS
     LOGICAL :: exists
+    INTEGER :: j
 
     IF (deck_state == c_ds_first) RETURN
 
@@ -74,6 +78,22 @@ SUBROUTINE qed_deck_finalise
     END IF
 
     IF (use_qed) need_random_state = .TRUE.
+
+    use_binary_collisions = use_LCS
+
+    IF (use_binary_collisions) THEN
+      DO j = 1, n_species
+        IF (species_list(j)%species_type == c_species_id_photon) THEN
+          species_list(j)%make_secondary_list = .TRUE.
+        END IF
+        IF (species_list(j)%species_type == c_species_id_electron) THEN
+          species_list(j)%make_secondary_list = .TRUE.
+        END IF
+        IF (species_list(j)%species_type == c_species_id_positron) THEN
+          species_list(j)%make_secondary_list = .TRUE.
+        END IF
+      END DO
+    END IF
 #else
     IF (use_qed) THEN
       IF (rank == 0) THEN
@@ -87,7 +107,6 @@ SUBROUTINE qed_deck_finalise
       CALL abort_code(c_err_pp_options_missing)
     END IF
 #endif
-
   END SUBROUTINE qed_deck_finalise
 
 
@@ -172,6 +191,16 @@ FUNCTION qed_block_handle_element(element, value) RESULT(errcode)
       RETURN
     END IF
 
+    IF(str_cmp(element, 'linear_compton_scattering')) THEN
+      use_LCS = as_logical_print(value, element, errcode)
+      RETURN
+    END IF
+
+    IF(str_cmp(element, 'LCS_differential_cross')) THEN
+      use_LCS_diff = as_logical_print(value, element, errcode)
+      RETURN
+    END IF
+
     errcode = c_err_unknown_element
 #endif
 

epoch1d/src/epoch1d.F90

@@ -216,7 +216,8 @@ PROGRAM pic
       ! .FALSE. this time to use load balancing threshold
       IF (use_balance) CALL balance_workload(.FALSE.)
       CALL push_particles
-      IF (use_particle_lists) THEN
+
+      IF (use_particle_lists .OR. use_binary_collisions) THEN
         ! Check whether this is a step with collisions or collisional ionisation
         collision_step = (MODULO(step, coll_n_step) == coll_n_step - 1) &
           .AND. use_collisions
@@ -227,7 +228,8 @@ PROGRAM pic
 
         ! After this line, the particles can be accessed on a cell by cell basis
         ! Using the particle_species%secondary_list property
-        IF (collision_step .OR. coll_ion_step .OR. recombine_step) THEN
+        IF (collision_step .OR. coll_ion_step .OR. recombine_step &
+            .OR. use_binary_collisions) THEN
           CALL reorder_particles_to_grid
         END IF
 
@@ -242,7 +244,13 @@ PROGRAM pic
 
         IF (recombine_step) CALL run_recombination
 
-        IF (collision_step .OR. coll_ion_step .OR. recombine_step) THEN
+#ifdef PHOTONS
+        IF (use_binary_collisions) THEN
+          CALL do_binary_collisions
+        END IF
+#endif
+        IF (collision_step .OR. coll_ion_step .OR. recombine_step &
+            .OR. use_binary_collisions) THEN
           CALL reattach_particles_to_mainlist
         END IF
       END IF