PROGRAM POLAR1D
C======================================================================
C        AN FEM SOLVER FOR INDUCTANCE POLAR COORDINATE PROBLEM
C                   EQUATION: D/DR(R(DV/DR))+Q=0; 
C      IBTYPE(1)=2 BV(1)=0 NEUMANN BOUNDARY CONDITION AT R=0
C IBTYPE(2)=1 BV(2)=PRESCRIBES DIRECHLET BOUNDARY CONDITION AT R=END
C      *************TRI-DIAGONAL MATRIX SOLVER******************
C                      FEB 11, 2022 EIJI FUKUMORI
C  X(I) = RADIUS OF POLAR COORDINATE
C  CHARGE(I) = CHARGE SUCH AS CURRENT INTO THE CENTER WIRE OF COAXIAL
C  IBTYPE(1) = BOUNDARY CONDITION AT R=0        2 FOR NEUMANN
C  IBTYPE(2) = BOUNDARY CONDITION AT R=RADIUSB  1 FOR DICHILET
C  BV(I) = BOUNDARY VALUES
C  NODEX(I,J) = REPRESENTS ELEMENT CONFIGURATION
C  RHS(I) = RIGHT HAND SIDE OF SIMULTANEOUS LINEAR EQUATIONS
C  A(I,J) = MATRIX OF SIMULTANEOUS LINEAR EQUATIONS
C======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER ( ND=2, MXE=10000, MXN=MXE+1,NBW=2*(ND-1)+1 )
      DIMENSION NODEX(MXE,ND),CHARGE(MXE),X(MXN),A(MXN,NBW),RHS(MXN),
     * IBTYPE(2), BV(2)
C======================================================================
C (1) READING OF DATA
      CALL INPUT ( MXE,MXN,ND,NE,NNODE,NODEX,CHARGE,X,IBTYPE,BV )
C======================================================================
C (2) CONSTRUCTION OF FEM-MATRIX EQUATION
      CALL MATRIX(MXE,MXN,ND,NBW,NE,NNODE,NODEX,CHARGE,X,A,RHS)
C======================================================================
C (3) IMPLEMENTATION OF BOUNDARY CONDITION
      CALL FORM ( MXN, NBW, NNODE, A, RHS, IBTYPE, BV )
C======================================================================
C (4) SOLVE FOR UNKNOWN VARIABLES
      CALL SYSTEM ( MXN, NBW, NNODE, A, RHS )
C======================================================================
C (5) PRINTING RESULTS
      OPEN ( 1,FILE='POLAR1D.SOL', STATUS='UNKNOWN')
      DO I = 1 , NNODE
      WRITE(1,*) I, X(I), RHS(I)
      END DO
      CLOSE (1)
      STOP' NORMAL TERMINATION'
      END
C
C
      SUBROUTINE INPUT (MXE,MXN,ND,NE,NNODE,NODEX,CHARGE,X,IBTYPE,BV )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),CHARGE(MXE),X(MXN),IBTYPE(2), BV(2)
      OPEN ( 1,FILE='POLAR.DAT', STATUS='OLD')
      READ(1,*) NE
      DO I = 1 , NE
      READ(1,*) IEL, (NODEX(IEL,J),J=1,ND),CHARGE(IEL)
      END DO
      NNODE = NE + 1
      DO I = 1 , NNODE
      READ(1,*) NODE, X(NODE)
      END DO
      READ(1,*) IBTYPE(1), BV(1)
      READ(1,*) IBTYPE(2), BV(2)
      CLOSE (1)
      RETURN
      END
C
C
      SUBROUTINE MATRIX ( MXE,MXN,ND,NBW,NE,NNODE,NODEX,
     *                    CHARGE,X,A,RHS )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),CHARGE(MXE),X(MXN),A(MXN,NBW),
     * RHS(MXN)
      DO I = 1 , NNODE
      A(I,1) = 0.D0
      A(I,2) = 0.D0
      A(I,3) = 0.D0
      RHS(I) = 0.D0
      END DO
      DO IEL = 1 , NE
      I = NODEX(IEL,1)
      J = NODEX(IEL,2)
      RL = X(J) - X(I)
      R1 = X(I)
      R2 = X(J)
      RSQ = R2**2 - R1**2
      RQB = R2**3 - R1**3
      C11 = (RQB/3.D0-R2*RSQ/2.D0)/RL**2
      C12 = -C11
      C21 = (R1*RSQ/2.D0-RQB/3.D0)/RL**2
      C22 = -C21
      B11 = RQB/(3.D0*RL**2)
      B12 = -B11
      B21 =  B12
      B22 =  B11
      Q1 = ( R2*RSQ/(2.D0*RL) - RQB/(3.D0*RL))*CHARGE(IEL)
      Q2 = (-R1*RSQ/(2.D0*RL) + RQB/(3.D0*RL))*CHARGE(IEL)
      A(I,2) = A(I,2) + C11 + B11
      A(I,3) = A(I,3) + C12 + B12
      A(J,1) = A(J,1) + C21 + B21
      A(J,2) = A(J,2) + C22 + B22
      RHS(I) = RHS(I) + Q1
      RHS(J) = RHS(J) + Q2
      END DO
      RETURN
      END
C
C
      SUBROUTINE FORM ( MXN, NBW, NNODE, A, RHS, IBTYPE, BV )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION IBTYPE(2),BV(2),A(MXN,NBW), RHS(MXN)
      IF ( IBTYPE(1) .EQ. 1 ) THEN
      A(1,2) = 1.D0
      A(1,3) = 0.D0
      RHS(1) = BV(1)
      RHS(2) = RHS(2) - BV(1)*A(2,1)
      A(2,1) = 0.D0
      ELSE
      RHS(1) = RHS(1) - BV(1)
      END IF
      IF ( IBTYPE(2) .EQ. 1 ) THEN
      A(NNODE,2) = 1.D0
      A(NNODE,1) = 0.D0
      RHS(NNODE) = BV(2)
      RHS(NNODE-1) = RHS(NNODE-1) - BV(2)*A(NNODE-1,3)
      A(NNODE-1,3) = 0.D0
      ELSE
      RHS(NNODE) = RHS(NNODE) - BV(2)
      END IF
      RETURN
      END
C
C
      SUBROUTINE SYSTEM ( MXN, NBW, NNODE, A , B )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION A(MXN,NBW) , B(MXN)
      B(1) = B(1) / A(1,2)
      A(1,2) = A(1,3) / A(1,2)
      DO I = 2 , NNODE
      P = A(I,2) - A(I,1) * A(I-1,2)
      A(I,2) = A(I,3) / P
      B(I) = ( B(I) - A(I,1)*B(I-1) ) / P
      END DO
C------ BACK SUBSTITUTION ----
      DO I = NNODE-1, 1,-1
      B(I) = B(I) - A(I,2) * B(I+1)
      END DO
      RETURN
      END