PROGRAM COSSININTE
C*************************************************************************
C----- INTEGRATION OF (1-COS(T))/T + J*SIN(T)/T FROM 0 TO PI
C                  27 APRIL 2023 EIJI FUKUMORI
C*************************************************************************
      PARAMETER ( INTEPT=6 )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION SAI(INTEPT),W(INTEPT)
      COMPLEX*16 S
C---------------------- BASIC PARAMETERS -------------------------
      PI = 4.D0*DATAN(1.D0)
      NE=40
      CALL GRULE ( INTEPT , SAI , W )
C*****************************************************************
C------------- FILE MANAGEMENT AND PRINTING HEADER ----------------
      OPEN ( 1, FILE="COSSININTE.DAT", STATUS='UNKNOWN')
      WRITE (1,*)'NE INPUT-IMPEDANCE'
C--------------------- CALCULATION OF IMPEDANCE -------------------
      X = 2.D0*PI
      CALL INTE( X, NE, INTEPT,SAI,W, S )
      WRITE (1,*) NE, S
      CLOSE (1)
      STOP 'NORMAL TERMINATION (1-COS(T))/T+JSIN(T)/T'
      END
C
C
      SUBROUTINE INTE( X, NE, INTEPT,SAI,W, S ) 
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION SAI(INTEPT),W(INTEPT)
      COMPLEX*16  CJ, S, SE
C================ CONSTANTS ======================
      CJ = DCMPLX(0.D0,1.D0)
      DT = X/NE
C-------- INTEGRATION STARTS HERE
      S = DCMPLX(0.D0,0.D0)
      DO IEL = 1 , NE
       T1 = DT*(IEL-1)
       T3 = DT*IEL
       T2 = ( T1 + T3 ) / 2.D0
       HS = ( T3 - T1 ) / 2.D0
C-------- GAUSS-LEGENDRE INTEGRATION FOR EACH SEGMENT
       SE = DCMPLX(0.D0,0.D0)
        DO K = 1 , INTEPT
          T = HS*SAI(K) + T2
          SE = SE + 1.D0/T*( 1.D0-DCOS(T) + CJ*DSIN(T)  )*W(K)
        END DO
          SE = SE*HS
          S = S + SE
      END DO
      S = 30.D0*S
      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