! 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 <http://www.gnu.org/licenses/>.
! 
! SACAMOS uses the EISPACK library (in /SRC/EISPACK). EISPACK is subject to 
! the GNU Lesser General Public License. A copy of the GNU Lesser General Public 
! License version can be found in the file GNU_LGPL in the root of EISPACK 
! (/SRC/EISPACK ) or at <http://www.gnu.org/licenses/>.
! 
! The University of Nottingham can be contacted at: ggiemr@nottingham.ac.uk
!
! File Contents:
! SUBROUTINE 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