!
! 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 .
!
! 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 .
!
! 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