PROGRAM SETLINQ16
C======================================================================
C                      INPUT DATA SET-UP PROGRAM
C                                 FOR
C                       BOUNDARY ELEMENT METHOD
C                              APPLIED TO
C              HELMHOLTZ EQUATION (HAMONIC WAVE EQUATION )
C                 ELEMENT TYPE: LINEAR ELEMENT
C    ********** WAVE AROUND CYLINDER IN INIFINITE DOMAIN ***********
C                PROGRAMMED BY EIJI FUKUMORI 2025 JAN 9
C======================================================================
      PARAMETER (MXE=1000, MXN=MXE, MXI=20, ND=2)
      DIMENSION NODEX(MXE,ND),IELTYPE(MXE),X(MXN),Y(MXN),
     *          XI(MXI),YI(MXI)
      COMPLEX*16 BV(MXE,ND)
C-------  RL = WAVE LENGTH,  DIAMETER = DIAMETER OF CYLINDER 
C------- NDIV = NUMBER OF ELEMENTS ON THE CYLINDER
C
C------- BASIC PARAMETERS
      NDIV=32
      RL = 1.D0
      DIAMETER = 1.6D0
C------- CONSTANTS
      PI = 4.D0*DATAN ( 1.D0 )
      OMEGA = 2.D0*PI/RL
      NE = NDIV
      NNODE = NE
C---------- NODAL COORDINATES ( CLOCKWISE )
      DTHETA = -2.D0*PI/NDIV
      R = DIAMETER/2.D0
      DO I = 1 , NDIV
      THETA = DTHETA*(I-1)
      X(I) = R*COS(THETA)
      Y(I) = R*SIN(THETA)
      END DO
C--------- BOUNDARY VALUES OF U AND V  ALSO NODEX(I,J)
      DO I = 1 , NE
      NODE1 = I
      NODE2 = I + 1
      IF ( I .EQ. NE ) NODE2 = 1
      NODEX(I,1) = NODE1
      NODEX(I,2) = NODE2
      DX = X(NODE2) - X(NODE1)
      DY = Y(NODE2) - Y(NODE1)
      DS = SQRT ( DX*DX + DY*DY )
      U1 =  OMEGA * (DY/DS) * SIN(OMEGA*X(NODE1))
      U2 =  OMEGA * (DY/DS) * SIN(OMEGA*X(NODE2))
      V1 = -OMEGA * (DY/DS) * COS(OMEGA*X(NODE1))
      V2 = -OMEGA * (DY/DS) * COS(OMEGA*X(NODE2))
      BV(I,1) = CMPLX( U1 , V1 )
      BV(I,2) = CMPLX( U2 , V2 )
      IELTYPE(I) = 2
      END DO
C-------- INTERNAL POINTS
      NIPT = 5
      DO I = 1 , 5
      XI(I) = I
      YI(I) = 0.
      END DO
C-------- FILE CREATION
      OPEN ( 1, FILE='BEM1.DAT', STATUS='UNKNOWN' )
      WRITE (1,*) OMEGA
      WRITE (1,*) NE
      DO I = 1 , NE
      WRITE (1,*) I,(NODEX(I,J),J=1,ND),IELTYPE(I), 
     * ( DREAL(BV(I,J)), DIMAG(BV(I,J)), J=1,ND )
      END DO
      WRITE (1,*) NNODE
      DO I = 1 , NNODE
      WRITE (1,*) I, X(I), Y(I)
      END DO
      WRITE (1,*) NIPT
      IF ( NIPT .GE. 1 ) THEN
      DO I = 1  , NIPT
      WRITE (1,*) I, XI(I), YI(I)
      END DO
      END IF
      CLOSE (1)
C------------------ BOUNDARY ELEMENTS ---------------
      OPEN ( 1, FILE='BOUNDARYELEMENTS.OUT', STATUS='UNKNOWN' )
      DO I = 1 , NE
      WRITE (1,*) X(NODEX(I,1)), Y(NODEX(I,1))
      WRITE (1,*) X(NODEX(I,2)), Y(NODEX(I,2))
      WRITE (1,*)
      END DO
      CLOSE (1)
      STOP
      END