PROGRAM DIPOLEANT
C======================================================================
C-- This computes characteristic impedance(Z) of The Dipole Antenna
C-------------------- FEEDDING POINT: THE CENTER
C---- GAUSS LEGENDRE INTEGRATION METHOD TO EVALUATE COMPLEX POWER
C-------- NE = number of SEGMENTS IN integration interval
C-------- H = Length of LAMBDA/4 OF ANTENA Element
C-------- The integral eq. solving the Z is normalied by wave length
C-------- INTEPT = NUMBER OF SAMPLING POINTS IN EACH INTERVAL.
C-------- SAI    = SAMPLING COORDINATES
C-------- W      = WEIGHTING FACTORS AT THE COORDINATES
C-------- GAUSS LEGENDRE coefficient ARE OBTAINED FROM SUROUTINE GRULE.
C---------------- HALF INTEGATION ( LIMITS: -H TO 0 )----------------
C-------- EIJI FUKUMORI   MAR. 19, 2023
C======================================================================
      PARAMETER ( INTEPT=6 )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION SAI(INTEPT),W(INTEPT)
      COMPLEX*16  ZD
C======================================================================
C--------------------------- BASIC PARAMETERS -------------------------
                        PI = 4.D0*DATAN (1.D0)
                        NE=40
                        FLAMBDA = 1.0D0
                        H = FLAMBDA / 4.D0
                        FKBASE = 2.D0*PI/FLAMBDA
                        NFREQUENCY = 1.D0
                        FK = NFREQUENCY*FKBASE
                        DZ =  2.D0*H/NE
C======================================================================
      WRITE(*,*) "Number of divisions in the integration interval=",NE
      WRITE (*,*) "H (LENGTH OF LAMBDA/4 ANTENA ELEMENT) =", H
C======================================================================
      CALL GRULE ( INTEPT , SAI , W )
      OPEN ( 1, FILE='EXACTCOS.DAT', STATUS='UNKNOWN')
      WRITE (1,*) 'NF-COS FDPOINT FRACTION Z-REAL Z-IMAG'
      NFEEDING = 40
      FRACTION = 1.D0/NFEEDING
      DO I = 1 , NFEEDING
      FDPOINT = (I-1)*FRACTION*H
      FRAC = (I-1)*FRACTION
      CALL INTE( FDPOINT, FK, NE, H, INTEPT,SAI,W, ZD )
      ZD = 2.D0*ZD
      WRITE (1,*) NFREQUENCY, FDPOINT, FRAC, DREAL(ZD),DIMAG(ZD)
      END DO
      CLOSE (1)
C======================================================================
      STOP 'NORMAL TERMINATION COS'
      END
C
C
      SUBROUTINE INTE ( FDPOINT, FK, NE, H, INTEPT,SAI,W, SUM )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION SAI(INTEPT),W(INTEPT)
      COMPLEX*16  CJ, SUM, SE, TERM1, TERM2
C-------- INTEGRATION LIMITS -H TO 0
C-------- TOTAL INTEGRATION SPAN = H
      PI = 4.D0*ATAN(1.D0)
      DZ = H/NE
      CJ = DCMPLX(0.D0,1.D0)
      DSAI = DZ/2.D0
C
      SUM  = DCMPLX(0.D0,0.D0)
      DO IEL = 1 , NE
               Z1 = DZ*(IEL-1) - H
               Z3 = DZ*IEL - H
               Z2 = ( Z1 + Z3 ) / 2.D0
               SE = DCMPLX(0.D0,0D0)
           DO K = 1 , INTEPT
                 Z = DSAI*SAI(K) + Z2
                 R1 = DABS(H+Z)
                 R2 = DABS(H-Z)
                 TERM1 = CDEXP(-CJ*FK*R1)/R1
                 TERM2 = CDEXP(-CJ*FK*R2)/R2
                 SE = SE + (TERM1 + TERM2)*DABS(DCOS(FK*DABS(Z)))*W(K)
           END DO
             SE = SE*DSAI
             SUM = SUM + SE
      END DO
      SUM = CJ*30.D0*SUM/DCOS(FK*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
C--------------------------------
      NN = N / 2
      DO  I = 1 , NN
      J = N - I + 1
      SAI(J) = - SAI(I)
      WT(J) = WT(I)
      END DO
      WRITE(*,*)' GAUSS-LEGENDRE INTEGRATION COEFFICIENTS'
      WRITE(*,*)' Number of sampling points in each interval=',N
      DO I = 1 , N
      WRITE (*,*) 'I,SAI(I),W(I)=',I,SAI(I),WT(I)
      END DO
      RETURN
      END