PROGRAM DIPOLE
C*************************************************************************
C----- COMPUTAION OF DIRECT(OR SELF) IMPEDFANCE(Z) OF DIPOLE ANTENNA
C                       ****** ASSUMPTION ******
C----- FEEDING POINT AT ANY POINT BETWEEN -LAMBDA/4 AND +LAMBDA/4 --------
C--------  Radius of cross section of element (A)  = 0
C-"A" IS OMITED SINCE A is very small in comparison with lengh of anttenna
C-------- NE = number of segments where integration is performed
C-------- NFREQUENCY = NUMBER OF WAVES ON DIPOLE ANTENNA
C-------- NFEEDSEG = NUMBER OF FEEDING POINTS
C-------- WAVE POWER IS FED THROUGH EACH FEEDING POINT
C-------- H = (LAMBDA/4) / (LAMBDA), DIMENTIONLESS NUMBER
C-------- 27 APRIL 2023 EIJI FUKUMORI
C*************************************************************************
      PARAMETER ( INTEPT=6 )
      IMPLICIT REAL*8 ( A-H , O-Z )
      COMPLEX*16  S
      DIMENSION SAI(INTEPT),W(INTEPT)
C*****************************************************************
C---------------------- BASIC PARAMETERS -------------------------
      PI = 4.D0*DATAN(1.D0)
      NE=40
      FLAMBDA = 1.D0
      H = FLAMBDA / 4.D0
      FKBASE = 2.D0*PI/FLAMBDA
      NFREQUENCY = 1
      FK = NFREQUENCY*FKBASE
C*****************************************************************
      CALL GRULE ( INTEPT , SAI , W )
C*****************************************************************
C------------- FILE MANAGEMENT AND PRINTING HEADER ----------------
      OPEN ( 1, FILE="EXACTSIN.DAT", STATUS='UNKNOWN')
      WRITE (1,*) 'NF-SIN FDPOINT FRACTION Z-REAL Z-IMAG'
C--------------------- CALCULATION OF IMPEDANCE -------------------
      NFEEDING = 40
      FRACTION = 1.D0/NFEEDING
      DO I = 1 , NFEEDING
      FDPOINT = (I-1)*FRACTION*H
      FRAC = (I-1)*FRACTION
      CALL INTE( NE, FDPOINT, FK, H, INTEPT,SAI,W, S )
      WRITE (1,*) NFREQUENCY, FDPOINT, FRAC, DREAL(S),DIMAG(S)
      END DO
      CLOSE (1)
      STOP 'NORMAL TERMINATION SIN'
      END
C
C
      SUBROUTINE INTE( NE, FDPOINT, FK, H, INTEPT,SAI,W, S ) 
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION SAI(INTEPT),W(INTEPT)
      COMPLEX*16  CJ, S, SE, T1, T2, T3
C================ CONSTANTS ======================
      C1 = 2.D0*DCOS(FK*H)
      CJ = DCMPLX(0.D0,1.D0)
      DZ = 2.D0*H/NE
C-------- INTEGRATION STARTS HERE
      S  = DCMPLX(0.D0,0.D0)
      DO IEL = 1 , NE
       Z1 = DZ*(IEL-1) - H
       Z3 = DZ*IEL - H
       Z2 = ( Z1 + Z3 ) / 2.D0
       HS = ( Z3 - Z1 ) / 2.D0
       SE = DCMPLX(0.D0,0D0)
C-------- GAUSS-LEGENDRE INTEGRATION FOR EACH SEGMENT
        DO K = 1 , INTEPT
          Z = HS*SAI(K) + Z2
          R0 = DABS(Z    )
          R1 = DABS(Z-H  )
          R2 = DABS(Z+H  )
          T1 =     CDEXP(-CJ*FK*R1)/R1
          T2 =     CDEXP(-CJ*FK*R2)/R2
          T3 = -C1*CDEXP(-CJ*FK*R0)/R0
          SE = SE + (T1+T2+T3)*DSIN(FK*(H-DABS(Z)))*W(K)
        END DO
          SE = SE*HS
          S = S + SE
      END DO
      S = CJ*30.D0*S/(DSIN(FK*(H+FDPOINT)))**2
      RETURN
      END
C
C
      SUBROUTINE GRULE ( N , SAI , WT )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION SAI(N) , WT(N)
      IF ( N .LE. 0 ) STOP'N<=0'
      IF ( N .GT. 6 ) STOP'N>6'
      IF ( N .EQ. 1 ) THEN
      SAI(1) = 0.D0
      WT(1) = 2.D0
      RETURN
      END IF
      IF ( N .EQ. 2 ) THEN
      SAI(1) = DSQRT(3.D0)/3.D0
      WT(1) = 1.D0
      END IF
      IF ( N .EQ. 3 ) THEN
      SAI(1) = DSQRT(3.D0/5.D0)
      SAI(2) = 0.D0
      WT(1) = 5.D0/ 9.D0
      WT(2) = 8.D0/ 9.D0
      END IF
      IF ( N .EQ. 4 ) THEN
      SAI(1) = 0.33998104358485D0
      SAI(2) = 0.86113631159405D0
        WT(1) = 0.65214515486255D0
        WT(2) = 0.34785484513745D0
      END IF
      IF ( N .EQ. 5 ) THEN
      SAI(1) = 0.90617984593866D0
      SAI(2) = 0.53846931010568D0
      SAI(3) = 0.D0
        WT(1) = 0.23692688505619D0
        WT(2) = 0.47862867049937D0
        WT(3) = 5.12D0 / 9.D0
      END IF
      IF ( N .EQ. 6 ) THEN
      SAI(1) = 0.23861918608319D0
      SAI(2) = 0.66120938646626D0
      SAI(3) = 0.93246951420315D0
        WT(1) = 0.46791393457269D0
        WT(2) = 0.36076157304814D0
        WT(3) = 0.17132449237917D0
      END IF
      NN = N / 2
      DO  I = 1 , NN
      J = N - I + 1
      SAI(J) = - SAI(I)
      WT(J) = WT(I)
      END DO
      RETURN
      END