c_simps_nonlin_func Function

function c_simps_nonlin_func(f,x) result(int)
  interface
     function f(x)
       real(8)    :: x
       complex(8) :: f
     end function f
  end interface
  real(8),dimension(:)             :: x
  real(8)                          :: dh
  integer                          :: L,M,i
  real(8),dimension(:),allocatable :: xx,wt,dx
  complex(8)                       :: int,int1,int2,int3
  L  = size(x)
  !
  int=0.d0
  allocate(dx(L+1))
  M=L-1
  dx=0d0
  forall(i=1:M)dx(i)=x(i+1)-x(i)
  int1=0d0;int2=0d0;int3=0d0
  i=0
  do while(i < L)
     if( (dx(i)==dx(i+1)) .AND. (dx(i)==dx(i+2)) )then !Simpson's 3/8 rule
        int1 = int1 + ( 3*dx(i)*( f(x(i)) + &
             3*( f(x(i+1)) + f(x(i+2)) ) + f(x(i+3)) ))/8
        i=i+3
     elseif(dx(i)==dx(i+1))then !Simpson's 1/3 rule
        int2 = int2+( 2*dx(i)*( f(x(i)) + &
             4*f(x(i+1)) + f(x(i+2))) )/6
        i=i+2
     elseif(dx(i)/=dx(i+1)) then !trapezoidal rule
        int3 = int3 + dx(i)*( f(x(i))+f(x(i+1)) )/2.d0
        i = i + 1
     endif
  enddo
  int = int1+int2+int3
end function c_simps_nonlin_func