1 files changed, 99 insertions, 0 deletions
diff --git a/Dragon/src/LIBLEX.f b/Dragon/src/LIBLEX.f
new file mode 100644
index 0000000..f7dfa30
--- /dev/null
+++ b/Dragon/src/LIBLEX.f
@@ -0,0 +1,99 @@
+*DECK LIBLEX
+      SUBROUTINE LIBLEX(NELE,PNTE,ELMT,NOTX,TERP)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute factors for Lagrangian 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): A. Hebert and G. Marleau
+*
+*Parameters: input
+* NELE    number of elements in table.
+* PNTE    extrapolation point.
+* ELMT    values of elements in table.
+* NOTX    order of interpolation:
+*         >0 order NOTX in square root of PNTE;
+*         <0 order -NOTX in PNTE.
+*
+*Parameters: output
+* TERP    extrapolation factor.
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER NELE,NOTX
+      REAL PNTE,ELMT(NELE)
+      DOUBLE PRECISION TERP(NELE)
+*----
+*  LOCAL VARIABLES
+*----
+      DOUBLE PRECISION FACTOR,TERM,SQPNTE
+*
+      I1=1
+      I2=NELE
+      NORD=ABS(NOTX)
+      TERP(:NELE)=0.0D0
+      DO 100 I=1,NELE
+        IF(PNTE.EQ.ELMT(I)) THEN
+           TERP(I)=1.0D0
+           RETURN
+        ELSE IF(PNTE.GT.ELMT(I)) THEN
+           I1=MAX(I-(NORD-1)/2,I1)
+        ELSE IF(PNTE.LT.ELMT(I)) THEN
+           I2=MIN(I+(NORD-1)/2,I2)
+        ENDIF
+ 100  CONTINUE
+*
+      IF(I1.EQ.1) THEN
+        I2=MIN(NELE,NORD+1)
+      ELSE IF(I2.EQ.NELE) THEN
+        I1=MAX(1,NELE-NORD)
+      ENDIF
+      FACTOR=1.0D0
+      IF(NOTX.LT.0) THEN
+        SQPNTE=DBLE(PNTE)
+        DO 110 I=I1,I2
+          IF(ABS(PNTE-ELMT(I)).LE.1.0E-5*ABS(PNTE)) THEN
+            TERP(I)=1.0D0
+            RETURN
+          ENDIF
+          FACTOR=FACTOR*( SQPNTE-DBLE(ELMT(I)) )
+ 110    CONTINUE
+        DO 120 I=I1,I2
+          TERM=FACTOR/( SQPNTE-DBLE(ELMT(I)) )
+          DO 130 J=I1,I2
+            IF(I.NE.J) TERM=TERM/
+     >               ( DBLE(ELMT(I))-DBLE(ELMT(J)) )
+ 130      CONTINUE
+          TERP(I)=TERM
+ 120    CONTINUE
+      ELSE
+        SQPNTE=DBLE(SQRT(PNTE))
+        DO 160 I=I1,I2
+          IF(ABS(PNTE-ELMT(I)).LE.1.0E-5*ABS(PNTE)) THEN
+            TERP(I)=1.0D0
+            RETURN
+          ENDIF
+          FACTOR=FACTOR*( SQPNTE-DBLE(SQRT(ELMT(I))) )
+ 160    CONTINUE
+        DO 170 I=I1,I2
+          TERM=FACTOR/( SQPNTE-DBLE(SQRT(ELMT(I))) )
+          DO 180 J=I1,I2
+            IF(I.NE.J) TERM=TERM/
+     >               ( DBLE(SQRT(ELMT(I)))-DBLE(SQRT(ELMT(J))) )
+ 180      CONTINUE
+          TERP(I)=TERM
+ 170    CONTINUE
+      ENDIF
+      RETURN
+      END