PROGRAM SETU
C=======================================================================
C      DATA GENERATING PROGRAM FOR BEM8LINU.FOR
C   PROJECT: DOUBLE-LET   PROJECT NAME: SKIN CUURENT
C    DOMAIN: INFINITE
C    ELEMENT: FOUR-NODED ISOPARAMETRIC ELEMENT
C    DOMAIN DISCRETIZATION: UNEVEN ELEMENTS WITH VERTICAL SCAN
C     EIJI FUKUMORI DECEMBER 29, 1993 
C=======================================================================
      PARAMETER ( ND=2,MXE=500,MXN=MXE, MXI=100 )
C=======================================================================
      DIMENSION NODEX(MXE,ND),XCOORD(MXN), YCOORD(MXN)
      DIMENSION XI(MXI), YI(MXI), IELTYPE(MXE),BV(MXE,ND)
C=======================================================================
      NC=16
      FMU = 1.
C=======================================================================
      PI = 4.* ATAN( 1.)
      RADIUS = 1.
      XC1 = 0.
      YC1 = 2.
      WRITE(*,*) 'TYPE IN HALF OF D.  I.E. (D/2)'
      READ (*,*) YC1
      IF ( YC1 .LT. 2*RADIUS ) STOP 'YC1>2.'
      XC2 = XC1
      YC2 = -YC1
C=======================================================================
C-----     RADIUS OF WELL = 1. I.E., A=1
      SEG = NC
      HSEG = SEG/2.
      DO I = 1 , NC
      XCOORD(I) = COS ( -I*PI/HSEG )
      YCOORD(I) = SIN ( -I*PI/HSEG )
      END DO
      DO I = 1 , NC
      XCOORD(I+NC) = XCOORD(I) + XC2
      YCOORD(I+NC) = YCOORD(I) + YC2
      XCOORD(I) = XCOORD(I) + XC1
      YCOORD(I) = YCOORD(I) + YC1
      END DO
      NNODE = NC*2
C=======================================================================
C                    ELEMENT CREATION
C       DISCHARGE (Q) = 1; SKIN CUURENT PER SKIN LENGTH = 1/(2PI)
      DO I = 1 , NC
      J = I+1
      IF ( I .EQ. NC ) J = 1
      NODEX(I,1) = I
      NODEX(I,2) = J
      IELTYPE(I) = 2
      BV(I,1) = 1./(2.*PI)
      BV(I,2) = 1./(2.*PI)
      END DO
      DO I = 1+NC , NC*2
      J = I+1
      IF ( I .EQ. NC*2 ) J = 1+NC
      NODEX(I,1) = I
      NODEX(I,2) = J
      IELTYPE(I) = 2
      BV(I,1) = -1./(2.*PI)
      BV(I,2) = -1./(2.*PI)
      END DO
      NE = NC*2
C=======================================================================
C--------INTERNAL POINTS
      NIP = 0
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = 0.0001
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = (YC1-RADIUS)/4.
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = (YC1-RADIUS)/2.
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = (YC1-RADIUS)/4.*3.
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = YC1 - RADIUS/4.*3.
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = YC1 - RADIUS/4.
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = YC1 + RADIUS/4.
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = YC1 + RADIUS/4.*3.
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = YC1 + RADIUS*1. + RADIUS/4.
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = YC1 + RADIUS*2.
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = YC1 + RADIUS*3.
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = YC1 + RADIUS*4.
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = YC1 + RADIUS*5.
C
      DO I = 3 , 30
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = YC1 + RADIUS*I*2
      END DO
C
      NNN = NIP
      DO I = 1 , NNN
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = -YI(I)
      END DO
C
C
      DO I = 1 , NIP-1
      YMIN = YI(I)
      K = I
      DO J = I+1 , NIP
      IF ( YI(J) .LT. YMIN ) THEN
      YMIN = YI(J)
      K = J
      END IF
      END DO
      YI(K) = YI(I)
      YI(I) = YMIN
      END DO
C=======================================================================
C                     MAKING DATA FILES
      OPEN ( 1, FILE='BEM1.DAT', STATUS='UNKNOWN' )
      WRITE (1,*) FMU
      WRITE (1,*) NE
      DO I = 1 , NE
      WRITE (1,100) I,(NODEX(I,J),J=1,ND),IELTYPE(I), (BV(I,J),J=1,ND)
      END DO
      NNODE = NE
      DO I = 1 , NNODE
      WRITE (1,200) I, XCOORD(I), YCOORD(I)
      END DO
C------ INTERNAL POINT
      WRITE(1,*) NIP
      IF ( NIP .GE. 1 ) THEN
      DO I = 1 , NIP
      WRITE(1,200) I, XI(I), YI(I)
      END DO
      END IF
      CLOSE (1)
  100 FORMAT ( 4I5,2E20.10 )
  200 FORMAT ( I5  , 2E20.10 )
      WRITE(*,*) 'NE=',NE
      WRITE(*,*) 'NNODE=',NNODE
      STOP
      END