Intensity using bessel function tending to infinity

Question

Making it short, my code is supposed to return a txt with my intensity values, instead, for all rs but 0, my intensity returns a value of +infinity. Don't know where my mistake is. This exercise is supposed to make us practice integration via Simpson's 1/3 method. All Bessels Jx txt values are working fine, the only problem relies within my intensity file. First code block is responsible for creating and filling Bessel Jx values in a txt file. Second part is responsible for creating and filling intensity values through a Bessel function (this is where the error is supposed to be, but i'm not sure). Third and fourth blocks are the Simpson 1/3 method and my Bessel function, respectively.

program intensidade
implicit none

real,parameter::pi=acos(-1.),lambda=500e-9
real::k,r,kr,intensidade1
real,external::bessel,simpson13
real::i
integer::j

open(0,file='besselj0.txt')
open(1,file='besselj1.txt')
open(2,file='besselj2.txt')
open(3,file='intensidade.txt')

do j=0,2
 i=0
  do while (i<=20)
   write(j,*)i,simpson13(bessel,j,i,0.,pi)
 i=i+1
 enddo
enddo

close(0);
close(1);
close(2);


r=0
k=2*pi/lambda
kr=k*r*10e-6

do while (r<=1)
  if(r==0) then
   write(3,*)r,(1/2)**2
  else
   write(3,*)r,(simpson13(bessel,1,kr,0.,pi)/kr)**2
  endif
 r=r+0.1
enddo

close(3)
pause
end program intensidade


real function simpson13(funcao,m,x,a,b)
implicit none

real,external::funcao
real,intent(in)::a,b,x
integer,intent(in)::m
integer::i
real::h

h=(b-a)/1000
simpson13=funcao(m,x,a)-funcao(m,x,b)

  do i=1,499
   simpson13=simpson13+4*funcao(m,x,a+h*(2*i-1))+2*funcao(m,x,a+2*h*i)
  enddo
simpson13=(h/3)*simpson13

end function simpson13

real function bessel(m,x,teta)
implicit none

real,parameter::pi=acos(-1.)
real,intent(in)::x,teta
integer,intent(in)::m

bessel=cos(m*teta-x*sin(teta))/pi

end function bessel

Answer 1

The main error arises because kr is not redefined in each loop.

Other improvements

• style: align and indent your code
• use file units provided by the system, i.e. open (newunit=...)
• remove the pause command

Furthermore, your line write(..)r,(1/2)**2 uses integer arithmetic s.t. 1/2 yields zero and (1/2)**2 is zero as well.

The following is a possible way to rewrite your program

program intensidade
  implicit none

  real, parameter :: pi=acos(-1.0), lambda=500e-9
  real            :: k, r, kr
  real, external  :: bessel, simpson13
  integer         :: ir, funit, ix, m
  character(128)  :: fname

  do m = 0, 2
    write (fname, "(A7, I1, A4)") 'besselj', m, '.txt'
    open (newunit=funit, file=fname)
    do ix = 0, 20
      write (funit, *) ix, simpson13(bessel, m, real(ix), 0.0, pi)
    end do
    close (funit)
  end do

  open (newunit=funit, file='intensidade.txt')
  r = 0
  k = 2*pi/lambda
  write (funit, *) r, (0.5)**2
  do ir = 1, 10
    r  = ir/10.0
    kr = k*r*10e-6
    write (funit, *) r, (simpson13(bessel, 1, kr, 0.0, pi)/kr)**2
  end do

  close (funit)
end program