!integration matrix test
! this is a routine for integration functions with an example at the end that integrates a matrix whose elements are functions
!It's written in fortran 2003

module intmat
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none

  ! Wrapper for a procedure pointer so we can make an array of them
  type ptr_wrapper
     procedure(f), nopass, pointer :: func
  end type ptr_wrapper


  ! Use an interface for the function rather than "extern" functions
  ! Change the interface to whatever you want. You could take multiple reals, integers, etc
  abstract interface
    pure function f(x1)
       import
       real(dp), intent(in) :: x1
       real(dp) :: f
     end function f
  end interface

! The intent(in attribute of argument x1 means that x1 cannot be changed inside the function. Note that the return type of func needs to be declared. If this is omitted, somme compilers will not compile

!When declaring variables inside functions and subroutines that need to be passed in or out, intent may be added to the declaration.

!intent(in) means that the variable value can enter, but not be changed

!intent(out) means the variable is set inside the procedure and sent back to the main program with any initial values ignored.

!intent(inout) means that the variable comes in with a value and leaves with a value (default).



contains

   subroutine DGS0(x1,x2,fx,S)
! Integra por gausseanos. LLama a la función DGS8.
! x1 y x2 son los límites, fx es la función y S el resultado
    real(dp), intent(in) :: x1
    real(dp), intent(in) :: x2
    procedure(f), pointer :: fx
    real(dp):: u, v, c, sl, sg, sp, sa, sf
	real(dp), intent(out) :: S
    real(dp) :: eps = 1.0E-6

	S = 0d0
	U = x1
1	V = x2
	IF ( U .LT. x2 ) THEN
		SF = DGS8 (u,v,fx)
2		C  = (U+V)/2d0
        SL = DGS8 (u,c,fx)
        SG = DGS8 (c,v,fx)
        SP = SL+SG
        SA = DABS(SP-SF)/EPS
        IF ( SA.GE.1.D0 ) THEN
  	      V  = C
          SF = SL
        GOTO 2
        END IF
  	      U = V
  	      S = S+SP
  	      GOTO 1
    END IF
end subroutine DGS0

function DGS8(x1,x2,fx) result (d)
! rutina que integra. Es llamada por DGS0. 
! x1 y x2 son los límites, fx la función, d el resultado
	implicit none
	real(dp), intent(in) :: x1
	real(dp), intent(in) :: x2
	procedure(f), pointer :: fx
!	real (dp), intent(out) :: d
	real(dp) :: d
	real(dp) :: c,h,x,s
	c = (x1+x2)/2
	h = (x2-x1)/2
    X = .960289856497536D0*H
    S = .101228536290376D0*(fx(C+X)+fx(C-X))
    X = .796666477413627D0*H
    S = S + .222381034453374D0*(fx(C+X)+fx(C-X))
    X = .525532409916329D0*H
    S = S + .313706645877887D0*(fx(C+X)+fx(C-X))
    X = .183434642495650D0*H
    S = S + .362683783378362D0*(fx(C+X)+fx(C-X))
    d = S * H
end function DGS8


  pure function test_func(x1)
    real(dp), intent(in) :: x1
    real(dp) test_func

    test_func = 2d0 * x1
  end function test_func

  pure function test_func2(x1)
    real(dp), intent(in) :: x1
    real(dp) test_func2

    test_func2 = x1**2
  end function test_func2


end module intmat

program integral_matrix
  use, intrinsic :: iso_fortran_env, only: dp => real64
  use intmat
  implicit none

  integer i, j
  real(dp) :: dx, S
!  real(dp) :: pi2 = acos(0d0)

  type(ptr_wrapper) :: w(2,2)

  !Point to the function that you want each index to point to
  w(1,1)%func => test_func
  w(1,2)%func => test_func2
  ! Pointing to the generic function wasn't working so I had to point to the specific double value trig functions
  w(2,1)%func => dcos 
  w(2,2)%func => dsin  

  do i = 1, 2
     do j = 1, 2
        call DGS0(0d0,1.5d0, w(i,j)%func,S)
		dx = S
        print *, dx
     end do
  end do

print *, 'Done :)'

end program integral_matrix

