!
! This file is part of SACAMOS, State of the Art CAble MOdels for Spice. 
! It was developed by the University of Nottingham and the Netherlands Aerospace 
! Centre (NLR) for ESA under contract number 4000112765/14/NL/HK.
! 
! Copyright (C) 2016-2018 University of Nottingham
! 
! SACAMOS is free software: you can redistribute it and/or modify it under the 
! terms of the GNU General Public License as published by the Free Software 
! Foundation, either version 3 of the License, or (at your option) any later 
! version.
! 
! SACAMOS is distributed in the hope that it will be useful, but 
! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 
! or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
! for more details.
! 
! A copy of the GNU General Public License version 3 can be found in the 
! file GNU_GPL_v3 in the root or at <http://www.gnu.org/licenses/>.
! 
! SACAMOS uses the EISPACK library (in /SRC/EISPACK). EISPACK is subject to 
! the GNU Lesser General Public License. A copy of the GNU Lesser General Public 
! License version can be found in the file GNU_LGPL in the root of EISPACK 
! (/SRC/EISPACK ) or at <http://www.gnu.org/licenses/>.
! 
! The University of Nottingham can be contacted at: ggiemr@nottingham.ac.uk
!
!
! File Contents:
! SUBROUTINE frequency_domain_MTL_solution
! SUBROUTINE frequency_domain_MTL_solution_V
!
! NAME
!     frequency_domain_MTL_solution
!
! AUTHORS
!     Chris Smartt
!
! DESCRIPTION
!     This subroutine implements the analytic solution for analysis of
!     multi-conductor transmission lines with resistive terminations
!     and voltages soruces at a single frequency. 
!     The subroutine returns the voltage of the specified conductor at 
!     the specified end of the transmission line
!
!     The comments in this file make reference to the project theory manual.
!     
! COMMENTS
!     Revised comments related to the theory docuemnt are incomplete - the theory
!     document needs to be completed.
!
! HISTORY
!
!     started 7/12/2015 CJS: STAGE_1 developments
!     12/1/2016        CJS: Use fortran intrinsic functions for matrix algebra
!     19/1/2016        CJS: Include comments which refer to the project theory document.
!     22/6/2016        CJS: Include incident field excitation. Note that this initial implementation 
!                           does NOT take proper account of shielded cables.
!     28/6/2016        CJS: Include a ground plane in the incident field excitation 
!     15/7/2016        CJS: Start to include solution for incident field excitation of shielded cables
!     7/3/2017         CJS: Add resistance and voltage source onto the reference coonductor 
!     8/5/2017         CJS: Include references to Theory_Manual
!
SUBROUTINE frequency_domain_MTL_solution(dim,Z_domain,Y_domain,MV,MVI,MI,MII, &
                                         Eamplitude,Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz,xcoord,ycoord,  &
                                         ground_plane_present,ground_plane_x,ground_plane_y,ground_plane_theta, & 
                                         length,Vs1,Z1,Vs2,Z2,is_shielded,f, &
                                         output_end,output_conductor,output_conductor_ref,Vout)

USE type_specifications
USE general_module
USE constants
USE cable_module
USE cable_bundle_module
USE spice_cable_bundle_module
USE maths

IMPLICIT NONE

! variables passed to the subroutine

integer,intent(IN)         :: dim                ! dimension of matrix system

complex(dp),intent(IN)     :: Z_domain(dim,dim)  ! domain based impedance matrix
complex(dp),intent(IN)     :: Y_domain(dim,dim)  ! domain based admittance matrix
complex(dp),intent(IN)     :: MV(dim,dim)        ! domain voltage decomposition matrix
complex(dp),intent(IN)     :: MVI(dim,dim)       ! inverse domain voltage decomposition matrix
complex(dp),intent(IN)     :: MI(dim,dim)        ! domain current decomposition matrix
complex(dp),intent(IN)     :: MII(dim,dim)       ! inverse domain current decomposition matrix

complex(dp),intent(IN) :: Eamplitude                        ! incident field amplitude
real(dp),intent(IN)    :: Ex                                ! Ex component of incident field
real(dp),intent(IN)    :: Ey                                ! Ey component of incident field
real(dp),intent(IN)    :: Ez                                ! Ez component of incident field
real(dp),intent(IN)    :: Hx                                ! Hx component of incident field
real(dp),intent(IN)    :: Hy                                ! Hy component of incident field
real(dp),intent(IN)    :: Hz                                ! Hz component of incident field
real(dp),intent(IN)    :: kx                                ! x component of incident field propagation vector
real(dp),intent(IN)    :: ky                                ! y component of incident field propagation vector
real(dp),intent(IN)    :: kz                                ! z component of incident field propagation vector

real(dp),intent(IN)    :: xcoord(dim+1)                     ! list of conductor x coordinates in bundle cross section
real(dp),intent(IN)    :: ycoord(dim+1)                     ! list of conductor x coordinates in bundle cross section

logical,intent(IN)     :: ground_plane_present              ! flag indicating the presence of a ground plane 
real(dp),intent(IN)    :: ground_plane_x,ground_plane_y     ! input: x and y coordinates of a point on the ground plane
real(dp),intent(IN)    :: ground_plane_theta                ! input: angle of the ground plane from the x axis

real(dp),intent(IN)    :: length                            ! length of bundle (m)

complex(dp),intent(IN)     :: Vs1(dim)                    ! list of voltage sources in end 1 of transmission line termination circuit
complex(dp),intent(IN)     :: Z1(dim,dim)                 ! impedance matrix for end 1 of transmission line termination circuit
complex(dp),intent(IN)     :: Vs2(dim)                    ! list of voltage sources in end 2 of transmission line termination circuit
complex(dp),intent(IN)     :: Z2(dim,dim)                 ! impedance matrix for end 2 of transmission line termination circuit

logical,intent(IN)         :: is_shielded(dim+1)            ! flag to indicate shielded conductors (i.e. those not illuminated by the incident field)

real(dp),intent(IN)        :: f                             ! frequency

integer,intent(IN)         :: output_end                    ! end of transmission line for conductor voltage output
integer,intent(IN)         :: output_conductor              ! conductor number for conductor voltage output
integer,intent(IN)         :: output_conductor_ref          ! conductor number for conductor voltage output reference
complex(dp),intent(OUT)    :: Vout                          ! conductor voltage output to be returned

! local variables

complex(dp)     :: Z(dim,dim)     ! glabal based impedance matrix
complex(dp)     :: Y(dim,dim)     ! glabal based admittance matrix

complex(dp)     :: YZ(dim,dim)    ! product of Y and Z matrices

complex(dp)     :: TV(dim,dim)     ! modal decomposition matrix            [Z][Y]=[TV][GAMMA_SQR][TVI]
complex(dp)     :: TVI(dim,dim)    ! inverse modal decomposition matrix 

complex(dp)     :: TI(dim,dim)     ! modal decomposition matrix            [Y][Z]=[TI][GAMMA_SQR][TII]
complex(dp)     :: TII(dim,dim)    ! inverse modal decomposition matrix 

complex(dp)     :: GAMMA_SQR(dim)  ! diagonal matrix elements in YZ/ ZY diagonalisation 

complex(dp)     :: Zm(dim,dim)     ! modal characteristic impedance matrix
complex(dp)     :: Ym(dim,dim)     ! modal characteristic admittance matrix
complex(dp)     :: Zmd(dim)        ! modal characteristic impedance list
complex(dp)     :: Ymd(dim)        ! modal characteristic impedance list

complex(dp)     :: GAMMA_C(dim)    ! complex square root of GAMMA_SQR
real(dp)        :: gamma_r(dim)    ! real part of the complex square root of GAMMA_SQR

complex(dp)     :: ZC(dim,dim)     ! Characteristic impedance matrix
complex(dp)     :: YC(dim,dim)     ! Characteristic admittance matrix

complex(dp)     :: Exp_p_gamma_l(dim,dim)  ! propagation matrix for modes in the +z direction
complex(dp)     :: Exp_m_gamma_l(dim,dim)  ! propagation matrix for modes in the -z direction

! Temporary matrices used in the matrix solution of the transmission line equations with 
! termination conditions appplied.

complex(dp)     :: D(dim,dim)
complex(dp)     :: sqrtDI(dim,dim)

complex(dp)     :: T1(dim,dim)
complex(dp)     :: T2(dim,dim)

complex(dp)     :: MC11(dim,dim)
complex(dp)     :: MC12(dim,dim)
complex(dp)     :: MC21(dim,dim)
complex(dp)     :: MC22(dim,dim)
complex(dp)     :: TC1(dim,dim)
complex(dp)     :: TC2(dim,dim)
complex(dp)     :: TC3(dim,dim)
complex(dp)     :: MC22I(dim,dim)

complex(dp)     :: TM1(dim,dim)

! Temporary vectors used in the matrix solution of the transmission line equations with 
! termination conditions appplied.
complex(dp)     :: VSC(dim)
complex(dp)     :: VLC(dim)
complex(dp)     :: Imp(dim)
complex(dp)     :: Imm(dim)

complex(dp)     :: TV1(dim)
complex(dp)     :: TV2(dim)
complex(dp)     :: TV3(dim)

! Conductor voltages at ends 1 and 2 referred to the reference conductor voltage at that end
complex(dp)     :: Vend1(dim)
complex(dp)     :: Vend2(dim)

complex(dp)     :: Vout_ref  ! voltage on the output reference conductor

! incident field excitation sources

complex(dp)     :: VFT(dim)
complex(dp)     :: IFT(dim)

real(dp)        :: w  ! angular frequency

! loop variables
integer :: row,col
integer :: i

! integer error indicator for the matrix inverse 
integer :: ierr

! START

! angular frequency
  w=2d0*pi*f                 

! calculate the global impedance matrix (Theory_Manual_Eqn 2.4,2.5) from the domain based impedance matrix
! and the domain decomposition matrices using Theory_Manual_Eqn 3.1 and 3.4 and the 
! definition of the voltages/ currents for the domain based impedance matrix 
! (V_domain)=[Z_domain] (I_domain), (V_domain)=[MV] (V_global), (I_domain)=[MI] )I_global) so
! [MV] (V_global)=[Z_domain] (I_domain)=[MI] )I_global)
! thus (V_global)=[MV^-1][Z_domain] [MI] (I_global)

  TM1=MATMUL(MVI,Z_domain)
  Z=MATMUL(TM1,MI)
    

! calculate the global admittance matrix  (Theory_Manual_Eqn 2.4,2.5) from the domain based admittance matrix
! and the domain decomposition matrices using Theory_Manual_Eqn 3.1 and 3.4 and the 
! definition of the voltages/ currents for the domain based admittance matrix 
! (I_domain)=[Y_domain] (V_domain), (V_domain)=[MV] (V_global), (I_domain)=[MI] )I_global) so
! [MI] (I_global)=[Y_domain] (V_domain)=[MV] )V_global)
! thus (I_global)=[MI^-1][Y_domain][MV] (V_global)

  TM1=MATMUL(MII,Y_domain)
  Y=MATMUL(TM1,MV)
  
! perform a modal decomposition on the YZ product Theory_Manual_Eqn 2.33, 2.36, 2.41 (characteristic impedance)

  CALL modal_decomposition_global(dim,Z_domain,Y_domain,MV,MVI,MI,MII,                                      &
                                  Y,Z,TI,TII,TV,TVI,GAMMA_C,GAMMA_SQR,gamma_r,D,sqrtDI,ZC,YC,Zm,Ym,Zmd,Ymd)
  
! We have assembled all the matrices required in the analysis so we can now solve for the termination voltages
  
! calculate diagonal modal propagation matrices for the modes as in Theory_Manual_Eqn 2.34,2.38
  Exp_p_gamma_l(:,:)=(0d0,0d0)
  Exp_m_gamma_l(:,:)=(0d0,0d0)
  do row=1,dim
     Exp_p_gamma_l(row,row)=exp( GAMMA_C(row)*length)
     Exp_m_gamma_l(row,row)=exp(-GAMMA_C(row)*length)
  end do
  
! Calculate the sources due to the incident field excitation on all conductors Theory_Manual_Eqn 2.60

  CALL calculate_lumped_incident_field_sources(xcoord,ycoord,is_shielded,Eamplitude,Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz,          &
                                               ground_plane_present,    &
                                               length,f,VFT,IFT,dim,TI,TII,Y,Z,ZC,YC,GAMMA_C)

! end1 voltage source
  VSC(:)=Vs1(:)
  
! end2 voltage source
! We add the incident field terms  [Z2](IFT) - (VFT)  at the load end. Theory_Manual_Eqn 2.62
  TV1=matmul(Z2,IFT)
  VLC(:)=Vs2(:)+TV1(:)-VFT(:)

! Fill the LHS matrix elements in Theory_Manual_Eqn 2.43, 2.44

! M11
  TC1(:,:)=Zc(:,:)+Z1(:,:)
  MC11=matmul(TC1,TI)

! M12
  TC1(:,:)=Zc(:,:)-Z1(:,:)
  MC12=matmul(TC1,TI)

! M21
  TC1(:,:)=Zc(:,:)-Z2(:,:)
  TC2=matmul(TC1,TI)
  MC21=matmul(TC2,Exp_m_gamma_l)

! M22
  TC1(:,:)=Zc(:,:)+Z2(:,:)
  TC2=matmul(TC1,TI)
  MC22=matmul(TC2,Exp_p_gamma_l)
  
! Solve the matrix system for Imp intially  Theory_Manual_Eqn 2.46
  if(verbose) write(*,*)'Invert MC22'
  ierr=0   ! set ierr=0 on input to matrix inverse to cause the program to stop if we have a singular matrix
  CALL cinvert_Gauss_Jordan(MC22,dim,MC22I,dim,ierr) 
  if(verbose) write(*,*)'Done: invert MC22'
  
  TC1=matmul(MC12,MC22I)  ! note Keep TC1 for use later
  TC2=matmul(TC1,MC21)

  TC3(:,:)=MC11(:,:)-TC2(:,:)
  
  if(verbose) write(*,*)'Invert TC3'
  ierr=0   ! set ierr=0 on input to matrix inverse to cause the program to stop if we have a singular matrix
  CALL cinvert_Gauss_Jordan(TC3,dim,TC2,dim,ierr) 
  if(verbose) write(*,*)'Done: invert TC3'
  
  TV1=matmul(TC1,VLC)
  
  TV2(:)=VSC(:)-TV1(:)
  
  Imp=matmul(TC2,TV2)
  
! next solve for Imm    Theory_Manual_Eqn 2.47
  TV1=matmul(MC21,Imp)
  
  TV2(:)=VLC(:)-TV1(:)
  Imm=matmul(MC22I,TV2)
  
! now solve for the source end voltages  Theory_Manual_Eqn  2.40 at z=0

  TC1=matmul(ZC,TI) 
  
  TV2(:)=Imm(:)+Imp(:)
  Vend1=matmul(TC1,TV2)
  
  
! now solve for the load end voltages  Theory_Manual_Eqn  2.40 at z=L
  TV1=matmul(Exp_m_gamma_l,Imp)
  TV2=matmul(Exp_p_gamma_l,Imm)
  
  TV3(:)=TV1(:)+TV2(:)
  Vend2=matmul(TC1,TV3)
  
! at this point Vend2 may include lumped sources due to the
! incident field excitation so we must remove these ( Theory_Manual_Eqn  2.62  )
  Vend2(:)=Vend2(:)+VFT(:)     
  
! From the conductor voltages, calculate the output voltage requested which may be the 
! voltage between any two conductors.

  Vout_ref=(0d0,0d0)    ! assume the output reference conductor is the transmission line reference conductor for now
  
  if (output_end.EQ.1) then
  
    if (output_conductor_ref.LE.dim) Vout_ref=Vend1(output_conductor_ref)
    Vout=Vend1(output_conductor)-Vout_ref
  else
  
    if (output_conductor_ref.LE.dim) Vout_ref=Vend2(output_conductor_ref)
    Vout=Vend2(output_conductor)-Vout_ref
    
  end if
  
! END OF THE CALCULATION, OPTIONAL OUTPUT FOR CHECKING FOLLOWS...
  if (.NOT.verbose) RETURN
  
  write(*,*)'YZ'

  YZ=matmul(Y,Z)
  CALL write_cmatrix(YZ,dim,0)
  
  write(*,*)'GAMMA_SQR'
  do i=1,dim
    write(*,*)i,GAMMA_SQR(i)
  end do
  write(*,*)'TI'
  CALL write_cmatrix(TI,dim,0)
  
  write(*,*)'Gamma_r'
  do row=1,dim
    write(*,*)row,gamma_r(row)
  end do

  write(*,*)'Mode velocities'
  do row=1,dim
    write(*,*)row,w/(gamma_r(row))
  end do

  write(*,*)'ZC'
  CALL write_cmatrix(ZC,dim,0)
  write(*,*)'YC'
  CALL write_cmatrix(YC,dim,0)
  
  write(*,*)'Z1:'
  CALL write_cmatrix_re(Z1,dim,0)

  write(*,*)'Vs1'
  do row=1,dim
    write(*,*)real(Vs1(row))
  end do
  
  write(*,*)'Z2:'
  CALL write_cmatrix_re(Z2,dim,0)

  write(*,*)'Vs2'
  do row=1,dim
    write(*,*)real(Vs2(row))
  end do
  
  write(*,*)'ZC:'
  CALL write_cmatrix(ZC,dim,0)
  
  write(*,*)'TI:'
  CALL write_cmatrix_re(TI,dim,0)
  
  write(*,*)'TII:'
  CALL write_cmatrix_re(TII,dim,0)
  
  write(*,*)'TV:'
  CALL write_cmatrix_re(TV,dim,0)
  
  write(*,*)'TVI:'
  CALL write_cmatrix_re(TVI,dim,0)

  write(*,*)'gamma'
  do row=1,dim
    write(*,*)GAMMA_C(row)
  end do
  
  write(*,*)'Check modal decomposition YZ=TI GAMMA_SQR TII'
  
  do row=1,dim
    do col=1,dim
      TC1(row,col)=TII(row,col)*GAMMA_C(row)*GAMMA_C(row)
    end do
  end do
  TC2=matmul(TI,TC1)
 
  write(*,*)'YZ'
  
  CALL write_cmatrix(YZ,dim,0)
  write(*,*)'TI GAMMA TII'
  CALL write_cmatrix(TC2,dim,0)
   
  write(*,*)'Zm:'
  CALL write_cmatrix(Zm,dim,0)
   
  write(*,*)'Ym:'
  CALL write_cmatrix(Ym,dim,0)
  
  write(*,*)'transmission line length'
  write(*,*)length,' (m)'
  
  write(*,*)'Mode Transmission Delay'
  do row=1,dim
    write(*,*)abs( length*sqrt(Zmd(row)*Ymd(row)/(-w*w)) )
  end do
  
  write(*,*)'Mode Transmission Delay 2'
  do row=1,dim
    write(*,*)length*GAMMA_C(row)/(j*w)
  end do
  
  write(*,*)'exp j gamma L:'
  CALL write_cmatrix(Exp_p_gamma_l,dim,0)
  
  write(*,*)'exp-j gamma L:'
  CALL write_cmatrix(Exp_m_gamma_l,dim,0)
  
  write(*,*)'VSc:'
  do row=1,dim
    write(*,*)VSC(row)
  end do
  
  write(*,*)'VLc:'
  do row=1,dim
    write(*,*)VLC(row)
  end do
    
  RETURN

END SUBROUTINE frequency_domain_MTL_solution
!
! NAME
!     frequency_domain_MTL_solution_V
!
! AUTHORS
!     Chris Smartt
!
! DESCRIPTION
!     This subroutine implements the analytic solution for analysis of
!     multi-conductor transmission lines with resistive terminations
!     and voltages soruces at a single frequency. 
!     The subroutine returns the voltage of the specified conductor at 
!     the specified end of the transmission line
!
!     The comments in this file make reference to the project theory manual.
!
!     This solution differs from the solution in frequency_domain_MTL_solution
!     in that we diagonalise the ZY product and base the calculations on modal voltages.
!     
! COMMENTS
!     This is used as a check only and is not the default solution for the
!     validation test cases, the theory is not in the Theory Manual however it
!     closely resembles the solution based on modal currents impemented in 
!     SUBROUTINE frequency_domain_MTL_solution
!
! HISTORY
!
!     started 7/12/2015 CJS: STAGE_1 developments
!     12/1/2016        CJS: Use fortran intrinsic functions for matrix algebra
!     19/1/2016        CJS: Include comments which refer to the project theory document.
!     22/6/2016        CJS: Include incident field excitation. Note that this initial implementation 
!                           does NOT take proper account of shielded cables.
!     28/6/2016        CJS: Include a ground plane in the incident field excitation 
!     15/7/2016        CJS: Start to include solution for incident field excitation of shielded cables
!     7/3/2017         CJS: Add resistance and voltage source onto the reference coonductor 
!     21/4/2017        CJS: Adapt the original frequency_domain_MTL_solution to work with modal voltages
!                           as an independent check of the solution
!
SUBROUTINE frequency_domain_MTL_solution_V(dim,Z_domain,Y_domain,MV,MVI,MI,MII, &
                                         Eamplitude,Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz,xcoord,ycoord,  &
                                         ground_plane_present,ground_plane_x,ground_plane_y,ground_plane_theta, & 
                                         length,Vs1,Z1,Vs2,Z2,is_shielded,f, &
                                         output_end,output_conductor,output_conductor_ref,Vout)

USE type_specifications
USE general_module
USE constants
USE cable_module
USE cable_bundle_module
USE spice_cable_bundle_module
USE maths

IMPLICIT NONE

! variables passed to the subroutine

integer,intent(IN)         :: dim                ! dimension of matrix system

complex(dp),intent(IN)     :: Z_domain(dim,dim)  ! domain based impedance matrix
complex(dp),intent(IN)     :: Y_domain(dim,dim)  ! domain based admittance matrix
complex(dp),intent(IN)     :: MV(dim,dim)        ! domain voltage decomposition matrix
complex(dp),intent(IN)     :: MVI(dim,dim)       ! inverse domain voltage decomposition matrix
complex(dp),intent(IN)     :: MI(dim,dim)        ! domain current decomposition matrix
complex(dp),intent(IN)     :: MII(dim,dim)       ! inverse domain current decomposition matrix

complex(dp),intent(IN) :: Eamplitude                        ! incident field amplitude
real(dp),intent(IN)    :: Ex                                ! Ex component of incident field
real(dp),intent(IN)    :: Ey                                ! Ey component of incident field
real(dp),intent(IN)    :: Ez                                ! Ez component of incident field
real(dp),intent(IN)    :: Hx                                ! Hx component of incident field
real(dp),intent(IN)    :: Hy                                ! Hy component of incident field
real(dp),intent(IN)    :: Hz                                ! Hz component of incident field
real(dp),intent(IN)    :: kx                                ! x component of incident field propagation vector
real(dp),intent(IN)    :: ky                                ! y component of incident field propagation vector
real(dp),intent(IN)    :: kz                                ! z component of incident field propagation vector

real(dp),intent(IN)    :: xcoord(dim+1)                     ! list of conductor x coordinates in bundle cross section
real(dp),intent(IN)    :: ycoord(dim+1)                     ! list of conductor x coordinates in bundle cross section

logical,intent(IN)     :: ground_plane_present              ! flag indicating the presence of a ground plane 
real(dp),intent(IN)    :: ground_plane_x,ground_plane_y     ! input: x and y coordinates of a point on the ground plane
real(dp),intent(IN)    :: ground_plane_theta                ! input: angle of the ground plane from the x axis

real(dp),intent(IN)    :: length                            ! length of bundle (m)

complex(dp),intent(IN)     :: Vs1(dim)                    ! list of voltage sources in end 1 of transmission line termination circuit
complex(dp),intent(IN)     :: Z1(dim,dim)                 ! impedance matrix for end 1 of transmission line termination circuit
complex(dp),intent(IN)     :: Vs2(dim)                    ! list of voltage sources in end 2 of transmission line termination circuit
complex(dp),intent(IN)     :: Z2(dim,dim)                 ! impedance matrix for end 2 of transmission line termination circuit

logical,intent(IN)         :: is_shielded(dim+1)            ! flag to indicate shielded conductors (i.e. those not illuminated by the incident field)

real(dp),intent(IN)        :: f                             ! frequency

integer,intent(IN)         :: output_end                    ! end of transmission line for conductor voltage output
integer,intent(IN)         :: output_conductor              ! conductor number for conductor voltage output
integer,intent(IN)         :: output_conductor_ref          ! conductor number for conductor voltage output reference
complex(dp),intent(OUT)    :: Vout                          ! conductor voltage output to be returned

! local variables

complex(dp)     :: Z(dim,dim)     ! glabal based impedance matrix
complex(dp)     :: Y(dim,dim)     ! glabal based admittance matrix

complex(dp)     :: ZY(dim,dim)     ! product of Z and Y matrices

complex(dp)     :: TV(dim,dim)     ! modal decomposition matrix            [Z][Y]=[TV][GAMMA_SQR][TVI]
complex(dp)     :: TVI(dim,dim)    ! inverse modal decomposition matrix 

complex(dp)     :: TI(dim,dim)     ! modal decomposition matrix            [Y][Z]=[TI][GAMMA_SQR][TII]
complex(dp)     :: TII(dim,dim)    ! inverse modal decomposition matrix 

complex(dp)     :: GAMMA_SQR(dim)  ! diagonal matrix elements in YZ/ ZY diagonalisation 

complex(dp)     :: Zm(dim,dim)     ! modal characteristic impedance matrix
complex(dp)     :: Ym(dim,dim)     ! modal characteristic admittance matrix
complex(dp)     :: Zmd(dim)        ! modal characteristic impedance list
complex(dp)     :: Ymd(dim)        ! modal characteristic impedance list

complex(dp)     :: GAMMA_C(dim)    ! complex square root of GAMMA_SQR
real(dp)        :: gamma_r(dim)    ! real part of the complex square root of GAMMA_SQR

complex(dp)     :: ZC(dim,dim)     ! Characteristic impedance matrix
complex(dp)     :: YC(dim,dim)     ! Characteristic admittance matrix

complex(dp)     :: Exp_p_gamma_l(dim,dim)  ! propagation matrix for modes in the +z direction
complex(dp)     :: Exp_m_gamma_l(dim,dim)  ! propagation matrix for modes in the -z direction

! Temporary matrices used in the matrix solution of the transmission line equations with 
! termination conditions appplied.

complex(dp)     :: D(dim,dim)
complex(dp)     :: sqrtDI(dim,dim)

complex(dp)     :: T1(dim,dim)
complex(dp)     :: T2(dim,dim)

complex(dp)     :: MC11(dim,dim)
complex(dp)     :: MC12(dim,dim)
complex(dp)     :: MC21(dim,dim)
complex(dp)     :: MC22(dim,dim)
complex(dp)     :: TC1(dim,dim)
complex(dp)     :: TC2(dim,dim)
complex(dp)     :: TC3(dim,dim)
complex(dp)     :: MC22I(dim,dim)

complex(dp)     :: TM1(dim,dim)

! Temporary vectors used in the matrix solution of the transmission line equations with 
! termination conditions appplied.
complex(dp)     :: VSC(dim)
complex(dp)     :: VLC(dim)
complex(dp)     :: Vmp(dim)
complex(dp)     :: Vmm(dim)

complex(dp)     :: TV1(dim)
complex(dp)     :: TV2(dim)
complex(dp)     :: TV3(dim)

! Conductor voltages at ends 1 and 2 referred to the reference conductor voltage at that end
complex(dp)     :: Vend1(dim)
complex(dp)     :: Vend2(dim)

complex(dp)     :: Vout_ref  ! voltage on the output reference conductor

! incident field excitation sources

complex(dp)     :: VFT(dim)
complex(dp)     :: IFT(dim)

real(dp)        :: w  ! angular frequency

! loop variables
integer :: row,col
integer :: i

! integer error indicator for the matrix inverse 
integer :: ierr

! START

! angular frequency
  w=2d0*pi*f                 

! calculate the global impedance matrix from the domain based impedance matrix
! and the domain decomposition matrices
  TM1=MATMUL(MVI,Z_domain)
  Z=MATMUL(TM1,MI)
    

! calculate the global admittance matrix from the domain based admittance matrix
! and the domain decomposition matrices
  TM1=MATMUL(MII,Y_domain)
  Y=MATMUL(TM1,MV)
  
! perform a modal decomposition on the ZY product

  CALL modal_decomposition_global_ZY(dim,Z_domain,Y_domain,MV,MVI,MI,MII,                                      &
                                  Y,Z,TI,TII,TV,TVI,GAMMA_C,GAMMA_SQR,gamma_r,D,sqrtDI,ZC,YC,Zm,Ym,Zmd,Ymd)
  
! We have assembled all the matrices required in the analysis so we can now solve for the termination voltages
      
! calculate diagonal modal propagation matrices for the modes
  Exp_p_gamma_l(:,:)=(0d0,0d0)
  Exp_m_gamma_l(:,:)=(0d0,0d0)
  do row=1,dim
     Exp_p_gamma_l(row,row)=exp( GAMMA_C(row)*length)
     Exp_m_gamma_l(row,row)=exp(-GAMMA_C(row)*length)
  end do
  
! Calculate the sources due to the incident field excitation on all conductors

  CALL calculate_lumped_incident_field_sources(xcoord,ycoord,is_shielded,Eamplitude,Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz,          &
                                               ground_plane_present,    &
                                               length,f,VFT,IFT,dim,TI,TII,Y,Z,ZC,YC,GAMMA_C)

! end1 voltage source
  VSC(:)=Vs1(:)
  
! end2 voltage source
! We add the incident field terms  [Z2](IFT) - (VFT)  at the load end
  TV1=matmul(Z2,IFT)
  VLC(:)=Vs2(:)+TV1(:)-VFT(:)

! Fill the LHS matrix elements in eqn

! M11
  TC1=matmul(Z1,Yc)
  TC2=matmul(TC1,TV)
  MC11(:,:)=TV(:,:)+TC2(:,:)

! M12
  TC1=matmul(Z1,Yc)
  TC2=matmul(TC1,TV)
  MC12(:,:)=TV(:,:)-TC2(:,:)

! M21
  TC1=matmul(Tv,Exp_m_gamma_l)
  TC2=matmul(Yc,TC1)
  TC3=matmul(Z2,TC2)
  MC21(:,:)=TC1(:,:)-TC3(:,:)

! M22
  TC1=matmul(Tv,Exp_p_gamma_l)
  TC2=matmul(Yc,TC1)
  TC3=matmul(Z2,TC2)
  MC22(:,:)=TC1(:,:)+TC3(:,:)
  
! Solve the matrix system for Vmp intially
  if(verbose) write(*,*)'Invert MC22'
  ierr=0   ! set ierr=0 on input to matrix inverse to cause the program to stop if we have a singular matrix
  CALL cinvert_Gauss_Jordan(MC22,dim,MC22I,dim,ierr) 
  if(verbose) write(*,*)'Done: invert MC22'
  
  TC1=matmul(MC12,MC22I)  ! note Keep TC1 for use later
  TC2=matmul(TC1,MC21)

  TC3(:,:)=MC11(:,:)-TC2(:,:)
  
  if(verbose) write(*,*)'Invert TC3'
  ierr=0   ! set ierr=0 on input to matrix inverse to cause the program to stop if we have a singular matrix
  CALL cinvert_Gauss_Jordan(TC3,dim,TC2,dim,ierr) 
  if(verbose) write(*,*)'Done: invert TC3'
  
  TV1=matmul(TC1,VLC)
  
  TV2(:)=VSC(:)-TV1(:)
  
  Vmp=matmul(TC2,TV2)
  
! next solve for Vmm
  TV1=matmul(MC21,Vmp)
  
  TV2(:)=VLC(:)-TV1(:)
  Vmm=matmul(MC22I,TV2)
  
! now solve for the source end voltages
  
  TV2(:)=Vmm(:)+Vmp(:)
  Vend1=matmul(TV,TV2)
  
  TV1=matmul(Exp_m_gamma_l,Vmp)
  TV2=matmul(Exp_p_gamma_l,Vmm)
  TV3(:)=TV1(:)+TV2(:)
  Vend2=matmul(TV,TV3)
  
! at this point Vend2 may include lumped sources due to the
! incident field excitation so we must remove these
  Vend2(:)=Vend2(:)+VFT(:) 
  
  Vout_ref=(0d0,0d0)    ! assume the output reference conductor is the transmission line reference conductor for now
  
  if (output_end.EQ.1) then
  
    if (output_conductor_ref.LE.dim) Vout_ref=Vend1(output_conductor_ref)
    Vout=Vend1(output_conductor)-Vout_ref
  else
  
    if (output_conductor_ref.LE.dim) Vout_ref=Vend2(output_conductor_ref)
    Vout=Vend2(output_conductor)-Vout_ref
    
  end if
    
  RETURN

END SUBROUTINE frequency_domain_MTL_solution_V