module Quanc8FMM
!!=============================================================================
!! Procedura Quanc8 
!! "Computer Methods for Mathematical  Computations",
!! by George E. Forsythe, Michael A. Malcolm, and Cleve B. Moler
!! Prentice-Hall, 1977.
!!-----------------------------------------------------------------------------
!! Prepis pro Fortran 90:
!! Last updated: 28 February 2005 by Ralph Carmichael ralph@pdas.com 
!! Public Domain Aeronautical Software 
!! P.O. Box 1438 Santa Cruz CA 95061 USA 
!!-----------------------------------------------------------------------------
!! Drobne upravy pro g95
!! Jan Cely, UFKL,Prirodovedecka fakulta MU, 270307,130208
!!=============================================================================
!!
!! SUBROUTINE Quanc8 (Fun,a,b,abserr,relerr,hodnota,errest,nofun,flag)
!!
!!  Estimate the integral of FUN(X) from A to B to a user provided 
!!  tolerance. This is an automatic adaptive routine based on the 8-panel 
!!  Newton-Cotes rule. 
!!	VSTUP:
!!	Fun		= integrand deklarovany takto : FUNCTION Fun(x) RESULT(f) 
!!	a,b		= integracni meze
!!	abserr	= maximum pro absolutni chybu (>=0)
!!	relerr	= hranice pro relativni chybu (>=0)
!!	VYSTUP:
!!	hodnota	= vypocteny odhad integralu
!!	errest	= odhad velikosti skutecne absoltni chyby
!!	flag	= informace o spolehlivosti vysledku
!!			  0.0 = vysledek pravdepodobne vyhovuje zadanym pozadavkum na chybu
!!			  XXX.YYY = cast zakladniho intervalu, ktera zustala nezpracovana,
!!			  kdyz se procedura priblizila k limitni hodnote pro "nofun"
!!			  Zadano v procedure jako INTEGER,PARAMETER :: nomax = 5000
!!	nofun	= pocet vypoctu Fun pri behu procedury
!!=============================================================================
!

implicit none

integer,parameter :: SP = selected_real_kind(4)
integer,parameter :: DP = selected_real_kind(8)

character(len=20),parameter,public :: Quanc8FMMVer = "1.0" ! 7.4.2008

private

public :: Quanc8

INTERFACE Quanc8
	MODULE PROCEDURE Quanc8Single,Quanc8Double
END INTERFACE

contains

SUBROUTINE Quanc8Single (Fun,a,b,abserr,relerr,hodnota,errest,nofun,flag)
! -----------------------------------------------------------------------------
! PURPOSE - Estimate the integral of FUN(X) from A to B to a user provided 
!  tolerance. This is an automatic adaptive routine based on the 8-panel 
!  Newton-Cotes rule. 

  	INTERFACE
		FUNCTION Fun(x) RESULT(f)  ! integrand function subprogram
			IMPLICIT NONE
			integer,parameter :: SP = selected_real_kind(4)
			REAL(SP),INTENT(IN)::x
			REAL(SP):: f
		END FUNCTION Fun
    END INTERFACE
   
   REAL(SP),INTENT(IN):: a        !  lower limit of integration
   REAL(SP),INTENT(IN):: b        !  upper limit of integration
   REAL(SP),INTENT(IN):: abserr   !  absolute error tolerance (>0)
   REAL(SP),INTENT(IN):: relerr   !  relative error tolerance (>0)

   REAL(SP),INTENT(OUT):: hodnota   !  approx. value of integral
   REAL(SP),INTENT(OUT):: errest   !  estimate of actual error
   INTEGER,INTENT(OUT):: nofun     !  no. of function evaluations
   REAL(SP),INTENT(OUT):: flag     !  reliability indicator. (=0 is O.K.)

   REAL(SP),PARAMETER:: W0=3956.0_SP/14175.0_SP
   REAL(SP),PARAMETER:: W1=23552.0_SP/14175.0_SP
   REAL(SP),PARAMETER:: W2=-3712.0_SP/14175.0_SP
   REAL(SP),PARAMETER:: W3=41984.0_SP/14175.0_SP
   REAL(SP),PARAMETER:: W4=-18160.0_SP/14175.0_SP

  REAL(SP):: area
  REAL(SP):: cor11  
  REAL(SP):: esterr
  REAL(SP), DIMENSION(0:16):: f,x
  REAL(SP), DIMENSION(8,30):: fsave,xsave
  INTEGER:: i,j
  INTEGER:: lev,nim
  INTEGER:: LEVMAX = 30
  INTEGER,PARAMETER :: LEVMIN = 1
  INTEGER,PARAMETER:: LEVOUT = 6
  INTEGER:: NOFIN     !...Trouble if NOFUN reaches NOFIN
  INTEGER,PARAMETER:: NOMAX = 5000
  REAL(SP):: qprev,qnow,qdiff,qleft
  REAL(SP):: stone,step
  REAL(SP):: tolerr
  REAL(SP),DIMENSION(31):: qright
  REAL(SP),PARAMETER:: ZERO=0.0, HALF=0.5, ONE=1.0, SIXTEEN=16.0
!------------------------------------------------------------------------------

!..... S t a g e   1  (General initialization)..........................
   NOFIN=NOMAX-8*(LEVMAX-LEVOUT+2**(LEVOUT+1))

!...initialize running sums to zero...
  area=ZERO
  cor11=ZERO
  flag=ZERO   ! flag is REAL
  hodnota=ZERO
  errest=ZERO
  nofun=0
  IF (a==b) RETURN
  
!..... S T A G E   2   (initialization for first interval)..............
  lev=0
  nim=1
  x(0)=a
  x(16)=b
  qprev=ZERO
  stone=(b-a)/SIXTEEN
  x(8)= HALF*(x(0)+x(16))
  x(4)= HALF*(x(0)+x(8))
  x(12)=HALF*(x(8)+x(16))
  x(2)= HALF*(x(0)+x(4))
  x(6)= HALF*(x(4)+x(8))
  x(10)=HALF*(x(8)+x(12))
  x(14)=HALF*(x(12)+x(16))

  DO j=0,16,2
    f(j)=Fun(x(j))
  END DO
  nofun=9

!..... S T A G E   3   (central calculation)............................
!...requires qprev,x0,x2,x4,....x(16),f0,f2,f4,....f16.
!...calculates x1,x3,....x15,f1,f3,....f15,qleft,qright,qnow,qdiff,area

N30: DO
	DO j=1,15,2  ! keeps coming back here until successful
      x(j)=HALF*(x(j-1)+x(j+1))
      f(j)=Fun(x(j))
    END DO
    nofun=nofun+8
    step=(x(16)-x(0))/SIXTEEN
    qleft=(w0*(f(0)+f(8))+w1*(f(1)+f(7))+w2*(f(2)+f(6))+   &
                           w3*(f(3)+f(5))+w4*f(4))*step
    qright(lev+1)=(w0*(f(8)+f(16))+w1*(f(9)+f(15))+w2*(f(10)+f(14))+   &
                    w3*(f(11)+f(13))+w4*f(12))*step
    qnow=qleft+qright(lev+1)
    qdiff=qnow-qprev
    area=area+qdiff

!..... S T A G E   4   (interval convergence test) .....................
    esterr=ABS(qdiff)/1023.0
    tolerr=MAX(abserr,relerr*ABS(area))*(step/stone)
    IF (lev < LEVMIN) GO TO 50
    
    
    IF (lev >= LEVMAX) THEN
      flag=flag+ONE
      goto 70
    END IF  
    IF (nofun > NOFIN) THEN      !..... S T A G E   6  (trouble section) 
      nofin=nofin+nofin  ! Number of function values is about to exceed limit
      levmax=levout
      flag=flag+(b-x(0))/(b-a)
      goto 70
    END IF
    IF (esterr <= TOLERR) goto 70

!..... S T A G E   5   (no convergence).................................
!...locate next interval...
50  continue
	nim=nim+nim
    lev=lev+1

    fsave(1:8,lev)=f(9:16)   ! store right hand elements for future use
    xsave(1:8,lev)=x(9:16)
!    DO i=1,8
!      fsave(i,lev)=f(i+8)
!      xsave(i,lev)=x(i+8)
!    END DO

    qprev=qleft
    DO I=1,8   ! assemble left hand elements for immediate use
      f(18-2*i)=f(9-i)
      x(18-2*i)=x(9-i)
    END DO
	goto 30

!..... S T A G E   7   (integral converged) ............................
!...add contributions into running sums...
70 continue
   hodnota=hodnota+qnow
   errest=errest+esterr
   cor11=cor11+qdiff/1023.0

  DO                                ! 72 loop in text 
    IF (nim == 2*(nim/2)) EXIT   
    nim=nim/2
    lev=lev-1
  END DO
  nim=nim+1
  
  IF (lev <= 0) GO TO 80   ! looks like success!
  qprev=qright(lev)   ! assemble elements required for next interval
  x(0)=x(16)
  f(0)=f(16)
  DO i=1,8
    f(2*i)=fsave(i,lev)
    x(2*i)=xsave(i,lev)
  END DO
 30 continue
 end do N30    ! puvodne goto 30

!..... S T A G E   8  (Finalize and return).............................
80 continue 
  hodnota=hodnota+cor11
  IF (errest /= ZERO) THEN
    DO  
      IF (ABS(hodnota)+errest > ABS(hodnota) ) EXIT
      errest=errest+errest  ! make sure ERREST is not < roundoff level ....
    END DO
  END IF

  RETURN
END Subroutine QuanC8Single   ! ---------------------------------------------------



SUBROUTINE Quanc8Double (Fun,a,b,abserr,relerr,hodnota,errest,nofun,flag)
! -----------------------------------------------------------------------------
! PURPOSE - Estimate the integral of FUN(X) from A to B to a user provided 
!  tolerance. This is an automatic adaptive routine based on the 8-panel 
!  Newton-Cotes rule. 

  	INTERFACE
		FUNCTION Fun(x) result(f)  ! integrand function subprogram
			IMPLICIT NONE
			integer,parameter :: DP = selected_real_kind(8)
			REAL(DP),INTENT(IN)::x
			REAL(DP):: f
		END FUNCTION Fun
	END INTERFACE
   
	REAL(DP),INTENT(IN):: a        !  lower limit of integration
	REAL(DP),INTENT(IN):: b        !  upper limit of integration
	REAL(DP),INTENT(IN):: abserr   !  absolute error tolerance (>0)
	REAL(DP),INTENT(IN):: relerr   !  relative error tolerance (>0)
    REAL(DP),INTENT(OUT):: hodnota  !  approx. value of integral
	REAL(DP),INTENT(OUT):: errest  !  estimate of actual error
	INTEGER,INTENT(OUT):: nofun     	!  no. of function evaluations
	REAL(DP),INTENT(OUT):: flag    ! reliability indicator. (=0 is O.K.)

	REAL(DP),PARAMETER:: W0=3956.0_DP/14175.0_DP
	REAL(DP),PARAMETER:: W1=23552.0_DP/14175.0_DP
	REAL(DP),PARAMETER:: W2=-3712.0_DP/14175.0_DP
	REAL(DP),PARAMETER:: W3=41984.0_DP/14175.0_DP
	REAL(DP),PARAMETER:: W4=-18160.0_DP/14175.0_DP

	REAL(DP):: area
	REAL(DP):: cor11  
	REAL(DP):: esterr
	REAL(DP), DIMENSION(0:16):: f,x
	REAL(DP), DIMENSION(8,30):: fsave,xsave
	INTEGER:: i,j
	INTEGER:: lev,nim
	INTEGER:: LEVMAX = 30
	INTEGER,PARAMETER :: LEVMIN = 1
	INTEGER,PARAMETER:: LEVOUT = 6
	INTEGER:: NOFIN     !...Trouble if NOFUN reaches NOFIN
	INTEGER,PARAMETER:: NOMAX = 5000
	REAL(DP):: qprev,qnow,qdiff,qleft
	REAL(DP):: stone,step
	REAL(DP):: tolerr
	REAL(DP),DIMENSION(31):: qright
	REAL(DP),PARAMETER:: ZERO=0.0, HALF=0.5, ONE=1.0, SIXTEEN=16.0
!------------------------------------------------------------------------------

!..... S t a g e   1  (General initialization)..........................
   NOFIN=NOMAX-8*(LEVMAX-LEVOUT+2**(LEVOUT+1))

!...initialize running sums to zero...
  area=ZERO
  cor11=ZERO
  flag=ZERO   ! flag is REAL
  hodnota=ZERO
  errest=ZERO
  nofun=0
  IF (a==b) RETURN
  
!..... S T A G E   2   (initialization for first interval)..............
  lev=0
  nim=1
  x(0)=a
  x(16)=b
  qprev=ZERO
  stone=(b-a)/SIXTEEN
  x(8)= HALF*(x(0)+x(16))
  x(4)= HALF*(x(0)+x(8))
  x(12)=HALF*(x(8)+x(16))
  x(2)= HALF*(x(0)+x(4))
  x(6)= HALF*(x(4)+x(8))
  x(10)=HALF*(x(8)+x(12))
  x(14)=HALF*(x(12)+x(16))

  DO j=0,16,2
    f(j)=Fun(x(j))
  END DO
  nofun=9

!..... S T A G E   3   (central calculation)............................
!...requires qprev,x0,x2,x4,....x(16),f0,f2,f4,....f16.
!...calculates x1,x3,....x15,f1,f3,....f15,qleft,qright,qnow,qdiff,area

30 Continue
N30: DO
	DO j=1,15,2  ! keeps coming back here until successful
      x(j)=HALF*(x(j-1)+x(j+1))
      f(j)=Fun(x(j))
    END DO
    nofun=nofun+8
    step=(x(16)-x(0))/SIXTEEN
    qleft=(w0*(f(0)+f(8))+w1*(f(1)+f(7))+w2*(f(2)+f(6))+   &
                           w3*(f(3)+f(5))+w4*f(4))*step
    qright(lev+1)=(w0*(f(8)+f(16))+w1*(f(9)+f(15))+w2*(f(10)+f(14))+   &
                    w3*(f(11)+f(13))+w4*f(12))*step
    qnow=qleft+qright(lev+1)
    qdiff=qnow-qprev
    area=area+qdiff

!..... S T A G E   4   (interval convergence test) .....................
    esterr=ABS(qdiff)/1023.0
    tolerr=MAX(abserr,relerr*ABS(area))*(step/stone)
    IF (lev < LEVMIN) GO TO 50
    
    
    IF (lev >= LEVMAX) THEN
      flag=flag+ONE
      GO TO 70
    END IF  
    IF (nofun > NOFIN) THEN      !..... S T A G E   6  (trouble section) 
      nofin=nofin+nofin  ! Number of function values is about to exceed limit
      levmax=levout
      flag=flag+(b-x(0))/(b-a)
      GO TO 70
    END IF
    IF (esterr <= TOLERR) GO TO 70

!..... S T A G E   5   (no convergence).................................
!...locate next interval...
50  Continue
	nim=nim+nim
    lev=lev+1

    fsave(1:8,lev)=f(9:16)   ! store right hand elements for future use
    xsave(1:8,lev)=x(9:16)
!    DO i=1,8
!      fsave(i,lev)=f(i+8)
!      xsave(i,lev)=x(i+8)
!    END DO

    qprev=qleft
    DO I=1,8   ! assemble left hand elements for immediate use
      f(18-2*i)=f(9-i)
      x(18-2*i)=x(9-i)
    END DO
    cycle N30


!..... S T A G E   7   (integral converged) ............................
!...add contributions into running sums...
70 Continue
   hodnota=hodnota+qnow
   errest=errest+esterr
   cor11=cor11+qdiff/1023.0

  DO                                ! 72 loop in text 
    IF (nim == 2*(nim/2)) EXIT   
    nim=nim/2
    lev=lev-1
  END DO
  nim=nim+1
  
  IF (lev <= 0) GO TO 80   ! looks like success!
  qprev=qright(lev)   ! assemble elements required for next interval
  x(0)=x(16)
  f(0)=f(16)
  DO i=1,8
    f(2*i)=fsave(i,lev)
    x(2*i)=xsave(i,lev)
  END DO
end do N30

!..... S T A G E   8  (Finalize and return).............................
80 Continue 
  hodnota=hodnota+cor11
  IF (errest /= ZERO) THEN
    DO  
      IF (ABS(hodnota)+errest > ABS(hodnota) ) EXIT
      errest=errest+errest  ! make sure ERREST is not < roundoff level ....
    END DO
  END IF

  RETURN
END Subroutine QuanC8Double  
! ---------------------------------------------------

end module Quanc8FMM
