diff --git a/src/misc.f90 b/src/misc.f90
index efc6fd8..38b17e7 100644
--- a/src/misc.f90
+++ b/src/misc.f90
@@ -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
@@ -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 = '___________________________________________________________________________'
@@ -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
diff --git a/test/test_misc.f90 b/test/test_misc.f90
index ac254e0..f886fa5 100644
--- a/test/test_misc.f90
+++ b/test/test_misc.f90
@@ -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