Resolve "Implement analytic inelastic scattering laws for Draglibs in module LIB:"

6 files changed, 386 insertions, 106 deletions

diff --git a/Dragon/data/uo2_kec1_ecco1962_light.x2m b/Dragon/data/uo2_kec1_ecco1962_light.x2m
index 4153470..4141c5a 100644
--- a/Dragon/data/uo2_kec1_ecco1962_light.x2m
+++ b/Dragon/data/uo2_kec1_ecco1962_light.x2m
@@ -1,6 +1,6 @@
 *----
 *
-*  TEST CASE uo2_kec1_ecco1962_light
+*  TEST CASE uo2_kec1_ecco1962_light_v5p1
 *  UO2 ROWLAND'S BENCHMARK
 *  DISTRIBUTED SELF-SHIELDING
 *  1962-GROUP ENDF/B-VIII DRAGLIB
@@ -83,7 +83,7 @@ LIBRARY := LIB: ::
 LIBRARY2 := USS: LIBRARY TRACK_SS :: EDIT 1 TRAN PASS 3 GRMIN 52 ;
 CP := ASM: LIBRARY2 TRACK TF_EXC :: EDIT 1 PIJ ;
 CALC := FLU: CP LIBRARY2 TRACK :: EDIT 1 TYPE K ;
-assertS CALC :: K-EFFECTIVE 1 1.400612 ;
+assertS CALC :: K-EFFECTIVE 1 1.394335 ;
 *
 * Tone's method
 LIBRARY CP CALC := DELETE: LIBRARY CP CALC ;
@@ -117,7 +117,7 @@ LIBRARY := LIB: ::
 LIBRARY := TONE: LIBRARY TRACK_SS :: EDIT 1 TRAN MXIT 3 SPH GRMIN 52 ;
 CP := ASM: LIBRARY TRACK TF_EXC :: EDIT 1 PIJ ;
 CALC := FLU: CP LIBRARY TRACK :: EDIT 1 TYPE K ;
-assertS CALC :: K-EFFECTIVE 1 1.400571 ;
+assertS CALC :: K-EFFECTIVE 1 1.394283 ;
 *
 * Resonance spectrum expansion method
 LIBRARY LIBRARY2 CP CALC := DELETE: LIBRARY LIBRARY2 CP CALC ;
@@ -152,7 +152,7 @@ LIBRARY := LIB: ::
 LIBRARY2 := USS: LIBRARY TRACK_SS :: EDIT 1 TRAN PASS 3 GRMIN 52 ;
 CP := ASM: LIBRARY2 TRACK TF_EXC :: EDIT 1 PIJ ;
 CALC := FLU: CP LIBRARY2 TRACK :: EDIT 1 TYPE K ;
-assertS CALC :: K-EFFECTIVE 1 1.401772 ;
+assertS CALC :: K-EFFECTIVE 1 1.395505 ;
 *
 * Autosecol UFG method
 LIBRARY LIBRARY2 CP CALC := DELETE: LIBRARY LIBRARY2 CP CALC ;
@@ -184,9 +184,9 @@ LIBRARY := LIB: ::
     H1      = H1     6.6988E-2
     O16     = O16    3.3494E-2
  ;
-LIBRARY2 := AUTO: LIBRARY TRACK_SS :: EDIT 1 GRMIN 52 MAXT 60000 SEED 123456 ;
+LIBRARY2 := AUTO: LIBRARY TRACK_SS :: EDIT 3 GRMIN 52 MAXT 60000 SEED 123456 ;
 CP := ASM: LIBRARY2 TRACK TF_EXC :: EDIT 1 PIJ ;
 CALC := FLU: CP LIBRARY2 TRACK :: EDIT 1 TYPE K ;
-assertS CALC :: K-EFFECTIVE 1 1.401040 ;
-ECHO "test uo2_kec1_ecco1962_light completed" ;
+assertS CALC :: K-EFFECTIVE 1 1.400297 ;
+ECHO "test uo2_kec1_ecco1962_light_v5p1 completed" ;
 QUIT "LIST" .
diff --git a/Dragon/src/LIBDRA.f b/Dragon/src/LIBDRA.f
index 9d23805..51438a9 100644
--- a/Dragon/src/LIBDRA.f
+++ b/Dragon/src/LIBDRA.f
@@ -87,7 +87,7 @@
       ALLOCATE(NFS(NGRO),ITYPRO(NL))
       ALLOCATE(AWR(NBISO),DELTA(NGRO),SIGS(NGRO,NL),SCAT(NGRO,NGRO,NL),
      1 TOTAL(NGRO),SIGF(NGRO,0:MAXDEL),CHI(NGRO,0:MAXDEL),
-     2 SADD(NGRO,NED),GOLD(NGRO),ZNPHI(NGRO))
+     2 SADD(NGRO,NED),ENER(NGRO+1),GOLD(NGRO),ZNPHI(NGRO))
       ALLOCATE(LSCAT(NL),LADD(NED))
 *----
 *  RECOVER THE GROUP STRUCTURE.
@@ -112,7 +112,6 @@
          WRITE (IOUT,'(40H LIBDRA: NUMBER OF ISOTOPES IN MICROLIB=,I6)')
      1   NBISO
       ENDIF
-      ALLOCATE(ENER(NGRO+1))
       CALL LCMLEN(IPDRL,'ENERGY',LENGT,ITYLCM)
       LENGT=LENGT-1
       IF(LENGT.NE.NGRO) CALL XABORT('LIBDRA: INVALID GROUP STRUCTURE.')
@@ -128,7 +127,6 @@
       ENDIF
       CALL LCMPUT(IPLIB,'ENERGY',NGRO+1,2,ENER)
       CALL LCMPUT(IPLIB,'DELTAU',NGRO,2,DELTA)
-      DEALLOCATE(ENER)
       CALL LCMLEN(IPDRL,'CHI-LIMITS',NBESP,ITYLCM)
       IF(NBESP.GT.0) THEN
          NBESP=NBESP-1
@@ -210,11 +208,11 @@
             NDEL=MAX(NDEL,NDEL0)
             IF(NDEL0.GT.MAXDEL) CALL XABORT('LIBDRA: MAXDEL OVERFLOW.')
             IF(NDEL0.GT.0) CALL LCMGET (IPDRL,'LAMBDA-D',ZLAMB)
-            CALL LIBDRB (IPDRL,NGRO,NL,NDEL0,NBESP,SN(1,IMX),SB(1,IMX),
-     1      NED,HVECT,DELTA,LBIN,NFS,EBIN,AWR(IMX),DELECC,IGECCO,IMPX,
-     2      NGF,NGFR,LSCAT,LSIGF,LADD,LGOLD,SIGS(1,1),SCAT(1,1,1),TOTAL,
-     3      ZNPHI,SIGF(1,1),CHI(1,1),CHI4G(1,1),SADD(1,1),GOLD(1),
-     4      BIN(1))
+            CALL LIBDRB (IPDRL,NGRO,NL,NDEL0,NBESP,ENER,SN(1,IMX),
+     1      SB(1,IMX),NED,HVECT,DELTA,LBIN,NFS,EBIN,AWR(IMX),DELECC,
+     2      IGECCO,IMPX,NGF,NGFR,LSCAT,LSIGF,LADD,LGOLD,SIGS(1,1),
+     3      SCAT(1,1,1),TOTAL,ZNPHI,SIGF(1,1),CHI(1,1),CHI4G(1,1),
+     4      SADD(1,1),GOLD(1),BIN(1))
          ELSE
 *----
 *  PERFORM TEMPERATURE LAGRANGIAN INTERPOLATION (ORDER ABS(NOTX)).
@@ -263,10 +261,11 @@
      >      1P,E12.4,18H KELVIN. FACTOR = ,E12.4)') TEMP(ITM),TERPM
             WRITE (CD,'(I4.4)') ITM
             CALL LCMSIX (IPDRL,'SUBTMP'//CD,1)
-            CALL LIBDRB (IPDRL,NGRO,NL,NDEL0,NBESP,SN(1,IMX),SB(1,IMX),
-     1      NED,HVECT,DELTA,LBIN,NFS,EBIN,AWR(IMX),DELECC,IGECCO,IMPX,
-     2      NGF,NGFR,LSCAT,LSIGF,LADD,LGOLD,SIGS2(1),SCAT2(1),TOTAL2,
-     3      ZNPHI2,SIGF2(1),CHI2(1),CHI4G2(1),SADD2(1),GOLD2(1),BIN2(1))
+            CALL LIBDRB (IPDRL,NGRO,NL,NDEL0,NBESP,ENER,SN(1,IMX),
+     1      SB(1,IMX),NED,HVECT,DELTA,LBIN,NFS,EBIN,AWR(IMX),DELECC,
+     2      IGECCO,IMPX,NGF,NGFR,LSCAT,LSIGF,LADD,LGOLD,SIGS2(1),
+     3      SCAT2(1),TOTAL2,ZNPHI2,SIGF2(1),CHI2(1),CHI4G2(1),SADD2(1),
+     4      GOLD2(1),BIN2(1))
             CALL LCMSIX (IPDRL,' ',2)
             DO 130 IG=1,NGRO
               TOTAL(IG)=TOTAL(IG)+TERPM*TOTAL2(IG)
@@ -403,8 +402,8 @@
 *  SCRATCH STORAGE DEALLOCATION
 *----
       DEALLOCATE(LADD,LSCAT)
-      DEALLOCATE(CHI4G,ZNPHI,GOLD,SADD,CHI,SIGF,TOTAL,SCAT,SIGS,DELTA,
-     1 AWR)
+      DEALLOCATE(CHI4G,ZNPHI,GOLD,ENER,SADD,CHI,SIGF,TOTAL,SCAT,SIGS,
+     1 DELTA,AWR)
       DEALLOCATE(ITYPRO,NFS)
       RETURN
 *
diff --git a/Dragon/src/LIBDRB.f b/Dragon/src/LIBDRB.f
index 600f737..2712363 100644
--- a/Dragon/src/LIBDRB.f
+++ b/Dragon/src/LIBDRB.f
@@ -1,7 +1,7 @@
 *DECK LIBDRB
-      SUBROUTINE LIBDRB (IPDRL,NGRO,NL,NDEL,NBESP,SN,SB,NED,HVECT,DELTA,
-     1 LBIN,NFS,BENER,AWR,DELECC,IGECCO,IMPX,NGF,NGFR,LSCAT,LSIGF,LADD,
-     2 LGOLD,SIGS,SCAT,TOTAL,ZNPHI,SIGF,CHI,CHI4G,SADD,GOLD,BIN)
+      SUBROUTINE LIBDRB (IPDRL,NGRO,NL,NDEL,NBESP,ENER,SN,SB,NED,HVECT,
+     1 DELTA,LBIN,NFS,BENER,AWR,DELECC,IGECCO,IMPX,NGF,NGFR,LSCAT,LSIGF,
+     2 LADD,LGOLD,SIGS,SCAT,TOTAL,ZNPHI,SIGF,CHI,CHI4G,SADD,GOLD,BIN)
 *
 *-----------------------------------------------------------------------
 *
@@ -25,13 +25,14 @@
 *         NL=1 or higher.
 * NDEL    number of delayed precursor groups.
 * NBESP   number of energy-dependent fission spectra.
+* ENER    energy limits of the coarse groups.
 * SN      dilution cross section in each energy group. A value of
 *         1.0E10 is used for infinite dilution.
 * SB      dilution cross section as used in Livolant and Jeanpierre
 *         normalization.
 * NED     number of extra vector edits.
 * HVECT   names of the extra vector edits.
-* DELTA   lethargy widths.
+* DELTA   lethargy widths of the coarse groups.
 * LBIN    number of fine groups.
 * NFS     number of fine groups per coarse group.
 * BENER   energy limits of the fine groups.
@@ -75,8 +76,8 @@
       CHARACTER*(*) HVECT(NED)
       TYPE(C_PTR) IPDRL
       INTEGER NGRO,NL,NDEL,NBESP,NED,LBIN,NFS(NGRO),IGECCO,IMPX,NGF,NGFR
-      REAL SN(NGRO),SB(NGRO),DELTA(NGRO),BENER(LBIN+1),AWR,DELECC,
-     1 SIGS(NGRO,NL),SCAT(NGRO,NGRO,NL),TOTAL(NGRO),ZNPHI(NGRO),
+      REAL ENER(NGRO+1),SN(NGRO),SB(NGRO),DELTA(NGRO),BENER(LBIN+1),AWR,
+     1 DELECC,SIGS(NGRO,NL),SCAT(NGRO,NGRO,NL),TOTAL(NGRO),ZNPHI(NGRO),
      2 SIGF(NGRO,0:NDEL),CHI(NGRO,0:NDEL),CHI4G(NGRO,NBESP),
      3 SADD(NGRO,NED),GOLD(NGRO),BIN(LBIN,3)
       LOGICAL LSCAT(NL),LSIGF,LADD(NED),LGOLD
@@ -84,7 +85,7 @@
 *  LOCAL VARIABLES
 *----
       CHARACTER CM*2,CD*4,HSMG*131,HNUSIG*12,HCHI*12,HTOTAL*5
-      PARAMETER (IOUT=6,MAXTRA=10000)
+      PARAMETER (IOUT=6)
       INTEGER KTOTLR,KSIGFR,KCHIR,KPHIR
       LOGICAL LPCAT
       DOUBLE PRECISION TMP,ZNGAR,SQD,SQ0,SQ1,SQ2,SQ3,FACT1,FACT2
@@ -93,8 +94,8 @@
 *  ALLOCATABLE ARRAYS
 *----
       INTEGER, ALLOCATABLE, DIMENSION(:) :: NJJ,IJJ,KADDR
-      REAL, ALLOCATABLE, DIMENSION(:) :: GAR,PRI,STIS,UUU,SSS
-      REAL, ALLOCATABLE, DIMENSION(:,:) :: TERP,SIGT
+      REAL, ALLOCATABLE, DIMENSION(:) :: GAR
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: TERP,SIGT,GAR2D
       LOGICAL, ALLOCATABLE, DIMENSION(:) :: LSDIL,LPDIL,LINF
 *----
 *  SCRATCH STORAGE ALLOCATION
@@ -148,9 +149,9 @@
       IF(LPCAT.AND.(IGECCO.EQ.0)) THEN
          CALL LCMGET(IPDRL,'NJJS'//CM,NJJ)
          CALL LCMGET(IPDRL,'IJJS'//CM,IJJ)
-         LENGT2=0
+         LENGT=0
          DO 20 I=1,NGRO
-         LENGT2=LENGT2+NJJ(I)
+         LENGT=LENGT+NJJ(I)
    20    CONTINUE
          GAR(:LENGT)=0.0
          CALL LCMGET(IPDRL,'SCAT'//CM,GAR)
@@ -163,41 +164,9 @@
    30    CONTINUE
    40    CONTINUE
       ELSE IF(LPCAT) THEN
-         ! on-flight elastic scattering kernel
-         CALL LCMGET(IPDRL,'NJJS'//CM,NJJ)
-         CALL LCMGET(IPDRL,'IJJS'//CM,IJJ)
-         ALLOCATE(PRI(MAXTRA),STIS(NGRO),UUU(NGRO),SSS(IGECCO))
-         CALL LIBPRI(MAXTRA,DELECC,AWR,0,IL,NPRI,PRI)
-         LENGT2=0
-         DO 50 I=1,NGRO
-         LENGT2=LENGT2+NJJ(I)
-   50    CONTINUE
-         GAR(:LENGT)=0.0
-         CALL LCMGET(IPDRL,'SCAT'//CM,GAR)
-         UUU(1)=DELTA(1)
-         DO 60 I=2,NGRO
-         UUU(I)=UUU(I-1)+DELTA(I)
-   60    CONTINUE
-         IGAR=0
-*        IG2 IS THE SECONDARY GROUP.
-         DO 90 IG2=1,NGRO
-         IF(IG2.LE.IGECCO) THEN
-           CALL LIBECT(MAXTRA,IG2,PRI,UUU,DELECC,DELTA,NPRI,1,MML,STIS)
-           IGAR=IGAR+NJJ(IG2)
-           SSS(IG2)=GAR(IGAR)
-           DO 70 I=1,MML
-           IG1=IG2-I+1
-           IF(IG1.LE.0) GO TO 90
-           SCAT(IG2,IG1,IL+1)=STIS(I)*SSS(IG1)
-   70      CONTINUE
-         ELSE
-           DO 80 IG1=IJJ(IG2),IJJ(IG2)-NJJ(IG2)+1,-1
-           IGAR=IGAR+1
-           SCAT(IG2,IG1,IL+1)=GAR(IGAR)
-   80      CONTINUE
-         ENDIF
-   90    CONTINUE
-         DEALLOCATE(SSS,UUU,STIS,PRI)
+         ! on-flight scattering kernel
+         CALL LIBECC(IPDRL,NGRO,IL,AWR,ENER,DELTA,DELECC,IGECCO,
+     1   SCAT(1,1,IL+1))
       ENDIF
       CALL LCMLEN(IPDRL,'SIGS'//CM,LENGT,ITYLCM)
       LSCAT(IL+1)=(LENGT.GT.0)
@@ -320,10 +289,10 @@
             WRITE (CM,'(I2.2)') IL
             CALL LCMLEN(IPDRL,'SCAT'//CM,LENGT,ITYLCM)
             IF(.NOT.LSDIL(IL+1))
-     >        LSDIL(IL+1)=(LENGT.GT.0).AND.LSCAT(IL+1)
+     1        LSDIL(IL+1)=(LENGT.GT.0).AND.LSCAT(IL+1)
             CALL LCMLEN(IPDRL,'SIGS'//CM,LENGT,ITYLCM)
             IF(.NOT.LPDIL(IL+1))
-     >        LPDIL(IL+1)=(LENGT.GT.0).AND.LSCAT(IL+1)
+     1        LPDIL(IL+1)=(LENGT.GT.0).AND.LSCAT(IL+1)
   220    CONTINUE
          CALL LCMLEN(IPDRL,HTOTAL,LENGT,ITYLCM)
          KTOTLR=MAX(KTOTLR,LENGT)
@@ -562,44 +531,18 @@
   480       CONTINUE
   490       CONTINUE
          ELSE IF(LSDIL(IL+1)) THEN
-            ! on-flight elastic scattering kernel
-            CALL LCMGET(IPDRL,'NJJS'//CM,NJJ)
-            CALL LCMGET(IPDRL,'IJJS'//CM,IJJ)
-            ALLOCATE(PRI(MAXTRA),STIS(NGRO),UUU(NGRO),SSS(IGECCO))
-            CALL LIBPRI(MAXTRA,DELECC,AWR,0,IL,NPRI,PRI)
-            LENGT2=0
-            DO 500 I=1,NGRO
-            LENGT2=LENGT2+NJJ(I)
-  500       CONTINUE
-            GAR(:LENGT)=0.0
-            CALL LCMGET(IPDRL,'SCAT'//CM,GAR)
-            UUU(1)=DELTA(1)
-            DO 510 I=2,NGRO
-            UUU(I)=UUU(I-1)+DELTA(I)
-  510       CONTINUE
-            IGAR=0
-*           IG2 IS THE SECONDARY GROUP.
+            ! on-flight scattering kernel
+            ALLOCATE(GAR2D(NGRO,NGRO))
+            CALL LIBECC(IPDRL,NGRO,IL,AWR,ENER,DELTA,DELECC,IGECCO,
+     1      GAR2D)
+            DO 550 IG1=1,NGRO
+            FNTRP=TERP(IDIL,IG1)
             DO 540 IG2=1,NGRO
-            IF(IG2.LE.IGECCO) THEN
-              CALL LIBECT(MAXTRA,IG2,PRI,UUU,DELECC,DELTA,NPRI,1,MML,
-     1        STIS)
-              IGAR=IGAR+NJJ(IG2)
-              SSS(IG2)=GAR(IGAR)
-              DO 520 I=1,MML
-              IG1=IG2-I+1
-              IF(IG1.LE.0) GO TO 540
-              SCAT(IG2,IG1,IL+1)=SCAT(IG2,IG1,IL+1)+TERP(IDIL,IG1)*
-     1        STIS(I)*SSS(IG1)
-  520         CONTINUE
-            ELSE
-              DO 530 IG1=IJJ(IG2),IJJ(IG2)-NJJ(IG2)+1,-1
-              IGAR=IGAR+1
-              SCAT(IG2,IG1,IL+1)=SCAT(IG2,IG1,IL+1)+TERP(IDIL,IG1)*
-     1        GAR(IGAR)
-  530         CONTINUE
-            ENDIF
+            ! IG2 is the secondary group
+            SCAT(IG2,IG1,IL+1)=SCAT(IG2,IG1,IL+1)+FNTRP*GAR2D(IG2,IG1)
   540       CONTINUE
-            DEALLOCATE(SSS,UUU,STIS,PRI)
+  550       CONTINUE
+            DEALLOCATE(GAR2D)
          ENDIF
          IF(LPDIL(IL+1)) THEN
             GAR(:NGRO)=0.0
diff --git a/Dragon/src/LIBECC.f b/Dragon/src/LIBECC.f
new file mode 100644
index 0000000..b0cea7d
--- /dev/null
+++ b/Dragon/src/LIBECC.f
@@ -0,0 +1,148 @@
+*DECK LIBECC
+      SUBROUTINE LIBECC(IPDRL,NGRO,IL,AWR,ENER,DELTA,DELECC,IGECCO,
+     > SCAT)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Construct the scattering matrix using analytical scattering kernels.
+*
+*Copyright:
+* Copyright (C) 2025 Ecole Polytechnique de Montreal
+* This library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public
+* License as published by the Free Software Foundation; either
+* version 2.1 of the License, or (at your option) any later version
+*
+*Author(s): A. Hebert
+*
+*Parameters: input
+* IPDRL   pointer to the draglib (L_DRAGLIB signature).
+* NGRO    number of energy groups.
+* IL      Legendre order (=0: isotropic kernel in LAB).
+* AWR     mass ratio for current isotope.
+* ENER    energy limits of the coarse groups.
+* DELTA   lethargy widths of the coarse groups.
+* DELECC  lethargy widths of eccolib libraries.
+* IGECCO  number of equal-width lethargy groups with eccolib libraries.
+* IMPX    print flag.
+*
+*Parameters: output
+* SCAT    scattering matrix.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPDRL
+      INTEGER NGRO,IL,IGECCO
+      REAL AWR,ENER(NGRO+1),DELTA(NGRO),DELECC,SCAT(NGRO,NGRO)
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER(MAXEDI=47,MAXTRA=10000)
+      CHARACTER CM*2
+      CHARACTER(LEN=8), SAVE, DIMENSION(MAXEDI) :: NAMEDI=
+     >  (/ 'NELAS   ','N2N     ','N3N     ','NNP     ','N4N     ',
+     >     'NINEL001','NINEL002','NINEL003','NINEL004','NINEL005',
+     >     'NINEL006','NINEL007','NINEL008','NINEL009','NINEL010',
+     >     'NINEL011','NINEL012','NINEL013','NINEL014','NINEL015',
+     >     'NINEL016','NINEL017','NINEL018','NINEL019','NINEL020',
+     >     'NINEL021','NINEL022','NINEL023','NINEL024','NINEL025',
+     >     'NINEL026','NINEL027','NINEL028','NINEL029','NINEL030',
+     >     'NINEL031','NINEL032','NINEL033','NINEL034','NINEL035',
+     >     'NINEL036','NINEL037','NINEL038','NINEL039','NINEL040',
+     >     'NINEL041','NINEL   '/)
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: NJJ,IJJ
+      REAL, ALLOCATABLE, DIMENSION(:) :: GAR,PRI,STIS,UUU,QQ
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: SSS2
+      LOGICAL, ALLOCATABLE, DIMENSION(:) :: LPAR
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(NJJ(NGRO),IJJ(NGRO),GAR(NGRO*NGRO))
+      ALLOCATE(PRI(MAXTRA),STIS(NGRO),UUU(NGRO))
+      ALLOCATE(LPAR(MAXEDI),SSS2(NGRO,MAXEDI),QQ(MAXEDI))
+*----
+*  RECOVER CROSS SECTIONS CONTRIBUTING TO THE SCATTERING MATRIX
+*----
+      SSS2(:NGRO,:MAXEDI)=0.0
+      LPAR(:MAXEDI)=.FALSE.
+      QQ(:MAXEDI)=0.0
+      DO I=1,MAXEDI
+        CALL LCMLEN(IPDRL,NAMEDI(I),LENGT,ITYLCM)
+        IF(LENGT.GT.0) THEN
+          LPAR(I)=.TRUE.
+          CALL LCMGET(IPDRL,NAMEDI(I),SSS2(1,I))
+          DO IG1=1,NGRO
+            IF(NAMEDI(I).EQ.'N2N') THEN
+              SSS2(IG1,I)=2.0*SSS2(IG1,I)
+            ELSE IF(NAMEDI(I).EQ.'N3N') THEN
+              SSS2(IG1,I)=3.0*SSS2(IG1,I)
+            ELSE IF(NAMEDI(I).EQ.'N4N') THEN
+              SSS2(IG1,I)=4.0*SSS2(IG1,I)
+            ENDIF
+          ENDDO
+          DO IG1=NGRO,1,-1
+            IF(SSS2(IG1,I).NE.0.0) EXIT
+            QQ(I)=-ENER(IG1)
+          ENDDO
+        ENDIF
+      ENDDO
+*----
+*  CONSTRUCT THE SCATTERING MATRIX
+*----
+      WRITE (CM,'(I2.2)') IL
+      CALL LCMGET(IPDRL,'NJJS'//CM,NJJ)
+      CALL LCMGET(IPDRL,'IJJS'//CM,IJJ)
+      LENGT=0
+      DO IG1=1,NGRO
+        LENGT=LENGT+NJJ(IG1)
+      ENDDO
+      GAR(:LENGT)=0.0
+      CALL LCMGET(IPDRL,'SCAT'//CM,GAR)
+      UUU(1)=DELTA(1)
+      DO IG1=2,NGRO
+        UUU(IG1)=UUU(IG1-1)+DELTA(IG1)
+      ENDDO
+      IGAR=0
+      SCAT(:NGRO,:NGRO)=0.0
+      DO IG1=1,IGECCO
+        DO I=1,MAXEDI
+          IF(LPAR(I)) THEN
+            IF(NAMEDI(I).EQ.'NELAS') THEN
+              CALL LIBPRI(MAXTRA,DELECC,AWR,0,IL,NPRI,PRI)
+            ELSE
+              ! treshold reaction
+              IF(ENER(IG1).LE.-QQ(I)*(AWR+1.0)/AWR) CYCLE
+              CALL LIBPRQ(MAXTRA,DELECC,AWR,ENER(IG1),QQ(I),0,IL,
+     >        NPRI,PRI)
+            ENDIF
+            DO IPRI=1,NPRI
+              IG2=IG1+IPRI-1 ! IG2 <-- IG1
+              IF(IG2.GT.IGECCO) EXIT
+              SCAT(IG2,IG1)=SCAT(IG2,IG1)+PRI(IPRI)*SSS2(IG1,I)
+            ENDDO
+          ENDIF
+        ENDDO
+        IGAR=IGAR+NJJ(IG1)
+      ENDDO ! IG1
+      DO IG2=IGECCO+1,NGRO
+        DO IG1=IJJ(IG2),IJJ(IG2)-NJJ(IG2)+1,-1
+          IGAR=IGAR+1
+          SCAT(IG2,IG1)=GAR(IGAR)
+        ENDDO
+      ENDDO ! IG2
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(QQ,SSS2,LPAR)
+      DEALLOCATE(UUU,STIS,PRI)
+      DEALLOCATE(GAR,IJJ,NJJ)
+      RETURN
+      END
diff --git a/Dragon/src/LIBPRI.f b/Dragon/src/LIBPRI.f
index 9461406..0589a48 100644
--- a/Dragon/src/LIBPRI.f
+++ b/Dragon/src/LIBPRI.f
@@ -22,7 +22,7 @@
 * AWR     mass ratio for current isotope.
 * IALTER  type of approximation (=0: use exponentials; =1: use Taylor
 *         expansions).
-* IL      Legendre order (=0: isotropic kernel).
+* IL      Legendre order (=0: isotropic kernel in LAB).
 *
 *Parameters: output
 * N       exact dimension of array PRI.
diff --git a/Dragon/src/LIBPRQ.f b/Dragon/src/LIBPRQ.f
new file mode 100644
index 0000000..08d9b40
--- /dev/null
+++ b/Dragon/src/LIBPRQ.f
@@ -0,0 +1,190 @@
+*DECK LIBPRQ
+      SUBROUTINE LIBPRQ(MAXTRA,DELI,AWR,E0,Q,IALTER,IL,N,PRI)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute the PRI array for various Legendre orders using Gaussian
+* integration. Inelastic scattering case.
+*
+*Copyright:
+* Copyright (C) 2025 Ecole Polytechnique de Montreal
+* This library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public
+* License as published by the Free Software Foundation; either
+* version 2.1 of the License, or (at your option) any later version
+*
+*Author(s): A. Hebert
+*
+*Parameters: input
+* MAXTRA  allocated dimension of array PRI.
+* DELI    elementary lethargy width of the equi-width lethargy mesh.
+* AWR     mass ratio for current isotope.
+* E0      energy corresponding to the upper limit of primary group.
+* Q       Q-value (negative value) for an inelastic diffusion.
+* IALTER  type of approximation (=0: use exponentials; =1: use Taylor
+*         expansions).
+* IL      Legendre order (=0: isotropic kernel in LAB).
+*
+*Parameters: output
+* N       exact dimension of array PRI.
+* PRI     array containing the slowing-down probabilities defined on
+*         an equi-width lethargy mesh.
+*
+*Reference:
+* M. Grandotto-Biettoli, "AUTOSECOL, un calcul automatique de
+* l'auto-protection des resonances des isotopes lourds,"
+* Note CEA-N-1961, Commissariat a l'Energie Atomique, 1977. 
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER MAXTRA,IALTER,IL,N
+      REAL DELI,AWR,E0,Q,PRI(MAXTRA)
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER(NGPT=6,MAXNL=50)
+      DOUBLE PRECISION AWRB,ALP,T0,FACT,COEF0,GAM,UM,UMIN,UMAX
+      REAL UI(NGPT),WI(NGPT),UJ(NGPT),WJ(NGPT)
+      REAL POLY(0:MAXNL),CALC(0:MAXNL,0:2)
+      CHARACTER HSMG*131
+      ZMU(AWR,RGAM,U)=0.5*(AWR+1.0)*EXP(-0.5*U)-0.5*(RGAM/(AWR+1.0)+
+     1 AWR-1.0)*EXP(0.5*U)
+*----
+*  COMPUTE THE LEGENDRE POLYNOMIAL OF ORDER IL.
+*----
+      IF(IL.GT.MAXNL) CALL XABORT('LIBPRQ: IL OVERFLOW.')
+      IF(IL.EQ.0) THEN
+        POLY(0)=1.0
+      ELSE IF(IL.EQ.1) THEN
+        POLY(0)=0.0
+        POLY(1)=1.0
+      ELSE
+        CALC(0:MAXNL,0:1)=0.0
+        CALC(0,0)=1.0
+        CALC(1,1)=1.0
+        DO J=2,IL
+          DO I=0,IL
+            T0=-REAL(J-1)*CALC(I,MOD(J-2,3))
+            IF(I.GT.0) T0=T0+(2.0*REAL(J-1)+1.0)*CALC(I-1,MOD(J-1,3))
+            CALC(I,MOD(J,3))=REAL(T0)/REAL(J)
+          ENDDO
+        ENDDO
+        POLY(0:IL)=CALC(0:IL,MOD(IL,3))
+      ENDIF
+*
+      AWRB=AWR
+      IF(AWR.LT.1.0001) AWRB=1.0001
+      GAM=AWRB*(AWRB+1.D0)*Q/E0
+      UMIN=LOG((AWRB+1.D0)**2/(AWRB**2+1.D0+GAM+2.D0*SQRT(AWRB**2+GAM)))
+      IF(-GAM.GT.AWRB**2) CALL XABORT('LIBPRQ: NEGATIVE SQRT ARGUMENT.')
+      IF(UMIN.LT.DELI) THEN
+        ! pathological case where the treshold is negligible
+        CALL LIBPRI(MAXTRA,DELI,AWR,IALTER,IL,N,PRI)
+        RETURN
+      ENDIF
+      PRI(:MAXTRA)=0.0
+      ALP=((AWRB-1.D0)/(AWRB+1.D0))**2
+      CALL ALGPT(NGPT,0.0,DELI,UI,WI)
+      N=0
+      DO I=1,NGPT
+        ! primary group base point
+        WII=WI(I)
+        UII=UI(I)
+        EN=E0*EXP(-UI(I)*DELI)
+        GAM=AWRB*(AWRB+1.D0)*Q/EN
+        UM=AWRB**2+1.D0+GAM
+        UMIN=LOG((AWRB+1.D0)**2/(UM+2.D0*SQRT(AWRB**2+GAM)))
+        UMAX=LOG((AWRB+1.D0)**2/(UM-2.D0*SQRT(AWRB**2+GAM)))
+        NMIN=1+INT((UMIN-1.D-6)/DELI)
+        NMAX=1+INT((UMAX-1.D-6)/DELI)
+        IF(NMAX.GT.MAXTRA) THEN
+          WRITE(HSMG,'(25HLIBPRQ: MAXTRA OVERFLOW (,I8,2H >,I8,2H).)')
+     1    NMAX,MAXTRA
+          CALL XABORT(HSMG)
+        ENDIF
+        COEF0=AWRB/SQRT(AWRB**2+GAM)
+        WII=WII*REAL(COEF0)
+        IF(NMAX.EQ.NMIN) THEN
+          CALL ALGPT(NGPT,REAL(UMIN),REAL(UMAX),UJ,WJ)
+          DO J=1,NGPT
+            FACT=POLY(0)
+            T0=1.0D0
+            DO K=1,IL
+              T0=T0*ZMU(AWR,REAL(GAM),UJ(J)-UII)
+              FACT=FACT+POLY(K)*T0
+            ENDDO
+            IF(IALTER.EQ.0) THEN
+              PRI(NMIN)=PRI(NMIN)+WII*WJ(J)*EXP(UII-UJ(J))*REAL(FACT)
+            ELSE
+              PRI(NMIN)=PRI(NMIN)+WII*WJ(J)*REAL(FACT/(UMAX-UMIN))
+            ENDIF
+          ENDDO ! J
+        ELSE
+          CALL ALGPT(NGPT,REAL(UMIN),NMIN*DELI,UJ,WJ)
+          DO J=1,NGPT
+            FACT=POLY(0)
+            T0=1.0D0
+            DO K=1,IL
+              T0=T0*ZMU(AWR,REAL(GAM),UJ(J)-UII)
+              FACT=FACT+POLY(K)*T0
+            ENDDO
+            IF(IALTER.EQ.0) THEN
+              PRI(NMIN)=PRI(NMIN)+WII*WJ(J)*EXP(UII-UJ(J))*REAL(FACT)
+            ELSE
+              PRI(NMIN)=PRI(NMIN)+WII*WJ(J)*REAL(FACT/(UMAX-UMIN))
+            ENDIF
+          ENDDO ! J
+          DO N=NMIN+1,NMAX-1
+            CALL ALGPT(NGPT,(N-1)*DELI,N*DELI,UJ,WJ)
+            DO J=1,NGPT
+              FACT=POLY(0)
+              T0=1.0D0
+              DO K=1,IL
+                T0=T0*ZMU(AWR,REAL(GAM),UJ(J)-UII)
+                FACT=FACT+POLY(K)*T0
+              ENDDO
+              IF(IALTER.EQ.0) THEN
+                PRI(N)=PRI(N)+WII*WJ(J)*EXP(UII-UJ(J))*REAL(FACT)
+              ELSE
+                PRI(N)=PRI(N)+WII*WJ(J)*REAL(FACT/(UMAX-UMIN))
+              ENDIF
+            ENDDO ! J
+          ENDDO ! N
+          CALL ALGPT(NGPT,(NMAX-1)*DELI,REAL(UMAX),UJ,WJ)
+          DO J=1,NGPT
+            FACT=POLY(0)
+            T0=1.0D0
+            DO K=1,IL
+              T0=T0*ZMU(AWR,REAL(GAM),UJ(J)-UII)
+              FACT=FACT+POLY(K)*T0
+            ENDDO
+            IF(IALTER.EQ.0) THEN
+              PRI(NMAX)=PRI(NMAX)+WII*WJ(J)*EXP(UII-UJ(J))*REAL(FACT)
+            ELSE
+              PRI(NMAX)=PRI(NMAX)+WII*WJ(J)*REAL(FACT/(UMAX-UMIN))
+            ENDIF
+          ENDDO ! J
+        ENDIF
+        N=MAX(N,NMAX)
+      ENDDO ! I
+      IF(IALTER.EQ.0) THEN
+        PRI(:N)=PRI(:N)/DELI/REAL(1.0D0-ALP)
+      ELSE
+        PRI(:N)=PRI(:N)/DELI
+      ENDIF
+*----
+*  NORMALIZATION
+*----
+      IF(IL.EQ.0) THEN
+        FACT=0.0D0
+        DO I=1,N
+          FACT=FACT+PRI(I)
+        ENDDO
+        PRI(:N)=PRI(:N)/REAL(FACT)
+      ENDIF
+      RETURN
+      END