PROGRAM CIRCLEAREACOMP
C====================== CIRCLE AREA COMPUTATION =======================
      IMPLICIT REAL*8 ( A-H , O-Z )
      COMPLEX*16 Z, Z0
      R = 1.D0
      PI = 4.D0*DATAN(1.D0)
      OPEN ( 1, FILE="AREA.VSN", STATUS="UNKNOWN" )
C***********************************************************************
      DO NE = 2,1000
      DTH = 2.D0*PI/(2*NE)
      Z0 = R*CDEXP(DCMPLX(0.D0,DTH))
      AREA = 0.D0
      Z = DCMPLX(R,0.D0)
C
C=======================================================================
      DO IEL = 1 , NE
      X1 = DREAL(Z)
      Y1 = DIMAG(Z)
      Z = Z*Z0
      X2 = DREAL(Z)
      Y2 = DIMAG(Z)
      Z = Z*Z0
      X3 = DREAL(Z)
      Y3 = DIMAG(Z)
C
      DX = X3-X1
      DY = Y3-Y1
      D2X = X1-2.D0*X2+X3
      D2Y = Y1-2.D0*Y2+Y3
      SUB1 = (1.D0/3.D0)*(DX*D2Y+DY*D2X/2.D0) + X2*DY
      SUB2 = (1.D0/3.D0)*(DY*D2X+DX*D2Y/2.D0) + Y2*DX
      AREA = AREA + 0.5D0*(SUB1-SUB2)
      END DO
C=======================================================================
      WRITE (1,*) NE, AREA
      END DO
C***********************************************************************
      CLOSE (1)
      STOP 'TERMINATION'
      END