program f90_simpson
    ! Calculation program of pi by Simpson rule
    implicit none
    integer::n,i
    real(8)::a,b
    real(8)::h,s,xa,xb,xc,ya,yb,yc
    character(len=20)::dummy
    real(8),parameter::pi=atan(1.0D0)*4.0D0

    call getarg(1,dummy)
    read(dummy,*) n

    a=-1.0D0
    b= 1.0D0
    h=(b-a)/dble(n)

    xa=a
    s=0.0D0
    ya=FUNC(xa)
    do i=1,n
        xb=xa+h
        xc=0.5D0*(xa+xb)
        yb=FUNC(xb)
        yc=FUNC(xc)
        s=s+(xb-xa)/6.0D0*(ya+4.0D0*yc+yb)
        xa=xb
        ya=yb
    end do
    write(6,'(i10,2(e20.12))') n,2.0D0*s,2.0D0*s-pi

stop

contains

    real(8) function FUNC(x)
        real(8),intent(in)::x
        real(8)::d
        FUNC=0.0D0
        d=1.0D0-x*x
        if(0.0D0<=d)FUNC=sqrt(d)
    end function FUNC

end program f90_simpson
