DE積分 振動型積分 cos版

\int_0^\infty f(x)\cos\omg x dx
f(x)=-exp(x)

module exp_lib
implicit none
real(8) pi

contains

subroutine calf(f,x)
real(8),intent(in) :: x
real(8),intent(out) :: f
f=exp(-x)
end subroutine calf

subroutine de2_int(sum,em,omg,n,calf)
integer,intent(in) :: n
real(8),intent(in) :: em,omg
real(8),intent(out) :: sum
real(8) h,t,e6,phi,dphi,x,f
integer i
h=pi/(em*omg)
t=0.d0+h*0.5d0
!!$ e6=1.d0
!!$ phi=1.d0/6.d0
!!$ dphi=1.d0/2.d0
e6=exp(-6.d0*sinh(t))
phi=t/(1.d0-e6)
dphi=(1.d0+(-6.d0*t*cosh(t)-1.d0)*e6)/(1.d0-e6)**2
x=em*phi
call calf(f,x)
sum=sum+f*cos(omg*x)*dphi
do i=1,n
t=i*h+h*0.5d0
e6=exp(-6.d0*sinh(t))
phi=t/(1.d0-e6)
dphi=(1.d0+(-6.d0*t*cosh(t)-1.d0)*e6)/(1.d0-e6)**2
x=em*phi
call calf(f,x)
sum=sum+f*cos(omg*x)*dphi
t=-t+h*0.5d0
e6=exp(-6.d0*sinh(t))
phi=t/(1.d0-e6)
dphi=(1.d0+(-6.d0*t*cosh(t)-1.d0)*e6)/(1.d0-e6)**2
x=em*phi
call calf(f,x)
sum=sum+f*cos(omg*x)*dphi
end do
sum=sum*em*h
end subroutine de2_int

end module exp_lib

program exp
use exp_lib
implicit none
external calf
integer n
real(8) em,omg,sum
pi=acos(-1.d0)
sum=0.d0
em=100.d0
omg=1.d0
n=75
call de2_int(sum,em,omg,n,calf)
write(*,*) sum
end program exp