!
! This file is part of SACAMOS, State of the Art CAble MOdels in 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-2017 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 calc_incident_field_components
! SUBROUTINE calculate_incident_field_sources
! SUBROUTINE calculate_ZT_incident_field_sources
! SUBROUTINE calc_incident_field_FD_values
! SUBROUTINE calc_incident_field_FD_values_GP
! SUBROUTINE calculate_lumped_incident_field_sources
!
! NAME
! calc_incident_field_components
!
! AUTHORS
! Chris Smartt
!
! DESCRIPTION
! This subroutine calculates the cartesian components of the incident field
! excitation and the propagation (k) vector from the inputs which are in
! a spherical coordinate system ! #FIGURE_REFERENCE_REQUIRED
!
! COMMENTS
!
!
! HISTORY
!
! started 13/6/2016 CJS: STAGE_6 developments
! 8/5/2017 CJS: Include references to Theory_Manual
!
SUBROUTINE calc_incident_field_components(Ktheta,Kphi,Etheta,Ephi, &
Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz)
USE type_specifications
USE general_module
USE constants
USE eispack
USE maths
IMPLICIT NONE
! variables passed to the subroutine
real(dp),intent(IN) :: ktheta ! theta component of k vector
real(dp),intent(IN) :: kphi ! phi component of k vector
real(dp),intent(IN) :: etheta ! theta component of E field polarisation vector
real(dp),intent(IN) :: ephi ! phi component of E field polarisation vector
! Cartesian field components
real(dp),intent(OUT) :: Ex
real(dp),intent(OUT) :: Ey
real(dp),intent(OUT) :: Ez
real(dp),intent(OUT) :: Hx
real(dp),intent(OUT) :: Hy
real(dp),intent(OUT) :: Hz
! Cartesian propagation vector components
real(dp),intent(OUT) :: kx
real(dp),intent(OUT) :: ky
real(dp),intent(OUT) :: kz
! local variables
real(dp) :: r,theta,phi
real(dp) :: x,y,z
real(dp) :: dxdt,dxdp
real(dp) :: dydt,dydp
real(dp) :: dzdt,dzdp
real(dp) :: vx,vy,vz
real(dp) :: norm
! START
! calculate unit length incident wave vector Theory_Manual_Eqn 3.135
r=1d0
theta=Ktheta
phi=Kphi
kx=r*sin(theta)*cos(phi)
ky=r*sin(theta)*sin(phi)
kz=r*cos(theta)
! Calculate cartesian unit vectors in the theta and phi directions. Theory_Manual_Eqn 3.137, 3.138
dxdt= r*cos(theta)*cos(phi)
dxdp=-r*sin(phi)
dydt= r*cos(theta)*sin(phi)
dydp= r*cos(phi)
dzdt=-r*sin(theta)
dzdp=0d0
! Calculate the cartesian E field vector from the sum of Etheta and Ephi contributions. Theory_Manual_Eqn 3.139
vx=Etheta*dxdt+Ephi*dxdp
vy=Etheta*dydt+Ephi*dydp
vz=Etheta*dzdt+Ephi*dzdp
! Normalise the E field vector to a length of 1
norm=sqrt(vx*vx+vy*vy+vz*vz)
if (norm.eq.0d0) then
if (verbose) then
write(*,*)'Error calculating incident field data'
write(*,8000)'E vector(theta)=',dxdt,dydt,dzdt
write(*,8000)'E vector(phi)=',dxdp,dydp,dzdp
write(*,8000)'E vector=',Etheta,Ephi
end if
run_status='ERROR in calc_incident_field_components: Zero length incident E field vector'
CALL write_program_status()
STOP 1
end if
Ex=vx/norm
Ey=vy/norm
Ez=vz/norm
! H = (K x E)/Z0. Theory_Manual_Eqn 3.142
Hx=(ky*Ez-kz*Ey)/Z0
Hy=(kz*Ex-kx*Ez)/Z0
Hz=(kx*Ey-ky*Ex)/Z0
if (verbose) then
write(*,*)'Incident_field_components:'
write(*,8000)'K=',kx,ky,kz
write(*,8000)'E=',Ex,Ey,Ez
write(*,8000)'H=',Hx,Hy,Hz
8000 format(A,3E16.6)
end if
RETURN
END SUBROUTINE calc_incident_field_components
!
! NAME
! calculate_incident_field_sources
!
! DESCRIPTION
! calculate constants required to determine the controlled sources
! for Spice models of transient incident field illumination using the
! method of characterisitcs based solution
!
!
! COMMENTS
! Refer to Theory_Manual_Section 3.7
! Also discussed in more depth in C.R. Paul, "Analysis of Multiconductor Transmission Lines"
! 1st edition, section 7.3.2
!
!
! HISTORY
!
! started 16/06/2016 CJS
! 28/6/2016 CJS :Include the presence of a ground plane if required.
! 8/5/2017 CJS: Include references to Theory_Manual
!
!
SUBROUTINE calculate_incident_field_sources(xcoord,ycoord,Eamplitude,Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz, &
ground_plane_present,ground_plane_x,ground_plane_y,ground_plane_theta, &
length,n_modes,TI,Tdi,Tdz,alpha_0,alpha_L)
USE type_specifications
USE constants
USE general_module
IMPLICIT NONE
! variables passed to subroutine
integer,intent(IN) :: n_modes ! input: number of modes in the illuminated domain
real(dp),intent(IN) :: xcoord(n_modes+1),ycoord(n_modes+1) ! input: centre coordinates of conductors
real(dp),intent(IN) :: Eamplitude,Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz ! input: incident field specification
logical ,intent(IN) :: ground_plane_present ! input: 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 ! input: length of transmission line
real(dp),intent(IN) :: TI(n_modes,n_modes) ! input: mode decomposition matrix
real(dp),intent(IN) :: Tdi(n_modes) ! input: mode delays
real(dp),intent(OUT) :: alpha_0(n_modes),alpha_L(n_modes) ! output: constants for controlled sources at z=0 and z=L
real(dp),intent(OUT) :: Tdz ! output: incident field delay time
! local variables
integer :: n_conductors
real(dp) :: xk,yk
real(dp) :: rvx,rvy,rvz
real(dp) :: Tdxyk(n_modes)
real(dp) :: vi
integer :: i,k
! START
n_conductors=n_modes+1
alpha_0(:)=0d0
alpha_L(:)=0d0
if (Eamplitude.eq.0d0) then ! there is no incident field so return zero constants
return
end if
! calculate 1/velocity in each direction (avoids divide by zeros using reciprocal values)
! obtained from Theory_Manual_Eqn 3.136
rvx=kx/c0
rvy=ky/c0
rvz=kz/c0
! Time delay in Z direction
Tdz=length*rvz
! loop over all conductors except the reference conductor
do k=1,n_conductors-1
! Assume that the final (reference) conductor lies at the origin of the xy plane
! and calculate the vector xk,yk from the reference conductor to conductor k.
! If there is a ground plane then calculate the vector to the closest point on the ground plane
! i.e. in the y direction since the ground plane is along the x axis.
if (ground_plane_present) then
xk=0d0
else
xk=xcoord(k)-xcoord(n_modes+1)
end if
yk=ycoord(k)-ycoord(n_modes+1)
! work out the delay propagating from this conductor to the reference
Tdxyk(k)=xk*rvx+yk*rvy
end do
! loop over modes, i
do i=1,n_modes
! loop over all conductors except the reference conductor
do k=1,n_conductors-1
! Again, assume that the final (reference) conductor lies at the origin of the xy plane
! and calculate the vector xk,yk from the reference conductor to conductor k.
! If there is a ground plane then calculate the vector to the closest point on the ground plane
! i.e. in the y direction since the ground plane is along the x axis.
if (ground_plane_present) then
! only the height over the ground plane matters so set x coordinate to zero and y is the height over ground
xk=0d0
else
xk=xcoord(k)-xcoord(n_modes+1)
end if
yk=ycoord(k)-ycoord(n_modes+1)
! note alpha_0,alpha_L are used to calculate the response of mode i due to excitations on all conductors, k
! hence the summation with over k with coefficients TI(k,i)
! Theory_Manual_Eqn 3.167
alpha_0(i)=alpha_0(i)+( Ez*Tdxyk(k)*length - (Ex*xk+Ey*yk)*(Tdi(i)+Tdz) )*TI(k,i) ! note, TI(k,i) = TI(conductor,mode)
alpha_L(i)=alpha_L(i)+( Ez*Tdxyk(k)*length + (Ex*xk+Ey*yk)*(Tdi(i)-Tdz) )*TI(k,i)
end do ! next conductor
end do ! next mode
if (ground_plane_present) then ! See Theory_Manual_Section 3.7.2
! the alpha terms are doubled
alpha_0(:)=2d0*alpha_0(:)
alpha_L(:)=2d0*alpha_L(:)
end if
if (.NOT.verbose) RETURN
write(*,*)'calculate_incident_field_sources'
write(*,*)'ex=',ex
write(*,*)'ey=',ey
write(*,*)'ez=',ez
write(*,*)'kx=',kx
write(*,*)'ky=',ky
write(*,*)'kz=',kz
write(*,*)'rvx=',rvx
write(*,*)'rvy=',rvy
write(*,*)'rvz=',rvz
write(*,*)'Tdz=',Tdz
RETURN
END SUBROUTINE calculate_incident_field_sources
!
! NAME
! calculate_ZT_incident_field_sources
!
! DESCRIPTION
! calculate the controlled sources for Spice modelling of transient incident field illumination
! of shielded cables (Xie's model), method of characterisitcs based solution
!
!
! COMMENTS
! Refer to Theory_Manual_Section 3.8
!
!
! HISTORY
!
! started 3/10/2016 CJS
! 8/5/2017 CJS: Include references to Theory_Manual
!
!
SUBROUTINE calculate_ZT_incident_field_sources(xcoord,ycoord,Eamplitude,Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz, &
ground_plane_present,ground_plane_x,ground_plane_y,ground_plane_theta, &
length,n_modes,shield_conductor,TI,Tdi,Tdz,c1,c2,c1b)
USE type_specifications
USE constants
USE general_module
IMPLICIT NONE
! variables passed to subroutine
integer,intent(IN) :: n_modes ! input: number of modes in the illuminated domain
real(dp),intent(IN) :: xcoord(n_modes+1),ycoord(n_modes+1) ! input: centre coordinates of conductors
real(dp),intent(IN) :: Eamplitude,Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz ! input: incident field specification
logical,intent(IN) :: ground_plane_present ! input: 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 ! input: length of transmission line
integer,intent(IN) :: shield_conductor ! input: the number of the shield conductor
real(dp),intent(IN) :: TI(n_modes,n_modes) ! input: mode decomposition matrix
real(dp),intent(IN) :: Tdi(n_modes) ! input: mode delays
real(dp),intent(OUT) :: c1(n_modes),c2(n_modes),c1b(n_modes) ! output: constants for controlled sources in shielded domains
real(dp),intent(OUT) :: Tdz ! output: incident field delay time
! local variables
integer :: n_conductors
real(dp) :: xk,yk
real(dp) :: rvx,rvy,rvz
real(dp) :: Tx(n_modes),dk,kw_rx,kw_ry,kw_r,vw_r
real(dp) :: vi
real(dp) :: Tsum
! loop variables
integer :: i,k
! START
if(verbose) write(*,*)'CALLED: calculate_ZT_incident_field_sources'
n_conductors=n_modes+1
if(verbose) then
write(*,*)'n_modes=',n_modes
write(*,*)'shield_conductor=',shield_conductor
end if
! Zero the coefficients initially
c1(:)=0d0
c2(:)=0d0
c1b(:)=0d0
if (Eamplitude.eq.0d0) then ! there is no incident field so return zero constants
return
end if
! calculate 1/velocity in each direction (avoids divide by zeros using reciprocal values)
! obtained from Theory_Manual_Eqn 3.136
rvx=kx/c0
rvy=ky/c0
rvz=kz/c0
Tdz=length*rvz
! loop over all conductors except the reference conductor
do k=1,n_conductors-1
! Assume that the final (reference) conductor lies at the origin of the xy plane
! and calculate the vector xk,yk from the reference conductor to conductor k.
! If there is a ground plane then calculate the vector to the closest point on the ground plane
! i.e. in the y direction since the ground plane is along the x axis.
if (ground_plane_present) then
! only the height over the ground plane matters so set x coordinate to zero and y is the height over ground
xk=0d0
else
xk=xcoord(k)-xcoord(n_modes+1)
end if
yk=ycoord(k)-ycoord(n_modes+1)
! distance from this conductor to the reference conductor
dk=sqrt(xk*xk+yk*yk)
! unit vector from this conductor to the reference conductor
kw_rx=xk/dk
kw_ry=yk/dk
! project the k vector along this direction
kw_r=kx*kw_rx+ky*kw_ry
! reciprocal velocity in this direction
vw_r=kw_r/c0
! Tx definition in Xie paper - see also Theory_Manual_Section 3.8 comment in the text below equation 3.198
! where here we have generalised to allow for general conductor offsets in the x-y plane of the bundle
Tx(k)=length*vw_r
end do
! Calculate C1(i), C2(i) , C1b(i) where i is the source domain mode.
! We only have contributions relating to the shield conductor whose transfer impedance we are
! including in the model.
! loop over modes, i in the source domain
do i=1,n_modes
! loop over all conductors except the reference conductor
do k=1,n_conductors-1
! Again, assume that the final (reference) conductor lies at the origin of the xy plane
! and calculate the vector xk,yk from the reference conductor to conductor k.
! If there is a ground plane then calculate the vector to the closest point on the ground plane
! i.e. in the y direction since the ground plane is along the x axis.
if (ground_plane_present) then
xk=0d0
else
xk=xcoord(k)-xcoord(n_modes+1)
end if
yk=ycoord(k)-ycoord(n_modes+1)
dk=sqrt(xk*xk+yk*yk)
! note c1,c2 are used to calculate the response of mode i due to excitations on all conductors, k
! hence the summation with over k with coefficients TI(k,i)
! See the two coefficients for the delayed incident field in Theory_Manual_Eqn 3.195 to 3.198
c1(i)=c1(i)+dk*length*(-Ez*Tx(k) +(Ex*xk/dk+Ey*yk/dk)*(Tdz-Tdi(i)) )*TI(k,i)
c2(i)=c2(i)+dk*length*(-Ez*Tx(k) +(Ex*xk/dk+Ey*yk/dk)*(Tdz+Tdi(i)) )*TI(k,i)
! C1b is for the special case when Tdz=Tdi, see the last paragraph in Theory_Manual_Section 3.8.3
c1b(i)=c1b(i)+dk*length*( Ex*xk/dk+Ey*yk/dk )*TI(k,i)
end do ! next conductor
end do ! next mode
if (ground_plane_present) then ! See Theory_Manual_Section 3.7.2
! the alpha terms are doubled
c1(:)=2d0*c1(:)
c2(:)=2d0*c2(:)
c1b(:)=2d0*c1b(:)
end if
if (.NOT.verbose) RETURN
write(*,*)'calculate_incident_field_sources'
write(*,*)'ex=',ex
write(*,*)'ey=',ey
write(*,*)'ez=',ez
write(*,*)'kx=',kx
write(*,*)'ky=',ky
write(*,*)'kz=',kz
write(*,*)'rvx=',rvx
write(*,*)'rvy=',rvy
write(*,*)'rvz=',rvz
write(*,*)'Tdz=',Tdz
! loop over conductors
do k=1,n_modes
xk=xcoord(k)-xcoord(n_modes+1)
yk=ycoord(k)-ycoord(n_modes+1)
! distance from this conductor to the reference conductor
dk=sqrt(xk*xk+yk*yk)
! unit vector from this conductor to the reference conductor
kw_rx=xk/dk
kw_ry=yk/dk
! project the k vector along this direction
kw_r=kx*kw_rx+ky*kw_ry
write(*,*)'conductor ',k
write(*,*)'xk',xk
write(*,*)'yk',yk
write(*,*)'dk=',dk
write(*,*)'kw_rx=',kw_rx
write(*,*)'kw_ry=',kw_ry
write(*,*)'kw_r=',kw_r
write(*,*)'Tx=',Tx(k)
end do
! loop over modes, i
do i=1,n_modes
write(*,*)'mode ',i
write(*,*)'Tdi=',Tdi(i)
write(*,*)'Tdi(i)+Tdz=',Tdi(i)+Tdz
write(*,*)'c1=',c1(i)
write(*,*)'c2=',c2(i)
write(*,*)'c1b=',c1b(i)
end do
RETURN
END SUBROUTINE calculate_ZT_incident_field_sources
!
! NAME
! calc_incident_field_values
!
! DESCRIPTION
! given coordinates in space and frequency calculate the incident E and H field component values
!
! COMMENTS
!
!
! HISTORY
!
! started 29/04/2015 CJS
! 8/5/2017 CJS: Include references to Theory_Manual
!
!
SUBROUTINE calc_incident_field_FD_values(x,y,z,f,Eamplitude,Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz, &
Ex_value,Ey_value,Ez_value,Hx_value,Hy_value,Hz_value)
USE type_specifications
USE constants
IMPLICIT NONE
complex(dp),intent(IN) :: Eamplitude ! E field amplitude
real(dp),intent(IN) :: Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz ! Cartesian components of field and propgation vector
real(dp),intent(IN) :: x,y,z ! coordinates in 3D space of field evaluation point
real(dp),intent(IN) :: f ! frequency
complex(dp),intent(OUT) :: Ex_value,Ey_value,Ez_value,Hx_value,Hy_value,Hz_value ! output cartesian field values
! local variables
real(dp) :: w ! angular frequency
complex(dp) :: phase ! phase at field evaluation point
! START
! Theory_Manual_Eqn 3.140
w=2d0*pi*f
phase=exp( -j*(w/c0)*( kx*x+ky*y+kz*z ) )
Ex_value=Eamplitude*Ex*phase
Ey_value=Eamplitude*Ey*phase
Ez_value=Eamplitude*Ez*phase
! Theory_Manual_Eqn 3.140 and 3.142
Hx_value=Eamplitude*Hx*phase
Hy_value=Eamplitude*Hy*phase
Hz_value=Eamplitude*Hz*phase
RETURN
END SUBROUTINE calc_incident_field_FD_values
!
! NAME
! calc_incident_field_FD_values_GP
!
! DESCRIPTION
! given coordinates in space and frequency calculate the incident E and H field component values
! with a ground plane present
! The ground plane is assumed to lie along the x axis
!
! COMMENTS
!
!
! HISTORY
!
! started 28/06/2016 CJS
!
!
SUBROUTINE calc_incident_field_FD_values_GP(x,y,z,f,Eamplitude,Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz, &
Ex_value,Ey_value,Ez_value,Hx_value,Hy_value,Hz_value)
USE type_specifications
USE constants
IMPLICIT NONE
complex(dp),intent(IN) :: Eamplitude ! E field amplitude
real(dp),intent(IN) :: Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz ! Cartesian components of field and propgation vector
real(dp),intent(IN) :: x,y,z ! coordinates in 3D space of field evaluation point
real(dp),intent(IN) :: f ! frequency
complex(dp),intent(OUT) :: Ex_value,Ey_value,Ez_value,Hx_value,Hy_value,Hz_value ! output cartesian field values
! local variables
real(dp) :: w ! angular frequency
complex(dp) :: phase ! phase at field evaluation point
! START
w=2d0*pi*f
! Theory_Manual_Eqn 3.140 with the addition of a reflected plane wave from the ground plane (y=0 plane)
! in which ky_reflected=-ky_incident and
! Ex_reflected=-Ex_incident, Ez_reflected=-Ez_incident, Hy_reflected=-Hy_incident
! Incident plane wave
phase=exp( -j*(w/c0)*( kx*x+ky*y+kz*z ) )
Ex_value=Eamplitude*Ex*phase
Ey_value=Eamplitude*Ey*phase
Ez_value=Eamplitude*Ez*phase
Hx_value=Eamplitude*Hx*phase
Hy_value=Eamplitude*Hy*phase
Hz_value=Eamplitude*Hz*phase
! Add reflected plane wave
phase=exp( -j*(w/c0)*(kx*x-ky*y+kz*z ) )
Ex_value=Ex_value-Eamplitude*Ex*phase
Ey_value=Ey_value+Eamplitude*Ey*phase
Ez_value=Ez_value-Eamplitude*Ez*phase
Hx_value=Hx_value+Eamplitude*Hx*phase
Hy_value=Hy_value-Eamplitude*Hy*phase
Hz_value=Hz_value+Eamplitude*Hz*phase
RETURN
END SUBROUTINE calc_incident_field_FD_values_GP
!
! NAME
! calculate_lumped_incident_field_sources
!
! DESCRIPTION
! calculate the lumped sources at z=l due to incident field excitation
!
!
! COMMENTS
! The code is based on Theory_Manual_Section 2.4
! Detailed discussions and derivations are found in
! C.R. Paul, "Analysis of Multiconductor Transmission Lines" 1st edition, section 7.2.1.
!
! The sources are calculated by numerical integration
!
! HISTORY
!
! started 29/04/2015 CJS ESA Spice cable modelling, Phase 1
! 8/5/2017 CJS: Include references to Theory_Manual
!
!
SUBROUTINE 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, &
T,TI,Y,Z,ZC,YC,GAMMA_C)
USE type_specifications
USE constants
USE maths
IMPLICIT NONE
integer,intent(IN) :: dim ! matrix system dimension
real(dp),intent(IN) :: xcoord(dim+1),ycoord(dim+1) ! note that we need the reference conductor position in space
logical,intent(IN) :: is_shielded(dim) ! flag to indicated shielded conductors
complex(dp),intent(IN) :: Eamplitude ! incident field amplitude
real(dp),intent(IN) :: Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz ! incident field cartesian field components and propagation vector
real(dp),intent(IN) :: length ! transmission line length
real(dp),intent(IN) :: f ! frequency
logical,intent(IN) :: ground_plane_present ! flag to indicate the presence of a ground plane
complex(dp),intent(IN) :: T(dim,dim) ! modal decomposition matrix
complex(dp),intent(IN) :: TI(dim,dim) ! modal decomposition matrix (inverse of [T])
complex(dp),intent(IN) :: Y(dim,dim) ! global admittance matrix
complex(dp),intent(IN) :: Z(dim,dim) ! global impedance matrix
complex(dp),intent(IN) :: ZC(dim,dim) ! characteristic impedance matrix
complex(dp),intent(IN) :: YC(dim,dim) ! characteristic admittance matrix
complex(dp),intent(IN) :: GAMMA_C(dim) ! modal propagation constants
complex(dp),intent(OUT) :: VFT(dim),IFT(dim) ! incident field voltage and current sources
! local variables
real(dp) :: xp,yp,zp
real(dp) :: vx,vy,norm
real(dp) :: w,lambda
real(dp) :: dz ! space step for numerical integration
integer :: nz,zloop ! loop variables for integration along the transmission line
! Cartesian field values on a conductor
complex(dp) :: Ex_c,Ey_c,Ez_c
complex(dp) :: Hx_c,Hy_c,Hz_c
! Cartesian field values on the reference conductor
complex(dp) :: Ex_r,Ey_r,Ez_r
complex(dp) :: Hx_r,Hy_r,Hz_r
! Temporary matrices for calculating the chain parameter matrices
complex(dp) :: M1(dim,dim),YI(dim,dim)
complex(dp) :: sinh_gamma_z(dim,dim)
complex(dp) :: cosh_gamma_z(dim,dim)
complex(dp) :: TSTI(dim,dim),TCTI(dim,dim)
! chain parameter matrix elements
complex(dp) :: PHI_11(dim,dim)
complex(dp) :: PHI_12(dim,dim)
complex(dp) :: PHI_21(dim,dim)
complex(dp) :: PHI_22(dim,dim)
!temporary vectors used in the calculations
complex(dp) :: VF1(dim),VF2(dim),VF3(dim),VF(dim),IF1(dim),IF(dim)
complex(dp) :: TV1(dim),TV2(dim),TV3(dim),TV4(dim)
! error code for matrix inverse
integer :: ierr
! loop variables
integer :: conductor
integer :: i
! START
IFT(:)=(0d0,0d0)
VFT(:)=(0d0,0d0)
if (Eamplitude.eq.0d0) then ! there is no incident field so return with zero value source terms
return
end if
! calculate angular frequency(rad/s) and wavelength(m)
w=2d0*pi*f
lambda=c0/f
ierr=0 ! exit on error with the matrix inverse
CALL cinvert_Gauss_Jordan(Y,dim,YI,dim,ierr)
dz=lambda/20d0 ! space step for numerical integration
nz=NINT(length/dz)
! ensure a minimum number of integration points
if (nz.LT.10) nz=10
dz=length/nz
! calculate source integrals
! loop over the transmission line length
do zloop=1,nz
zp=(zloop-1)*dz+dz/2d0 ! z coordinate at the current point. We assume the integrand is uniform over intervals, dz.
! Calculate the chain parameter matrices for the current z value Theory_Manual_Eqns 2.58
sinh_gamma_z(:,:)=(0d0,0d0)
cosh_gamma_z(:,:)=(0d0,0d0)
do i=1,dim
cosh_gamma_z(i,i)=0.5d0*( exp( GAMMA_C(i)*(length-zp) )+exp( -GAMMA_C(i)*(length-zp) ) )
sinh_gamma_z(i,i)=0.5d0*( exp( GAMMA_C(i)*(length-zp) )-exp( -GAMMA_C(i)*(length-zp) ) )
end do
M1=matmul(T,cosh_gamma_z)
TCTI=matmul(M1,TI)
M1=matmul(T,sinh_gamma_z)
TSTI=matmul(M1,TI)
! Phi11
M1=matmul(YI,TCTI)
PHI_11=matmul(M1,Y)
! Phi12
M1=matmul(ZC,TSTI)
PHI_12=-M1
! Phi21
M1=matmul(TSTI,YC)
PHI_21=-M1
! Phi22
PHI_22=TCTI
! loop over conductors (other than the reference)
do conductor=1,dim
! calculate the cartesian field values at the current conductor and the reference conductor
if (.NOT.ground_plane_present) then
! get the field values at the reference conductor i.e. the (dim+1)th conductor
xp=xcoord(dim+1)
yp=ycoord(dim+1)
CALL calc_incident_field_FD_values(xp,yp,zp,f,Eamplitude,Ex,Ey,Ez,Hx,Hy,Hz, &
kx,ky,kz,Ex_r,Ey_r,Ez_r,Hx_r,Hy_r,Hz_r)
! get the field values at the conductor
xp=xcoord(conductor)
yp=ycoord(conductor)
CALL calc_incident_field_FD_values(xp,yp,zp,f,Eamplitude,Ex,Ey,Ez,Hx,Hy,Hz, &
kx,ky,kz,Ex_c,Ey_c,Ez_c,Hx_c,Hy_c,Hz_c)
! calculate a vector from the reference to the conductor
vx=xp-xcoord(dim+1)
vy=yp-ycoord(dim+1)
else
! get the field values at the reference conductor i.e. the (dim+1)th conductor
! under the conductor we are interested in
xp=xcoord(dim+1) ! The 'reference' is at the origin when a ground plane is present
yp=ycoord(dim+1)
CALL calc_incident_field_FD_values_GP(xp,yp,zp,f,Eamplitude,Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz, &
Ex_r,Ey_r,Ez_r,Hx_r,Hy_r,Hz_r)
! get the field values at the conductor
xp=xcoord(conductor)
yp=ycoord(conductor)
CALL calc_incident_field_FD_values_GP(xp,yp,zp,f,Eamplitude,Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz, &
Ex_c,Ey_c,Ez_c,Hx_c,Hy_c,Hz_c)
! calculate a vector from the reference to the conductor
vx=xp-xcoord(dim+1) ! Made the vector consistent with the field evaluation points 16/8/2016
vy=yp-ycoord(dim+1) ! The 'reference' is at the origin when a ground plane is present
end if
! the first part is due to the z component of the E field
VF1(conductor)=(Ez_c-Ez_r) ! C. Paul, equation 7.24(e), Ez contribution ! Theory_Manual_Eqn 2.53(second term) using 2.54
! the second part is due to the z derivative of the incident field integrated from the reference to conductor
! note the field is assumed to be linearly varying between its endpoint values in the integration
VF2(conductor)=j*(w/c0)*kz*(0.5d0*(Ex_c+Ex_r)*vx+0.5d0*(Ey_c+Ey_r)*vy) ! C. Paul, equation 7.24(e), Et contribution, Theory_Manual_Eqn 2.53, first term
! An alternative is the time derivative of the B field integrated from reference to conductor
! VF3(conductor)=j*w*mu0*(-0.5d0*(Hx_c+Hx_r)*vy+0.5d0*(Hy_c+Hy_r)*vx) ! C. Paul, equation 7.24(c), B.dA contribution (H field form)
! the IF source term is due to the incident field integrated from the reference to conductor
! note the field is assumed to be linearly varying between its endpoint values in the integration
! this needs to be pre-multiplied by G+jwC but this happens later...
! C. Paul, equation 7.24(d), Et contribution, Theory_Manual_Eqn 2.55 using a linear interpolation of the field
! between the reference conductor and the conductor k.
IF1(conductor)=-(0.5d0*(Ex_c+Ex_r)*vx+0.5d0*(Ey_c+Ey_r)*vy)
end do ! next conductor
VF(:)=VF1(:)+VF2(:) ! E field form
! VF(:)=VF3(:) ! H field form
! multiply IF1 by Y, the MTL admittance matrix (G+jwC) ! Theory_Manual_Eqn 2.56 transformed to the time domain
IF=matmul(Y,IF1)
! Special process for shielded conductors - set IF to zero
do conductor=1,dim
if (is_shielded(conductor)) then
! set the current sources to zero, note that voltage sources must remain even for shielded conductors.
! Note also that provided shielded cables are centred on their respective shield then this step should
! not be necessary due to the structure of the [Y] matrix and (IF)=[Y](IF1)
IF(conductor)=(0d0,0d0)
end if
end do ! next conductor
! add the contribution to the integration in Theory_Manual_Eqns 2.60
! C. Paul, equation 7.29b,7.29c
TV1=matmul(PHI_11,VF)
TV2=matmul(PHI_12,IF)
TV3=matmul(PHI_21,VF)
TV4=matmul(PHI_22,IF)
! Add contribution to the integral over z assuming constant value of integrand over this section, dz
do conductor=1,dim
VFT(conductor)=VFT(conductor)+dz*(TV1(conductor)+TV2(conductor))
IFT(conductor)=IFT(conductor)+dz*(TV3(conductor)+TV4(conductor))
end do
end do ! next z value
END SUBROUTINE calculate_lumped_incident_field_sources