      PROGRAM STRMFC
C=======================================================================
C     COMPUTATION OF STREAM FUNCTION
C     EIJI FUKUMORI NOVEMBER 28, 1995 
C=======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER ( N=11, M=41 )
      DIMENSION XCOORD(M), YCOORD(N), FAI(N)
      COMMON /BASEDAT/ MXN
      COMMON /BASEDATA/ C, U, A, D
	COMMON /ANGLE/ PI
C=======================================================================
      Q = 20.
      D = 10.
      U=1.
      A=5.
C---------
      PI = 4.*DATAN( 1.D0 )
      C = Q/(2.*PI)
      MXN = 1000
C======================================================================= 
      DIV = (N-1)/D
      DO I = 1, N
      YCOORD(I) = (I-1)/DIV
      END DO
      DO I = 1, M
      XCOORD(I) = (M-I)/DIV
      END DO
C=======================================================================
C                     MAKING DATA FILES
C---------- 'PROJECT'.JNK
      IR = 1
      OPEN ( IR, FILE='EXACT.DAT', STATUS = 'UNKNOWN' )
      WRITE (IR,'(7H X//Y->, 21E15.6)') ( YCOORD(I),I=1,N )
      DO I = 1,M
      DO J = 1, N
      FAI(J) =  STREAM( XCOORD(I), YCOORD(J) )
      END DO
      WRITE (IR,'(22E15.6)') XCOORD(I), FAI
      END DO
      CLOSE (1)
      STOP
      END
C
C
      FUNCTION STREAM( X, Y )
      IMPLICIT REAL * 8 ( A-H , O-Z )
      COMMON /BASEDAT/ MXN
      COMMON /BASEDATA/ C, U, A, D
      SUM0 = DDATAN(Y,(X+A))-DDATAN(Y,(X-A))
      SUM1 = 0.
      SUM2 = 0.
      DO M = 1 , MXN
      SUM1 = SUM1 + DDATAN((Y+2.*D*M),(X+A))-DDATAN((Y+2.*D*M),(X-A))
      SUM2 = SUM2 + DDATAN((Y-2.*D*M),(X+A))-DDATAN((Y-2.*D*M),(X-A))
      END DO
      STREAM = U*Y+C* ( SUM0 + SUM1 + SUM2 )
      RETURN
      END
C
C
      FUNCTION DDATAN ( Y,X )
      IMPLICIT REAL * 8 ( A-H , O-Z )
      COMMON /ANGLE/ PI
      IF ( (X .EQ. 0.) .AND. (Y.EQ. 0.) ) THEN
      DDATAN = PI
      RETURN
      END IF
      IF ( (X .EQ. 0.) .AND. (Y.GT.0.) ) THEN
      DDATAN = PI/2.
      RETURN
      END IF
      IF ( (X .EQ. 0.) .AND. (Y.LT.0.) ) THEN
      DDATAN = 1.5*PI
      RETURN
      END IF
      ANG = DATAN( Y/X )
      IF ( (X .GT. 0. ) .AND. (Y.GE. 0.) ) THEN
      DDATAN = ANG
      RETURN
      END IF
      IF ( (X .GT. 0. ) .AND. (Y.LT. 0.) ) THEN
      DDATAN = 2.*PI - DABS(ANG)
      RETURN
      END IF
      IF ( (X .LT. 0.) .AND. (Y.GE. 0.) ) THEN
      DDATAN = PI - DABS(ANG)
      RETURN
      END IF
      IF ( (X .LT. 0.) .AND. (Y.LT. 0.) ) THEN
      DDATAN = PI + DABS(ANG)
      RETURN
      END IF
      RETURN
      END