PROGRAM DIPOLENUM
C***********************************************************************
C----- COMPUTAION OF INPUT IMPEDANCE OF DIPOE ANTENNA THROUGH
C         NUMERICAL INTEGRATION OF VECTOR POTENTIAL AND ELECTRIC FIELD
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 antenna
C             ********* VARIABLES IN THIS PROGRAM *********
C-------- INTEPT = INTEGRATION SAMPLING POINT OF GAUSS-LEGEDRE METHOD
C-------- SAI(INTEPT) = DIMENSION LESS COORDINATES OF THE ABOVE METHOD
C-------- W(INTEPT) = WEIGHTING FACTORS OF THE ABOVE METHOD
C-------- NE = NUMBER OF SEGMENTS
C-------- FLAMBDA = WAVE LENGTH
C-------- H = QUATER WAVE LENGTH
C-------- Z0 = INPUT IMPEDANCE AT FEEDING POINT OF DIPOLE ANTENNA
C-------- FK = 2PI/FLAMBDA
C-------- AZ(I) = VECTOR POTENTIAL
C-------- D2AZ = SECOND DERIVATIVE OF VECTOR POTENTIAL AZ
C-------- EZ = ELECTRICAL FIELD
C-------- CURRENT = EITHER DSIN (FK*(H-DABS(Z))) OR DCOS (FK*Z )
C-------- ZD = INPUT IMPEDANCE AT THE FEEDING POINT OF DIPOLE ANTENNA
C-------- 27 APRIL 2023 EIJI FUKUMORI
C***********************************************************************
      PARAMETER ( NE=100, INTEPT=6 )
      IMPLICIT REAL*8 ( A-H , O-Z )
      COMPLEX*16  AZ(NE), ZD
      DIMENSION SAI(INTEPT),W(INTEPT), OBPOINT(NE)
C---------------------------- BASIC CONSTANTS --------------------------
      PI = 4.D0*DATAN(1.D0)
      CALL GRULE ( INTEPT , SAI , W )
C--------------------------- BASIC PARAMETERS --------------------------
      FLAMBDA = 1.0D0
      H = FLAMBDA / 4.D0
      FKBASE = 2.D0*PI/FLAMBDA
      NFREQUENCY = 2
      FK = NFREQUENCY*FKBASE
C***********************************************************************
      OPEN ( 1, FILE="NUMERICAL.DAT", STATUS='UNKNOWN')
C***********************************************************************
      WRITE (1,*) 'NF-NUMERICAL 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(PI,FDPOINT,NE,OBPOINT,AZ,FK,H,INTEPT,SAI,W,ZD)
      WRITE (1,*) NFREQUENCY, FDPOINT, FRAC, DREAL(ZD),DIMAG(ZD)
      END DO
C***********************************************************************
      CLOSE (1)
C***********************************************************************
      STOP 'NORMAL TERMINATION NUMERICAL'
      END
C
C
      SUBROUTINE INTE(PI,FDPOINT,NE,OBPOINT,AZ,FK,H,INTEPT,SAI,W,ZD)
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION SAI(INTEPT),W(INTEPT),OBPOINT(NE)
      COMPLEX*16  CJ, AZ(NE), D2AZ, EZ,  ZD
      DATA SPL, PERMEA / 299792458.D0, 1.25663706212D-6 /
C***********************************************************************
C-------- COMPUTATION OF AZ BY GAUSS-LEGENDRE INTEGRATION METHOD
C----- SPL = SPEED OF LIGHT
C----- PERMEA = PERMEABILITY IN VACUUM
C----- CJ = UNIT IMAGINARY NUMBER
C----- C2 = PERMEABILITY / (4PI)
C----- OMEGA = angular frequency
C----- ANTLENG = LENGTH OF ANTENNA
C----- DZ = NUMERICAL INTEGRATION SEGMENT
C----- HS = JACOBIAN MATRIX IN COORDINATE TRANSFORMATION
C----- IOB = INDEX OF OBSERVATION POINT
C----- OBPOINT(IOB) = COORDINATE OF OBSERVATION POINT
C----- IEL = INDEX OF INTEGRATION SEGMENT
C----- DZ = WIDTH OF INTEGRATION SEGMENT
C----- Z = COORDINATE OF SOURCE POINT
C----- Z2 = SOURCE POINT COORDINATE OF EACH INTEGRATION SEGMENT
C----- R = LENGTH BETWEEN SOURCE POINT AND OBSERVATION POINT
C-------- ZD = INPUT IMPEDANCE AT THE FEEDING POINT OF DIPOLE ANTENNA
C***********************************************************************
      CJ = DCMPLX(0.D0,1.D0)
      C2 = PERMEA / (4.D0*PI)
      OMEGA = 2.D0*PI*SPL/(4.D0*H)
      ANTLENG = 2.D0*H
      DZ = ANTLENG/NE
      HS = DZ/2.D0
C-------- DZ*IOB-H-HS IS THE CENTER COORDINATE OF EACH SEGMENT
C-------------- Z2 IS ALSO THE CENTER COORDINATE OF EACH SEGMENT
      DO IOB = 1 , NE
           OBPOINT(IOB) = DZ*IOB-H-HS
           AZ(IOB) = DCMPLX(0.D0,0.D0)
         DO IEL = 1 , NE
               Z2 = DZ*IEL-H-HS
          DO K = 1 , INTEPT
           Z = HS*SAI(K) + Z2
           R  = DABS(Z - OBPOINT(IOB) )
           AZ(IOB)=AZ(IOB)+CURRENT(FK,H,Z)*C2*CDEXP(-CJ*FK*R)/R*W(K)*HS
          END DO
         END DO
      END DO
C-------- COMPUTATION OF THE SECOND DERIVATIVE OF AZ BY FDM AND
C-------- COMPUTATION OF E AND ANTENNA IMPEDANCE AT THE FEEDING POINT
      ZD = DCMPLX(0.D0,0.D0)
      DO I = 1 , NE
      EZ = -CJ * OMEGA * ( AZ(I) + D2AZ(I,NE,DZ,AZ)/FK**2 )
      ZD = ZD - EZ*CURRENT(FK,H,OBPOINT(I))*DZ
      END DO
      ZD = ZD/CURRENT(FK,H,FDPOINT)**2
      RETURN
      END
C
C
      FUNCTION CURRENT ( FK, H, Z )
      IMPLICIT REAL*8 ( A-H , O-Z )
CCCCCCCC      CURRENT = DCOS(FK*Z)
      CURRENT = DABS ( DSIN (FK*(H-DABS(Z))) )
      RETURN
      END
C
C
      FUNCTION D2AZ ( I, NE, DZ, AZ )
      IMPLICIT REAL*8 ( A-H , O-Z )
C-------- AZ(NE) = VECTOR POTENTIAL
C-------- D2AZ = SECOND DERIVATIVE OF AZ BY FDM
      COMPLEX*16  AZ(NE), D2AZ
      IF ( I .EQ. 1 ) THEN
      D2AZ=(2.D0*AZ(I)-5.D0*AZ(I+1)+4.D0*AZ(I+2)-AZ(I+3))/DZ**2
      RETURN
      END IF
      IF ( I .EQ. NE ) THEN
      D2AZ=(-AZ(I-3)+4.D0*AZ(I-2)-5.D0*AZ(I-1)+2.D0*AZ(I))/DZ**2
      RETURN
      END IF
      D2AZ = ( AZ(I-1) + AZ(I+1) - 2.D0*AZ(I) ) /DZ**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