PROGRAM SETU
C=======================================================================
C      DATA GENERATING PROGRAM FOR BEM8LINU.FOR
C   PROJECT: DOUBLE-LET   PROJECT NAME: NEUMANN BOUNDARY CONDITION
C    DOMAIN: INFINITE
C    ELEMENT: TWO-NODED LINEAR BOUNDARY  ELEMENT
C     EIJI FUKUMORI DECEMBER 29, 1993 
C=======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      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.D0
C=======================================================================
      PI = 4.* DATAN( 1.D0)
      RADIUS = 1.0D0
      XC1 = 0.D0
      WRITE (*,*) 'RADIUS = ', RADIUS
      WRITE(*,*) 'TYPE IN GAP BETWEEN WIRE ELEMENTS'
      READ (*,*) GAP
      YC1 = RADIUS + GAP/2.D0
      IF ( GAP .LE. 0.D0 ) STOP 'A VALUE OF GAP MUST BE MORE THAN ZERO'
      XC2 =  XC1
      YC2 = -YC1
C=======================================================================
C-----     NORDAL COORDINATES CREATION
      SEG = NC
      HSEG = SEG/2.D0
      DO I = 1 , NC
      XCOORD(I) = RADIUS*COS ( -I*PI/HSEG )
      YCOORD(I) = RADIUS*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----------------  CIRCUMFERENCE = 2*PI*RADIUS
C  MAGNETIC FLUX DENSITY ON SURFACE = ELECTRIC CURRENT/(CIRCUMFERENCE)
C-------- B OR D INTO DOMAIN     ------> DIRECTIN = -1   ---------------
C-------- B OR D OUT FROM DOMAIN ------> DIRECTIN = +1   ---------------
      B = 1.D0/(2.D0*PI*RADIUS)
      DIRECTIN = -1.D0
      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) = B * DIRECTIN 
      BV(I,2) = B * DIRECTIN
      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) = -B * DIRECTIN
      BV(I,2) = -B * DIRECTIN
      END DO
      NE = NC*2
C=======================================================================
C--------INTERNAL POINTS
      NIP = 0
      NIP = NIP + 1
      XI(NIP) = 0.D0
      YI(NIP) = 0.0001D0
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = (YC1-RADIUS)/4.D0
      NIP = NIP + 1
      XI(NIP) = 0.D0
      YI(NIP) = (YC1-RADIUS)/2.D0
      NIP = NIP + 1
      XI(NIP) = 0.D0
      YI(NIP) = (YC1-RADIUS)/4.D0*3.D0
      NIP = NIP + 1
      XI(NIP) = 0.D0
      YI(NIP) = YC1 - RADIUS/4.D0*3.D0
      NIP = NIP + 1
      XI(NIP) = 0.D0
      YI(NIP) = YC1 - RADIUS/4.D0
      NIP = NIP + 1
      XI(NIP) = 0.D0
      YI(NIP) = YC1 + RADIUS/4.D0
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = YC1 + RADIUS/4.D0*3.D0
      NIP = NIP + 1
      XI(NIP) = 0.D0
      YI(NIP) = YC1 + RADIUS*1.D0 + RADIUS/4.D0
      NIP = NIP + 1
      XI(NIP) = 0.D0
      YI(NIP) = YC1 + RADIUS*2.D0
      NIP = NIP + 1
      XI(NIP) = 0.D0
      YI(NIP) = YC1 + RADIUS*3.D0
      NIP = NIP + 1
      XI(NIP) = 0.D0
      YI(NIP) = YC1 + RADIUS*4.D0
      NIP = NIP + 1
      XI(NIP) = 0.D0
      YI(NIP) = YC1 + RADIUS*5.D0
C
      DO I = 3 , 30
      NIP = NIP + 1
      XI(NIP) = 0.
      YI(NIP) = YC1 + RADIUS*I*2.D0
      END DO
C
      NNN = NIP
      DO I = 1 , NNN
      NIP = NIP + 1
      XI(NIP) = 0.D0
      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,*) 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,*) I, XCOORD(I), YCOORD(I)
      END DO
C------ INTERNAL POINT
      WRITE(1,*) NIP
      IF ( NIP .GE. 1 ) THEN
      DO I = 1 , NIP
      WRITE(1,*) I, XI(I), YI(I)
      END DO
      END IF
      CLOSE (1)
      WRITE(*,*) 'NE=',NE
      WRITE(*,*) 'NNODE=',NNODE
C------------- BOUNDARY configuration
      OPEN ( 2, FILE='BOUNDARY.DAT', STATUS='UNKNOWN' )
      DO IEL = 1 , NE
      WRITE (2,*) XCOORD(NODEX(IEL,1)), YCOORD(NODEX(IEL,1))
      WRITE (2,*) XCOORD(NODEX(IEL,2)), YCOORD(NODEX(IEL,2))
      WRITE (2,*)
      END DO
      CLOSE (2)
      STOP
      END