! 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 write_transfer_impedance_circuit
!
! NAME
! write_transfer_impedance_circuit
!
! AUTHORS
! Chris Smartt
!
! DESCRIPTION
! This code writes the circuit components required to implement the tranfer impdedance coupling model
! The Spice circuit model is seen in Theory_Manual_Figures 3.12, 3.13 and 3.14
!
! INPUTS REQUIRED
! 1. source domain number of modes, domain based shield conductor number, modal decomposition matrices, mode delays and characteristic impedances
! 2. victim domain number of modes, domain based shield conductor number, modal decomposition matrices and mode delays
! 3. victim domain propagation correction filter functions for each mode
! 4. transfer impedance model
! 5. bundle length
!
! OUTPUTS
! The components required to implement the model are written to
! the subcircuit file
!
! COMMENTS
! Write the circuit elements required for the transfer impedance model
!
! Need to check that we are not including more components than necessary e.g. we have source domain delay terms for each source AND victim mode
!
! HISTORY
!
! STAGE 5 developments started 17/05/2016 CJS. Single mode source and victim only to start with
! 16/6/2016 CJS write the s-domain transfer function using the subroutine write_s_domain_controlled_voltage_source
! 26/3/2016 CJS Fix error for purely resistive transfer impedance functions - the final filter function had order -1.
! 9/5/2017 CJS Document software with reference to Theory_Manual
!
SUBROUTINE write_transfer_impedance_circuit(n_source_domain_modes,n_victim_domain_modes, &
Vv_end1_node,Vv_end2_node,Vv_ref_end1_node,Vv_ref_end2_node, &
Vs_minus_node,Vs_plus_node, &
T_source,Z_source,T_victim, &
source_domain_shield_conductor, &
victim_domain_shield_conductor, &
length, &
TVI_source,TI_source, &
TVI_victim,TI_victim, &
Hpv_filter,ZT_filter, &
next_free_node,vref_node,ZT_model,source_domain,victim_domain )
USE type_specifications
USE general_module
USE constants
USE spice_cable_bundle_module
USE filter_module
USE maths
IMPLICIT NONE
! variables passed to the subroutine
integer,intent(IN) :: n_source_domain_modes ! number of modes in the source domain (and hence array dimensions)
integer,intent(IN) :: n_victim_domain_modes ! number of modes in the source domain (and hence array dimensions)
integer,intent(IN) :: Vv_end1_node(1:n_victim_domain_modes) ! first nodes for victim domain source terms, end 1
integer,intent(IN) :: Vv_end2_node(1:n_victim_domain_modes) ! first nodes for victim domain source terms, end 2
integer,intent(IN) :: Vv_ref_end1_node(1:n_victim_domain_modes) ! second nodes for victim domain source terms, end 1
integer,intent(IN) :: Vv_ref_end2_node(1:n_victim_domain_modes) ! second nodes for victim domain source terms, end 2
integer,intent(IN) :: Vs_minus_node(1:n_source_domain_modes) ! nodes for source domain characteristic variables in -z direction
integer,intent(IN) :: Vs_plus_node(1:n_source_domain_modes) ! nodes for source domain characteristic variables in +z direction
real(dp),intent(IN) :: T_source(1:n_source_domain_modes) ! array of source domain mode propagation delay times
real(dp),intent(IN) :: Z_source(1:n_source_domain_modes) ! array of source domain mode characteristic impedances
real(dp),intent(IN) :: T_victim(1:n_victim_domain_modes) ! array of victim domain mode propagation delay times
real(dp),intent(IN) :: TVI_source(1:n_source_domain_modes,1:n_source_domain_modes) ! Source domain Inverse of Voltage modal decomposition matrix
real(dp),intent(IN) :: TI_source(1:n_source_domain_modes,1:n_source_domain_modes) ! Source domain Current modal decomposition matrix
real(dp),intent(IN) :: TVI_victim(1:n_victim_domain_modes,1:n_victim_domain_modes) ! Victim domain Inverse of Voltage modal decomposition matrix
real(dp),intent(IN) :: TI_victim(1:n_victim_domain_modes,1:n_victim_domain_modes) ! Victim domain Current modal decomposition matrix
integer,intent(IN) :: source_domain_shield_conductor ! shield conductor number in source domain
integer,intent(IN) :: victim_domain_shield_conductor ! shield conductor number in victim domain
real(dp),intent(IN) :: length ! bundle length (m)
TYPE(Sfilter),intent(IN) :: Hpv_filter(1:n_victim_domain_modes) ! Victim domain mode propagation correction filter functions
TYPE(Sfilter),intent(IN) :: ZT_filter ! Transfer impedance filter function for the shield
integer,intent(INOUT) :: next_free_node ! next free spice subcircuit node number
integer,intent(IN) :: vref_node ! sub-circuit reference node number
integer,intent(IN) :: ZT_model ! Transfer impedance model number used for unique naming of components only
integer,intent(IN) :: source_domain ! Source domain number used for unique naming of components only
integer,intent(IN) :: victim_domain ! Victim domain number used for unique naming of components only
! local variables
! loop variables for source and victim modes
integer :: s_mode
integer :: v_mode
real(dp) :: Rlarge=1D10 ! large resistance
! Names for delay lines and associated source and load components
character(len=spice_name_length) :: delay_line_pz_Ts_name(1:n_source_domain_modes,1:n_victim_domain_modes)
character(len=spice_name_length) :: delay_line_pz_Tv_name(1:n_source_domain_modes,1:n_victim_domain_modes)
character(len=spice_name_length) :: delay_line_pz_TsPTv_name(1:n_source_domain_modes,1:n_victim_domain_modes)
character(len=spice_name_length) :: delay_line_mz_Ts_name(1:n_source_domain_modes,1:n_victim_domain_modes)
character(len=spice_name_length) :: delay_line_mz_Tv_name(1:n_source_domain_modes,1:n_victim_domain_modes)
character(len=spice_name_length) :: delay_line_mz_TsPTv_name(1:n_source_domain_modes,1:n_victim_domain_modes)
character(len=spice_name_length) :: delay_line_ZC_pz_Ts_name(1:n_source_domain_modes,1:n_victim_domain_modes)
character(len=spice_name_length) :: delay_line_ZC_pz_Tv_name(1:n_source_domain_modes,1:n_victim_domain_modes)
character(len=spice_name_length) :: delay_line_ZC_pz_TsPTv_name(1:n_source_domain_modes,1:n_victim_domain_modes)
character(len=spice_name_length) :: delay_line_ZC_mz_Ts_name(1:n_source_domain_modes,1:n_victim_domain_modes)
character(len=spice_name_length) :: delay_line_ZC_mz_Tv_name(1:n_source_domain_modes,1:n_victim_domain_modes)
character(len=spice_name_length) :: delay_line_ZC_mz_TsPTv_name(1:n_source_domain_modes,1:n_victim_domain_modes)
character(len=spice_name_length) :: delay_line_E1_pz_Ts_name(1:n_source_domain_modes,1:n_victim_domain_modes)
character(len=spice_name_length) :: delay_line_E1_pz_Tv_name(1:n_source_domain_modes,1:n_victim_domain_modes)
character(len=spice_name_length) :: delay_line_E1_pz_TsPTv_name(1:n_source_domain_modes,1:n_victim_domain_modes)
character(len=spice_name_length) :: delay_line_E1_mz_Ts_name(1:n_source_domain_modes,1:n_victim_domain_modes)
character(len=spice_name_length) :: delay_line_E1_mz_Tv_name(1:n_source_domain_modes,1:n_victim_domain_modes)
character(len=spice_name_length) :: delay_line_E1_mz_TsPTv_name(1:n_source_domain_modes,1:n_victim_domain_modes)
! names for the special case circuit for Tsource=Tvictim
character(len=spice_name_length) :: G_Vplus_derivative_name(1:n_source_domain_modes,1:n_victim_domain_modes)
character(len=spice_name_length) :: G_Vminus_derivative_name(1:n_source_domain_modes,1:n_victim_domain_modes)
character(len=spice_name_length) :: L_Vplus_derivative_name(1:n_source_domain_modes,1:n_victim_domain_modes)
character(len=spice_name_length) :: L_Vminus_derivative_name(1:n_source_domain_modes,1:n_victim_domain_modes)
character(len=spice_name_length) :: E_ZT_l_name(1:n_victim_domain_modes)
character(len=spice_name_length) :: E_ZT_s_name(1:n_victim_domain_modes)
! working strings
character(len=spice_name_length) :: name1,name2,name3
character(len=spice_name_length) :: ZT_string
! circuit to combine transfer impedance terms
character(len=spice_name_length) :: combine_delays_s_E_name(1:n_source_domain_modes,1:n_victim_domain_modes,1:4)
character(len=spice_name_length) :: R_combine_delays_s_name(1:n_victim_domain_modes)
character(len=spice_name_length) :: combine_delays_l_E_name(1:n_source_domain_modes,1:n_victim_domain_modes,1:4)
character(len=spice_name_length) :: R_combine_delays_l_name(1:n_victim_domain_modes)
real(dp) :: Evalue
! Nodes for delay lines
! delay line Ts in +z (Forward) direction
integer :: delay_line_pz_Ts_s_nodes(1:n_source_domain_modes,1:n_victim_domain_modes)
integer :: delay_line_pz_Ts_l_nodes(1:n_source_domain_modes,1:n_victim_domain_modes)
! delay line Tv in +z (Forward) direction
integer :: delay_line_pz_Tv_s_nodes(1:n_source_domain_modes,1:n_victim_domain_modes)
integer :: delay_line_pz_Tv_l_nodes(1:n_source_domain_modes,1:n_victim_domain_modes)
! delay line Ts+Tv in +z (Forward) direction
integer :: delay_line_pz_TsPTv_s_nodes(1:n_source_domain_modes,1:n_victim_domain_modes)
integer :: delay_line_pz_TsPTv_l_nodes(1:n_source_domain_modes,1:n_victim_domain_modes)
! delay line Ts in -z (Backward) direction
integer :: delay_line_mz_Ts_s_nodes(1:n_source_domain_modes,1:n_victim_domain_modes)
integer :: delay_line_mz_Ts_l_nodes(1:n_source_domain_modes,1:n_victim_domain_modes)
! delay line Tv in -z (Backward) direction
integer :: delay_line_mz_Tv_s_nodes(1:n_source_domain_modes,1:n_victim_domain_modes)
integer :: delay_line_mz_Tv_l_nodes(1:n_source_domain_modes,1:n_victim_domain_modes)
! delay line Ts+Tv in -z (Backward) direction
integer :: delay_line_mz_TsPTv_s_nodes(1:n_source_domain_modes,1:n_victim_domain_modes)
integer :: delay_line_mz_TsPTv_l_nodes(1:n_source_domain_modes,1:n_victim_domain_modes)
integer :: last_Hjw_source_s_node(1:n_source_domain_modes,1:n_victim_domain_modes)
integer :: last_Hjw_source_l_node(1:n_source_domain_modes,1:n_victim_domain_modes)
! nodes for the special case circuit for Tsource=Tvictim
integer :: Vplus_derivative_node(1:n_source_domain_modes,1:n_victim_domain_modes)
integer :: Vminus_derivative_node(1:n_source_domain_modes,1:n_victim_domain_modes)
integer :: combine_delays_s_Enode
integer :: combine_delays_l_Enode
! filter variables for the time integration of the Zt filter and the propagation correction filter
TYPE(Sfilter) :: integrator_filter
TYPE(Sfilter) :: integrate_ZT_filter
TYPE(Sfilter) :: Hp_integrate_ZT_filter
TYPE(Sfilter) :: Hp_integrate_ZT_filter_with_cancellation
real(dp) :: ZT_Rdc ! d.c. resistance of shield calculated from the transfer impedance
type(Sfilter) :: ZT_filter_minus_Rdc ! filter function for transfer impedance with d.c. resistance subtracted
real(dp) :: gain ! temporary variable used for filter gain
real(dp) :: Ts_minus_Tv ! difference between source and victim mode velocities
logical :: Ts_equal_Tv ! if source and victim mode velocities are the same then we have a different circuit topology
real(dp) :: TI_source_row(1:n_source_domain_modes) ! row extracted from source domain modal decomposition matrix
real(dp) :: TVI_victim_row(1:n_victim_domain_modes) ! row extracted from victim domain modal decomposition matrix
real(dp) :: PS_PV
integer :: first_combine_Zt_l_node ! nodes used in the part of the subcircuit in which controlled sources are connected in series.
integer :: first_combine_Zt_s_node
integer :: row,i ! temporary loop variables
integer :: aorder,border ! numerator and denominator orders for filter functions
! START
if (verbose) then
write(*,*)'CALLED write_transfer_impedance_circuit'
write(*,*)'TI_source'
CALL dwrite_matrix(TI_source,n_source_domain_modes,n_source_domain_modes,n_source_domain_modes,0)
write(*,*)'TVI_source'
CALL dwrite_matrix(TVI_source,n_source_domain_modes,n_source_domain_modes,n_source_domain_modes,0)
write(*,*)'TI_victim'
CALL dwrite_matrix(TI_victim,n_victim_domain_modes,n_victim_domain_modes,n_victim_domain_modes,0)
write(*,*)'TVI_victim'
CALL dwrite_matrix(TVI_victim,n_victim_domain_modes,n_victim_domain_modes,n_victim_domain_modes,0)
end if
! work out the contribution from each of the modes in the source domain to the shield current
! i.e. get the appropriate row of the source domain TI matrix Theory_Manual_Section 3.7.1
if (source_domain_shield_conductor.LE.n_source_domain_modes) then
! the shield conductor is not the reference conductor in the source domain so pull out the corresponding row of the TI matrix
! Here, TI_source_row(i)=P_s,i in Theory_Manual_Eqn 3.1.24, 3.1.26
TI_source_row(1:n_source_domain_modes)=TI_source(source_domain_shield_conductor,1:n_source_domain_modes)
else
! the shield is the reference conductor so the required row is -(sum all the rows)
! Here, TI_source_row(i)=P_s,i in Theory_Manual_Eqn 3.1.25, 3.1.26
TI_source_row(1:n_source_domain_modes)=0d0
do i=1,n_source_domain_modes
do row=1,n_source_domain_modes
TI_source_row(i)=TI_source_row(i)-TI_source(row,i)
end do
end do
end if
if (verbose) then
write(*,*)'TI_source_row:'
write(*,*)TI_source_row(1:n_source_domain_modes)
end if
! work out the contribution to each of the modes in the victim domain from the transfer impedance source term on the shield conductor
! i.e. get the appropriate row of the victim domain TVI matrix
if (victim_domain_shield_conductor.LE.n_victim_domain_modes) then
! the shield conductor is not the reference conductor in the victim domain so pull out the corresponding col of the TVI matrix
! Here, TVI_victim_row(i)=P_v,i in Theory_Manual_Eqn 3.1.28, 3.1.30
! sign error found 13/10/2016
TVI_victim_row(1:n_victim_domain_modes)=-TVI_victim(1:n_victim_domain_modes,victim_domain_shield_conductor)
else
! The shield is the reference conductor so (sum the columns of the TVI_victim_row matrix) gives the rows of the TVI_victim_row array
! Here, TVI_victim_row(i)=P_v,i in Theory_Manual_Eqn 3.1.29, 3.1.30
TVI_victim_row(1:n_victim_domain_modes)=0d0
do row=1,n_victim_domain_modes
do i=1,n_victim_domain_modes
TVI_victim_row(row)=TVI_victim_row(row)+TVI_victim(row,i)
end do
end do
end if
if (verbose) then
write(*,*)'TVI_victim_row:'
write(*,*)TVI_victim_row(1:n_victim_domain_modes)
end if
CALL write_spice_comment('START OF TRANSFER IMPEDANCE COUPLING MODELS')
! Each mode in the victim domain gets a contribution from each mode in the source domain
! so loop over all victim and source mode combinations
do v_mode=1,n_victim_domain_modes ! loop over victim domain modes
do s_mode=1,n_source_domain_modes ! loop over source domain modes
! create ZT_string which labels the transfer impedance model number plus the source mode and victim mode numbers
name1='ZT'
CALL add_integer_to_string(name1,ZT_model,name2)
name1=trim(name2)//'_sm_'
CALL add_integer_to_string(name1,s_mode,name2)
name1=trim(name2)//'_vm_'
CALL add_integer_to_string(name1,v_mode,ZT_string)
! we always need the T_source(s_mode)+T_victim(v_mode) delay lines so write these components now
! Set TsPTv delay line nodes
! delay line nodes for positive z propagation, source end
CALL create_new_node(delay_line_pz_TsPTv_s_nodes(s_mode,v_mode),next_free_node)
! delay line nodes for positive z propagation, load end
CALL create_new_node(delay_line_pz_TsPTv_l_nodes(s_mode,v_mode),next_free_node)
! delay line nodes for negative z propagation, source end
CALL create_new_node(delay_line_mz_TsPTv_s_nodes(s_mode,v_mode),next_free_node)
! delay line nodes for negative z propagation, load end
CALL create_new_node(delay_line_mz_TsPTv_l_nodes(s_mode,v_mode),next_free_node)
! Set TsPTv delay line component names
delay_line_pz_TsPTv_name(s_mode,v_mode)='T_pz_TsPTv_'//trim(ZT_string)
delay_line_mz_TsPTv_name(s_mode,v_mode)='T_mz_TsPTv_'//trim(ZT_string)
! mode impedance
delay_line_ZC_pz_TsPTv_name(s_mode,v_mode)='RZC_pz_TsPTv_'//trim(ZT_string)
delay_line_ZC_mz_TsPTv_name(s_mode,v_mode)='RZC_mz_TsPTv_'//trim(ZT_string)
! source terms
delay_line_E1_pz_TsPTv_name(s_mode,v_mode)='E1_pz_TsPTv_'//trim(ZT_string)
delay_line_E1_mz_TsPTv_name(s_mode,v_mode)='E1_mz_TsPTv_'//trim(ZT_string)
! Write TsPTv delay lines
CALL write_spice_comment('Delay lines for positive z propagation, T_source(s_mode)+T_victim(v_mode)') ! Theory_Manaul_Eqn 3.132, line 2, term2.
write(spice_model_file_unit,'(A30,4I6,A4,E16.6,A4,E16.6)')delay_line_pz_TsPTv_name(s_mode,v_mode),&
delay_line_pz_TsPTv_s_nodes(s_mode,v_mode),vref_node, &
delay_line_pz_TsPTv_l_nodes(s_mode,v_mode),vref_node, &
' Z0=',Z_source(s_mode),' TD=',T_source(s_mode)+T_victim(v_mode)
CALL write_spice_comment('Delay lines for negative z propagation, T_source(s_mode)+T_victim(v_mode)') ! Theory_Manaul_Eqn 3.131, line 3, term2.
write(spice_model_file_unit,'(A30,4I6,A4,E16.6,A4,E16.6)')delay_line_mz_TsPTv_name(s_mode,v_mode),&
delay_line_mz_TsPTv_s_nodes(s_mode,v_mode),vref_node, &
delay_line_mz_TsPTv_l_nodes(s_mode,v_mode),vref_node, &
' Z0=',Z_source(s_mode),' TD=',T_source(s_mode)+T_victim(v_mode)
! modal impedances on modal delay lines, TsPTv
CALL write_spice_comment('Modal impedances: T_source(s_mode)+T_victim(v_mode)')
write(spice_model_file_unit,'(A30,2I6,E16.6)')delay_line_ZC_pz_TsPTv_name(s_mode,v_mode), &
delay_line_pz_TsPTv_l_nodes(s_mode,v_mode),vref_node,Z_source(s_mode)
CALL write_spice_comment('Modal impedances: T_source(s_mode)+T_victim(v_mode)')
write(spice_model_file_unit,'(A30,2I6,E16.6)')delay_line_ZC_mz_TsPTv_name(s_mode,v_mode), &
delay_line_mz_TsPTv_l_nodes(s_mode,v_mode),vref_node,Z_source(s_mode)
! delay line controlled source for positive z propagation, TsPTv
CALL write_spice_comment('Delay line controlled sources for positive z propagation')
write(spice_model_file_unit,'(A30,4I6,E16.6)')delay_line_E1_pz_TsPTv_name(s_mode,v_mode),&
delay_line_pz_TsPTv_s_nodes(s_mode,v_mode),vref_node, &
Vs_plus_node(s_mode),vref_node,1.0
! delay line controlled source for negative z propagation, TsPTv
CALL write_spice_comment('Delay line controlled sources for negative z propagation')
write(spice_model_file_unit,'(A30,4I6,E16.6)')delay_line_E1_mz_TsPTv_name(s_mode,v_mode) &
,delay_line_mz_TsPTv_s_nodes(s_mode,v_mode),vref_node, &
Vs_minus_node(s_mode),vref_node,1.0
! End of tsptv delay lines
! Set Ts delay line node numbers
! delay line nodes for positive z propagation, source end
CALL create_new_node(delay_line_pz_Ts_s_nodes(s_mode,v_mode),next_free_node)
! delay line nodes for positive z propagation, load end
CALL create_new_node(delay_line_pz_Ts_l_nodes(s_mode,v_mode),next_free_node)
! delay line nodes for negative z propagation, source end
CALL create_new_node(delay_line_mz_Ts_s_nodes(s_mode,v_mode),next_free_node)
! delay line nodes for negative z propagation, load end
CALL create_new_node(delay_line_mz_Ts_l_nodes(s_mode,v_mode),next_free_node)
! Set Ts delay line component names
delay_line_pz_Ts_name(s_mode,v_mode)='T_pz_Ts_'//trim(ZT_string)
delay_line_mz_Ts_name(s_mode,v_mode)='T_mz_Ts_'//trim(ZT_string)
! mode impedance
delay_line_ZC_pz_Ts_name(s_mode,v_mode)='RZC_pz_Ts_'//trim(ZT_string)
delay_line_ZC_mz_Ts_name(s_mode,v_mode)='RZC_mz_Ts_'//trim(ZT_string)
! source terms
delay_line_E1_pz_Ts_name(s_mode,v_mode)='E1_pz_Ts_'//trim(ZT_string)
delay_line_E1_mz_Ts_name(s_mode,v_mode)='E1_mz_Ts_'//trim(ZT_string)
! Write T_source delay line components
CALL write_spice_comment('Delay lines for positive z propagation, T_source(s_mode) delay') ! Theory_Manaul_Eqn 3.131, line 2, term2.
write(spice_model_file_unit,'(A30,4I6,A4,E16.6,A4,E16.6)')delay_line_pz_Ts_name(s_mode,v_mode), &
delay_line_pz_Ts_s_nodes(s_mode,v_mode),vref_node, &
delay_line_pz_Ts_l_nodes(s_mode,v_mode),vref_node, &
' Z0=',Z_source(s_mode),' TD=',T_source(s_mode)
CALL write_spice_comment('Delay lines for negative z propagation, T_source(s_mode) delay') ! Theory_Manaul_Eqn 3.132, line 3, term2.
write(spice_model_file_unit,'(A30,4I6,A4,E16.6,A4,E16.6)')delay_line_mz_Ts_name(s_mode,v_mode), &
delay_line_mz_Ts_s_nodes(s_mode,v_mode),vref_node, &
delay_line_mz_Ts_l_nodes(s_mode,v_mode),vref_node, &
' Z0=',Z_source(s_mode),' TD=',T_source(s_mode)
! modal impedances on modal delay lines, Ts
CALL write_spice_comment('Modal impedances: for pz propagation, T_source(s_mode)')
write(spice_model_file_unit,'(A30,2I6,E16.6)')delay_line_ZC_pz_Ts_name(s_mode,v_mode), &
delay_line_pz_Ts_l_nodes(s_mode,v_mode),vref_node,Z_source(s_mode)
CALL write_spice_comment('Modal impedances: for mz propagation, T_source(s_mode)')
write(spice_model_file_unit,'(A30,2I6,E16.6)')delay_line_ZC_mz_Ts_name(s_mode,v_mode), &
delay_line_mz_Ts_l_nodes(s_mode,v_mode),vref_node,Z_source(s_mode)
! delay line controlled source for positive z propagation, Ts
CALL write_spice_comment('Delay line controlled sources for positive z propagation, Vs+(0,t)')
write(spice_model_file_unit,'(A30,4I6,E16.6)')delay_line_E1_pz_Ts_name(s_mode,v_mode),&
delay_line_pz_Ts_s_nodes(s_mode,v_mode),vref_node, &
Vs_plus_node(s_mode),vref_node,1.0
! delay line controlled source for negative z propagation, Ts
CALL write_spice_comment('Delay line controlled sources for negative z propagation Vs-(L,t)')
write(spice_model_file_unit,'(A30,4I6,E16.6)')delay_line_E1_mz_Ts_name(s_mode,v_mode), &
delay_line_mz_Ts_s_nodes(s_mode,v_mode),vref_node, &
Vs_minus_node(s_mode),vref_node,1.0
! End of T_source(s_mode) delay lines
! Test for the special case when the source and victim mode delays are the same (or very close)
! In this case we need to use a different model to avoid a singularity in the normal model
Ts_minus_Tv=T_source(s_mode)-T_victim(v_mode)
if (abs(Ts_minus_Tv).GT.ZT_min_delay) then
! The whole transfer impedance coupling circuit is implemented using delay lines
! as in equations Theory_Manaul_Eqn 3.131, 3.132
! Set Tv delay line nodes
! delay line nodes for positive z propagation, source end
CALL create_new_node(delay_line_pz_Tv_s_nodes(s_mode,v_mode),next_free_node)
! delay line nodes for positive z propagation, load end
CALL create_new_node(delay_line_pz_Tv_l_nodes(s_mode,v_mode),next_free_node)
! delay line nodes for negative z propagation, source end
CALL create_new_node(delay_line_mz_Tv_s_nodes(s_mode,v_mode),next_free_node)
! delay line nodes for negative z propagation, load end
CALL create_new_node(delay_line_mz_Tv_l_nodes(s_mode,v_mode),next_free_node)
! Set Tv delay line component names
delay_line_pz_Tv_name(s_mode,v_mode)='T_pz_Tv_'//trim(ZT_string)
delay_line_mz_Tv_name(s_mode,v_mode)='T_mz_Tv_'//trim(ZT_string)
! mode impedance
delay_line_ZC_pz_Tv_name(s_mode,v_mode)='RZC_pz_Tv_'//trim(ZT_string)
delay_line_ZC_mz_Tv_name(s_mode,v_mode)='RZC_mz_Tv_'//trim(ZT_string)
! source terms
delay_line_E1_pz_Tv_name(s_mode,v_mode)='E1_pz_Tv_'//trim(ZT_string)
delay_line_E1_mz_Tv_name(s_mode,v_mode)='E1_mz_Tv_'//trim(ZT_string)
! Write Tv delay lines
CALL write_spice_comment('Delay lines for positive z propagation, T_victim(v_mode)') ! Theory_Manaul_Eqn 3.131, line 2, term1
write(spice_model_file_unit,'(A30,4I6,A4,E16.6,A4,E16.6)')delay_line_pz_Tv_name(s_mode,v_mode), &
delay_line_pz_Tv_s_nodes(s_mode,v_mode),vref_node, &
delay_line_pz_Tv_l_nodes(s_mode,v_mode),vref_node, &
' Z0=',Z_source(s_mode),' TD=',T_victim(v_mode)
CALL write_spice_comment('Delay lines for negative z propagation, T_victim(v_mode)') ! Theory_Manaul_Eqn 3.132, line 3, term1
write(spice_model_file_unit,'(A30,4I6,A4,E16.6,A4,E16.6)')delay_line_mz_Tv_name(s_mode,v_mode), &
delay_line_mz_Tv_s_nodes(s_mode,v_mode),vref_node, &
delay_line_mz_Tv_l_nodes(s_mode,v_mode),vref_node, &
' Z0=',Z_source(s_mode),' TD=',T_victim(v_mode)
! modal impedances on modal delay lines, T_victim(v_mode)
CALL write_spice_comment('Modal impedances for pz propagation: T_victim(v_mode)')
write(spice_model_file_unit,'(A30,2I6,E16.6)')delay_line_ZC_pz_Tv_name(s_mode,v_mode), &
delay_line_pz_Tv_l_nodes(s_mode,v_mode),vref_node,Z_source(s_mode)
CALL write_spice_comment('Modal impedances for Mz propagation: T_victim(v_mode)')
write(spice_model_file_unit,'(A30,2I6,E16.6)')delay_line_ZC_mz_Tv_name(s_mode,v_mode), &
delay_line_mz_Tv_l_nodes(s_mode,v_mode),vref_node,Z_source(s_mode)
! delay line controlled source for positive z propagation, T_victim(v_mode)
CALL write_spice_comment('Delay line controlled sources for positive z propagation')
write(spice_model_file_unit,'(A30,4I6,E16.6)')delay_line_E1_pz_Tv_name(s_mode,v_mode) &
,delay_line_pz_Tv_s_nodes(s_mode,v_mode),vref_node, &
Vs_plus_node(s_mode),vref_node,1.0
! delay line controlled source for negative z propagation, T_victim(v_mode)
CALL write_spice_comment('Delay line controlled sources for negative z propagation')
write(spice_model_file_unit,'(A30,4I6,E16.6)')delay_line_E1_mz_Tv_name(s_mode,v_mode) &
,delay_line_mz_Tv_s_nodes(s_mode,v_mode),vref_node, &
Vs_minus_node(s_mode),vref_node,1.0
else
! **** The special case required time derivative circuits operating on the delayed source domain modes****
! as in Theory_Manaul_Section 3.7.2 Theory_Manaul_Eqn 3.133, 3.134
! new node for time derivative of +z travelling wave
CALL create_new_node(Vplus_derivative_node(s_mode,v_mode),next_free_node)
! new node for time derivative of -z travelling wave
CALL create_new_node(Vminus_derivative_node(s_mode,v_mode),next_free_node)
! ****** Names for the special case circuit for Tsource=Tvictim
G_Vplus_derivative_name(s_mode,v_mode)='G_Vp_ddt_'//trim(ZT_string)
G_Vminus_derivative_name(s_mode,v_mode)='G_Vm_ddt_'//trim(ZT_string)
L_Vplus_derivative_name(s_mode,v_mode)='L_Vp_ddt_'//trim(ZT_string)
L_Vminus_derivative_name(s_mode,v_mode)='L_Vm_ddt_'//trim(ZT_string)
! ***** Inductive circuit to calculate the time derivative of Vs+ *****
! See Theory_Manual_Figure 3.14
CALL write_spice_comment('Controlled source for derivative of positive z propagating voltage wave, Vs+(0,t-Ts)')
write(spice_model_file_unit,'(A30,4I6,E16.6)')G_Vplus_derivative_name(s_mode,v_mode), &
Vplus_derivative_node(s_mode,v_mode),vref_node, &
delay_line_pz_Ts_l_nodes(s_mode,v_mode),vref_node,1.0
CALL write_spice_comment('1H inductor for derivative of positive z propagating voltage wave, Vs+(0,t-Ts)')
write(spice_model_file_unit,'(A30,2I6,E16.6)')L_Vplus_derivative_name(s_mode,v_mode), &
Vplus_derivative_node(s_mode,v_mode),vref_node,1.0
! ***** Inductive circuit to calculate the time derivative of Vs- *****
! delay line controlled source for negative z propagation, Ts
! See Theory_Manual_Figure 3.14
CALL write_spice_comment('Controlled sources for derivative of negative z propagating voltage wave Vs-(L,t-Ts)')
write(spice_model_file_unit,'(A30,4I6,E16.6)')G_Vminus_derivative_name(s_mode,v_mode), &
Vminus_derivative_node(s_mode,v_mode),vref_node, &
delay_line_mz_Ts_l_nodes(s_mode,v_mode),vref_node,1.0
CALL write_spice_comment('1H inductor for derivative of negative z propagating voltage wave, Vs-(0,t-Ts)')
write(spice_model_file_unit,'(A30,2I6,E16.6)')L_Vminus_derivative_name(s_mode,v_mode), &
Vminus_derivative_node(s_mode,v_mode),vref_node,1.0
end if ! Special case Ts-Tv =0
end do ! next source mode
! The remaining part of the circuit combines all the contributions to the victim mode voltage source
! We create the nodes for the summation circuit as we go
do s_mode=1,n_source_domain_modes ! loop over source domain modes
! calculate the scaling factor originating from the modal decomposition matrices
! PS_PV is the P_s,iP_v,j term in Theory_Manual_Eqns 3.131, 3.132
PS_PV=TI_source_row(s_mode)*TVI_victim_row(v_mode)
if(verbose) write(*,*)'v_mode=',v_mode,' s_mode=',s_mode,' PS_PV=',PS_PV
! create ZT_string which labels the transfer impedance model number plus the source mode and victim mode numbers
name1='ZT'
CALL add_integer_to_string(name1,ZT_model,name2)
name1=trim(name2)//'_sm_'
CALL add_integer_to_string(name1,s_mode,name2)
name1=trim(name2)//'_vm_'
CALL add_integer_to_string(name1,v_mode,ZT_string)
! first add the contributions which are common to both forms of circuit
! Transfer impedance voltage source names
combine_delays_s_E_name(s_mode,v_mode,1)='E_zt_dsum_s_'//trim(ZT_string)//'_E1'
combine_delays_s_E_name(s_mode,v_mode,2)='E_zt_dsum_s_'//trim(ZT_string)//'_E2'
combine_delays_l_E_name(s_mode,v_mode,1)='E_zt_dsum_l_'//trim(ZT_string)//'_E1'
combine_delays_l_E_name(s_mode,v_mode,2)='E_zt_dsum_l_'//trim(ZT_string)//'_E2'
! START OF CIRCUIT TO COMBINE TRANSFER IMPEDANCE TERMS
if (s_mode.eq.1) then
first_combine_Zt_l_node=vref_node
first_combine_Zt_s_node=vref_node
else
first_combine_Zt_l_node=combine_delays_l_Enode
first_combine_Zt_s_node=combine_delays_s_Enode
end if
! Forward (pz) propagating modes: calculation of V_victim at z=L
CALL write_spice_comment('Circuit to combine transfer impedance terms')
! Vs_minus source, no delay: Theory_Manual_Eqn 3.132, line 2, term1.
Evalue=-length*PS_PV/(2d0*Z_source(s_mode)*(T_source(s_mode)+T_victim(v_mode)))
write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_l_E_name(s_mode,v_mode,1) &
,next_free_node ,first_combine_Zt_l_node &
,Vs_minus_node(s_mode),vref_node &
,Evalue
CALL create_new_node(combine_delays_l_Enode,next_free_node)
! Vs_minus source, delay=T_victim+T_source: Theory_Manual_Eqn 3.132, line 2, term2
Evalue=+length*PS_PV/(2d0*Z_source(s_mode)*(T_source(s_mode)+T_victim(v_mode)))
write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_l_E_name(s_mode,v_mode,2) &
,next_free_node,combine_delays_l_Enode &
,delay_line_mz_TsPTv_l_nodes(s_mode,v_mode),vref_node &
,Evalue
CALL create_new_node(combine_delays_l_Enode,next_free_node)
! Backward (mz) propagating modes: calculation of V_victim at z=0
! Vs_plus source, no delay: Theory_Manual_Eqn 3.131, line 3, term1.
Evalue=-length*PS_PV/(2d0*Z_source(s_mode)*(T_source(s_mode)+T_victim(v_mode)))
write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_s_E_name(s_mode,v_mode,1) &
,next_free_node,first_combine_Zt_s_node &
,Vs_plus_node(s_mode),vref_node &
,Evalue
CALL create_new_node(combine_delays_s_Enode,next_free_node)
! Vs_plus source, delay=T_source+T_victim: Theory_Manual_Eqn 3.131, line 3, term2.
Evalue=+length*PS_PV/(2d0*Z_source(s_mode)*(T_source(s_mode)+T_victim(v_mode)))
write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_s_E_name(s_mode,v_mode,2) &
,next_free_node,combine_delays_s_Enode &
,delay_line_pz_TsPTv_l_nodes(s_mode,v_mode),vref_node &
,Evalue
CALL create_new_node(combine_delays_s_Enode,next_free_node)
! Check for special case when T_source=T_victim
Ts_minus_Tv=T_source(s_mode)-T_victim(v_mode)
if (abs(Ts_minus_Tv).GT.ZT_min_delay) then
! normal firm based on delay lines, Theory_Manual_Eqns 3.131, 3.132
combine_delays_s_E_name(s_mode,v_mode,3)='E_zt_dsum_s_'//trim(ZT_string)//'_E3'
combine_delays_s_E_name(s_mode,v_mode,4)='E_zt_dsum_s_'//trim(ZT_string)//'_E4'
combine_delays_l_E_name(s_mode,v_mode,3)='E_zt_dsum_l_'//trim(ZT_string)//'_E3'
combine_delays_l_E_name(s_mode,v_mode,4)='E_zt_dsum_l_'//trim(ZT_string)//'_E4'
! we need to combiine contributions from the normal delay line circuit
! Vs_plus source, delay=T_victim: Theory_Manual_Eqn 3.131, line 2, term1.
Evalue=+length*PS_PV/(2d0*Z_source(s_mode)*(T_source(s_mode)-T_victim(v_mode)))
write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_l_E_name(s_mode,v_mode,3) &
,next_free_node,combine_delays_l_Enode &
,delay_line_pz_Tv_l_nodes(s_mode,v_mode),vref_node &
,Evalue
CALL create_new_node(combine_delays_l_Enode,next_free_node)
! Vs_plus source, delay=T_source: Theory_Manual_Eqn 3.131, line 2, term2
Evalue=-length*PS_PV/(2d0*Z_source(s_mode)*(T_source(s_mode)-T_victim(v_mode)))
write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_l_E_name(s_mode,v_mode,4) &
,next_free_node,combine_delays_l_Enode &
,delay_line_pz_Ts_l_nodes(s_mode,v_mode),vref_node &
,Evalue
CALL create_new_node(combine_delays_l_Enode,next_free_node)
! Vs_minus source, delay=T_victim: Theory_Manual_Eqn 3.132, line 3, term1
Evalue=+length*PS_PV/(2d0*Z_source(s_mode)*(T_source(s_mode)-T_victim(v_mode)))
write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_s_E_name(s_mode,v_mode,4) &
,next_free_node,combine_delays_s_Enode &
,delay_line_mz_Tv_l_nodes(s_mode,v_mode),vref_node &
,Evalue
CALL create_new_node(combine_delays_s_Enode,next_free_node)
! Vs_minus source, delay=T_source: Theory_Manual_Eqn 3.132, line 3, term2
Evalue=-length*PS_PV/(2d0*Z_source(s_mode)*(T_source(s_mode)-T_victim(v_mode)))
write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_s_E_name(s_mode,v_mode,3) &
,next_free_node,combine_delays_s_Enode &
,delay_line_mz_Ts_l_nodes(s_mode,v_mode),vref_node &
,Evalue
CALL create_new_node(combine_delays_s_Enode,next_free_node)
else
! we need to combine contributions from the delay lines and time derivative circuits
! as in Theory_Manual_Section 3.7.2, Theory_Manual_Eqns 3.133, 3.134
combine_delays_s_E_name(s_mode,v_mode,3)='E_zt_dsum_s_'//trim(ZT_string)//'_E3'
combine_delays_l_E_name(s_mode,v_mode,3)='E_zt_dsum_l_'//trim(ZT_string)//'_E3'
! Forward (pz) propagating modes: calculation of V_victim at z=L
! Vs_plus source, time derivative of Vs+: Theory_Manual_Eqn 3.131 line 2, terms 1 and 2 with with Ts=Tv
Evalue=-length*PS_PV/(2d0*Z_source(s_mode))
write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_l_E_name(s_mode,v_mode,3) &
,next_free_node,combine_delays_l_Enode &
,Vplus_derivative_node(s_mode,v_mode),vref_node &
,Evalue
CALL create_new_node(combine_delays_l_Enode,next_free_node)
! Backward (mz) propagating modes: calculation of V_victim at z=0
! Vs_minus source, time derivative of Vs- Theory_Manual_Eqn 3.132 line 3, terms 1 and 2 with with Ts=Tv
Evalue=-length*PS_PV/(2d0*Z_source(s_mode))
write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_s_E_name(s_mode,v_mode,3) &
,next_free_node,combine_delays_s_Enode &
,Vminus_derivative_node(s_mode,v_mode),vref_node &
,Evalue
CALL create_new_node(combine_delays_s_Enode,next_free_node)
end if ! special case, Ts=Tv
end do ! next source domain mode
! Large resistance to complete the circuit for the series voltage sources
R_combine_delays_l_name(v_mode)='R_zt_dsum_l_'//trim(ZT_string)
write(spice_model_file_unit,'(A30,2I6,E16.6)')R_combine_delays_l_name(v_mode),combine_delays_l_Enode,vref_node,Rlarge
! Large resistance to complete the circuit for the series voltage sources
R_combine_delays_s_name(v_mode)='R_zt_dsum_s_'//trim(ZT_string)
write(spice_model_file_unit,'(A30,2I6,E16.6)')R_combine_delays_s_name(v_mode),combine_delays_s_Enode,vref_node,Rlarge
! Write filter function for integral of transfer impedance filter with the propagation correction i.e. Hp(jw)*(ZT(jw)-ZT_dc)/jw
! See Theory_Manual_Section 3.7, Theory_Manual_Equation 3.118.
if (.NOT.high_freq_Zt_model) then
! The d.c. transfer impedance has been included on the conductor termination so we must remove it here
ZT_Rdc=ZT_filter%a%coeff(0)/ZT_filter%b%coeff(0)
ZT_filter_minus_Rdc=ZT_filter ! NOTE:assumes a%order>=b%order
do i=0,ZT_filter_minus_Rdc%b%order
ZT_filter_minus_Rdc%a%coeff(i)=ZT_filter_minus_Rdc%a%coeff(i)-ZT_Rdc*ZT_filter_minus_Rdc%b%coeff(i)
end do
! set up a filter function with transfer function 1/jw
integrator_filter=allocate_Sfilter(0,1)
integrator_filter%wnorm=1d0
integrator_filter%a%coeff(0)=1d0
integrator_filter%b%coeff(0)=0d0
integrator_filter%b%coeff(1)=1d0
integrate_ZT_filter=integrator_filter*ZT_filter_minus_Rdc
! multiply the mode propagation correction by the time integral function
! Note the order of multiplication... This keeps the wnormalisation from Hpvfilter in the result.
Hp_integrate_ZT_filter=integrate_ZT_filter*Hpv_filter(v_mode)
! This filter function now has a0=b0=0 so divide top and bottom by s to give the final filter
aorder=Hp_integrate_ZT_filter%a%order
border=Hp_integrate_ZT_filter%b%order
! if ZT is purely resistive i.e. Zt=rdc then aorder-1=-1 so set Hp_integrate_ZT_filter_with_cancellation to a zero filter
if (aorder.EQ.0) then
Hp_integrate_ZT_filter_with_cancellation=0d0
else
Hp_integrate_ZT_filter_with_cancellation=allocate_Sfilter(aorder-1,border-1)
Hp_integrate_ZT_filter_with_cancellation%wnorm=Hp_integrate_ZT_filter%wnorm
! numerator terms
do i=0,aorder-1
Hp_integrate_ZT_filter_with_cancellation%a%coeff(i)=Hp_integrate_ZT_filter%a%coeff(i+1)
end do
! denominator terms
do i=0,border-1
Hp_integrate_ZT_filter_with_cancellation%b%coeff(i)=Hp_integrate_ZT_filter%b%coeff(i+1)
end do
end if
else ! use the high_freq_Zt_model
! set up a filter function with transfer function 1/jw
integrator_filter=allocate_Sfilter(0,1)
integrator_filter%wnorm=1d0
integrator_filter%a%coeff(0)=1d0
integrator_filter%b%coeff(0)=0d0
integrator_filter%b%coeff(1)=1d0
integrate_ZT_filter=integrator_filter*ZT_filter
! multiply the mode propagation correction by the time integral function
! Note the order of multiplication... This keeps the wnormalisation from Hpvfilter in the result.
Hp_integrate_ZT_filter_with_cancellation=integrate_ZT_filter*Hpv_filter(v_mode)
end if
! Use new subroutines for writing s-domain transfer function sources here
! Theory_Manual_Eqn 3.132
CALL write_spice_comment('Transfer impedance sources, end 1')
E_ZT_s_name(v_mode)='ZT_s_'//trim(ZT_string)
CALL write_s_domain_controlled_voltage_source(E_ZT_s_name(v_mode), &
combine_delays_s_Enode,vref_node, &
Vv_end1_node(v_mode),Vv_ref_end1_node(v_mode), &
Hp_integrate_ZT_filter_with_cancellation,1d0, &
vref_node,next_free_node) ! note: gain set to 1.0
! Theory_Manual_Eqn 3.131
CALL write_spice_comment('Transfer impedance sources, end 2')
E_ZT_l_name(v_mode)='ZT_l_'//trim(ZT_string)
CALL write_s_domain_controlled_voltage_source(E_ZT_l_name(v_mode), &
combine_delays_l_Enode,vref_node, &
Vv_end2_node(v_mode),Vv_ref_end2_node(v_mode), &
Hp_integrate_ZT_filter_with_cancellation,1d0, &
vref_node,next_free_node) ! note: gain set to 1.0
! deallocate the temporary filter data
CALL deallocate_Sfilter(integrate_ZT_filter)
CALL deallocate_Sfilter(Hp_integrate_ZT_filter)
CALL deallocate_Sfilter(ZT_filter_minus_Rdc)
CALL deallocate_Sfilter(Hp_integrate_ZT_filter_with_cancellation)
end do ! next victim mode
RETURN
END SUBROUTINE write_transfer_impedance_circuit