Skip to content

Commit

Permalink
Fixed some indentation issues.
Browse files Browse the repository at this point in the history
  • Loading branch information
nakib committed Nov 13, 2024
1 parent 28f846f commit 0406afe
Show file tree
Hide file tree
Showing 2 changed files with 67 additions and 62 deletions.
79 changes: 42 additions & 37 deletions src/misc.f90
Original file line number Diff line number Diff line change
Expand Up @@ -1535,8 +1535,8 @@ subroutine interpolate(coarsemesh, refinement, f, q, interpolation)
end select

interpolation = aux
end subroutine interpolate

end subroutine interpolate
pure function Pade_coeffs(iomegas, us)
!! Evaluate eqs. A2 from the following article:
!! Solving the Eliashberg equations by means of N-point Pade' approximants
Expand Down Expand Up @@ -1641,11 +1641,11 @@ end subroutine invert_complex_square

subroutine subtitle(text)
!! Subroutine to print a subtitle.

character(len = *), intent(in) :: text
integer(i64) :: length
character(len = 75) :: string2print

length = len(text)

string2print = '___________________________________________________________________________'
Expand All @@ -1655,38 +1655,43 @@ subroutine subtitle(text)
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"

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 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) = -(fx(k + 2) - fx(k) + term2 + term3)/pi
end do
!! Does Hilbert tranform for a given function
!! 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(:)

! 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 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) = -(fx(k + 2) - fx(k) + term2 + term3)/pi
end do
end subroutine Hilbert_transform
end module misc
50 changes: 25 additions & 25 deletions test/test_misc.f90
Original file line number Diff line number Diff line change
Expand Up @@ -386,29 +386,29 @@ program test_misc

if(tests_all%get_status() .eqv. .false.) error stop -1

contains
! reference functions for the Hilbert transform test
! fx1 is an even function
pure elemental real(r64) function fx1(x)
real(r64), intent(in) :: x
fx1 = 1/(1.0_r64 + x**2)
end function fx1

! Hfx1 is actual hilbert transform of fx1, is an odd function
pure elemental real(r64) function hfx1(x)
real(r64), intent(in) :: x
hfx1 = x/(1.0_r64 + x**2)
end function hfx1

! fx2 is an odd function
pure elemental real(r64) function fx2(x)
real(r64), intent(in) :: x
fx2 = sin(x)/(1.0_r64 + x**2)
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 = (exp(-1.0_r64) - cos(x))/(1.0_r64 + x**2)
end function hfx2
contains
! Some reference functions and their Hilbert transforms:
pure elemental real(r64) function fx1(x)
real(r64), intent(in) :: x

fx1 = 1/(1.0_r64 + x**2)
end function fx1

pure elemental real(r64) function hfx1(x)
real(r64), intent(in) :: x

hfx1 = x/(1.0_r64 + x**2)
end function hfx1

pure elemental real(r64) function fx2(x)
real(r64), intent(in) :: x

fx2 = sin(x)/(1.0_r64 + x**2)
end function fx2

pure elemental real(r64) function hfx2(x)
real(r64), intent(in) :: x

hfx2 = (exp(-1.0_r64) - cos(x))/(1.0_r64 + x**2)
end function hfx2
end program test_misc

0 comments on commit 0406afe

Please sign in to comment.