Line data    Source code

       1             : MODULE m_smooth
       2             : !
       3             : !     the function f(x), defined on a linear mesh,
       4             : !     is smoothened by a gaussian
       5             : !
       6             : 
       7             :    INTERFACE smooth
       8             :       PROCEDURE smooth_r, smooth_c
       9             :    END INTERFACE
      10             : 
      11             :    CONTAINS
      12             : 
      13          89 :    SUBROUTINE smooth_r(e,f,sigma,n)
      14             : 
      15             :       INTEGER, INTENT(IN)    :: n
      16             :       REAL,    INTENT(INOUT) :: f(:)
      17             :       REAL,    INTENT(IN)    :: sigma , e(:)
      18             : 
      19          89 :       COMPLEX, ALLOCATABLE :: f_c(:)
      20             : 
      21         267 :       ALLOCATE(f_c(n))
      22      117747 :       f_c = f
      23          89 :       CALL smooth_c(e,f_c,sigma,n)
      24      117658 :       f = REAL(f_c)
      25             : 
      26          89 :    END SUBROUTINE smooth_r
      27             : 
      28          89 :    SUBROUTINE smooth_c(e,f,sigma,n)
      29             : 
      30             :    USE m_constants, ONLY: tpi_const
      31             :    IMPLICIT NONE
      32             : 
      33             : !    Arguments
      34             :    INTEGER, INTENT(IN)    :: n
      35             :    COMPLEX, INTENT(INOUT) :: f(:)
      36             :    REAL,    INTENT(IN)    :: sigma , e(:)
      37             : 
      38             : !    Locals
      39             :    REAL :: c , c2 , dx
      40             :    INTEGER :: i , j , j1 , j2 , m1, m
      41          89 :    COMPLEX :: f0(n)
      42             : 
      43          89 :    REAL, ALLOCATABLE :: ee(:)
      44             : 
      45          89 :    dx = e(2) - e(1)
      46          89 :    c = dx/(sigma*sqrt(tpi_const))
      47          89 :    c2 = -0.5*(dx/sigma)**2
      48             : 
      49          89 :    m = NINT(sqrt(log(1.0e-8/c)/c2))+1
      50         267 :    ALLOCATE ( ee(m) )
      51       10326 :    DO i = 1, m
      52       10326 :       ee(i) = c * exp(c2*(i-1)**2)
      53       10326 :       IF ( ee(i).LT.1.E-8 ) EXIT
      54             :    ENDDO
      55          89 :    m1=i-1
      56      117658 :    f0 = f
      57      117658 :    f = 0.
      58             : 
      59      117658 :    DO i = 1 , N
      60      117569 :       j1 = i - m1 + 1
      61      117569 :       IF ( j1.LT.1 ) j1 = 1
      62      117569 :       j2 = i + m1 - 1
      63             :       IF ( j2.GT.N ) j2 = N
      64    25860447 :       DO j = j1 , j2
      65    25860358 :          f(i) = f(i) + ee(IABS(j-i)+1)*f0(j)
      66             :       ENDDO
      67             :    ENDDO
      68          89 :    DEALLOCATE ( ee )
      69             : 
      70          89 :    END SUBROUTINE smooth_c
      71             : END MODULE m_smooth