CODEWITHSUNDEEP
fortran
si wang
si wang

Reputation: 11

how to express delta function in Fortran

In Fortran, to do the Monte Carlo integration, I want to times an extra factor which contains a Dirac Delta function like delta(1-x), how do I express delta function in Fortran? Is there anything like log(x)? or do I have to define a function? how?

Upvotes: 0

Views: 633

Answers (1)

Federico Perini
Federico Perini

Reputation: 1416

This should do the job:

program test_dirac
        use iso_fortran_env
        implicit none

        integer :: i
        real(real64) :: x

        do i=-10,10
           x = real(i,real64)
           print *, 'x=',x,' dirac=',dirac(x)
        end do

        contains


        elemental real(real64) function dirac(x)
            real(real64), intent(in) :: x
            real(real64) :: absx
            absx = abs(x)
            dirac = merge(1.0_real64/absx,0.0_real64,absx<epsilon(x))
        end function dirac

end program test_dirac

On gfortran-10, it returns:

 x=  -10.000000000000000       dirac=   0.0000000000000000
 x=  -9.0000000000000000       dirac=   0.0000000000000000
 x=  -8.0000000000000000       dirac=   0.0000000000000000
 x=  -7.0000000000000000       dirac=   0.0000000000000000
 x=  -6.0000000000000000       dirac=   0.0000000000000000
 x=  -5.0000000000000000       dirac=   0.0000000000000000
 x=  -4.0000000000000000       dirac=   0.0000000000000000
 x=  -3.0000000000000000       dirac=   0.0000000000000000
 x=  -2.0000000000000000       dirac=   0.0000000000000000
 x=  -1.0000000000000000       dirac=   0.0000000000000000
 x=   0.0000000000000000       dirac=                  Infinity
 x=   1.0000000000000000       dirac=   0.0000000000000000
 x=   2.0000000000000000       dirac=   0.0000000000000000
 x=   3.0000000000000000       dirac=   0.0000000000000000
 x=   4.0000000000000000       dirac=   0.0000000000000000
 x=   5.0000000000000000       dirac=   0.0000000000000000
 x=   6.0000000000000000       dirac=   0.0000000000000000
 x=   7.0000000000000000       dirac=   0.0000000000000000
 x=   8.0000000000000000       dirac=   0.0000000000000000
 x=   9.0000000000000000       dirac=   0.0000000000000000
 x=   10.000000000000000       dirac=   0.0000000000000000

Upvotes: 1

Related Questions