From 72ca4af0f99fbce82d143b99649669923b4b1161 Mon Sep 17 00:00:00 2001
From: Wojtek Kosior <kwojtus@protonmail.com>
Date: Wed, 3 Jul 2019 21:12:59 +0200
Subject: return NaN for wrong poly order argument and handle this error upper
 in the call stack

---
 src/main.f90        |  7 +++++++
 src/quadratures.f90 | 15 ++++++++-------
 2 files changed, 15 insertions(+), 7 deletions(-)

(limited to 'src')

diff --git a/src/main.f90 b/src/main.f90
index c0ecd48..e49603e 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -29,6 +29,7 @@ PROGRAM integrator
   USE quadratures
   USE functions
   USE iso_fortran_env, ONLY: error_unit
+  USE, intrinsic :: ieee_arithmetic
   IMPLICIT none
 
   real(kind=8) :: ibeg, iend, val
@@ -112,6 +113,12 @@ PROGRAM integrator
 
   val = numerical_int(ibeg, iend, fun, poly_order)
 
+  IF (ieee_is_nan(val)) THEN
+     errmsg = "invalid 5th argument (polynomial order)"
+     isok = .false.
+     GOTO 1
+  END IF
+  
   if (this_image() == 1) write (*,*) val
 
   1 if (this_image() == 1 .and. .not. isok) write(*,*) trim(errmsg)
diff --git a/src/quadratures.f90 b/src/quadratures.f90
index 06d4c6c..6fbb5cb 100644
--- a/src/quadratures.f90
+++ b/src/quadratures.f90
@@ -131,6 +131,8 @@ CONTAINS
   END FUNCTION integrate_generic
 
   FUNCTION newton_cotes(ibeg, iend, fun, p) result(val)
+    USE, intrinsic :: ieee_arithmetic
+
     IMPLICIT none
     real(kind=8), intent(in) :: ibeg
     real(kind=8), intent(in) :: iend
@@ -141,8 +143,6 @@ CONTAINS
     procedure(quadrature), pointer :: quad
     
     SELECT CASE (p)
-    CASE (:-1)
-       STOP "negative interpolationg polynomial order passed"
     CASE (0)
        quad => rectangle
     CASE (1)
@@ -150,7 +150,8 @@ CONTAINS
     CASE (2)
        quad => simpson_1_third
     CASE default
-       STOP "Newton-Cotes quadratures only implemented for order < 3"
+       val = ieee_value(val, ieee_quiet_nan)
+       RETURN
     END SELECT
 
     val = integrate_generic(ibeg, iend, fun, quad)
@@ -185,6 +186,8 @@ CONTAINS
   END FUNCTION simpson_1_third
 
   FUNCTION gauss(ibeg, iend, fun, p) result(val)
+    USE, intrinsic :: ieee_arithmetic
+
     IMPLICIT none
     real(kind=8), intent(in) :: ibeg
     real(kind=8), intent(in) :: iend
@@ -195,8 +198,6 @@ CONTAINS
     procedure(quadrature), pointer :: quad
     
     SELECT CASE (p)
-    CASE (:0)
-       STOP "non-positive Legendre polynomial order passed"
     CASE (1)
        quad => gauss_n1
     CASE (2)
@@ -204,8 +205,8 @@ CONTAINS
     CASE (3)
        quad => gauss_n3
     CASE default
-       STOP "Gauss quadratures only implemented with roots of" &
-            // " Legendgre polynomial of order <= 3"
+       val = ieee_value(val, ieee_quiet_nan)
+       RETURN
     END SELECT
 
     val = integrate_generic(ibeg, iend, fun, quad)
-- 
cgit v1.2.3