Line data    Source code

       1             : MODULE m_lorentzian_smooth
       2             : 
       3             :    USE m_constants
       4             : 
       5             :    IMPLICIT NONE
       6             : 
       7             :    !PARAMETER FOR LORENTZIAN SMOOTHING
       8             :    REAL,    PARAMETER :: cut              = 1e-8
       9             : 
      10             :    INTERFACE lorentzian_smooth
      11             :       PROCEDURE lorentzian_smooth_r, lorentzian_smooth_c
      12             :    END INTERFACE
      13             : 
      14             :    CONTAINS
      15             : 
      16           0 :    SUBROUTINE lorentzian_smooth_r(e,f,sigma,n)
      17             : 
      18             :       INTEGER, INTENT(IN)    :: n
      19             :       REAL,    INTENT(INOUT) :: f(:)
      20             :       REAL,    INTENT(IN)    :: sigma , e(:)
      21             : 
      22             :       COMPLEX, ALLOCATABLE :: f_c(:)
      23             : 
      24           0 :       f_c = f
      25           0 :       CALL lorentzian_smooth_c(e,f_c,sigma,n)
      26           0 :       f = REAL(f_c)
      27             : 
      28           0 :    END SUBROUTINE lorentzian_smooth_r
      29             : 
      30             :    !This is essentially smooth out of m_smooth but with a lorentzian distribution
      31         196 :    SUBROUTINE lorentzian_smooth_c(e,f,sigma,n)
      32             : 
      33             :       INTEGER, INTENT(IN)    :: n
      34             :       COMPLEX, INTENT(INOUT) :: f(:)
      35             :       REAL,    INTENT(IN)    :: sigma
      36             :       REAL,    INTENT(IN)    :: e(:)
      37             : 
      38             :       REAL :: dx
      39         196 :       COMPLEX :: f0(n)
      40         196 :       REAL :: ee(n)
      41             :       INTEGER :: ie , je , j1 , j2 , numPoints
      42             : 
      43     1058596 :       f0 = f
      44     1058596 :       f = 0.0
      45     1058596 :       ee = 0.0
      46         196 :       dx = e(2)-e(1)
      47             : 
      48     1058596 :       DO ie =1,  n
      49             : 
      50     1058400 :          ee(ie) = 1/pi_const * sigma/dx * 1.0/((ie-1)**2+(sigma/dx)**2)
      51             : 
      52     1058596 :          IF ( ee(ie).LT.cut ) EXIT
      53             :       ENDDO
      54         196 :       numPoints = ie - 1
      55             : 
      56     1058596 :       DO ie = 1, n
      57             : 
      58     1058400 :          j1 = ie - numPoints + 1
      59     1058400 :          j1 = MERGE(1,j1,j1.LT.1)
      60     1058400 :          j2 = ie + numPoints - 1
      61     1058400 :          j2 = MERGE(n,j2,j2.GT.n)
      62  5716418596 :          DO je = j1 , j2
      63  5716418400 :             f(ie) = f(ie) + ee(IABS(je-ie)+1)*f0(je)
      64             :          ENDDO
      65             : 
      66             :       ENDDO
      67             : 
      68         196 :    END SUBROUTINE lorentzian_smooth_c
      69             : 
      70             : END MODULE m_lorentzian_smooth