!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  PROGRAM TrapInt
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  (c) 2005, Christian B. Lang, Mathematische Methoden in der Physik
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Eine numerische Integration wird mit der erweiterten Trapezformel mit
!  sich verdoppelnder Stuetzstellen Anzahl durchgefuehrt. Das Ergebnis
!  der entsprechenden Simpson-Regel wird ebenfalls angeben.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

PROGRAM TrapInt

   implicit none

   real:: a,b,AlterWert,Integral
   integer:: MaximaleStufe,n

!  Die nachfolgenden Parameter koennten auch eingelesen und veraendert werden

   MaximaleStufe = 12
   a = 0
   b = 3.14159265358 

   write(*,*) 'Das Integral wird mit zunehmender Genauigkeit berechnet:'
   call Trapez(a, b, 1, Integral)
   write(*,'(a,i2,a,e15.8)') 'Stufe ',1,' Integral: ',Integral
   do n = 2,MaximaleStufe
     AlterWert = Integral
     call Trapez(a, b, n, Integral)
     write(*,'(a,i2,a,e15.8,a,e15.8)') 'Stufe ',n,' Integral: ',Integral,'  Simpsonwert: ',(4*Integral - AlterWert)/3
   end do

END PROGRAM TrapInt

REAL FUNCTION fun(x)

   implicit none

   real, intent(in)::x
!  Testfunktion (Polynom)
!   a=0, b=1               >> Integral=1/3=0.3333

!   fun=x*x

!  Testfunktion (Winkelfunktion)
!  a=0, b=3.14159265358   >> Integral=2}

   fun= SIN(x)

!  Testfunktion (Bogenlaenge Kreis)
!  a=0, b=1/Sqrt(2)=0.70711 >> Integral=pi/4=0.7854

!  fun= 1/Sqrt(1-x*x)

END FUNCTION fun

SUBROUTINE Trapez (a,b,nn,Integ)

!  Trapez-Integrationsformel
   implicit none

   real, intent(in)::a,b
   integer, intent(in)::nn
   real, intent(out)::Integ

   real:: del,x,s,fun
   integer:: np,j,n
 
   n=nn
   if(n<1)then
     write(*,*) 'Order too small, reset to n=2'
     n=2
   end if
   if(n==1)then
     Integ= (b-a)*(fun(a) + fun(b))/2
   else
     np=2**(n-2)
     del = (b-a)/np
     x = a + del/2
     s = 0
     do j = 1,np
        s = s + fun(x)
        x = x + del
     end do
!  Neuer Wert = (Alter Wert + Neue Beitraege) /2

   Integ = ((b-a)*s/np + Integ)/2
   end if
  
END SUBROUTINE Trapez
