c_simpson_dh_sample Function

function c_simpson_dh_sample(f,dh) result(sum)
  complex(8)           :: f(:)
  integer              :: n
  real(8)              :: dh
  complex(8)           :: sum,sum1,sum2,int1,int2
  integer              :: i,p,m
  complex(8),parameter :: zero=cmplx(0d0,0d0,8)
  N=size(f)
  sum=zero  
  if(N==1)return
  sum1=zero
  sum2=zero
  int1=zero
  int2=zero
  if(mod(n-1,2)==0)then                !if n-1 is even:
     do i=2,N-1,2
        sum1 = sum1 + f(i)
     enddo
     do i=3,N-2,2
        sum2 = sum2 + f(i)
     enddo
     sum = (f(1) + 4d0*sum1 + 2d0*sum2 + f(n))*dh/3d0
  else                        !if n-1 is odd, use Simpson's for N-3 slices + 3/8rule for the last
     if (N>=6) then
        do i=2,N-4,2
           sum1 = sum1 + f(i)
        enddo
        do i=3,N-5,2
           sum2 = sum2 + f(i)
        enddo
        int1 = (f(1) + 4d0*sum1 + 2d0*sum2 + f(n-3))*dh/3d0
     endif
     int2 = (f(n-3)+3d0*f(n-2)+3d0*f(n-1)+f(n))*dh*3d0/8d0
     sum  = int1 + int2
  end if
end function c_simpson_dh_sample