!
! 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_analysis_FH2
!
! NAME
! frequency_domain_analysis_FH2
!
! AUTHORS
! Chris Smartt
!
! DESCRIPTION
! This subroutine controls the analytic solution for the frequency domain analysis of
! multi-conductor transmission lines.
! The solution is obtained using the full dimension transmission line equations
! i.e. we are NOT using the weak form of transfer impedance coupling
! Note also that frequency dependent quantities are evaluated separately at
! each frequency of analysis, i.e. the frequency dependence of the solution is rigorous
! given only the frequency dependence of the dielectrics is modelled using impedance/ admittance
! matrices whose elements are rational frequency dependent filter functions.
!
! INPUTS:
! spice_bundle_model structure
! spice_validation_test structure
!
! OUTPUT
! analytic frequency domain termination voltage for the specified validation test case written to file
!
! COMMENTS
! STAGE_1: frequency independent parameter solution
! STAGE_2: multi-conductor solution
! STAGE_3: shielded cable solution
! STAGE_4: frequency dependent model
! STAGE_5: transfer impedance model
!
! HISTORY
!
! started 7/12/2015 CJS: STAGE_1 developments
! 24/03/2016 CJS: STAGE_3 developments -shielded cables
! 22/04/2016 CJS: STAGE_4 developments -frequency dependent model
! Include general conductor impedance model 12/05/2016 CJS
! Fix bug with conductor impedance contributions 12/05/2016 CJS
! 25/08/2016 CJS Include revised transfer impedance/ condcutor impedance model for shields
! 8/09/2016 CJS Correct the common mode/ differential mode loss terms for twisted pairs
! 13/10/2016 CJS Correct transfer impedance for multiple modes in external domain
! 7/3/2017 CJS: Add resistance and voltage source onto the reference coonductor
! 8/5/2017 CJS: Include references to Theory_Manual
!
! 26/10/2023 CJS: use FastHenry impedance matrix in analytic solution
!
SUBROUTINE frequency_domain_analysis_FH2(spice_bundle_model,spice_validation_test)
USE type_specifications
USE general_module
USE constants
USE cable_module
USE cable_bundle_module
USE spice_cable_bundle_module
USE maths
USE frequency_spec
IMPLICIT NONE
! variables passed to subroutine
TYPE(spice_model_specification_type),intent(IN):: spice_bundle_model ! Spice cable bundle model structure
TYPE(spice_validation_test_type),intent(IN) :: spice_validation_test ! Spice validation circuit structure
! local variables
real(dp) :: f,w ! frequency and angular frequency
integer :: frequency_loop ! frequency loop variable
integer :: dim ! dimension of matrix system to solve
! domain based impedance and admittance matrices
complex(dp),allocatable :: Z_domain(:,:)
complex(dp),allocatable :: Y_domain(:,:)
! domain based conductor impedance terms
complex(dp),allocatable ::Z_domain_conductor_impedance_correction(:,:)
! Vectors and matrices used in the frequency domain solution of the transmission line equations with termination conditions
complex(dp),allocatable :: Vs1(:)
complex(dp),allocatable :: Z1(:,:)
complex(dp),allocatable :: Vs2(:)
complex(dp),allocatable :: Z2(:,:)
complex(dp) :: Vout ! complex output voltage value
! domain transformation matrices
complex(dp),allocatable :: MI(:,:)
complex(dp),allocatable :: MII(:,:)
complex(dp),allocatable :: MV(:,:)
complex(dp),allocatable :: MVI(:,:)
! temporary working matrices
complex(dp),allocatable :: TM1(:,:)
! temporary variables
integer :: conductor,inner_domain,outer_domain
integer :: domain1,inner_domain1,outer_domain1
integer :: conductor1,reference_conductor1
integer :: domain_conductor1,domain_reference_conductor1
logical :: is_shield1
integer :: domain2,inner_domain2,outer_domain2
integer :: conductor2,reference_conductor2
integer :: domain_conductor2,domain_reference_conductor2
logical :: is_shield2
! conductor based impedance (loss) and transfer impedance model data
complex(dp) :: Zint_c ! conductor surface impedance
complex(dp) :: Zint_c_ref ! reference conductor surface impedance
real(dp) :: Rdc_c ! d.c. resistance of conductor
real(dp) :: Rdc_c_ref ! d.c. resistance of reference conductor
complex(dp) :: Zint_t ! conductor transfer impedance
complex(dp) :: Zint_t_ref ! reference conductor transfer impedance
real(dp) :: Rdc_t ! d.c. resistance of conductor (from transfer impedance)
real(dp) :: Rdc_t_ref ! d.c. resistance of reference conductor (from transfer impedance)
! complex amplitude of incident field
complex(dp) :: Einc
logical,allocatable :: is_shielded_flag(:) ! flags conductors which are not exposed to the incident field
integer :: shield_conductor ! temporary variable, shield conductor number for shielded conductors
real(dp),allocatable :: local_conductor_x_offset(:) ! x coordinate in bundle cross section of conductors
real(dp),allocatable :: local_conductor_y_offset(:) ! y coordinate in bundle cross section of conductors
integer :: n_conductors_outer_domain ! for shield conductors, the number of conductors in the domain outside the shield
integer :: shield_conductor_number_in_outer_domain ! for shield conductors, the conductor number in the domain outside the shield
! loop variables
integer :: row,col
integer :: i
integer :: ierr ! error code for matrix inverse calls
! FastHenry2 stuff...
character(LEN=256) :: line
character(LEN=256) :: freq_and_dim_string
integer :: local_line_length
integer :: n_freq_in
real(dp) :: freq
integer :: nrows,ncols,r,c
complex(dp),allocatable :: Zc_in(:,:,:)
real(dp),allocatable :: freq_in(:)
real(dp),allocatable :: values(:),re,im
integer :: loop,floop
! START
! Open output file
open(unit=analytic_soln_file_unit,file=trim(analytic_soln_filename))
! write the file header line
if (spice_validation_test%analysis_freq_spec%freq_range_type.EQ.'log') then
write(analytic_soln_file_unit,'(A)')log_freq_header
else if (spice_validation_test%analysis_freq_spec%freq_range_type.EQ.'lin') then
write(analytic_soln_file_unit,'(A)')lin_freq_header
end if
dim=spice_bundle_model%bundle%system_dimension
! allocate memory
ALLOCATE( Z_domain(dim,dim) )
ALLOCATE( Y_domain(dim,dim) )
ALLOCATE( Z_domain_conductor_impedance_correction(dim,dim) )
ALLOCATE( Vs1(dim) )
ALLOCATE( Z1(dim,dim) )
ALLOCATE( Vs2(dim) )
ALLOCATE( Z2(dim,dim) )
! domain transformation matrices
ALLOCATE( MI(dim,dim) )
ALLOCATE( MII(dim,dim) )
ALLOCATE( MV(dim,dim) )
ALLOCATE( MVI(dim,dim) )
! temporary working matrices
ALLOCATE( TM1(dim,dim) )
ALLOCATE( is_shielded_flag(1:dim+1) )
ALLOCATE( local_conductor_x_offset(1:dim+1) )
ALLOCATE( local_conductor_y_offset(1:dim+1) )
! loop over conductors to work out which are in shielded domains and which are in the external domain
! also get the position of the conductor in the bundle cross section for incident field excitation
do i=1,dim+1
if (spice_bundle_model%bundle%terminal_conductor_to_outer_domain(i).EQ.spice_bundle_model%bundle%tot_n_domains) then
is_shielded_flag(i)=.FALSE.
local_conductor_x_offset(i)=spice_bundle_model%bundle%conductor_x_offset(i)
local_conductor_y_offset(i)=spice_bundle_model%bundle%conductor_y_offset(i)
else
is_shielded_flag(i)=.TRUE.
! work out the conductor number of the shield
shield_conductor=spice_bundle_model%bundle%terminal_conductor_to_reference_terminal_conductor(i)
! shielded conductors pick up the coordinate of the shield for the purposes of incident field excitation
local_conductor_x_offset(i)=spice_bundle_model%bundle%conductor_x_offset(shield_conductor)
local_conductor_y_offset(i)=spice_bundle_model%bundle%conductor_y_offset(shield_conductor)
end if
end do
! build the termination specifications and convert to complex
Vs1(1:dim)=cmplx( spice_validation_test%Vs_end1(1:dim)-spice_validation_test%Vs_end1(dim+1) )
Vs2(1:dim)=cmplx( spice_validation_test%Vs_end2(1:dim)-spice_validation_test%Vs_end2(dim+1) )
Z1(:,:) =cmplx( spice_validation_test%R_end1(dim+1) )
Z2(:,:) =cmplx( spice_validation_test%R_end2(dim+1) )
do i=1,dim
Z1(i,i) =Z1(i,i)+cmplx( spice_validation_test%R_end1(i) )
Z2(i,i) =Z2(i,i)+cmplx( spice_validation_test%R_end2(i) )
end do
! Copy the domain transformation matrices and calculate the inverses
MI(:,:)=cmplx(spice_bundle_model%bundle%global_MI%mat(:,:))
MV(:,:)=cmplx(spice_bundle_model%bundle%global_MV%mat(:,:))
if (verbose) write(*,*)'Invert MI'
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(MI,dim,MII,dim,ierr)
if (verbose) then
write(*,*)'Transpose[MII]'
do row=1,dim
write(*,8000)(real(MII(col,row)),col=1,dim)
end do
write(*,*)'[MV]'
do row=1,dim
write(*,8000)(real(MV(row,col)),col=1,dim)
end do
8000 format(20F4.1)
end if ! verbose
if (verbose) write(*,*)'Invert MV'
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(MV,dim,MVI,dim,ierr)
if (verbose) then
write(*,*)
write(*,*)' Processing frequency dependent impedance file from FastHenry2'
write(*,*)
end if
OPEN(unit=fh2_output_file_unit,file='Zc.mat')
do loop=1,2
n_freq_in=0
10 CONTINUE
read(fh2_output_file_unit,'(A)',END=1000,ERR=1000)line
if (line(1:32).NE.'Impedance matrix for frequency =') GOTO 10
! if we get here then we have some impedance matrix data to read
! write(*,*)'Found impedance matrix:',trim(line)
! increase the number of frequencies for which we have an impedance matrix
n_freq_in=n_freq_in+1
local_line_length=len(trim(line))
freq_and_dim_string=line(33:local_line_length)
! write(*,*)'reading matrix for:',trim(freq_and_dim_string)
! replace the 'x' by a space then read the frequency, nrows and ncols
local_line_length=len(trim(freq_and_dim_string))
do i=1,local_line_length
if (freq_and_dim_string(i:i).EQ.'x') freq_and_dim_string(i:i)=' '
end do
! write(*,*)'data string:',trim(freq_and_dim_string)
read(freq_and_dim_string,*)freq,nrows,ncols
if (nrows.NE.ncols) then
write(*,*)'****** ERROR: nrows.NE.ncols ******'
end if
! write(*,*)'frequency=',freq,' n=',nrows
if (loop.EQ.2) then
freq_in(n_freq_in)=freq
end if
do r=1,nrows
read(fh2_output_file_unit,'(A)',END=1000,ERR=1000)line
if (loop.EQ.2) then
! replace the 'j's by a space then read the complex data
local_line_length=len(trim(line))
do i=1,local_line_length
if (line(i:i).EQ.'j') line(i:i)=' '
end do
read(line,*)(values(i),i=1,2*ncols)
do c=1,ncols
re=values(c*2-1)
im=values(c*2)
Zc_in(n_freq_in,r,c)=cmplx(re,im)
end do ! next column
end if
end do ! next row
GOTO 10 ! continue to read the file
! Jump here when the file has been read
1000 CONTINUE
if (loop.Eq.1) then
if (verbose) write(*,*)'Number of frequencies=',n_freq_in
ALLOCATE( freq_in(n_freq_in) )
ALLOCATE( Zc_in(n_freq_in,nrows,ncols) )
ALLOCATE( values(2*ncols) )
rewind(unit=fh2_output_file_unit)
end if
end do ! next file read loop
CLOSE(unit=fh2_output_file_unit)
! Check dimensions
if (nrows.NE.dim) then
write(*,*)'ERROR in frequency_domain_analysis_FH2: nrows.NE.dim'
write(*,*)'nrows from FH2 =',nrows,' dim=',dim
write(run_status,*)'ERROR in frequency_domain_analysis_FH2: dimension error:',dim,nrows
CALL write_program_status()
STOP 1
end if
! ************* LOOP OVER FREQUENCY ***************
! Loop over the specified frequencies
do frequency_loop=1,n_freq_in
! get the frequency and angular frequency values
f=freq_in(frequency_loop)
w=2d0*pi*f
! Use the global domain based L and C matrices and the domain voltage and current
! domain transformation matrices to calculate the impedance [Z] and admittance [Y] matrices
do row=1,dim
do col=1,dim
! Evaluate the cable impedance filter function
Z_domain(row,col)=Zc_in(frequency_loop,row,col)
! Evaluate the cable admittance filter function
Y_domain(row,col)=evaluate_Sfilter_frequency_response(spice_bundle_model%bundle%global_Y%sfilter_mat(row,col),f)
end do ! next col
end do ! next row
! Contribution to the impedance matrices from the conductor based impedance models is zero here
Z_domain_conductor_impedance_correction(1:dim,1:dim)=(0d0,0d0)
if (verbose) write(*,*)'Add transfer impedance contributions'
! add transfer impedance contributions *********** NOT SURE WHETHER WE NEED TO DO ANYTHING HERE *********
! loop over conductors looking for shields. Note include all conductors including the reference here
do conductor=1,dim+1
is_shield1=spice_bundle_model%bundle%terminal_conductor_is_shield_flag(conductor)
if (is_shield1) then
! add transfer impedance contributions to inner and outer domain conductors
inner_domain=spice_bundle_model%bundle%terminal_conductor_to_inner_domain(conductor)
outer_domain=spice_bundle_model%bundle%terminal_conductor_to_outer_domain(conductor)
CALL evaluate_conductor_impedance_model(spice_bundle_model%bundle%conductor_impedance(conductor), &
f,Zint_c,Rdc_c,Zint_t,Rdc_t)
! Check whether the shield is the reference conductor in the outer domain - the contributions
! are different if this is the case.
n_conductors_outer_domain=spice_bundle_model%bundle%n_conductors(outer_domain)
shield_conductor_number_in_outer_domain=spice_bundle_model%bundle%terminal_conductor_to_local_domain_conductor(conductor)
! number of conductors in a domain is spice_bundle_model%bundle%n_conductors(domain)
if (shield_conductor_number_in_outer_domain.NE.n_conductors_outer_domain) then
! loop over all conductors
do row=1,dim
! get the domain of row conductor
domain1=spice_bundle_model%bundle%terminal_conductor_to_outer_domain(row)
if (domain1.EQ.inner_domain) then
! The row conductor is in the inner shield domain and so gets a transfer impedance contribution from the shield conductor
! the shield couples these two domains so add the transfer impedance term - also include term to make the matrix symmetric
domain_conductor1=spice_bundle_model%bundle%terminal_conductor_to_global_domain_conductor(row)
domain_conductor2=spice_bundle_model%bundle%terminal_conductor_to_global_domain_conductor(conductor)
Z_domain_conductor_impedance_correction(domain_conductor1,domain_conductor2)= &
Z_domain_conductor_impedance_correction(domain_conductor1,domain_conductor2) -Zint_t
Z_domain_conductor_impedance_correction(domain_conductor2,domain_conductor1)= &
Z_domain_conductor_impedance_correction(domain_conductor2,domain_conductor1) -Zint_t
if (verbose) then
write(*,*)'Shield conductor',conductor,' inner domain',inner_domain,' outer domain',outer_domain
write(*,*)'row',row,' col',col,' row domain',domain1,' col domain',domain2
write(*,*)'Contribution to Zt(',domain_conductor1,domain_conductor2,')'
write(*,*)'Contribution to Zt(',domain_conductor2,domain_conductor1,')'
write(*,*)'Zt conductor:',-Zint_t
end if ! verbose
end if ! transfer impedance term required
end do ! next row conductor
else ! shield IS reference conductor in outer domain
! loop over all conductors
do row=1,dim
! get the domain of row conductor
domain1=spice_bundle_model%bundle%terminal_conductor_to_outer_domain(row)
if (domain1.EQ.inner_domain) then
! The row conductor is in the inner shield domain and so gets a transfer impedance contribution from the shield conductor
! the shield couples these two domains so add the transfer impedance term - also include term to make the matrix symmetric
domain_conductor1=spice_bundle_model%bundle%terminal_conductor_to_global_domain_conductor(row)
! As the shield conductor is the reference we need to find all the conductors contributing to the shield current
! note that the contribution is then -ve of the normal transfer impedance contribution as the currents are in the
! opposite direction
do col=1,dim
domain2=spice_bundle_model%bundle%terminal_conductor_to_outer_domain(col)
! Check the domain of the col conductor. If it is an outer domain conductor of the shield then it contributes
if (domain2.EQ.outer_domain) then
domain_conductor2=spice_bundle_model%bundle%terminal_conductor_to_global_domain_conductor(col)
Z_domain_conductor_impedance_correction(domain_conductor1,domain_conductor2)= &
Z_domain_conductor_impedance_correction(domain_conductor1,domain_conductor2) +Zint_t
Z_domain_conductor_impedance_correction(domain_conductor2,domain_conductor1)= &
Z_domain_conductor_impedance_correction(domain_conductor2,domain_conductor1) +Zint_t
if (verbose) then
write(*,*)'Shield conductor',conductor,' inner domain',inner_domain,' outer domain',outer_domain
write(*,*)'row',row,' col',col,' row domain',domain1,' col domain',domain2
write(*,*)'Contribution to Zt(',domain_conductor1,domain_conductor2,')'
write(*,*)'Contribution to Zt(',domain_conductor2,domain_conductor1,')'
write(*,*)'Zt conductor:',-Zint_t
end if ! verbose
end if ! transfer impedance term required for this col conductor
end do ! next condutor to check
end if ! transfer impedance term required for this row conductor
end do ! next row conductor
end if ! shield is/ is not reference conductor in outer domain
end if ! conductor is a shield
end do ! next conductor
! Add the conductor impedance contributions to the domain based impedance matrix
Z_domain(:,:)=Z_domain(:,:)+Z_domain_conductor_impedance_correction(:,:)
if (verbose) then
write(*,*)'[R_domain]=Re[Z_domain]'
do row=1,dim
write(*,8020)(real(Z_domain(row,col)),col=1,dim)
end do
end if ! verbose
if (verbose) then
write(*,*)'Im[Z_domain]'
do row=1,dim
write(*,8010)(aimag(Z_domain(row,col)),col=1,dim)
end do
write(*,*)'[R_domain]=Re[Z_domain]'
do row=1,dim
write(*,8020)(real(Z_domain(row,col)),col=1,dim)
end do
write(*,*)'Im[Y_domain]'
do row=1,dim
write(*,8010)(aimag(Y_domain(row,col)),col=1,dim)
end do
8010 format(20ES10.2)
8020 format(20F12.4)
end if ! verbose
! Get the incident field amplitude
Einc=cmplx(spice_bundle_model%Eamplitude)
! Solve the frequency domain multi-conductor transmission line equations with the specified termination circuit and
! incident field excitation, return the required conductor voltage in Vout.
if (.NOT.run_validation_test_Vbased) then
CALL frequency_domain_MTL_solution(dim,Z_domain,Y_domain,MV,MVI,MI,MII, &
Einc,spice_bundle_model%Ex,spice_bundle_model%Ey,spice_bundle_model%Ez, &
spice_bundle_model%Hx,spice_bundle_model%Hy,spice_bundle_model%Hz, &
spice_bundle_model%kx,spice_bundle_model%ky,spice_bundle_model%kz, &
local_conductor_x_offset, &
local_conductor_y_offset, &
spice_bundle_model%bundle%ground_plane_present, &
spice_bundle_model%bundle%ground_plane_x, &
spice_bundle_model%bundle%ground_plane_y, &
spice_bundle_model%bundle%ground_plane_theta, &
spice_bundle_model%length,Vs1,Z1,Vs2,Z2, &
is_shielded_flag, &
f,spice_validation_test%output_end,spice_validation_test%output_conductor, &
spice_validation_test%output_conductor_ref,Vout)
else
CALL frequency_domain_MTL_solution_V(dim,Z_domain,Y_domain,MV,MVI,MI,MII, &
Einc,spice_bundle_model%Ex,spice_bundle_model%Ey,spice_bundle_model%Ez, &
spice_bundle_model%Hx,spice_bundle_model%Hy,spice_bundle_model%Hz, &
spice_bundle_model%kx,spice_bundle_model%ky,spice_bundle_model%kz, &
local_conductor_x_offset, &
local_conductor_y_offset, &
spice_bundle_model%bundle%ground_plane_present, &
spice_bundle_model%bundle%ground_plane_x, &
spice_bundle_model%bundle%ground_plane_y, &
spice_bundle_model%bundle%ground_plane_theta, &
spice_bundle_model%length,Vs1,Z1,Vs2,Z2, &
is_shielded_flag, &
f,spice_validation_test%output_end,spice_validation_test%output_conductor, &
spice_validation_test%output_conductor_ref,Vout)
end if
! Output the result to file
if (spice_validation_test%output_type.EQ.'li') then
if (plot_real) then
write(analytic_soln_file_unit,*)f,real(Vout),aimag(Vout)
else
write(analytic_soln_file_unit,*)f,abs(Vout),atan2(aimag(Vout),real(Vout))
end if
else if (spice_validation_test%output_type.EQ.'dB') then
write(analytic_soln_file_unit,*)f,20d0*log10(abs(Vout))
end if ! output format (linear or dB)
end do ! next frequency in frequency loop
! Close output file
Close(unit=analytic_soln_file_unit)
! deallocate memory
DEALLOCATE( Z_domain )
DEALLOCATE( Y_domain )
DEALLOCATE( Z_domain_conductor_impedance_correction )
DEALLOCATE( Vs1 )
DEALLOCATE( Z1 )
DEALLOCATE( Vs2 )
DEALLOCATE( Z2 )
! domain transformation matrices
DEALLOCATE( MI )
DEALLOCATE( MII )
DEALLOCATE( MV )
DEALLOCATE( MVI )
DEALLOCATE( is_shielded_flag )
DEALLOCATE( local_conductor_x_offset )
DEALLOCATE( local_conductor_y_offset )
! temporary working matrices
DEALLOCATE( TM1 )
DEALLOCATE( freq_in )
DEALLOCATE( Zc_in )
DEALLOCATE( values )
RETURN
END SUBROUTINE frequency_domain_analysis_FH2