implicit real*8(a-h,o-z)
h=0.05d0
omg=1.d0
dafmax=0.d0
omgd=omg*sqrt(1.d0-h**2)
do i=0, 10000
t=real(i)/100.d0
 u=1-exp(-h*omg*t)*(omg/omgd*h*sin(omgd*t)+cos(omgd*t))
 print*, t, u
if(u>dafmax) dafmax=u
end do
pi=acos(-1.d0)
print *, "#       DAF_max=",dafmax
print *, "# 1+e^{-0.05π/√1-0.05^2}=", 1.d0+exp(-0.05d0*pi/sqrt(1-0.05**2))
end
