From 4d3d4edfcc4b38e5c805b91275dbedd72d1fdbc9 Mon Sep 17 00:00:00 2001 From: Dwaipayan Paul Date: Tue, 12 Nov 2024 17:27:24 +0100 Subject: [PATCH 1/4] Moved Hilbert transform to misc, and wrote unit test (after revert #162) --- src/misc.f90 | 47 ++++++++++++++++++++++++- src/screening.f90 | 50 ++------------------------- test/test_misc.f90 | 86 ++++++++++++++++++++++++++++++++++++++++++++-- 3 files changed, 132 insertions(+), 51 deletions(-) diff --git a/src/misc.f90 b/src/misc.f90 index 30254ed..563a527 100644 --- a/src/misc.f90 +++ b/src/misc.f90 @@ -22,7 +22,7 @@ module misc #endif use precision, only: r128, r64, i64 - use params, only: kB, twopi + use params, only: kB, twopi, pi implicit none @@ -1653,4 +1653,49 @@ subroutine subtitle(text) string2print(75 - length + 1 : 75) = text if(this_image() == 1) write(*,'(A75)') string2print end subroutine subtitle + + subroutine Hilbert_transform(fx, Hfx) + !! Does Hilbert tranform for a given function + !! Ref - EQ (4)3, R. Balito et. al. + !! - An algorithm for fast Hilbert transform of real + !! functions + !! + !! fx is the function + !! Hfx is the Hilbert transform of the function + !! + + real(r64), allocatable, intent(out) :: Hfx(:) + real(r64), intent(in) :: fx(:) + + ! Local variables + integer :: n, k, nfx + real(r64) :: term2, term3, b + + nfx = size(fx) + allocate(Hfx(nfx)) + + ! Hilbert function is zero at the edges + Hfx(1) = 0.0_r64 + Hfx(nfx) = 0.0_r64 + + ! Both Hfx and fx follow 1-indexing system unlike the paper + do k = 1, nfx - 2 ! Run over the internal points + term2 = 0.0_r64 ! 2nd term in Bilato Eq. 4 + term3 = 0.0_r64 ! 3rd term in Bilato Eq. 4 + + do n = 1, nfx - 2 - k ! Partial sum over internal points + b = log((n + 1.0_r64)/n) + term2 = term2 - (1.0_r64 - (n + 1.0_r64)*b)*fx(k + n + 1) + & + (1.0_r64 - n*b)*fx(k + n + 2) + end do + + do n = 1, k - 1 ! Partial sum over internal points + b = log((n + 1.0_r64)/n) + term3 = term3 + (1.0_r64 - (n + 1.0_r64)*b)*fx(k - n + 1) - & + (1.0_r64 - n*b)*fx(k - n) + end do + + Hfx(k + 1) = (-1.0_r64/pi)*(fx(k + 2) - fx(k) + term2 + term3) + end do + end subroutine Hilbert_transform end module misc diff --git a/src/screening.f90 b/src/screening.f90 index 4734fcb..2caaa40 100644 --- a/src/screening.f90 +++ b/src/screening.f90 @@ -7,7 +7,8 @@ module screening_module use numerics_module, only: numerics use misc, only: linspace, mux_vector, binsearch, Fermi, print_message, & compsimps, twonorm, write2file_rank2_real, write2file_rank1_real, & - distribute_points, sort, qdist, operator(.umklapp.), Bose + distribute_points, sort, qdist, operator(.umklapp.), Bose, & + Hilbert_transform use wannier_module, only: wannier use delta, only: delta_fn, get_delta_fn_pointer @@ -58,53 +59,6 @@ subroutine calculate_qTF(crys, el) end if end subroutine calculate_qTF - subroutine Hilbert_transform(f, Hf) - !! Does Hilbert tranform of spectral head of bare polarizability - !! Ref - EQ (4), R. Balito et. al. - !! - !! Hf H.f(x), the Hilbert transform - !! f(x) the function - - real(r64), intent(in) :: f(:) - real(r64), allocatable, intent(out) :: Hf(:) - - ! Locals - real(r64) :: term2, term3, term1, b - integer(i64) :: nxs, k, n - - !Number points on domain grid - nxs = size(f) - - allocate(Hf(nxs)) - - !Assume that f vanishes at the edges, and Hf also - Hf(1) = 0.0_r64 - Hf(nxs) = 0.0_r64 - - do k = 2, nxs - 1 !Run over the internal points - term2 = 0.0_r64 !2nd term in Bilato Eq. 4 - term3 = 0.0_r64 !3rd term in Bilato Eq. 4 - - do n = 2, nxs - 1 - k !Partial sum over internal points - b = log((n + 1.0_r64)/n) - - term2 = term2 - (1.0_r64 - (n + 1.0_r64)*b)*f(k + n) + & - (1.0_r64 - n*b)*f(k + n + 1) - end do - - do n = 2, k - 2 !Partial sum over internal points - b = log((n + 1.0_r64)/n) - - term3 = term3 + (1.0_r64 - (n + 1.0_r64)*b)*f(k - n) - & - (1.0_r64 - n*b)*f(k - n - 1) - end do - - term1 = f(k + 1) - f(k - 1) - - Hf(k) = (-1.0_r64/pi)*(term1 + term2 + term3) - end do - end subroutine Hilbert_transform - !!$ subroutine head_polarizability_real_3d_T(Reeps_T, Omegas, spec_eps_T, Hilbert_weights_T) !!$ !! Head of the bare real polarizability of the 3d Kohn-Sham system using !!$ !! Hilbert transform for a given set of temperature-dependent quantities. diff --git a/test/test_misc.f90 b/test/test_misc.f90 index 5cab0e4..e77acfb 100644 --- a/test/test_misc.f90 +++ b/test/test_misc.f90 @@ -7,12 +7,12 @@ program test_misc twonorm, binsearch, mux_vector, demux_vector, interpolate, coarse_grained, & unique, linspace, compsimps, mux_state, demux_state, demux_mesh, expm1, & Fermi, Bose, Pade_continued, precompute_interpolation_corners_and_weights, & - interpolate_using_precomputed, operator(.umklapp.), shrink + interpolate_using_precomputed, operator(.umklapp.), shrink, Hilbert_transform implicit none integer :: itest - integer, parameter :: num_tests = 28 + integer, parameter :: num_tests = 32 type(testify) :: test_array(num_tests), tests_all integer(i64) :: index, quotient, remainder, int_array(5), v1(3), v2(3), & v1_muxed, v2_muxed, ik, ik1, ik2, ik3, ib1, ib2, ib3, wvmesh(3), & @@ -23,6 +23,9 @@ program test_misc real_array(5), result, q1(3, 4), q2(3, 4), q3(3, 4) real(r64), allocatable :: integrand(:), domain(:), im_axis(:), real_func(:), & widc(:, :), f_coarse(:), f_interp(:), array_of_reals(:) + real(r64), allocatable :: hfx1_even(:), hfx1_odd(:), hfx2_even(:), hfx2_odd(:), & + ind_even(:), ind_odd(:), x_even(:), x_odd(:), xmin, xmax + integer(i64) :: n_even, n_odd print*, '<>' @@ -334,8 +337,87 @@ program test_misc call shrink(array_of_reals, 2_i64) call test_array(itest)%assert(array_of_reals, [1, 2]*1.0_r64) + ! Hilbert transform tests (H) + ! fx1 -> function 1, fx2 -> function 2 + ! - hfx1_even stores hilbert transform calculated for fx1, and for even number + ! - of points + xmin = -30.0 + xmax = 30.0 + n_even = 4000 + n_odd = 4001 + ! ind_even are indices to compare in case of even number of points + ! ind_odd are indices to compare in case of odd number of points + allocate(ind_even(6),ind_odd(5)) + ind_even = [801, 1201, 1601, 2001, 2401, 2801] + ind_odd = [889, 1333, 1777, 2221, 2665] + + itest = itest + 1 + test_array(itest) = testify("Hilbert transform: even function, even points") + allocate(x_even(n_even), hfx1_even(n_even)) + call linspace(x_even, xmin, xmax, n_even) + call Hilbert_transform(fx1_array(x_even), hfx1_even) + call test_array(itest)%assert(hfx1_even(ind_even), hfx1_array(x_even(ind_even)), & + tol = 2e-4_r64) + + itest = itest + 1 + test_array(itest) = testify("Hilbert transform: odd function, even points") + allocate(hfx2_even(n_even)) + call Hilbert_transform(fx2_array(x_even), hfx2_even) + call test_array(itest)%assert(hfx2_even(ind_even), hfx2_array(x_even(ind_even)), & + tol = 1e-4_r64) + + itest = itest + 1 + test_array(itest) = testify("Hilbert transform: even function, odd points") + allocate(x_odd(n_odd), hfx1_odd(n_odd)) + call linspace(x_odd, xmin, xmax, n_odd) + call Hilbert_transform(fx1_array(x_odd), hfx1_odd) + call test_array(itest)%assert(hfx1_odd(ind_odd), hfx1_array(x_odd(ind_odd)), & + tol = 4e-4_r64) + + itest = itest + 1 + test_array(itest) = testify("Hilbert transform: odd function, odd points") + allocate(hfx2_odd(n_odd)) + call Hilbert_transform(fx2_array(x_odd), hfx2_odd) + call test_array(itest)%assert(hfx2_odd(ind_odd), hfx2_array(x_odd(ind_odd)), & + tol = 1e-5_r64) + tests_all = testify(test_array) call tests_all%report if(tests_all%get_status() .eqv. .false.) error stop -1 + + contains + ! reference functions for the Hilbert transform test + ! fx1 is an even function + function fx1_array(x) result(fx1) + real(r64), intent(in) :: x(:) + real(r64), allocatable :: fx1(:) + allocate(fx1(size(x))) + fx1 = 1/(1 + x**2) + end function fx1_array + + ! Hfx1 is actual hilbert transform of fx1, is an odd function + function hfx1_array(x) result(hfx1) + real(r64), intent(in) :: x(:) + real(r64), allocatable :: hfx1(:) + allocate(hfx1(size(x))) + hfx1 = x/(1 + x**2) + end function hfx1_array + + ! fx2 is an odd function + function fx2_array(x) result(fx2) + real(r64), intent(in) :: x(:) + real(r64), allocatable :: fx2(:) + allocate(fx2(size(x))) + fx2 = sin(x)/(1 + x**2) + end function fx2_array + + ! Hfx2 is actual hilbert transform of fx2, is an even function + function hfx2_array(x) result(hfx2) + real(r64), intent(in) :: x(:) + real(r64), allocatable :: hfx2(:) + real(r64), parameter :: e = 2.718281 + allocate(hfx2(size(x))) + hfx2 = (1/e - cos(x))/(1 + x**2) + end function hfx2_array end program test_misc From eda1e54d5cde18d7e1f59309add982434f7c3ce3 Mon Sep 17 00:00:00 2001 From: Dwaipayan Paul Date: Wed, 13 Nov 2024 15:01:43 +0100 Subject: [PATCH 2/4] Resolved requested changes --- src/misc.f90 | 16 ++++------ test/test_misc.f90 | 73 +++++++++++++++++++++------------------------- 2 files changed, 40 insertions(+), 49 deletions(-) diff --git a/src/misc.f90 b/src/misc.f90 index 563a527..d6cab81 100644 --- a/src/misc.f90 +++ b/src/misc.f90 @@ -1657,15 +1657,12 @@ end subroutine subtitle subroutine Hilbert_transform(fx, Hfx) !! Does Hilbert tranform for a given function !! Ref - EQ (4)3, R. Balito et. al. - !! - An algorithm for fast Hilbert transform of real - !! functions - !! - !! fx is the function - !! Hfx is the Hilbert transform of the function + !! "An algorithm for fast Hilbert transform of real + !! functions" !! - real(r64), allocatable, intent(out) :: Hfx(:) real(r64), intent(in) :: fx(:) + real(r64), allocatable, intent(out) :: Hfx(:) ! Local variables integer :: n, k, nfx @@ -1677,8 +1674,8 @@ subroutine Hilbert_transform(fx, Hfx) ! Hilbert function is zero at the edges Hfx(1) = 0.0_r64 Hfx(nfx) = 0.0_r64 - - ! Both Hfx and fx follow 1-indexing system unlike the paper + ! Note: In the reference, we have N + 1 points and indexing is 0-based + ! whereas here we have N points and indexing is 1-based do k = 1, nfx - 2 ! Run over the internal points term2 = 0.0_r64 ! 2nd term in Bilato Eq. 4 term3 = 0.0_r64 ! 3rd term in Bilato Eq. 4 @@ -1688,14 +1685,13 @@ subroutine Hilbert_transform(fx, Hfx) term2 = term2 - (1.0_r64 - (n + 1.0_r64)*b)*fx(k + n + 1) + & (1.0_r64 - n*b)*fx(k + n + 2) end do - do n = 1, k - 1 ! Partial sum over internal points b = log((n + 1.0_r64)/n) term3 = term3 + (1.0_r64 - (n + 1.0_r64)*b)*fx(k - n + 1) - & (1.0_r64 - n*b)*fx(k - n) end do - Hfx(k + 1) = (-1.0_r64/pi)*(fx(k + 2) - fx(k) + term2 + term3) + Hfx(k + 1) = -(fx(k + 2) - fx(k) + term2 + term3)/pi end do end subroutine Hilbert_transform end module misc diff --git a/test/test_misc.f90 b/test/test_misc.f90 index e77acfb..788c168 100644 --- a/test/test_misc.f90 +++ b/test/test_misc.f90 @@ -339,8 +339,8 @@ program test_misc ! Hilbert transform tests (H) ! fx1 -> function 1, fx2 -> function 2 - ! - hfx1_even stores hilbert transform calculated for fx1, and for even number - ! - of points + ! hfx1_even stores hilbert transform calculated for fx1, and for even number + ! of points xmin = -30.0 xmax = 30.0 n_even = 4000 @@ -352,33 +352,33 @@ program test_misc ind_odd = [889, 1333, 1777, 2221, 2665] itest = itest + 1 - test_array(itest) = testify("Hilbert transform: even function, even points") + test_array(itest) = testify("Hilbert transform: f(x) = 1/(1 + x^2), even points") allocate(x_even(n_even), hfx1_even(n_even)) call linspace(x_even, xmin, xmax, n_even) - call Hilbert_transform(fx1_array(x_even), hfx1_even) - call test_array(itest)%assert(hfx1_even(ind_even), hfx1_array(x_even(ind_even)), & + call Hilbert_transform(fx1_el(x_even), hfx1_even) + call test_array(itest)%assert(hfx1_even(ind_even), hfx1_el(x_even(ind_even)), & tol = 2e-4_r64) itest = itest + 1 - test_array(itest) = testify("Hilbert transform: odd function, even points") + test_array(itest) = testify("Hilbert transform: f(x) = sin(x)/(1 + x^2), even points") allocate(hfx2_even(n_even)) - call Hilbert_transform(fx2_array(x_even), hfx2_even) - call test_array(itest)%assert(hfx2_even(ind_even), hfx2_array(x_even(ind_even)), & + call Hilbert_transform(fx2_el(x_even), hfx2_even) + call test_array(itest)%assert(hfx2_even(ind_even), hfx2_el(x_even(ind_even)), & tol = 1e-4_r64) itest = itest + 1 - test_array(itest) = testify("Hilbert transform: even function, odd points") + test_array(itest) = testify("Hilbert transform: f(x) = 1/(1 + x^2), odd points") allocate(x_odd(n_odd), hfx1_odd(n_odd)) call linspace(x_odd, xmin, xmax, n_odd) - call Hilbert_transform(fx1_array(x_odd), hfx1_odd) - call test_array(itest)%assert(hfx1_odd(ind_odd), hfx1_array(x_odd(ind_odd)), & + call Hilbert_transform(fx1_el(x_odd), hfx1_odd) + call test_array(itest)%assert(hfx1_odd(ind_odd), hfx1_el(x_odd(ind_odd)), & tol = 4e-4_r64) itest = itest + 1 - test_array(itest) = testify("Hilbert transform: odd function, odd points") + test_array(itest) = testify("Hilbert transform: f(x) = sin(x)/(1 + x^2), odd points") allocate(hfx2_odd(n_odd)) - call Hilbert_transform(fx2_array(x_odd), hfx2_odd) - call test_array(itest)%assert(hfx2_odd(ind_odd), hfx2_array(x_odd(ind_odd)), & + call Hilbert_transform(fx2_el(x_odd), hfx2_odd) + call test_array(itest)%assert(hfx2_odd(ind_odd), hfx2_el(x_odd(ind_odd)), & tol = 1e-5_r64) tests_all = testify(test_array) @@ -389,35 +389,30 @@ program test_misc contains ! reference functions for the Hilbert transform test ! fx1 is an even function - function fx1_array(x) result(fx1) - real(r64), intent(in) :: x(:) - real(r64), allocatable :: fx1(:) - allocate(fx1(size(x))) - fx1 = 1/(1 + x**2) - end function fx1_array + elemental function fx1_el(x) result(fx1) + real(r64), intent(in) :: x + real(r64) :: fx1 + fx1 = 1/(1.0_r64 + x**2) + end function fx1_el ! Hfx1 is actual hilbert transform of fx1, is an odd function - function hfx1_array(x) result(hfx1) - real(r64), intent(in) :: x(:) - real(r64), allocatable :: hfx1(:) - allocate(hfx1(size(x))) - hfx1 = x/(1 + x**2) - end function hfx1_array + elemental function hfx1_el(x) result(hfx1) + real(r64), intent(in) :: x + real(r64) :: hfx1 + hfx1 = x/(1.0_r64 + x**2) + end function hfx1_el ! fx2 is an odd function - function fx2_array(x) result(fx2) - real(r64), intent(in) :: x(:) - real(r64), allocatable :: fx2(:) - allocate(fx2(size(x))) - fx2 = sin(x)/(1 + x**2) - end function fx2_array + elemental function fx2_el(x) result(fx2) + real(r64), intent(in) :: x + real(r64) :: fx2 + fx2 = sin(x)/(1.0_r64 + x**2) + end function fx2_el ! Hfx2 is actual hilbert transform of fx2, is an even function - function hfx2_array(x) result(hfx2) - real(r64), intent(in) :: x(:) - real(r64), allocatable :: hfx2(:) - real(r64), parameter :: e = 2.718281 - allocate(hfx2(size(x))) - hfx2 = (1/e - cos(x))/(1 + x**2) - end function hfx2_array + elemental function hfx2_el(x) result(hfx2) + real(r64), intent(in) :: x + real(r64) :: hfx2 + hfx2 = (1/exp(1.0_r64) - cos(x))/(1.0_r64 + x**2) + end function hfx2_el end program test_misc From 09009f7e193fc08d9ca9579f60456f207c95b8b2 Mon Sep 17 00:00:00 2001 From: Dwaipayan Paul Date: Wed, 13 Nov 2024 15:50:26 +0100 Subject: [PATCH 3/4] Resolved second round of requested changes --- src/misc.f90 | 1 - test/test_misc.f90 | 36 ++++++++++++++++-------------------- 2 files changed, 16 insertions(+), 21 deletions(-) diff --git a/src/misc.f90 b/src/misc.f90 index d6cab81..8957ba0 100644 --- a/src/misc.f90 +++ b/src/misc.f90 @@ -1659,7 +1659,6 @@ subroutine Hilbert_transform(fx, Hfx) !! Ref - EQ (4)3, R. Balito et. al. !! "An algorithm for fast Hilbert transform of real !! functions" - !! real(r64), intent(in) :: fx(:) real(r64), allocatable, intent(out) :: Hfx(:) diff --git a/test/test_misc.f90 b/test/test_misc.f90 index 788c168..cd389c1 100644 --- a/test/test_misc.f90 +++ b/test/test_misc.f90 @@ -355,30 +355,30 @@ program test_misc test_array(itest) = testify("Hilbert transform: f(x) = 1/(1 + x^2), even points") allocate(x_even(n_even), hfx1_even(n_even)) call linspace(x_even, xmin, xmax, n_even) - call Hilbert_transform(fx1_el(x_even), hfx1_even) - call test_array(itest)%assert(hfx1_even(ind_even), hfx1_el(x_even(ind_even)), & + call Hilbert_transform(fx1(x_even), hfx1_even) + call test_array(itest)%assert(hfx1_even(ind_even), hfx1(x_even(ind_even)), & tol = 2e-4_r64) itest = itest + 1 test_array(itest) = testify("Hilbert transform: f(x) = sin(x)/(1 + x^2), even points") allocate(hfx2_even(n_even)) - call Hilbert_transform(fx2_el(x_even), hfx2_even) - call test_array(itest)%assert(hfx2_even(ind_even), hfx2_el(x_even(ind_even)), & + call Hilbert_transform(fx2(x_even), hfx2_even) + call test_array(itest)%assert(hfx2_even(ind_even), hfx2(x_even(ind_even)), & tol = 1e-4_r64) itest = itest + 1 test_array(itest) = testify("Hilbert transform: f(x) = 1/(1 + x^2), odd points") allocate(x_odd(n_odd), hfx1_odd(n_odd)) call linspace(x_odd, xmin, xmax, n_odd) - call Hilbert_transform(fx1_el(x_odd), hfx1_odd) - call test_array(itest)%assert(hfx1_odd(ind_odd), hfx1_el(x_odd(ind_odd)), & + call Hilbert_transform(fx1(x_odd), hfx1_odd) + call test_array(itest)%assert(hfx1_odd(ind_odd), hfx1(x_odd(ind_odd)), & tol = 4e-4_r64) itest = itest + 1 test_array(itest) = testify("Hilbert transform: f(x) = sin(x)/(1 + x^2), odd points") allocate(hfx2_odd(n_odd)) - call Hilbert_transform(fx2_el(x_odd), hfx2_odd) - call test_array(itest)%assert(hfx2_odd(ind_odd), hfx2_el(x_odd(ind_odd)), & + call Hilbert_transform(fx2(x_odd), hfx2_odd) + call test_array(itest)%assert(hfx2_odd(ind_odd), hfx2(x_odd(ind_odd)), & tol = 1e-5_r64) tests_all = testify(test_array) @@ -389,30 +389,26 @@ program test_misc contains ! reference functions for the Hilbert transform test ! fx1 is an even function - elemental function fx1_el(x) result(fx1) + pure elemental real(r64) function fx1(x) real(r64), intent(in) :: x - real(r64) :: fx1 fx1 = 1/(1.0_r64 + x**2) - end function fx1_el + end function fx1 ! Hfx1 is actual hilbert transform of fx1, is an odd function - elemental function hfx1_el(x) result(hfx1) + pure elemental real(r64) function hfx1(x) real(r64), intent(in) :: x - real(r64) :: hfx1 hfx1 = x/(1.0_r64 + x**2) - end function hfx1_el + end function hfx1 ! fx2 is an odd function - elemental function fx2_el(x) result(fx2) + pure elemental real(r64) function fx2(x) real(r64), intent(in) :: x - real(r64) :: fx2 fx2 = sin(x)/(1.0_r64 + x**2) - end function fx2_el + end function fx2 ! Hfx2 is actual hilbert transform of fx2, is an even function - elemental function hfx2_el(x) result(hfx2) + pure elemental real(r64) function hfx2(x) real(r64), intent(in) :: x - real(r64) :: hfx2 hfx2 = (1/exp(1.0_r64) - cos(x))/(1.0_r64 + x**2) - end function hfx2_el + end function hfx2 end program test_misc From b9caf4091dd375cca32651a01ae0e877fc92825f Mon Sep 17 00:00:00 2001 From: Dwaipayan Paul Date: Wed, 13 Nov 2024 16:12:37 +0100 Subject: [PATCH 4/4] Resolved third round of requested changes --- src/misc.f90 | 16 ++++++---------- test/test_misc.f90 | 2 +- 2 files changed, 7 insertions(+), 11 deletions(-) diff --git a/src/misc.f90 b/src/misc.f90 index 8957ba0..efc6fd8 100644 --- a/src/misc.f90 +++ b/src/misc.f90 @@ -1657,39 +1657,35 @@ end subroutine subtitle subroutine Hilbert_transform(fx, Hfx) !! Does Hilbert tranform for a given function !! Ref - EQ (4)3, R. Balito et. al. - !! "An algorithm for fast Hilbert transform of real - !! functions" + !! "An algorithm for fast Hilbert transform of real functions" real(r64), intent(in) :: fx(:) real(r64), allocatable, intent(out) :: Hfx(:) - ! Local variables integer :: n, k, nfx real(r64) :: term2, term3, b nfx = size(fx) allocate(Hfx(nfx)) - ! Hilbert function is zero at the edges Hfx(1) = 0.0_r64 Hfx(nfx) = 0.0_r64 ! Note: In the reference, we have N + 1 points and indexing is 0-based ! whereas here we have N points and indexing is 1-based - do k = 1, nfx - 2 ! Run over the internal points - term2 = 0.0_r64 ! 2nd term in Bilato Eq. 4 - term3 = 0.0_r64 ! 3rd term in Bilato Eq. 4 + do k = 1, nfx - 2 ! Run over the internal points + term2 = 0.0_r64 ! 2nd term in Bilato Eq. 4 + term3 = 0.0_r64 ! 3rd term in Bilato Eq. 4 - do n = 1, nfx - 2 - k ! Partial sum over internal points + do n = 1, nfx - 2 - k ! Partial sum over internal points b = log((n + 1.0_r64)/n) term2 = term2 - (1.0_r64 - (n + 1.0_r64)*b)*fx(k + n + 1) + & (1.0_r64 - n*b)*fx(k + n + 2) end do - do n = 1, k - 1 ! Partial sum over internal points + do n = 1, k - 1 ! Partial sum over internal points b = log((n + 1.0_r64)/n) term3 = term3 + (1.0_r64 - (n + 1.0_r64)*b)*fx(k - n + 1) - & (1.0_r64 - n*b)*fx(k - n) end do - Hfx(k + 1) = -(fx(k + 2) - fx(k) + term2 + term3)/pi end do end subroutine Hilbert_transform diff --git a/test/test_misc.f90 b/test/test_misc.f90 index cd389c1..ac254e0 100644 --- a/test/test_misc.f90 +++ b/test/test_misc.f90 @@ -409,6 +409,6 @@ end function fx2 ! Hfx2 is actual hilbert transform of fx2, is an even function pure elemental real(r64) function hfx2(x) real(r64), intent(in) :: x - hfx2 = (1/exp(1.0_r64) - cos(x))/(1.0_r64 + x**2) + hfx2 = (exp(-1.0_r64) - cos(x))/(1.0_r64 + x**2) end function hfx2 end program test_misc