From 7dfcc480ba1e19bd3232349fc733caef94034292 Mon Sep 17 00:00:00 2001
From: stainer_t <thomas.stainer@oecd-nea.org>
Date: Mon, 8 Sep 2025 13:48:49 +0200
Subject: Initial commit from Polytechnique Montreal

---
 Dragon/src/XDREXP.f | 106 ++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 106 insertions(+)
 create mode 100644 Dragon/src/XDREXP.f

(limited to 'Dragon/src/XDREXP.f')

diff --git a/Dragon/src/XDREXP.f b/Dragon/src/XDREXP.f
new file mode 100644
index 0000000..88ce471
--- /dev/null
+++ b/Dragon/src/XDREXP.f
@@ -0,0 +1,106 @@
+*DECK XDREXP
+      SUBROUTINE XDREXP(DX,NBX,PARAM,E00,E01,E10,E11)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Construct exponential tables for linear interpolation.
+*
+*Copyright:
+* Copyright (C) 2002 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): R. Roy, R. Le Tellier 
+*
+*Parameters: input
+* DX      step for tables (here, DX=0.02d0).
+* NBX     order of tables (here, NBX=7936).
+*
+*Parameters: output
+* PARAM   exponential table characteristics.
+* E00     exponential table.
+* E01     exponential table.
+* E10     exponential table.
+* E11     exponential table.
+*
+*Comments:
+*  Modified in order to tabulate (1-exp(-x))/x instead of (1-exp(-x))).
+*
+*-----------------------------------------------------------------------
+*
+      IMPLICIT NONE
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER          NBX
+      DOUBLE PRECISION DX
+      REAL             PARAM(3), E00(0:NBX), E01(0:NBX), E10(0:NBX),
+     >                 E11(0:NBX)
+*----
+*  LOCAL VARIABLES
+*----
+      INTEGER          I
+      DOUBLE PRECISION X0, X1, X2, X3, X4, PAS, XLIM, EX00, EX01,
+     >                 EX10, EX11
+      DOUBLE PRECISION DEPS, DREF, DZERO, DONE, TAUDMIN
+      PARAMETER      ( DEPS= 1.D-10, DREF= 1.D0/512.D0,
+     >                 DZERO= 0.D0, DONE=1.D0, TAUDMIN=2.D-2 )
+*----
+*  SPACE AND VALUES FOR EXPONENTIAL TABLES
+*----
+      INTEGER          MEX1
+      PARAMETER      ( MEX1=7936 )
+*
+      IF( NBX.NE.  MEX1 ) GO TO 97
+      IF( DX .GT. DREF+DEPS .OR. DX .LT. DREF-DEPS  ) GO TO 98
+      PAS=   DONE/DX
+      XLIM=  DBLE(NBX)*DX
+*----
+*  WE CONSTRUCT THE LINEAR TABLES USING ACCURATE *EXP* VALUES
+*----
+      X0=  DZERO
+      EX00= DZERO
+      EX10= DONE
+      DO 20 I= 0, NBX-1
+*
+*        STORE CONSTANT VALUE:
+         X1= X0 + DX
+         IF (X1.LE.TAUDMIN) THEN
+            X2=0.5D0*X1
+            X3=X1/3.D0
+            X4=0.5D0*X2
+            EX11=DONE-X2*(DONE-X3*(DONE-X4))
+         ELSE
+            EX11=(DONE - EXP(-X1))/X1
+         ENDIF
+         EX01=(DONE - EXP(-X1))
+*
+*        STORE STEP AND CONSTANT VALUES:
+         X2=(EX01-EX00)/(X1-X0)
+         X3=(EX11-EX10)/(X1-X0)
+         E01(I)= REAL(X2)
+         E00(I)= REAL(EX00-X0*X2)
+         E11(I)= REAL(X3)
+         E10(I)= REAL(EX10-X0*X3)
+         X0= X1
+         EX00= EX01
+         EX10= EX11
+   20 CONTINUE
+      PARAM(1)=REAL(PAS)
+      PARAM(2)=REAL(DX)
+      PARAM(3)=REAL(XLIM)
+      E10(NBX)= REAL(1.D0/XLIM)
+      E11(NBX)= 0.0
+      E00(NBX)= 1.0
+      E01(NBX)= 0.0
+      RETURN
+*----
+*  ERROR SECTION
+*----
+   97 CALL XABORT('XDREXP: EXP LINEAR TABLES HAVE MORE THAN 7936 ELEM'
+     > //'ENTS')
+   98 CALL XABORT('XDREXP: EXP LINEAR TABLES HAVE A STEP OF 1/512')
+      END
-- 
cgit v1.2.3