*----
*  TEST CASE TCM78
*
*  Coupled transport of photons, electrons and positrons in water
*  Cross-sections library produced by Radiant package
*  10 MeV normally incident electron source
*  Convergence of the energy deposition profile
*
*  Ref: Bienvenue, Charles, et al. "Toward Highly Accurate 
*  Multigroup Coupled Photon-Electron-Positron Cross-Sections
*  for the Boltzmann Fokker-Planck Equation." Journal of
*  Computational Physics (2025): 113740.
*
*----

*----
*  Define STRUCTURES and MODULES used
*----

LOGICAL IS_CONV ; 
INTEGER NPART NGEN_MAX IPART IGEN NX IX ;
DOUBLE EDEP EDEP1 EPSILON EPSILON1 EPSILON_EDEP ;
MODULE GEO: FMAC: MAC: SNT: ASM: FLU: PSOUR: OUT: DELETE: END: GREP: HEAT: ABORT: ;
PROCEDURE assertS_low ;
LINKED_LIST GEOM ETRACK TRACK GTRACK EMACRO GMACRO PMACRO MACRO ESOUR SOUR FLUX FLUX1 FLUX2 SYSTEM
     RR ERR GRR PRR DEPOS DEPOS1 ;
SEQ_ASCII FMAC_M :: FILE './radiant_fmac.txt' ;

*----
*  Set and track the geometry
*----

EVALUATE NX := 40 ;
GEOM := GEO: :: CAR1D 1
     X- VOID  X+ VOID
     MIX 1
     MESHX 0.0 10.0
     SPLITX <<NX>>
     ;
ETRACK := SNT: GEOM ::
     MAXR 300
     NSHT 
     SCHM 2 DIAM 0     ! Linear discontinuous Galerkin for streaming term
     ESCHM 2 EDIAM 0   ! Linear discontinuous Galerkin for CSD term
     QUAD 2 GQ 1       ! Gauss-Legendre with Galerkin quadrature method
     SN 8 SCAT 8 
     EPSI 1.E-5
     MAXI 1000 
     BFPG 
     ;
GTRACK := ETRACK ;
GTRACK := SNT: GTRACK GEOM :: BTE ;

*----
*  Extract the cross-sections
*----

! Electrons
EMACRO := FMAC: FMAC_M :: PARTICLE B ;
MACRO := EMACRO ;
EMACRO := MAC: EMACRO MACRO :: MIX 1 1 OLDL ;
MACRO := DELETE: MACRO ;

! Photons
GMACRO := FMAC: FMAC_M :: PARTICLE G ;
MACRO := GMACRO ;
GMACRO := MAC: GMACRO MACRO :: MIX 1 1 OLDL ;
MACRO := DELETE: MACRO ;

! Positrons
PMACRO := FMAC: FMAC_M :: PARTICLE C ;
MACRO := PMACRO ;
PMACRO := MAC: PMACRO MACRO :: MIX 1 1 OLDL ;
MACRO := DELETE: MACRO ;

*----
*  Fixed external electron source definition
*----

ESOUR := PSOUR: EMACRO ETRACK GEOM ::
     MONO 1
     X-
     INTG 1 1.0
     DIR 1.0 0.0 0.0
     ;

*----
*  Coupled transport calculations
*----

EVALUATE NPART := 3 ;             ! Number of particles
EVALUATE NGEN_MAX := 10 ;         ! Maximum number of particle generations
EVALUATE EPSILON_EDEP := 1.D-5 ;  ! Convergence criterion on energy deposition profile
EVALUATE IS_CONV := $False_L ;

EVALUATE IGEN := 1 ;
WHILE IGEN NGEN_MAX <= IS_CONV NOT * DO
     EVALUATE IPART := 1 ;
     WHILE IPART NPART <= DO

          ECHO " " ;
          ECHO " COUPLED TRANSPORT INDEXES -------" ;
          ECHO " IPART: " IPART ;
          ECHO " IGEN: " IGEN ;
          ECHO "----------------------------------" ;
          ECHO " " ;

          *----
          *  Set the cross-sections and tracking
          *----

          IF IPART 1 = THEN
               MACRO := EMACRO ;
               TRACK := ETRACK ;
          ELSEIF IPART 2 = THEN
               MACRO := GMACRO ;
               TRACK := GTRACK ;
          ELSEIF IPART 3 = THEN
               MACRO := PMACRO ;
               TRACK := ETRACK ;
          ENDIF ;

          *----
          *  Set the sources of particles
          *----

          IF IPART 1 = THEN
               IF IGEN 1 = THEN
                    SOUR := ESOUR ;
               ELSE
                    SOUR := PSOUR: MACRO TRACK ETRACK FLUX1 :: PARTICLE C ;
                    SOUR := PSOUR: MACRO TRACK GTRACK FLUX2 SOUR :: PARTICLE G ;
               ENDIF ;
          ELSEIF IPART 2 = THEN
               SOUR := PSOUR: MACRO TRACK ETRACK FLUX1 :: PARTICLE B ;
               IF IGEN 1 > THEN
                    SOUR := PSOUR: MACRO TRACK ETRACK FLUX2 SOUR :: PARTICLE C ;
               ENDIF ;
          ELSEIF IPART 3 = THEN
               SOUR := PSOUR: MACRO TRACK GTRACK FLUX1 :: PARTICLE G ;
               SOUR := PSOUR: MACRO TRACK ETRACK FLUX2 SOUR :: PARTICLE B ;
          ENDIF ;

          *----
          *  Compute the particle flux
          *----

          SYSTEM := ASM: MACRO TRACK :: ARM ;
          FLUX := FLU: MACRO TRACK SYSTEM SOUR :: TYPE S EXTE 100 ;
          RR := OUT: FLUX TRACK MACRO GEOM :: COND INTG NONE ;
          IF IPART 1 = THEN
               ERR := RR ;
          ELSEIF IPART 2 = THEN
               GRR := RR ;
          ELSEIF IPART 3 = THEN
               PRR := RR ;
          ENDIF ;

          *----
          *  Compute the energy deposition
          *----

          IF IPART 3 = IGEN 1 = * THEN
               DEPOS := HEAT: ERR GRR PRR :: EDIT 1 NORM 1.0 ;
               ERR GRR PRR := DELETE: ERR GRR PRR ;
          ELSEIF IPART 3 = THEN
               DEPOS1 := DEPOS ;
               DEPOS := HEAT: DEPOS ERR GRR PRR :: EDIT 1 NORM 1.0 ;
               ERR GRR PRR := DELETE: ERR GRR PRR ;
          ENDIF ;

          *----
          * Verify convergence of the energy deposition profile
          *----

          IF IPART 3 = IGEN 1 = NOT * THEN
               EVALUATE IX := 1 ;
               EVALUATE EPSILON := 1.D-10 ;
               WHILE IX NX <= DO
                    GREP: DEPOS :: GETVAL 'EDEPOS_TOT' <<IX>> >>EDEP<< ;
                    GREP: DEPOS1 :: GETVAL 'EDEPOS_TOT' <<IX>> >>EDEP1<< ;
                    IF EDEP EPSILON_EDEP > THEN
                         EVALUATE EPSILON1 := EDEP EDEP1 - EDEP / ;
                    ELSE
                         EVALUATE EPSILON1 := EDEP EDEP1 - EPSILON_EDEP / ;
                    ENDIF ;
                    IF EPSILON1 EPSILON > THEN
                         EVALUATE EPSILON := EPSILON1 ;
                    ENDIF ;
                    EVALUATE IX := IX 1 + ;
               ENDWHILE ;
               IF EPSILON EPSILON_EDEP <= THEN
                    EVALUATE IS_CONV := $True_L ;
               ENDIF ;
               DEPOS1 := DELETE: DEPOS1 ;
          ENDIF ;

          *----
          *  Set the fluxes
          *----

          IF IPART 1 = IGEN 1 = * THEN
               FLUX1 := FLUX ; 
          ELSEIF IPART 2 = IGEN 1 = * THEN
               FLUX2 := FLUX1 ;
               FLUX1 := DELETE: FLUX1 ;
               FLUX1 := FLUX ;
          ELSE
               FLUX2 := DELETE: FLUX2 ;
               FLUX2 := FLUX1 ;
               FLUX1 := DELETE: FLUX1 ;
               FLUX1 := FLUX ;
          ENDIF ;

          SOUR FLUX SYSTEM MACRO TRACK RR := DELETE: SOUR FLUX SYSTEM MACRO TRACK RR ;
          EVALUATE IPART := IPART 1 + ;

     ENDWHILE ;
     EVALUATE IGEN := IGEN 1 + ;
ENDWHILE ;

*----
*  Asserts
*----
assertS_low DEPOS :: 'EDEPOS_TOT'  1 1.909428E+00  1.0E-4 ;
assertS_low DEPOS :: 'EDEPOS_TOT'  3 2.100268E+00  1.0E-4 ;
assertS_low DEPOS :: 'EDEPOS_TOT' 10 2.191387E+00  1.0E-4 ;
assertS_low DEPOS :: 'EDEPOS_TOT' 30 1.893547E-01  1.0E-4 ;
assertS_low DEPOS :: 'CDEPOS_TOT'  1 -1.488792E-01 5.0E-4 ;
assertS_low DEPOS :: 'CDEPOS_TOT'  3 -4.382314E-03 5.0E-2 ;
assertS_low DEPOS :: 'CDEPOS_TOT' 10 1.891187E-01  5.0E-4 ;
assertS_low DEPOS :: 'CDEPOS_TOT' 30 4.212810E-02  5.0E-4 ;
ECHO "test TCM78 completed" ;
END: ;