PROGRAM SETNEUM4
C=======================================================================
C    ********************** HYBRID ELEMENT ************************
C              DATA GENERATING PROGRAM FOR BEM8QUDQ.FOR
C                   PROJECT:COAXIAL CABLE ANALYSIS
C   COAXIAL NEUMANN BOUNDARY ON INNER ALONG WITH  DIRICHLET ON OUTER
C              INNER SHAPE: CIRCLE   OUTER SHAPE: RECTANGULAR--
C************************** QUARTER DOMAIN *****************************
C                            DOMAIN: FINITE
C             ELEMENT: 2-NODED LINEAR WITH HYBRID CONSTANT
C                     EIJI FUKUMORI MARCH 19, 2022
C=======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER ( ND=2,MXE=50000,MXN=MXE, MXI=10000 )
C=======================================================================
      DIMENSION NODEX(MXE,ND),XCOORD(MXN), YCOORD(MXN)
      DIMENSION XI(MXI), YI(MXI), IELTYPE(MXE),BV(MXE,ND)
      COMPLEX*16 BASEVEC , UNITVEC 
C=======================================================================
C--------------OUTER BOUNDARY IS DEVIDED INTO 8 SEGMENTS
C---------------- ELEMENTS IN EACH SEGMENT "NEOUTER".
      NEOUTER = 4
C=======================================================================
      PI = 4.* DATAN(1.D0)
      RADIUSA =  6.35D0
      RADIUSB =  9.5D0
      RADIUSC = 12.0D0
      IF ( RADIUSB .LE. RADIUSA ) STOP 'RADIUSB MUST BE >  RADIUSA.'
      IF ( RADIUSC .LT. RADIUSB ) STOP 'RADIUSC MUST BE >= RADIUSB.'
      WRITE (*,*) 'RADIUS-A(INNER) = ', RADIUSA
      WRITE (*,*) 'RADIUS-B(OUTER-VERTICA  L) = ', RADIUSB
      WRITE (*,*) 'RADIUS-C(OUTER-HORIZONTAL) = ', RADIUSC
C=======================================================================
C----------------------- BOUNDARY CONDITIONS
C------ ON OUTER BOUNDARY: DIRICHLET ----> POTENTIAL(POUTER)    = 0
C------ ON INNER BOUNDARY: NEUMANN ------> POTENTIAL DERIVATIVE = BINNER
C---------- B INTO DOMAIN     -----> THE VALUE IS NEGATIVE
C---------- B OUT FROM DOMAIN -----> THE VALUE IS POSITIVE
       Q = 1.D0
       POUTER =  0.D0
       BINNER = -Q/(2.D0*PI*RADIUSA)
       BVN = 0.D0
       WRITE(*,*) 'NEUMANN BOUNDARY VALUE =', BINNER
C=======================================================================
C------------------- NORDAL COORDINATES CREATION
C------------------------- OUTER ELEMENTS -------------
      DX = RADIUSC/NEOUTER
      DXX = (RADIUSC-RADIUSA)/NEOUTER
      DY = RADIUSB/NEOUTER
      DYY = (RADIUSB-RADIUSA)/NEOUTER
C----------- FIRST SEGMENT STARTS FROM (RADIUSC,0) TO (RADIUSC,RADIUSB)
      NODE = 0
      DO I = 1 , NEOUTER
      NODE = NODE + 1
      XCOORD (NODE) =  RADIUSC
      YCOORD (NODE) =  DY * (I-1)
      END DO
C---------- SECOND SEGMENT STARTS FROM (RADIUSC,RADIUSB) TO (0,RADIUSB)
      DO I = 1 , NEOUTER
      NODE = NODE + 1
      XCOORD (NODE) =  RADIUSC - DX*(I-1)
      YCOORD (NODE) =  RADIUSB
      END DO
C----------- THIRD SEGMENT STARTS FROM (0,RADIUSB)  TO (0,0) ----------
      DO I = 1 , NEOUTER
      NODE = NODE + 1
      XCOORD (NODE) =  0.D0
      YCOORD (NODE) =   RADIUSB - DYY*(I-1)
      END DO
C------------ 4TH SEGMENT ALONG INNER CONDUCTOR SURFACE -------------
      BASEVEC = DCMPLX( 0.D0, RADIUSA )
      UNITANG = -(PI/2.D0) / NEOUTER
      UNITVEC =  DCMPLX( DCOS(UNITANG), DSIN(UNITANG) )
      DO I = 1 , NEOUTER
      NODE = NODE + 1
      XCOORD (NODE) =  DREAL( BASEVEC )
      YCOORD (NODE) =  DIMAG( BASEVEC )
      BASEVEC = BASEVEC * UNITVEC 
      END DO
C--------- 5TH SEGMENT STARTS FROM (RADIUSA,0) TO (RADIUSC,0) ----------
      DO I = 1 , NEOUTER
      NODE = NODE + 1
      XCOORD (NODE) =   RADIUSA+DXX*(I-1)
      YCOORD (NODE) =  0.D0
      END DO
C=======================================================================
C                    ELEMENT CREATION
C------- SEGMENT 1 AND 2: OUTER ELENENTS: TWO SEGMENTS ARE DIRICLET-
      DO I = 1 , 2*NEOUTER
      II = I
      JJ = II + 1
      NODEX(I,1) = II
      NODEX(I,2) = JJ
      BV(I,1) = POUTER
      BV(I,2) = POUTER
      IELTYPE(I) = 1
      END DO
C-------------- SEGMENT 3: NATURAL BOUNDARY: VERTICAL SECTION
      DO I = 2*NEOUTER+1 , 3*NEOUTER
      II = I
      JJ = II + 1
      NODEX(I,1) = II
      NODEX(I,2) = JJ
      BV(I,1) = BVN
      BV(I,2) = BVN
      IELTYPE(I) = 2
      END DO
C---------------- SEGMENT 4: INNER CONDUCTOR SURFACE SECTION IS NEUMANN------
      DO I = 3*NEOUTER+1 , 4*NEOUTER
      II = I
      JJ = II + 1
      NODEX(I,1) = II
      NODEX(I,2) = JJ
      BV(I,1) = BINNER
      BV(I,2) = BINNER
      IELTYPE(I) = 2
      END DO
C-------------- SEGMENT 5: NATURAL BOUNDARY: HORIZONTAL SECTION
      DO I = 4*NEOUTER+1 , 5*NEOUTER
      II = I
      JJ = II + 1
      IF ( I .EQ. 5*NEOUTER ) JJ = 1
      NODEX(I,1) = II
      NODEX(I,2) = JJ
      BV(I,1) = BVN
      BV(I,2) = BVN
      IELTYPE(I) = 2
      END DO
C=======================================================================
      NE = 5*NEOUTER
      NNODE = NODE
C=======================================================================
C--------INTERNAL POINTS
      NIP = 20
      REND = DSQRT(RADIUSB**2+RADIUSC**2)
      SINANG = RADIUSB/REND
      COSANG = RADIUSC/REND
      DR = (REND - RADIUSA) / (NIP+1)
      DO I = 1 , NIP
      R = RADIUSA + I*DR
      YI(I) = R*SINANG
      XI(I) = R*COSANG
      END DO
C=======================================================================
C                     MAKING DATA FILES
      OPEN ( 1, FILE='BEMHYB.DAT', STATUS='UNKNOWN' )
      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
      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