add gauss quadratures, generic way (untestedgit add quadratures.f90)

1 files changed, 72 insertions, 0 deletions

diff --git a/src/quadratures.f90 b/src/quadratures.f90
index 43e6ac6..cc279ed 100644
--- a/src/quadratures.f90
+++ b/src/quadratures.f90
@@ -65,6 +65,33 @@ CONTAINS
     val = integrate_generic(ibeg, iend, fun, quad)
   END FUNCTION newton_cotes
   
+  FUNCTION gauss(ibeg, iend, fun, p) result(val)
+    IMPLICIT none
+    real(kind=8), intent(in) :: ibeg
+    real(kind=8), intent(in) :: iend
+    procedure(funint) :: fun
+    integer(kind=4), intent(in) :: p
+    real(kind=8) :: val
+    
+    procedure(quadrature), pointer :: quad
+    
+    SELECT CASE (p)
+    CASE (:0)
+       STOP "non-positive Legendre polynomial order passed"
+    CASE (1)
+       quad => gauss_n1
+    CASE (2)
+       quad => gauss_n2
+    CASE (3)
+       quad => gauss_n3
+    CASE default
+       STOP "Gauss quadratures only implemented with roots of" &
+            // " Legendgre polynomial of order <= 3"
+    END SELECT
+
+    val = integrate_generic(ibeg, iend, fun, quad)
+  END FUNCTION gauss
+  
   FUNCTION integrate_generic(ibeg, iend, fun, quad) result(val)
     IMPLICIT none
     real(kind=8), intent(in) :: ibeg
@@ -140,4 +167,49 @@ CONTAINS
          (fun(qbeg) + 4 * fun ((qbeg + qend) * 0.5) + fun(qend))
   END FUNCTION simpson_1_third
 
+  FUNCTION gauss_generic(mid, halfwidth, fun, points, weights)    &
+       result(val)
+    IMPLICIT none
+    real(kind=8), intent(in) :: mid, halfwidth, points(:), weights(:)
+    procedure(funint) :: fun
+    real(kind=8) :: val
+    
+    integer(kind=4) :: i
+
+    val = halfwidth * sum(weights * &
+         [(fun(points(i) * halfwidth + mid), i = 1, size(points))])
+  END FUNCTION gauss_generic
+  
+  FUNCTION gauss_n1(qbeg, qend, fun) result(val)
+    IMPLICIT none
+    real(kind=8), intent(in) :: qbeg, qend
+    procedure(funint) :: fun
+    real(kind=8) :: val, weights(1) = [2], points(1) = [0]
+
+    val = gauss_generic((qbeg + qend) * 0.5, (qend - qbeg) * 0.5, &
+         fun, points, weights)
+  END FUNCTION gauss_n1
+  
+  FUNCTION gauss_n2(qbeg, qend, fun) result(val)
+    IMPLICIT none
+    real(kind=8), intent(in) :: qbeg, qend
+    procedure(funint) :: fun
+    real(kind=8) :: val, weights(2) = [1, 1],                     &
+         points(2) = [1 / sqrt(3.0), -1 / sqrt(3.0)]
+         
+    val = gauss_generic((qbeg + qend) * 0.5, (qend - qbeg) * 0.5, &
+         fun, points, weights)
+  END FUNCTION gauss_n2
+  
+  FUNCTION gauss_n3(qbeg, qend, fun) result(val)
+    IMPLICIT none
+    real(kind=8), intent(in) :: qbeg, qend
+    procedure(funint) :: fun
+    real(kind=8) :: val, weights(3) = [8 / 9.0, 5 / 9.0, 5 / 9.0],&
+         points(3) = [0.0, sqrt(3 / 5.0), -sqrt(3 / 5.0)]
+
+    val = gauss_generic((qbeg + qend) * 0.5, (qend - qbeg) * 0.5, &
+         fun, points, weights)
+  END FUNCTION gauss_n3
+
 END MODULE quadratures