Hilbert transform #157

dpaulzc · 2024-10-25T10:40:18Z

Implemented Hilbert transform using the numerical method from R. Balito et. al. Eq(4).

nakib

Dear Dwaipayan,

Thanks for the pull request. I left some comments. When possible, please take a look.

Best,
Nakib

nakib · 2024-10-25T10:47:56Z

src/screening.f90

@@ -127,6 +127,65 @@ subroutine calculate_Hilbert_weights(w_disc, w_cont, zeroplus, Hilbert_weights)
    end do
  end subroutine calculate_Hilbert_weights

+   subroutine hilbert_transform(Reeps_T, Omegas, spec_eps_T)
+    !! EQ (4), R. Balito et. al. 


Add meaningful comment.

nakib · 2024-10-25T10:48:09Z

src/screening.f90

@@ -127,6 +127,65 @@ subroutine calculate_Hilbert_weights(w_disc, w_cont, zeroplus, Hilbert_weights)
    end do
  end subroutine calculate_Hilbert_weights

+   subroutine hilbert_transform(Reeps_T, Omegas, spec_eps_T)


Write Hilbert instead of hilbert.

nakib · 2024-10-25T10:49:21Z

src/screening.f90

+    real(r64) :: in1, in2, fin, bl !, domegas
+    integer(i64) :: nx, kh, np


Use meaningful names for variables unless you are staying consistent with the notation used in the reference paper's equation(s).

nakib · 2024-10-25T10:50:10Z

src/screening.f90

+        in1=0.0_r64
+        in2=0.0_r64


Use space on both sides of =, here and elsewhere.

nakib · 2024-10-25T10:51:06Z

src/screening.f90

+        in1=0.0_r64
+        in2=0.0_r64
+        do np=1,nx-1-kh
+            bl = log(real(np+1.0_r64)/real(np))


Add space around =, +, and -, but not around * and /.

nakib · 2024-10-25T10:51:30Z

src/screening.f90

+        in2=0.0_r64
+        do np=1,nx-1-kh
+            bl = log(real(np+1.0_r64)/real(np))
+            if ( (kh+np).lt.nx) then


Use < instead of .lt..

nakib · 2024-10-25T10:51:55Z

src/screening.f90

+
+        do np=1,kh-1
+           bl = log(real(np+1.0_r64)/real(np))
+           if ((kh-np).gt.1) then


Use > instead of .gt..

nakib · 2024-10-25T10:53:46Z

src/screening.f90

+        if (kh.ge.2 .and. kh.le.(nx-1)) then
+           fin = spec_eps_T(kh+1)-spec_eps_T(kh-1)
+         else
+            fin = 1e-12_r64
+         end if


Use proper indentation. For example, the else and end if should start at the same level of indentation as if. In a decent editor (like emacs), selecting a region and pressing TAB fixes all indentation issues. Perhaps your editor is also decent like that...

nakib · 2024-10-25T10:54:05Z

src/screening.f90

+            fin = 1e-12_r64
+         end if
+        Reeps_T(kh) = (-1.0_r64/pi)*(fin+in1+in2)
+


Remove this empty space.

nakib · 2024-10-25T10:54:16Z

src/screening.f90

+        Reeps_T(kh) = (-1.0_r64/pi)*(fin+in1+in2)
+
+    end do
+


Remove this empty space.

nakib · 2024-10-25T11:00:33Z

src/screening.f90

+        in1=0.0_r64
+        in2=0.0_r64
+        do np=1,nx-1-kh
+            bl = log(real(np+1.0_r64)/real(np))


If you write np + 1.0_r64, do you have to wrap further with real? What does real do? Does it respect the precision you are seeking? Also, does the np in the denominator have to be wrapped by real since writing the numerator as (np + 1.0_r64) already makes it a r64 type number.

nakib · 2024-10-25T11:02:36Z

src/screening.f90

+        in2=0.0_r64
+        do np=1,nx-1-kh
+            bl = log(real(np+1.0_r64)/real(np))
+            if ( (kh+np).lt.nx) then


It would be nice to add comments explaining the reasoning for these conditionals. You should do this unless it is absolutely obvious from the context.

nakib · 2024-10-25T11:02:54Z

src/screening.f90

+
+    end do
+
+  end subroutine hilbert_transform 


hilbert -> Hilbert

nakib · 2024-10-25T11:23:01Z

src/screening.f90

+            end if
+        end do
+
+        if (kh.ge.2 .and. kh.le.(nx-1)) then


Write >= and <= instead of .ge. and .le., respectively.

nakib · 2024-10-25T13:59:57Z

src/screening.f90

+        in2=0.0_r64
+        do np=1,nx-1-kh
+            bl = log(real(np+1.0_r64)/real(np))
+            if ( (kh+np).lt.nx) then


Write if((kh + np) < nx) instead of if ( (kh+np).lt.nx). Fix similar issues elsewhere. There should not be randomly inserted spaces. We should try to stick to a consistent style.

nakib · 2024-10-25T14:02:29Z

src/screening.f90

+            else 
+              !print *,"Hilbert transform: Boundary breached",kh+np,nx
+              in1 = in1 - (1.0_r64-(np+1.0_r64)*bl)*spec_eps_T(kh+np) + &
+                      (1.0_r64-np*bl)*1e-12_r64


What is this magic 1e-12_r64? If it is supposed to be a small number, define above as a parameter along with a reference to the equation where it appears. The same number appears below again.

nakib

Thanks for the changes. A few more minor comments.

Best,
Nakib

nakib · 2024-10-25T14:36:09Z

src/screening.f90

+        end do
+
+        ! To avoid edge error
+        if (kh >= 2 .and. kh <= (nOmegas-1)) then


Space missing around -. But I'll fix it after merging.

nakib · 2024-10-25T14:36:50Z

src/screening.f90

+        do nh = 1, kh - 1
+           bh = log((nh + 1.0_r64)/nh)
+           ! To avoid edge error in the lower limit
+           if ((kh - nh) > 1) then


Why is there a space after if? I don't think this done anywhere else in this code base.

nakib · 2024-10-25T14:38:25Z

src/screening.f90

+        ! For outside the boundary
+            diff = smalln
+        end if
+        Reeps_T(kh) = (-1.0_r64/pi)*(diff + sum1 + sum2)


It is sufficient to write -1.0_r64/pi*(...). The brackets are unnecessary. Fortran understands the usual order of arithmetic.

dpaulzc added 2 commits October 25, 2024 11:47

Added Hilbert transform

9f0aba8

Changed app/elphbolt back to previous state

88e9e71

nakib requested changes Oct 25, 2024

View reviewed changes

nakib reviewed Oct 25, 2024

View reviewed changes

Addressed feedback to the pull request

89a9fbb

nakib reviewed Oct 25, 2024

View reviewed changes

nakib merged commit 439d37f into prot Oct 25, 2024
2 checks passed

dpaulzc deleted the paul branch November 6, 2024 17:24

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Hilbert transform #157

Hilbert transform #157

dpaulzc commented Oct 25, 2024

nakib left a comment

nakib Oct 25, 2024

nakib Oct 25, 2024

nakib Oct 25, 2024

nakib Oct 25, 2024

nakib Oct 25, 2024

nakib Oct 25, 2024

nakib Oct 25, 2024

nakib Oct 25, 2024 •

edited

Loading

nakib Oct 25, 2024

nakib Oct 25, 2024

nakib Oct 25, 2024

nakib Oct 25, 2024

nakib Oct 25, 2024

nakib Oct 25, 2024

nakib Oct 25, 2024

nakib Oct 25, 2024

nakib left a comment

nakib Oct 25, 2024

nakib Oct 25, 2024

nakib Oct 25, 2024

		real(r64) :: in1, in2, fin, bl !, domegas
		integer(i64) :: nx, kh, np

Hilbert transform #157

Hilbert transform #157

Conversation

dpaulzc commented Oct 25, 2024

nakib left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nakib Oct 25, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nakib left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nakib Oct 25, 2024 •

edited

Loading