Merge branch '16-implement-analytic-inelastic-scattering-laws-for-draglibs-in-module-lib' into 'main'

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

See merge request dragon/5.1!29

1 files changed, 190 insertions, 0 deletions

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