! Copyright 2019 Wojciech Kosior
  
! This is free and unencumbered software released into the public domain.

! Anyone is free to copy, modify, publish, use, compile, sell, or
! distribute this software, either in source code form or as a compiled
! binary, for any purpose, commercial or non-commercial, and by any
! means.

! In jurisdictions that recognize copyright laws, the author or authors
! of this software dedicate any and all copyright interest in the
! software to the public domain. We make this dedication for the benefit
! of the public at large and to the detriment of our heirs and
! successors. We intend this dedication to be an overt act of
! relinquishment in perpetuity of all present and future rights to this
! software under copyright law.

! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
! EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
! MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
! IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR
! OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
! ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
! OTHER DEALINGS IN THE SOFTWARE.

! For more information, please refer to <http://unlicense.org/>

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
  integer(kind=4) :: poly_order
  procedure(integrate), pointer :: numerical_int
  procedure(funint), pointer :: fun
  procedure(analytical_integral), pointer :: analytical_int
  character(100) :: arg, errmsg
  logical :: isok = .true.

  IF (command_argument_count() < 4) THEN
     errmsg = "at least 4 arguments required"
     isok = .false.
     GOTO 1
  END IF

  call get_command_argument(2, arg)

  IF      (arg == "exp")  THEN
     fun            => my_exp
     analytical_int => my_exp_int
  ELSE IF (arg == "sin")  THEN
     fun            => my_sin
     analytical_int => my_sin_int
  ELSE IF (arg == "poly") THEN
     fun            => my_poly
     analytical_int => my_poly_int
  ELSE
     errmsg =  "bad second argument (should be 'exp', 'sin' or 'poly')"
     isok = .false.
     GOTO 1
  END IF

  
  call get_command_argument(3, arg)

  read (arg,*) ibeg
  
  call get_command_argument(4, arg)

  read (arg,*) iend

  
  call get_command_argument(1, arg)

  IF      (arg == "gauss")        THEN
     numerical_int => gauss
  ELSE IF (arg == "newton-cotes") THEN
     numerical_int => newton_cotes
  ELSE IF (arg == "analytical")   THEN
     if (this_image() == 1) write (*,*) analytical_int(ibeg, iend)
     GOTO 1
  ELSE
     errmsg = "bad first argument (should be 'gauss',"            &
          // " 'newton-cotes' or 'analytical')"
     isok = .false.
     GOTO 1
  END IF

  
  IF (command_argument_count() < 5) THEN
     errmsg = "5th argument (polynomial order) required"          &
          // " for numerical integration"
     isok = .false.
     GOTO 1
  END IF
  
  call get_command_argument(5, arg)
  read (arg,*) poly_order

  
  IF (command_argument_count() > 5) THEN
     call get_command_argument(6, arg)
     read (arg,*) subintervals
     IF (subintervals < 1) THEN
        errmsg = "subintervals number must be positive"
        isok = .false.
        GOTO 1
     END IF
  END IF

  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 (*,'(es21.14," ",es21.14)')        &
       val, abs(val - analytical_int(ibeg, iend))

  1 if (this_image() == 1 .and. .not. isok) write(*,*) trim(errmsg)
  
CONTAINS

END PROGRAM integrator