module para
   double precision,dimension(1:8)::x,y
end module para

program integral
   use para
   implicit none
   integer,parameter::m=8
   double precision exact,aint,aerr,rerr
   integer n
   integer j

   exact=datan(1.0d0)*4.0d0/2.0d0
   n=4
   do j=1,8
      n=n*2
      call trapezoid(n,aint)
      aerr=dabs(aint-exact)
      rerr=aerr/exact
      x(j)=dble(n)
      y(j)=rerr
      write(*,100) n,aint,aerr
   end do
   100 format(1x,'n=',1x,i3,1x,'integral=',1x,f7.5,1x,'error=',1x,f7.5,1x)

   call init
   call viewport(0.2, 0.2, 0.8, 0.8)
   call xyworld(0.0, 0.0, 4.0, 4.0)
   call frame
   call plotd(m)
   call fin
   stop
   end

subroutine trapezoid(n,ans)
   implicit none
   integer n
   double precision f,x,dx,ans
   integer i

   f(x)=sqrt(1.0d0-x*x)

   dx=2.0d0/dble(n)
   x=-1.0d0
   ans=f(x)
   do i=1,n-1
      x=dble(i)*dx-1.0d0
      ans=ans+2.0d0*f(x)
   end do
   x=1.0
   ans=(ans+f(x))*dx/2.0d0
   return
   end

subroutine plotd(n)
   use para
   implicit none
   integer n
   integer i
   do i=1,n
      x(i)=log10(x(i))
      y(i)=log10(y(i)*100000.0d0)
      write(*,*) x(i),y(i)
   end do
   call linewidth(1.5)
   call linety(1)
   call setgray(0.0)
   call plot(real(x(1)), real(y(1)), 3)
   do i=2,n
      call plot(real(x(i)),real(y(i)),2)
   end do
   call stroke
   call circ(real(x(1)),real(y(1)),0.03)
   call stroke
   do i=2,n
      call circ(real(x(i)),real(y(i)),0.03)
      call stroke
   end do
   return
   end

subroutine frame
   implicit none
   integer ix,nx,iy,ny,ixx,iyy
   double precision x,y,dx,dy
   call linewidth(1.0)
   call rect(0.0, 0.0, 4.0, 4.0)
   nx=4
   do ix=0,nx-1
      call setgray(0.3)
      x =dble(ix)
      dy=4.0
      if(ix/=0)then
         call plot(real(x),0.0,3)
         call plot(real(x),real(dy),2)
         call stroke
      end if
      do ixx=2,9
         call setgray(0.8)
         x =dble(ix)+log10(dble(ixx))
         call plot(real(x),0.0,3)
         call plot(real(x),real(dy),2)
         call stroke
      end do
   end do
   ny=4
   do iy=0,ny-1
      call setgray(0.5)
      y =dble(iy)
      dx=4.0
      if(iy/=0)then
         call plot(0.0,real(y),3)
         call plot(real(dx),real(y),2)
         call stroke
      end if
      do iyy=2,9
         call setgray(0.8)
         y =dble(iy)+log10(dble(iyy))
         call plot(0.0,real(y),3)
         call plot(real(dx),real(y),2)
         call stroke
      end do
   end do
   call setgray(0.0)
   call textx(3.5, 0.0,'n')
   call textx(0.0, 0.02, '1')
   call textx(1.0, 0.02, '10')
   call textx(2.0, 0.02, '100')
   call textx(3.0, 0.02, '1000')
   call textx(4.0, 0.02, '10000')
   call texty(0.0, 3.5, 'eps')
   call texty(0.0, 0.0, '0.00001')
   call texty(0.0, 1.0, '0.0001')
   call texty(0.0, 2.0, '0.001')
   call texty(0.0, 3.0, '0.01')
   call texty(0.0, 4.0, '0.1')
   return
   end