      program besselplots ! bernd berg, Nov 20 2016.
c subroutines adapted from numerical recipes.
c potts time series (ts) with systematic metropolis updating.
c n    number of point in the plot
c xmax maximum values of x.
      parameter(n=512,xmax=10.0d00,iuo=6,iud1=11,iud=12)
c
      write(iuo,'(/," program bessel0_plots.f:",/)')
      write(iuo,'(" gnuplot: n, xmax =",i9,f10.1,/)') n,xmax
      open(iud1,file='bessj0.d',form='formatted',status='unknown')
      dx=xmax/n
      do i=0,n
        x=i*dx
        y0=bessj0(x)
        write(iud1,'(2f15.6)') x,y0
      enddo
      close(iud1)
c
      x01=2.405 ! From Table 2.1 of Lecture Notes.
      y01=bessj0(x01)
      write(iuo,'(" x01,y01 =",2f15.6,3x,"(Test of zero.)")') x01,y01
c
      s01=sinh(x01)
      write(iuo,'(" x01,s01 =",2f15.6,/)') x01,s01
      y0=bessj0(zero)
      write(iuo,'(" zero,y0 =",2f15.6)') zero,y0
      x01h=x01/2
      s=sinh(x01h)
      write(iuo,'(" x01/2,s =",2f15.6)') x01h,s
      rsh=s/s01 
      write(iuo,'(" Potential at rho=0: s/s01 =",1f12.6,/)') rsh
c
      y=bessj0(x01h)
      write(iuo,'(" x01/2,y =",2f15.6)') x01h,y
      rsh=rsh*y
      write(iuo,'(" Potential at rho=x01/2: s/s01=",1f9.6,/)') rsh
c
      stop "bessel0_plots: all done, gnuplot plotbessj.plt."
      end

      include 'bessj0.f'
