complex function hintegrate(x0, x1, y0, y1, z0, z1, n, parx, pary, parz, m, t, p, f, gaussk, func, g, dgx, dgy, dgz)
    use functions
    real(8) :: x0, x1, y0, y1, z0, z1, n
    real(8) :: parx, pary, parz

    real(8) :: m, t, p, f, gaussk

    real(8) :: stepx, x, px, xc, dx
    real(8) :: stepy, y, py, yc, dy
    real(8) :: stepz, z, pz, zc, dz
    real(8) :: wx, wy, wz

    real(8) :: sinc, func, g, dgx, dgy, dgz

    complex I, fk

    stepx = (x1 - x0) / n
    stepy = (y1 - y0) / n
    stepz = (z1 - z0) / n

    px = x0
    x = px + stepx

    py = y0
    y = py + stepy

    pz = z0
    z = pz + stepz

    I = 0
    do while (x < x1)
        do while(y < y1)
            do while(z < z1)
                dx = x - px
                dy = y - py
                dz = z - pz

                xc = (px + x) / 2
                yc = (py + y) / 2
                zc = (pz + z) / 2

                wx = sinc(dgx(m, t, p, f, parx, pary, parz, xc, yc, zc)*dx)
                wy = sinc(dgy(m, t, p, f, parx, pary, parz, xc, yc, zc)*dy)
                wz = sinc(dgz(m, t, p, f, parx, pary, parz, xc, yc, zc)*dz)

                fk = func(xc, yc, zc, gaussk) * cdexp((CMPLX(0,1) / h) * g(m, t, p, f, parx, pary, parz, xc, yc, zc))

                I = I + dx * dy * dz * f * wx * wy * wz

                pz = z
                z = pz + stepz

                end do

            py = y
            y = py + stepy
            end do

        px = x
        x = px + stepx
        end do

    hintegrate = I
    return
end function hintegrate