!
! 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 <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_ZT_incident_field_circuit
!
! NAME
!     write_ZT_incident_field_circuit
!
! AUTHORS
!     Chris Smartt
!
! DESCRIPTION
!     This code writes the circuit components required to implement the incident field excitation model
!     The Spice model is described in Theory_Manual_Section 3.9 
!     and the Spice circuit is seen in Theory_Manual_Figures 3.19, and 3.20.
!
!     INPUTS REQUIRED
!     1. source (external) 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, mode delays and characteristic impedances
!     3. victim domain propagation correction filter functions for each mode
!     4. Incident field source function nodes and incident field delay in z, Tz
!     5. transfer impedance model
!     6. bundle length
!     7. Precalculated constants for the model
!
!     OUTPUTS
!     The components required to implement the model are written to
!     the subcircuit file
!     
! COMMENTS
!     Write the circuit elements required for the incident field excitation model
!
!     The excitation function comes from the voltage between two subcircuit termination nodes.
!
!     We will need to deal with the cases when Tz is less than zero and when Tz is equal to zero.
!
! HISTORY
!
!     STAGE 9 developments started 27/06/2016 CJS.
!     9/5/2017 CJS Document software with reference to Theory_Manual 
!     15/10/2017 CJS Include source scaling to keep the voltages and currents in a sensible range
!     20/10/2017 CJS call subroutine to write delay lines so we can choose to use LTRA or T elements
!     16/11/2017 CJS Include network synthesis process to replace s-domain transfer functions
!
SUBROUTINE write_ZT_incident_field_circuit(ZT_model,source_domain,victim_domain,n_source_domain_modes,n_victim_domain_modes, &
                                           Vv_end1_node, &
                                           Vv_end2_node, &
                                           Vv_end1_ref_node, &
                                           Vv_end2_ref_node, &
                                           Einc_node1,Einc_node2, &
                                           Vplus_node,Vminus_node, &
                                           TVI_source,TI_source, &
                                           T_source, &
                                           Zc_source, &
                                           TVI_victim,TI_victim, &
                                           source_domain_shield_conductor, &
                                           victim_domain_shield_conductor, &
                                           T_victim, &
                                           Zc_victim, &
                                           ZT_filter, &
                                           Hpv_filter, &
                                           Tz, &
                                           c1,c2,c1b, &
                                           length, &
                                           next_free_node,vref_node)
                                                                                     
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) :: ZT_model                        ! transfer impedance number
integer,intent(IN) :: source_domain                   ! source domain number
integer,intent(IN) :: victim_domain                   ! victim domain number
integer,intent(IN) :: n_source_domain_modes           ! number of modes in the source domain 
integer,intent(IN) :: n_victim_domain_modes           ! number of modes in the victim domain

integer,intent(IN) :: Vv_end1_node(1:n_victim_domain_modes)  ! first node for end 1 source terms in the victim domain
integer,intent(IN) :: Vv_end2_node(1:n_victim_domain_modes)  ! first node for end 1 source terms in the victim domain

integer,intent(IN) :: Vv_end1_ref_node(1:n_victim_domain_modes)  ! second(reference) node for end 1 source terms in the victim domain
integer,intent(IN) :: Vv_end2_ref_node(1:n_victim_domain_modes)  ! second(reference) node for end 2 source terms in the victim domain

integer,intent(IN) :: Einc_node1  ! first node for incident field excitation function (voltage)
integer,intent(IN) :: Einc_node2  ! second node for incident field excitation function (voltage)

integer,intent(IN) :: Vplus_node(1:n_source_domain_modes)   ! node for +z travelling characteristic in source domain
integer,intent(IN) :: Vminus_node(1:n_source_domain_modes)  ! node for -z travelling characteristic in source domain

real(dp),intent(IN) :: TVI_source(1:n_source_domain_modes,1:n_source_domain_modes)  ! mode decomposition TVI matrix, source domain
real(dp),intent(IN) :: TI_source(1:n_source_domain_modes,1:n_source_domain_modes)   ! mode decomposition TI matrix, source domain

real(dp),intent(IN) :: T_source(1:n_source_domain_modes)                            ! mode delay, source domain

real(dp),intent(IN) :: Zc_source(1:n_source_domain_modes)                           ! mode impedance, source domain

real(dp),intent(IN) :: TVI_victim(1:n_victim_domain_modes,1:n_victim_domain_modes)  ! mode decomposition TVI matrix, victim domain
real(dp),intent(IN) :: TI_victim(1:n_victim_domain_modes,1:n_victim_domain_modes)   ! mode decomposition TI matrix, victim domain

integer,intent(IN)  :: source_domain_shield_conductor                               ! conductor number of source domain shield conductor 
integer,intent(IN)  :: victim_domain_shield_conductor                               ! conductor number of victim domain shield conductor 

real(dp),intent(IN) :: T_victim(1:n_victim_domain_modes)                            ! mode delay, victim domain

real(dp),intent(IN) :: Zc_victim(1:n_victim_domain_modes)                           ! mode impedance, victim domain

TYPE(Sfilter),intent(IN) :: ZT_filter                                               ! Transfer impedance function

TYPE(Sfilter),intent(IN) :: Hpv_filter(1:n_victim_domain_modes)                     ! propagation correction, victim domain

real(dp),intent(IN) :: Tz                                                           ! incident field excitation delay time

real(dp),intent(IN) :: c1(1:n_source_domain_modes)                                  ! Constant 1 realted to incident field specification
real(dp),intent(IN) :: c2(1:n_source_domain_modes)                                  ! Constant 2 realted to incident field specification
real(dp),intent(IN) :: c1b(1:n_source_domain_modes)                                 ! Constant 1b realted to incident field specification for Tz=Ts

real(dp),intent(IN) :: length                                                       ! transmission line length

integer,intent(INOUT):: next_free_node    ! next free node - used to keep track for creation of new nodes in the model

integer,intent(IN) :: vref_node         ! reference node for the multi-conductor model

! local variables

! loop variables for source and victim modes
integer :: v_mode
integer :: s_mode

real(dp) :: Z0_delay=50d0         ! transmission line impedance for pure delay lines for which the delay only, not the characteristic impedance matters

! Names for delay lines

character(len=spice_name_length) :: delay_line_Tz_name

character(len=spice_name_length) :: RTz_name

character(len=spice_name_length) :: delay_line_Tv_name(1:n_victim_domain_modes)
character(len=spice_name_length) :: delay_line_TzPTv_name(1:n_victim_domain_modes)

character(len=spice_name_length) :: delay_line_Ts_name(1:n_source_domain_modes)
character(len=spice_name_length) :: delay_line_TzPTs_name(1:n_source_domain_modes)

character(len=spice_name_length) :: delay_line_TvPTs_name(1:n_victim_domain_modes,1:n_source_domain_modes)
character(len=spice_name_length) :: delay_line_TzPTvPTs_name(1:n_victim_domain_modes,1:n_source_domain_modes)

! Names for delay line termination impedances

character(len=spice_name_length) :: delay_line_ZC_Tz_name

character(len=spice_name_length) :: delay_line_ZC_Tv_name(1:n_victim_domain_modes)
character(len=spice_name_length) :: delay_line_ZC_TzPTv_name(1:n_victim_domain_modes)

character(len=spice_name_length) :: delay_line_ZC_Ts_name(1:n_source_domain_modes)
character(len=spice_name_length) :: delay_line_ZC_TzPTs_name(1:n_source_domain_modes)

character(len=spice_name_length) :: delay_line_ZC_TvPTs_name(1:n_victim_domain_modes,1:n_source_domain_modes)
character(len=spice_name_length) :: delay_line_ZC_TzPTvPTs_name(1:n_victim_domain_modes,1:n_source_domain_modes)

! Names for delay line sources

character(len=spice_name_length) :: delay_line_E1_Tz_name

character(len=spice_name_length) :: delay_line_E1_Tv_name(1:n_victim_domain_modes)
character(len=spice_name_length) :: delay_line_E1_TzPTv_name(1:n_victim_domain_modes)

character(len=spice_name_length) :: delay_line_E1_Ts_name(1:n_source_domain_modes)
character(len=spice_name_length) :: delay_line_E1_TzPTs_name(1:n_source_domain_modes)

character(len=spice_name_length) :: delay_line_E1_TvPTs_name(1:n_victim_domain_modes,1:n_source_domain_modes)
character(len=spice_name_length) :: delay_line_E1_TzPTvPTs_name(1:n_victim_domain_modes,1:n_source_domain_modes)

! Nodes for delay lines 

! delay line Tz 
integer :: delay_line_Tz_s_nodes
integer :: delay_line_Tz_l_nodes

! delay line Tv 
integer :: delay_line_Tv_s_nodes(1:n_victim_domain_modes)
integer :: delay_line_Tv_l_nodes(1:n_victim_domain_modes)

! delay line Tz+Tv 
integer :: delay_line_TzPTv_s_nodes(1:n_victim_domain_modes)
integer :: delay_line_TzPTv_l_nodes(1:n_victim_domain_modes)

! delay line Ts 
integer :: delay_line_Ts_s_nodes(1:n_source_domain_modes)
integer :: delay_line_Ts_l_nodes(1:n_source_domain_modes)

! delay line Tz+Ts 
integer :: delay_line_TzPTs_s_nodes(1:n_source_domain_modes)
integer :: delay_line_TzPTs_l_nodes(1:n_source_domain_modes)

! delay line Tv+Ts 
integer :: delay_line_TvPTs_s_nodes(1:n_victim_domain_modes,1:n_source_domain_modes)
integer :: delay_line_TvPTs_l_nodes(1:n_victim_domain_modes,1:n_source_domain_modes)

! delay line Tz+Tv+Ts
integer :: delay_line_TzPTvPTs_s_nodes(1:n_victim_domain_modes,1:n_source_domain_modes)
integer :: delay_line_TzPTvPTs_l_nodes(1:n_victim_domain_modes,1:n_source_domain_modes)

! circuit to combine incident field terms

character(len=spice_name_length) :: combine_delays_s_E_name(1:n_victim_domain_modes,1:8)
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_victim_domain_modes,1:8)
character(len=spice_name_length) :: R_combine_delays_l_name(1:n_victim_domain_modes)

real(dp)       :: Evalue

real(dp)       :: Const_Fs1,Const_Fs2,Const_Fs3,Const_Fs4  ! temporary constants used in calculating controlled source values

integer :: last_Hjw_source_s_node(1:n_victim_domain_modes)
integer :: last_Hjw_source_l_node(1:n_victim_domain_modes)

integer :: combine_delays_s_Enode
integer :: combine_delays_l_Enode

real(dp) :: gain      ! temporary varaible for gain of s-domain transfer functions

! temporary variables for string creation
character(len=spice_name_length) :: name1
character(len=spice_name_length) :: name2
character(len=spice_name_length) :: ZT_Einc_tag
character(len=spice_name_length) :: ZT_Einc_string

! loop variable
integer :: i

! 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
type(Sfilter) :: ZT_filter_minus_Rdc

integer :: row
integer :: aorder,border

character(len=spice_name_length) :: E_EincZT_l_name(1:n_victim_domain_modes)
character(len=spice_name_length) :: E_EincZT_s_name(1:n_victim_domain_modes)

character(len=spice_name_length) :: E_EincZT_filter_name(1:n_victim_domain_modes)

integer  :: first_node

real(dp) :: TI_source_row(1:n_source_domain_modes)
real(dp) :: TVI_victim_row(1:n_victim_domain_modes)
real(dp) :: PS_PV

! flags for special cases

logical :: Tv_equals_Ts(1:n_victim_domain_modes,1:n_source_domain_modes)
logical :: Tv_equals_Tz(1:n_victim_domain_modes)
logical :: Ts_equals_Tz(1:n_source_domain_modes)

! additional component names and nodes for special cases

character(len=spice_name_length) :: Inductor_Tv_equals_Ts_name1(1:n_victim_domain_modes,1:n_source_domain_modes)
character(len=spice_name_length) :: Inductor_Tv_equals_Ts_name2(1:n_victim_domain_modes,1:n_source_domain_modes)

character(len=spice_name_length) :: G_Tv_equals_Ts_name1(1:n_victim_domain_modes,1:n_source_domain_modes)
character(len=spice_name_length) :: G_Tv_equals_Ts_name2(1:n_victim_domain_modes,1:n_source_domain_modes)

integer :: Tv_equals_Ts_nodes1(1:n_victim_domain_modes,1:n_source_domain_modes)
integer :: Tv_equals_Ts_nodes2(1:n_victim_domain_modes,1:n_source_domain_modes)

character(len=spice_name_length) :: Inductor_Tv_equals_Tz_name(1:n_victim_domain_modes)

character(len=spice_name_length) :: G_Tv_equals_Tz_name(1:n_victim_domain_modes)

integer :: Tv_equals_Tz_nodes(1:n_victim_domain_modes)

real(dp) :: source_scale,Lsource_scale,max_E1     ! scaling factor for controlled sources to ensure that 
                                    ! voltage/ currents stay in a sensible range

! START
if (verbose) then
  write(*,*)'CALLED write_ZT_incient_field_excitation_circuit'
  write(*,*)
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.9.2, 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.199
   
  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.199

  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,  Theory_Manual_Section 3.9.2, 3.7.1

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

  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

! 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

  do s_mode=1,n_source_domain_modes

! calculate the scaling factor originating from the modal decomposition matrices  
    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
    
  end do
  
end do
! write delay lines for the incident field function

CALL write_spice_comment('START OF INCIDENT FIELD EXCITATION COUPLED THROUGH TRANSFER IMPEDANCE TERMS') 
write(spice_model_file_unit,'(A,I6)')'* Victim domain: ',victim_domain
write(spice_model_file_unit,'(A)')'*'

! ****** Tz DELAY LINE ****** Theory_Manual_Figure 3.20

! Set Tz delay line node numbers

! delay line nodes for positive z propagation, source end
next_free_node=next_free_node+1
delay_line_Tz_s_nodes=next_free_node

! delay line nodes for positive z propagation, load end
next_free_node=next_free_node+1
delay_line_Tz_l_nodes=next_free_node

! Label all components with ZT_EINC_+ZT_model
        
name1='ZTEi_'
CALL add_integer_to_string(name1,ZT_model,ZT_Einc_tag)

! Set Tz delay line component names
! mode impedance
  
delay_line_ZC_Tz_name='RZC_Tz_'//trim(ZT_Einc_tag)

! source terms
  
delay_line_E1_Tz_name='E1_Tz_'//trim(ZT_Einc_tag)
  
! Write Tz delay line components

if (abs(Tz).GT.Einc_min_delay) then

  delay_line_Tz_name='T_Tz_'//trim(ZT_Einc_tag)

! Implement Tz as a delay line
  
  CALL write_spice_comment('Incident field delay lines, Tz delay')
  
!  write(spice_model_file_unit,'(A30,4I6,A4,E16.6,A4,E16.6)')delay_line_Tz_name,  &
!                                           delay_line_Tz_s_nodes,vref_node, &
!                                           delay_line_Tz_l_nodes,vref_node, &
!                                           ' Z0=',Z0_delay,' TD=',Tz
  CALL write_delay_line(delay_line_Tz_name,  &
                        delay_line_Tz_s_nodes,vref_node, &
                        delay_line_Tz_l_nodes,vref_node, &
                        Z0_delay,Tz,length)
  
else

! Tz is zero so replace the delay line with a very small resistance 
! modal impedances on modal delay lines, Tz

  RTz_name='R_Tz_'//trim(ZT_Einc_tag)
  
  CALL write_spice_comment('Incident field delay lines, Special case for Tz=0, delay replaced by small resistance')
  
  write(spice_model_file_unit,'(A30,2I6,E16.6)')RTz_name,  &
                                           delay_line_Tz_s_nodes,delay_line_Tz_l_nodes,Rsmall

end if

CALL write_spice_comment('Matched impedance: Tz delay')
  
write(spice_model_file_unit,'(A30,2I6,E16.6)')delay_line_ZC_Tz_name, &
                             delay_line_Tz_l_nodes,vref_node,Z0_delay

! delay line controlled source for positive z propagation, Ts
  
CALL write_spice_comment('Delay line controlled sources Tz delay')

write(spice_model_file_unit,'(A30,4I6,E16.6)')delay_line_E1_Tz_name,&
                            delay_line_Tz_s_nodes,vref_node, &
                            Einc_node1,Einc_node2,1.0 
 
! End of Tz delay line

do v_mode=1,n_victim_domain_modes ! loop over victim domain modes
  
! create ZT_Einc_string which labels the transfer impedance model number plus the source mode and victim mode numbers
  name1=trim(ZT_Einc_tag)//'_vm_'
  CALL add_integer_to_string(name1,v_mode,ZT_Einc_string)

! ****** Tv DELAY LINE ****** Theory_Manual_Figure 3.20

! Set Tv delay line nodes

! delay line nodes, source end
  next_free_node=next_free_node+1
  delay_line_Tv_s_nodes(v_mode)=next_free_node

! delay line nodes, load end
  next_free_node=next_free_node+1
  delay_line_Tv_l_nodes(v_mode)=next_free_node

! Set Tv delay line component names

  delay_line_Tv_name(v_mode)='T_Tv_'//trim(ZT_Einc_string)
  
! mode impedance
  
  delay_line_ZC_Tv_name(v_mode)='RZC_Tv_'//trim(ZT_Einc_string)

! source terms
  
  delay_line_E1_Tv_name(v_mode)='E1_Tv_'//trim(ZT_Einc_string)  
    
! Write Tv delay lines
  
  CALL write_spice_comment('Delay line T_victim(v_mode)')
  
!  write(spice_model_file_unit,'(A30,4I6,A4,E16.6,A4,E16.6)')delay_line_Tv_name(v_mode),&
!                                       delay_line_Tv_s_nodes(v_mode),vref_node, &
!                                       delay_line_Tv_l_nodes(v_mode),vref_node, &
!                                       ' Z0=',Z0_delay,' TD=',T_victim(v_mode)
  CALL write_delay_line(delay_line_Tv_name(v_mode),&
                        delay_line_Tv_s_nodes(v_mode),vref_node, &
                        delay_line_Tv_l_nodes(v_mode),vref_node, &
                        Z0_delay,T_victim(v_mode),length)

! modal impedances on modal delay lines, Tv
  
  CALL write_spice_comment('Matched impedance:  T_victim(v_mode) delay')
  
  write(spice_model_file_unit,'(A30,2I6,E16.6)')delay_line_ZC_Tv_name(v_mode),        &
                               delay_line_Tv_l_nodes(v_mode),vref_node,Z0_delay

! delay line controlled source, Tv
  
  CALL write_spice_comment('Incident field delay line controlled source:  T_victim(v_mode) delay')
  
  write(spice_model_file_unit,'(A30,4I6,E16.6)')delay_line_E1_Tv_name(v_mode),&
                               delay_line_Tv_s_nodes(v_mode),vref_node, &
                               Einc_node1,Einc_node2,1.0

! End of Tv delay lines      

! ****** Tz + Tv DELAY LINE ****** Theory_Manual_Figure 3.20

! we always need the Tz+T_victim(v_mode) delay lines so write these components now

! Set TzPTv delay line nodes

! delay line nodes, source end
  next_free_node=next_free_node+1
  delay_line_TzPTv_s_nodes(v_mode)=next_free_node

! delay line nodes, load end
  next_free_node=next_free_node+1
  delay_line_TzPTv_l_nodes(v_mode)=next_free_node

! Set TzPTv delay line component names

  delay_line_TzPTv_name(v_mode)='T_TzPTv_'//trim(ZT_Einc_string)
  
! mode impedance
  
  delay_line_ZC_TzPTv_name(v_mode)='RZC_TzPTv_'//trim(ZT_Einc_string)

! source terms
  
  delay_line_E1_TzPTv_name(v_mode)='E1_TzPTv_'//trim(ZT_Einc_string)  
    
! Write TzPTv delay lines
  
  CALL write_spice_comment('Delay line Tz+T_victim(v_mode)')
  
!  write(spice_model_file_unit,'(A30,4I6,A4,E16.6,A4,E16.6)')delay_line_TzPTv_name(v_mode),&
!                                       delay_line_TzPTv_s_nodes(v_mode),vref_node, &
!                                       delay_line_TzPTv_l_nodes(v_mode),vref_node, &
!                                       ' Z0=',Z0_delay,' TD=',Tz+T_victim(v_mode)
  CALL write_delay_line(delay_line_TzPTv_name(v_mode),&
                        delay_line_TzPTv_s_nodes(v_mode),vref_node, &
                        delay_line_TzPTv_l_nodes(v_mode),vref_node, &
                        Z0_delay,Tz+T_victim(v_mode),length)

! modal impedances on modal delay lines, TzPTv
  
  CALL write_spice_comment('Matched impedance:  Tz+T_victim(v_mode) delay')
  
  write(spice_model_file_unit,'(A30,2I6,E16.6)')delay_line_ZC_TzPTv_name(v_mode),        &
                               delay_line_TzPTv_l_nodes(v_mode),vref_node,Z0_delay

! delay line controlled source, TzPTv
  
  CALL write_spice_comment('Incident field delay line controlled source:  Tz+T_victim(v_mode) delay')
  
  write(spice_model_file_unit,'(A30,4I6,E16.6)')delay_line_E1_TzPTv_name(v_mode),&
                               delay_line_TzPTv_s_nodes(v_mode),vref_node, &
                               Einc_node1,Einc_node2,1.0

! End of TzPTv delay lines      
                               
end do ! next victim domain mode

do s_mode=1,n_source_domain_modes ! loop over source domain modes
  
! create ZT_Einc_string which labels the transfer impedance model number plus the source mode and source mode numbers
  name1=trim(ZT_Einc_tag)//'_sm_'
  CALL add_integer_to_string(name1,s_mode,ZT_Einc_string)

! ****** Ts DELAY LINE ****** Theory_Manual_Figure 3.20

! we always need the T_source(s_mode) delay lines so write these components now

! Set Ts delay line nodes

! delay line nodes, source end
  next_free_node=next_free_node+1
  delay_line_Ts_s_nodes(s_mode)=next_free_node

! delay line nodes, load end
  next_free_node=next_free_node+1
  delay_line_Ts_l_nodes(s_mode)=next_free_node

! Set Ts delay line component names

  delay_line_Ts_name(s_mode)='T_Ts_'//trim(ZT_Einc_string)
  
! mode impedance
  
  delay_line_ZC_Ts_name(s_mode)='RZC_Ts_'//trim(ZT_Einc_string)

! source terms
  
  delay_line_E1_Ts_name(s_mode)='E1_Ts_'//trim(ZT_Einc_string)  
    
! Write Ts delay lines
  
  CALL write_spice_comment('Delay line T_source(s_mode)')
  
!  write(spice_model_file_unit,'(A30,4I6,A4,E16.6,A4,E16.6)')delay_line_Ts_name(s_mode),&
!                                       delay_line_Ts_s_nodes(s_mode),vref_node, &
!                                       delay_line_Ts_l_nodes(s_mode),vref_node, &
!                                       ' Z0=',Z0_delay,' TD=',T_source(s_mode)
  CALL write_delay_line(delay_line_Ts_name(s_mode),&
                        delay_line_Ts_s_nodes(s_mode),vref_node, &
                        delay_line_Ts_l_nodes(s_mode),vref_node, &
                        Z0_delay,T_source(s_mode),length)

! modal impedances on modal delay lines, Ts
  
  CALL write_spice_comment('Matched impedance:  T_source(s_mode) delay')
  
  write(spice_model_file_unit,'(A30,2I6,E16.6)')delay_line_ZC_Ts_name(s_mode),        &
                               delay_line_Ts_l_nodes(s_mode),vref_node,Z0_delay

! delay line controlled source, Ts
  
  CALL write_spice_comment('Incident field delay line controlled source:  T_source(s_mode) delay')
  
  write(spice_model_file_unit,'(A30,4I6,E16.6)')delay_line_E1_Ts_name(s_mode),&
                               delay_line_Ts_s_nodes(s_mode),vref_node, &
                               Einc_node1,Einc_node2,1.0

! End of Ts delay lines 
     
! ****** Tz + Ts DELAY LINE ******  Theory_Manual_Figure 3.20
! Set TzPTs delay line nodes 

! delay line nodes, source end
  next_free_node=next_free_node+1
  delay_line_TzPTs_s_nodes(s_mode)=next_free_node

! delay line nodes, load end
  next_free_node=next_free_node+1
  delay_line_TzPTs_l_nodes(s_mode)=next_free_node

! Set TzPTs delay line component names

  delay_line_TzPTs_name(s_mode)='T_TzPTs_'//trim(ZT_Einc_string)
  
! mode impedance
  
  delay_line_ZC_TzPTs_name(s_mode)='RZC_TzPTs_'//trim(ZT_Einc_string)

! source terms
  
  delay_line_E1_TzPTs_name(s_mode)='E1_TzPTs_'//trim(ZT_Einc_string)  
    
! Write TzPTs delay lines
  
  CALL write_spice_comment('Delay line Tz+T_source(s_mode)')
  
!  write(spice_model_file_unit,'(A30,4I6,A4,E16.6,A4,E16.6)')delay_line_TzPTs_name(s_mode),&
!                                       delay_line_TzPTs_s_nodes(s_mode),vref_node, &
!                                       delay_line_TzPTs_l_nodes(s_mode),vref_node, &
!                                       ' Z0=',Z0_delay,' TD=',Tz+T_source(s_mode)
  CALL write_delay_line(delay_line_TzPTs_name(s_mode),&
                        delay_line_TzPTs_s_nodes(s_mode),vref_node, &
                        delay_line_TzPTs_l_nodes(s_mode),vref_node, &
                        Z0_delay,Tz+T_source(s_mode),length)

! modal impedances on modal delay lines, TzPTs
  
  CALL write_spice_comment('Matched impedance:  Tz+T_source(s_mode) delay')
  
  write(spice_model_file_unit,'(A30,2I6,E16.6)')delay_line_ZC_TzPTs_name(s_mode),        &
                               delay_line_TzPTs_l_nodes(s_mode),vref_node,Z0_delay

! delay line controlled source, TzPTs
  
  CALL write_spice_comment('Incident field delay line controlled source:  Tz+T_source(s_mode) delay')
  
  write(spice_model_file_unit,'(A30,4I6,E16.6)')delay_line_E1_TzPTs_name(s_mode),&
                               delay_line_TzPTs_s_nodes(s_mode),vref_node, &
                               Einc_node1,Einc_node2,1.0

! End of TzPTs delay lines      
                               
end do ! next source domain mode


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_Einc_string which labels the transfer impedance model number plus the source mode and victim mode numbers
    name1=trim(ZT_Einc_tag)//'_vm_'
    CALL add_integer_to_string(name1,v_mode,name2)
    name1=trim(name2)//'_sm_'
    CALL add_integer_to_string(name1,s_mode,ZT_Einc_string)

! ****** Tv + Ts DELAY LINE ******  Theory_Manual_Figure 3.20

! Set TvPTs delay line nodes

! delay line nodes, source end
    next_free_node=next_free_node+1
    delay_line_TvPTs_s_nodes(v_mode,s_mode)=next_free_node

! delay line nodes, load end
    next_free_node=next_free_node+1
    delay_line_TvPTs_l_nodes(v_mode,s_mode)=next_free_node

! Set TvPTs delay line component names

    delay_line_TvPTs_name(v_mode,s_mode)='T_TvPTs_'//trim(ZT_Einc_string)
  
! mode impedance
  
    delay_line_ZC_TvPTs_name(v_mode,s_mode)='RZC_TvPTs_'//trim(ZT_Einc_string)

! source terms
  
    delay_line_E1_TvPTs_name(v_mode,s_mode)='E1_TvPTs_'//trim(ZT_Einc_string)  
    
! Write TvPTs delay lines
  
    CALL write_spice_comment('Delay line T_source(s_mode)+T_victim(v_mode)')  
  
!    write(spice_model_file_unit,'(A30,4I6,A4,E16.6,A4,E16.6)')delay_line_TvPTs_name(v_mode,s_mode),&
!                                         delay_line_TvPTs_s_nodes(v_mode,s_mode),vref_node, &
!                                         delay_line_TvPTs_l_nodes(v_mode,s_mode),vref_node, &
!                                         ' Z0=',Z0_delay,' TD=',T_source(s_mode)+T_victim(v_mode)
    CALL write_delay_line(delay_line_TvPTs_name(v_mode,s_mode),&
                          delay_line_TvPTs_s_nodes(v_mode,s_mode),vref_node, &
                          delay_line_TvPTs_l_nodes(v_mode,s_mode),vref_node, &
                          Z0_delay,T_source(s_mode)+T_victim(v_mode),length)

! modal impedances on modal delay lines, TvPTs
  
    CALL write_spice_comment('Matched impedance:  T_source(s_mode)+T_victim(v_mode) delay')
  
    write(spice_model_file_unit,'(A30,2I6,E16.6)')delay_line_ZC_TvPTs_name(v_mode,s_mode),        &
                                 delay_line_TvPTs_l_nodes(v_mode,s_mode),vref_node,Z0_delay

! delay line controlled source, TvPTs
  
    CALL write_spice_comment('Incident field delay line controlled source:  T_source(s_mode)+T_victim(v_mode) delay')
   
    write(spice_model_file_unit,'(A30,4I6,E16.6)')delay_line_E1_TvPTs_name(v_mode,s_mode),&
                                 delay_line_TvPTs_s_nodes(v_mode,s_mode),vref_node, &
                                 Einc_node1,Einc_node2,1.0

! End of TvPTs delay lines      

! ****** Tz + Ts + Tv DELAY LINE ******  Theory_Manual_Figure 3.20

! Set TzPTvPTs delay line nodes

! delay line nodes, source end
    next_free_node=next_free_node+1
    delay_line_TzPTvPTs_s_nodes(v_mode,s_mode)=next_free_node

! delay line nodes, load end
    next_free_node=next_free_node+1
    delay_line_TzPTvPTs_l_nodes(v_mode,s_mode)=next_free_node

! Set TzPTvPTs delay line component names

    delay_line_TzPTvPTs_name(v_mode,s_mode)='T_TzPTvPTs_'//trim(ZT_Einc_string)
  
! mode impedance
  
    delay_line_ZC_TzPTvPTs_name(v_mode,s_mode)='RZC_TzPTvPTs_'//trim(ZT_Einc_string)

! source terms
  
    delay_line_E1_TzPTvPTs_name(v_mode,s_mode)='E1_TzPTvPTs_'//trim(ZT_Einc_string)  
    
! Write TzPTvPTs delay lines
  
    CALL write_spice_comment('Delay line Tz+T_source(s_mode)+T_victim(v_mode)')
  
!    write(spice_model_file_unit,'(A30,4I6,A4,E16.6,A4,E16.6)')delay_line_TzPTvPTs_name(v_mode,s_mode),&
!                                         delay_line_TzPTvPTs_s_nodes(v_mode,s_mode),vref_node, &
!                                         delay_line_TzPTvPTs_l_nodes(v_mode,s_mode),vref_node, &
!                                         ' Z0=',Z0_delay,' TD=',Tz+T_source(s_mode)+T_victim(v_mode)
    CALL write_delay_line(delay_line_TzPTvPTs_name(v_mode,s_mode),&
                          delay_line_TzPTvPTs_s_nodes(v_mode,s_mode),vref_node, &
                          delay_line_TzPTvPTs_l_nodes(v_mode,s_mode),vref_node, &
                          Z0_delay,Tz+T_source(s_mode)+T_victim(v_mode),length)
 
! modal impedances on modal delay lines, TzPTvPTs
  
    CALL write_spice_comment('Matched impedance:  Tz+T_source(s_mode)+T_victim(v_mode) delay')
  
    write(spice_model_file_unit,'(A30,2I6,E16.6)')delay_line_ZC_TzPTvPTs_name(v_mode,s_mode),        &
                                 delay_line_TzPTvPTs_l_nodes(v_mode,s_mode),vref_node,Z0_delay

! delay line controlled source, TzPTvPTs
  
    CALL write_spice_comment('Incident field delay line controlled source:  Tz+T_source(s_mode)+T_victim(v_mode) delay')
   
    write(spice_model_file_unit,'(A30,4I6,E16.6)')delay_line_E1_TzPTvPTs_name(v_mode,s_mode),&
                                 delay_line_TzPTvPTs_s_nodes(v_mode,s_mode),vref_node, &
                                 Einc_node1,Einc_node2,1.0

! End of TzPTvPTs delay lines      


  end do ! next source mode
  
end do ! next victim mode

! Get an estimate of the controlled source values so that we can rescale them to sensible values
! just use the Fs2 and Fs4 terms to avoid the special cases here
max_E1=0d0

do v_mode=1,n_victim_domain_modes

  do s_mode=1,n_source_domain_modes
  
    Const_Fs2=c2(s_mode)/(2d0*Zc_source(s_mode)*(Tz+T_source(s_mode)))  ! initial sign from substitution in Xie equation 12a with 9, 13 and 20
  
    Const_Fs4=+c2(s_mode)/(2d0*Zc_source(s_mode)*(Tz+T_source(s_mode)))  ! initial sign from substitution in equation 12b with 9, 13 and 20

    Evalue=abs( Const_Fs2/(Tz+T_victim(v_mode)) )
    max_E1=max(max_E1,Evalue)
    
    Evalue=abs( Const_Fs4/(T_victim(v_mode)+T_source(s_mode)) )
    max_E1=max(max_E1,Evalue)
   
  end do
  
end do

! we can now scale the sources by max_E1 to reduce the source voltage values to order 1V
source_scale=max_E1
Lsource_scale=c0

! Check for special cases and set up the appropriate time derivative circuits as required

! TV_EQUALS_TZ. See Theory_Manual_Section 3.9.3 and Theory_Manual_Section 3.7.2 for the implementation details of these special cases

do v_mode=1,n_victim_domain_modes ! loop over victim domain modes

  if (abs(Tz-T_victim(v_mode)).GT.Einc_min_delay) then
  
    Tv_equals_Tz(v_mode)=.FALSE.
    
  else
  
    Tv_equals_Tz(v_mode)=.TRUE.
! create a new node for the time derivative circuit for Einc(t-Tv) required    
 
    next_free_node=next_free_node+1
    Tv_equals_Tz_nodes(v_mode)=next_free_node

! create names for the time derivative circuit elements      .

! create ZT_Einc_string which labels the transfer impedance model number plus the source mode and victim mode numbers
    name1=trim(ZT_Einc_tag)//'_vm_'
    CALL add_integer_to_string(name1,v_mode,ZT_Einc_string)

    Inductor_Tv_equals_Tz_name(v_mode)='L_Tv_EQ_Tz_'//trim(ZT_Einc_string)
                                                
    G_Tv_equals_Tz_name(v_mode)='G_Tv_EQ_Tz_'//trim(ZT_Einc_string)
     
! ***** Inductive circuit to calculate the time derivative of delayed Einc *****
  
    CALL write_spice_comment('Special case: Tv_equals_Tz, Controlled source for derivative of delayed Einc i.e. Einc(0,t-Tv)')

    write(spice_model_file_unit,'(A30,4I6,E16.6)')G_Tv_equals_Tz_name(v_mode), &
                                                  Tv_equals_Tz_nodes(v_mode),vref_node, &
                                                  delay_line_Tv_l_nodes(v_mode),vref_node,1.0 
                    
    CALL write_spice_comment('1H inductor for derivative  of delayed Einc i.e. Einc(0,t-Tv)')
   
    write(spice_model_file_unit,'(A30,2I6,E16.6)')Inductor_Tv_equals_Tz_name(v_mode), &
                                                  Tv_equals_Tz_nodes(v_mode),vref_node,1.0 
  end if
  
end do ! next v_mode

! TS_EQUALS_TZ. See Theory_Manual_Section 3.9.3 and Theory_Manual_Section 3.7.2 for the implementation details of these special cases

do s_mode=1,n_source_domain_modes ! loop over source domain modes
  if (abs(Tz-T_source(s_mode)).GT.Einc_min_delay) then
  
    Ts_equals_Tz(s_mode)=.FALSE.
    
  else
  
    Ts_equals_Tz(s_mode)=.TRUE.

  end if

! no additional circuit elements are required in this case

end do ! next s_mode

! TV_EQUALS_TS. See Theory_Manual_Section 3.9.3 and Theory_Manual_Section 3.7.2 for the implementation details of these special cases

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
  
    if (abs(T_source(s_mode)-T_victim(v_mode)).GT.Einc_min_delay) then
  
      Tv_equals_Ts(v_mode,s_mode)=.FALSE.
    
    else
  
      Tv_equals_Ts(v_mode,s_mode)=.TRUE.

! create a new node for the time derivative circuit for Einc(t-Tz-Tv) required     
      next_free_node=next_free_node+1
      Tv_equals_Ts_nodes1(v_mode,s_mode)=next_free_node

! create names for the tie derivative circuit elements      .
! create ZT_Einc_string which labels the transfer impedance model number plus the source mode and victim mode numbers
      name1=trim(ZT_Einc_tag)//'_vm_'
      CALL add_integer_to_string(name1,v_mode,name2)
      name1=trim(name2)//'_sm_'
      CALL add_integer_to_string(name1,s_mode,ZT_Einc_string)

      Inductor_Tv_equals_Ts_name1(v_mode,s_mode)='L1_Tv_EQ_Ts_'//trim(ZT_Einc_string)
                                                
      G_Tv_equals_Ts_name1(v_mode,s_mode)='G1_Tv_EQ_Ts_'//trim(ZT_Einc_string)
     
! ***** Inductive circuit to calculate the time derivative of delayed Einc *****
  
      CALL write_spice_comment('Special case Tv_equals_Ts: Controlled source for derivative of delayed Einc, Einc(0,t-(Tz+Ts))')

      write(spice_model_file_unit,'(A30,4I6,E16.6)')G_Tv_equals_Ts_name1(v_mode,s_mode), &
                                                    Tv_equals_Ts_nodes1(v_mode,s_mode),vref_node, &
                                                    delay_line_TzPTv_l_nodes(v_mode),vref_node,1.0 
                      
      CALL write_spice_comment('Inductor for derivative  of delayed Einc i.e. Einc(0,t-(Tz+Ts))')
  
      write(spice_model_file_unit,'(A30,2I6,E16.6)')Inductor_Tv_equals_Ts_name1(v_mode,s_mode), &
                                                    Tv_equals_Ts_nodes1(v_mode,s_mode),vref_node,1D0/Lsource_scale

! create a new node for the time derivative circuit for Einc(t-Tv) required     
      next_free_node=next_free_node+1
      Tv_equals_Ts_nodes2(v_mode,s_mode)=next_free_node

! create names for the time derivative circuit elements      .

      Inductor_Tv_equals_Ts_name2(v_mode,s_mode)='L2_Tv_EQ_Ts_'//trim(ZT_Einc_string)
                                                
      G_Tv_equals_Ts_name2(v_mode,s_mode)='G2_Tv_EQ_Ts_'//trim(ZT_Einc_string)
     
! ***** Inductive circuit to calculate the time derivative of delayed Einc *****
  
      CALL write_spice_comment('Controlled source for derivative of delayed Einc i.e. Einc(0,t-Tv)') 

      write(spice_model_file_unit,'(A30,4I6,E16.6)')G_Tv_equals_Ts_name2(v_mode,s_mode), &
                                                    Tv_equals_Ts_nodes2(v_mode,s_mode),vref_node, &
                                                    delay_line_Tv_l_nodes(v_mode),vref_node,1.0 
                      
      CALL write_spice_comment('Inductor for derivative  of delayed Einc i.e. Einc(0,t-Tv)')
    
      write(spice_model_file_unit,'(A30,2I6,E16.6)')Inductor_Tv_equals_Ts_name2(v_mode,s_mode), &
                                                    Tv_equals_Ts_nodes2(v_mode,s_mode),vref_node,1D0/Lsource_scale
                                                    
    end if ! Tsource=Tvictim
    
  end do ! next source mode
  
end do ! next victim mode

! START OF CIRCUIT TO COMBINE INCIDENT FIELD EXCITATION TERMS
! note that constants c1, c2, c1b, are calculated in SUBROUTINE calculate_ZT_incident_field_sources
! and include the cable bundle length.


do v_mode=1,n_victim_domain_modes

  do s_mode=1,n_source_domain_modes

! calculate the scaling factor originating from the modal decomposition matrices  
    PS_PV=TI_source_row(s_mode)*TVI_victim_row(v_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
    
! create ZT_Einc_string which labels the transfer impedance model number plus the source mode and victim mode numbers
    name1=trim(ZT_Einc_tag)//'_v_'
    CALL add_integer_to_string(name1,v_mode,name2)
    name1=trim(name2)//'_s_'
    CALL add_integer_to_string(name1,s_mode,ZT_Einc_string)

! Incident field excitation voltage source names      
! We have 8 terms for Fs1 and Fs2, equations 20a and b in Xie
! Theory_Manaul_Eqn 3.197

    name1='Ezt_dsum_s_'//trim(ZT_Einc_string)//'_E'
    do i=1,8
      CALL add_integer_to_string(name1,i,combine_delays_s_E_name(v_mode,i))
    end do
 
!  calculation of V_victim at z=0
   
    CALL write_spice_comment('Circuit to combine incident field excitation terms for V_victim at end 1')  

! Coefficients of delayed incident field function from Theory_Manual_Eqn  3.197, 3.198 substituted into 3.195
! These are the constants multiplying the incident field function in the curly brackets in 3.197,3.198 when substituted into 3.195
  
    if (.NOT.Ts_equals_Tz(s_mode)) then
      Const_Fs1=+c1(s_mode)/(2d0*Zc_source(s_mode)*(Tz-T_source(s_mode)))  ! initial sign from substitution in Xie equation 12a with 9, 13 and 20
    else
! special case for Tz=T_source(s_mode)
      write(spice_model_file_unit,'(A)')'* Using Tz=Ts model '
      Const_Fs1=+c1b(s_mode)/(2d0*Zc_source(s_mode))
    end if
    Const_Fs2=c2(s_mode)/(2d0*Zc_source(s_mode)*(Tz+T_source(s_mode)))  ! initial sign from substitution in Xie equation 12a with 9, 13 and 20
      
    if (s_mode.EQ.1) then
      first_node=vref_node
    else
      first_node=combine_delays_s_Enode
    end if

! Einc, no delay:  Xie equation 20a, term 1, Theory_Manual_Equation 3.197, term1 substituted into 3.195
    next_free_node=next_free_node+1
    Evalue=-Const_Fs1*PS_PV/(T_victim(v_mode)+T_source(s_mode))
    Evalue=Evalue/source_scale
    write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_s_E_name(v_mode,1) &
                                ,next_free_node,first_node &
                                ,Einc_node1,Einc_node2 &
                                ,Evalue
  combine_delays_s_Enode=next_free_node
    
! Einc, delay=T_victim+T_source: Xie equation 20a,   term 2, Theory_Manual_Equation 3.197, term2 substituted into 3.195
    next_free_node=next_free_node+1
    Evalue=+Const_Fs1*PS_PV/(T_victim(v_mode)+T_source(s_mode))
    Evalue=Evalue/source_scale
    write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_s_E_name(v_mode,2) &
                                 ,next_free_node,combine_delays_s_Enode &
                                 ,delay_line_TvPTs_l_nodes(v_mode,s_mode),vref_node &
                                 ,Evalue 
  combine_delays_s_Enode=next_free_node
    
! Einc, no delay: Xie equation 20a,   term 3, Theory_Manual_Equation 3.197, term3 substituted into 3.195
    next_free_node=next_free_node+1
    Evalue=+Const_Fs1*PS_PV/(Tz+T_victim(v_mode))
    Evalue=Evalue/source_scale
    write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_s_E_name(v_mode,3) &
                                 ,next_free_node,combine_delays_s_Enode &
                                 ,Einc_node1,Einc_node2 &
                                 ,Evalue 
  combine_delays_s_Enode=next_free_node
    
! Einc, delay=Tz+T_victim: Xie equation 20a, Theory_Manual_Equation 3.197, term4 substituted into 3.195
    next_free_node=next_free_node+1
    Evalue=-Const_Fs1*PS_PV/(Tz+T_victim(v_mode))
    Evalue=Evalue/source_scale
    write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_s_E_name(v_mode,4) &
                                 ,next_free_node,combine_delays_s_Enode &
                                 ,delay_line_TzPTv_l_nodes(v_mode),vref_node &
                                 ,Evalue 
    combine_delays_s_Enode=next_free_node
    
! Einc, no delay:  Xie equation 20b, term 1, Theory_Manual_Equation 3.198, term1 substituted into 3.195
    next_free_node=next_free_node+1
    Evalue=-Const_Fs2*PS_PV/(Tz+T_victim(v_mode))
    Evalue=Evalue/source_scale
    write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_s_E_name(v_mode,5) &
                                ,next_free_node,combine_delays_s_Enode &
                                ,Einc_node1,Einc_node2 &
                                ,Evalue
    combine_delays_s_Enode=next_free_node
    
! Einc, delay=Tz+T_victim: Xie equation 20b,   term 2, Theory_Manual_Equation 3.198, term2 substituted into 3.195
    next_free_node=next_free_node+1
    Evalue=+Const_Fs2*PS_PV/(Tz+T_victim(v_mode))
    Evalue=Evalue/source_scale
    write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_s_E_name(v_mode,6) &
                                 ,next_free_node,combine_delays_s_Enode &
                                 ,delay_line_TzPTv_l_nodes(v_mode),vref_node &
                                 ,Evalue 
    combine_delays_s_Enode=next_free_node
          
    if ( .NOT.Tv_equals_Ts(v_mode,s_mode) ) then
! normal delay line implementation

! Einc, delay=Tz+T_victim: Xie equation 20b,   term 3, Theory_Manual_Equation 3.198, term3 substituted into 3.195
      next_free_node=next_free_node+1
      Evalue=-Const_Fs2*PS_PV/(T_victim(v_mode)-T_source(s_mode))
      Evalue=Evalue/source_scale
      write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_s_E_name(v_mode,7) &
                                   ,next_free_node,combine_delays_s_Enode &
                                   ,delay_line_TzPTv_l_nodes(v_mode),vref_node &
                                   ,Evalue 
    combine_delays_s_Enode=next_free_node
    
! Einc, delay=Tz+T_source: Xie equation 20b,   term 4, Theory_Manual_Equation 3.198, term4 substituted into 3.195
      next_free_node=next_free_node+1
      Evalue=+Const_Fs2*PS_PV/(T_victim(v_mode)-T_source(s_mode))
      Evalue=Evalue/source_scale
      write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_s_E_name(v_mode,8) &
                                   ,next_free_node,combine_delays_s_Enode &
                                   ,delay_line_TzPTs_l_nodes(s_mode),vref_node &
                                   ,Evalue 
      combine_delays_s_Enode=next_free_node
    
    else
! special case using derivative circuit to implement Theory_Manual_Equation 3.198, terms 3,4 substituted into 3.195

      next_free_node=next_free_node+1
      Evalue=-Const_Fs2*PS_PV
      Evalue=Evalue*Lsource_scale/source_scale
      write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_s_E_name(v_mode,8) &
                                   ,next_free_node,combine_delays_s_Enode &
                                   ,Tv_equals_Ts_nodes1(v_mode,s_mode),vref_node &
                                   ,Evalue 
      combine_delays_s_Enode=next_free_node

    end if  ! Tvictim=Tsource special case
    
  end do ! next s_mode
  
  CALL write_spice_comment('Large resistance to complete the circuit for the series voltage sources, end 1')
   
! Resistance to complete the circuit for the series voltage sources
  R_combine_delays_s_name(v_mode)='R_Einc_dsum_s_'//trim(ZT_Einc_string)
  write(spice_model_file_unit,'(A30,2I6,E16.6)')R_combine_delays_s_name(v_mode),combine_delays_s_Enode,vref_node,Rcombine_sources 

  do s_mode=1,n_source_domain_modes

! calculate the scaling factor originating from the modal decomposition matrices  
    PS_PV=TI_source_row(s_mode)*TVI_victim_row(v_mode)

! Incident field excitation voltage source names      
! We have 8 terms for Fs3 and Fs4, equations 20c and d in Xie
    
! create ZT_Einc_string which labels the transfer impedance model number plus the source mode and victim mode numbers
    name1=trim(ZT_Einc_tag)//'_v_'
    CALL add_integer_to_string(name1,v_mode,name2)
    name1=trim(name2)//'_s_'
    CALL add_integer_to_string(name1,s_mode,ZT_Einc_string)

    name1='E_zt_dsum_l_'//trim(ZT_Einc_string)//'_E'
    do i=1,8
      CALL add_integer_to_string(name1,i,combine_delays_l_E_name(v_mode,i))
    end do
 
 ! START OF CIRCUIT TO COMBINE INCIDENT FIELD EXCITATION TERMS

!  calculation of V_victim at z=L
  
    CALL write_spice_comment('Circuit to combine incident field excitation terms for V_victim at end 2')
  
! Coefficients of delayed incident field function from Theory_Manual_Eqn  3.199, 3.200 substituted into  3.196
! These are the constants multiplying the incident field function in the curly brackets in 3.199,3.200 when substituted into 3.196
    
    if (.NOT.Ts_equals_Tz(s_mode)) then
      Const_Fs3=c1(s_mode)/(2d0*Zc_source(s_mode)*(Tz-T_source(s_mode)))  ! initial sign from substitution in equation 12b with 9, 13 and 20
    else
! special case for Tz=T_source(s_mode)
      write(spice_model_file_unit,'(A)')'* Using Tz=Ts model '
      Const_Fs3=c1b(s_mode)/(2d0*Zc_source(s_mode))
    end if
  
    Const_Fs4=+c2(s_mode)/(2d0*Zc_source(s_mode)*(Tz+T_source(s_mode)))  ! initial sign from substitution in equation 12b with 9, 13 and 20
      
    if (s_mode.EQ.1) then
      first_node=vref_node
    else
      first_node=combine_delays_l_Enode
    end if

! Einc, delay: Tz Xie equation 20c, term 1, Theory_Manual_Equation 3.199, term1 substituted into 3.196
    if ( .NOT.Tv_equals_Tz(v_mode) ) then
! normal delay line implementation

      next_free_node=next_free_node+1
      Evalue=+Const_Fs3*PS_PV/(Tz-T_victim(v_mode))
      Evalue=Evalue/source_scale
      write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_l_E_name(v_mode,1) &
                                  ,next_free_node,first_node &
                                  ,delay_line_Tz_l_nodes,vref_node &
                                  ,Evalue
      combine_delays_l_Enode=next_free_node
    
! Einc, delay=T_victim: Xie equation 20c,   term 2, Theory_Manual_Equation 3.199, term2 substituted into 3.196
      next_free_node=next_free_node+1
      Evalue=-Const_Fs3*PS_PV/(Tz-T_victim(v_mode))
      Evalue=Evalue/source_scale
      write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_l_E_name(v_mode,2) &
                                   ,next_free_node,combine_delays_l_Enode &
                                   ,delay_line_Tv_l_nodes(v_mode),vref_node &
                                   ,Evalue 
      combine_delays_l_Enode=next_free_node
      
    else
    
! special case using derivative circuit, Theory_Manual_Equation 3.199, terms 1 and 2 substituted into 3.196

      next_free_node=next_free_node+1
      Evalue=+Const_Fs3*PS_PV
      Evalue=Evalue*Lsource_scale/source_scale
      write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_l_E_name(v_mode,2) &
                                   ,next_free_node,first_node &
                                   ,Tv_equals_Tz_nodes(v_mode),vref_node &
                                   ,Evalue 
      combine_delays_l_Enode=next_free_node
   
    end if
          
    if ( .NOT.Tv_equals_Ts(v_mode,s_mode) ) then
! normal delay line implementation
    
! Einc, delay=T_victim Xie equation 20c,   term 3, Theory_Manual_Equation 3.199, term3 substituted into 3.196
      next_free_node=next_free_node+1
      Evalue=-Const_Fs3*PS_PV/(T_victim(v_mode)-T_source(s_mode))
      Evalue=Evalue/source_scale
      write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_l_E_name(v_mode,3) &
                                   ,next_free_node,combine_delays_l_Enode &
                                   ,delay_line_Tv_l_nodes(v_mode),vref_node &
                                   ,Evalue 
    combine_delays_l_Enode=next_free_node
    
! Einc, delay=T_source: Xie equation 20c,   term 4, Theory_Manual_Equation 3.199, term4 substituted into 3.196
      next_free_node=next_free_node+1
      Evalue=+Const_Fs3*PS_PV/(T_victim(v_mode)-T_source(s_mode))
      Evalue=Evalue/source_scale
      write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_l_E_name(v_mode,4) &
                                   ,next_free_node,combine_delays_l_Enode &
                                   ,delay_line_Ts_l_nodes(s_mode),vref_node &
                                   ,Evalue 
      combine_delays_l_Enode=next_free_node
    
    else
! special case using derivative circuit, Theory_Manual_Equation 3.199, terms 3 and 4 substituted into 3.196

      next_free_node=next_free_node+1
      Evalue=-Const_Fs3*PS_PV
      Evalue=Evalue*Lsource_scale/source_scale
      write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_l_E_name(v_mode,4) &
                                   ,next_free_node,combine_delays_l_Enode &
                                   ,Tv_equals_Ts_nodes2(v_mode,s_mode),vref_node &
                                   ,Evalue 
      combine_delays_l_Enode=next_free_node

    end if  ! Tvictim=Tsource special case
    
    if ( .NOT.Tv_equals_Tz(v_mode) ) then
! normal delay line implementation

! Einc, delay: Tz Xie equation 20d, term 1, Theory_Manual_Equation 3.200, term 1 substituted into 3.196
      next_free_node=next_free_node+1
      Evalue=-Const_Fs4*PS_PV/(Tz-T_victim(v_mode))
      Evalue=Evalue/source_scale
      write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_l_E_name(v_mode,5) &
                                  ,next_free_node,combine_delays_l_Enode   &
                                  ,delay_line_Tz_l_nodes,vref_node &
                                  ,Evalue
      combine_delays_l_Enode=next_free_node
    
! Einc, delay=Tz+T_victim: Xie equation 20d,   term 2, Theory_Manual_Equation 3.200, term 2 substituted into 3.196
      next_free_node=next_free_node+1
      Evalue=+Const_Fs4*PS_PV/(Tz-T_victim(v_mode))
      Evalue=Evalue/source_scale
      write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_l_E_name(v_mode,6) &
                                   ,next_free_node,combine_delays_l_Enode &
                                   ,delay_line_Tv_l_nodes(v_mode),vref_node &
                                   ,Evalue 
      combine_delays_l_Enode=next_free_node
    
    else
 ! special case using derivative circuit, Theory_Manual_Equation 3.200, terms 1 and 2  substituted into 3.196

      next_free_node=next_free_node+1
      Evalue=-Const_Fs4*PS_PV
      Evalue=Evalue*Lsource_scale/source_scale
      write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_l_E_name(v_mode,6) &
                                   ,next_free_node,combine_delays_l_Enode &
                                   ,Tv_equals_Tz_nodes(v_mode),first_node &
                                   ,Evalue 
      combine_delays_l_Enode=next_free_node

    end if
    
! Einc, delay: Tz  Xie equation 20d,   term 3, Theory_Manual_Equation 3.200, term 3 substituted into 3.196
    next_free_node=next_free_node+1
    Evalue=-Const_Fs4*PS_PV/(T_victim(v_mode)+T_source(s_mode))
    Evalue=Evalue/source_scale
    write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_l_E_name(v_mode,7) &
                                 ,next_free_node,combine_delays_l_Enode &
                                 ,delay_line_Tz_l_nodes,vref_node &
                                 ,Evalue 
  combine_delays_l_Enode=next_free_node
    
! Einc, delay=Tz+T_source: Xie equation 20d,   term 4, Theory_Manual_Equation 3.200, term 4 substituted into 3.196
    next_free_node=next_free_node+1
    Evalue=Const_Fs4*PS_PV/(T_victim(v_mode)+T_source(s_mode))
    Evalue=Evalue/source_scale
    write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_l_E_name(v_mode,8) &
                                 ,next_free_node,combine_delays_l_Enode &
                                 ,delay_line_TzPTvPTs_l_nodes(v_mode,s_mode),vref_node &
                                 ,Evalue 
    combine_delays_l_Enode=next_free_node
    
  end do ! next s_mode
  
  CALL write_spice_comment('Large resistance to complete the circuit for the series voltage sources, load end')
  
! Resistance to complete the circuit for the series voltage sources
  R_combine_delays_l_name(v_mode)='R_Einc_dsum_l_'//trim(ZT_Einc_string)
  write(spice_model_file_unit,'(A30,2I6,E16.6)')R_combine_delays_l_name(v_mode),combine_delays_l_Enode,vref_node,Rcombine_sources

! write the controlled sources with s-domain transfer functions and propagation corrections for each mode
  
! Write filter function for integral of transfer impedance filter with the propagation correction i.e. Hp(jw)*(ZT(jw)-ZT_dc)/jw
! Theory_Manual_Section 3.9 text below Theory_Manual_Eqn 3.196

! The d.c. transfer impedance has been included on the conductor termination so we must remove it here

  if (.NOT.high_freq_Zt_model) then
  
    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
  
    Hp_integrate_ZT_filter_with_cancellation=allocate_Sfilter(aorder-1,border-1)
    write(*,*)'DONE Set  Hp_integrate_ZT_filter_with_cancellation '
  
    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
    
  else  ! don't remove Rdc.
    
    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

  CALL write_spice_comment('Transfer impedance sources, end 1')
  E_EincZT_s_name(v_mode)='ZT_s_'//trim(ZT_Einc_string)
  CALL write_s_domain_controlled_voltage_source(E_EincZT_s_name(v_mode),                           &
                                                combine_delays_s_Enode,vref_node,              &
                                                Vv_end1_node(v_mode),Vv_end1_ref_node(v_mode), &
                                                Hp_integrate_ZT_filter_with_cancellation,source_scale,  &
                                                vref_node,next_free_node)         ! note: gain set to 1.0 but with source_scale applied


  CALL write_spice_comment('Transfer impedance sources, end 2')
  E_EincZT_l_name(v_mode)='ZT_l_'//trim(ZT_Einc_string)
  CALL write_s_domain_controlled_voltage_source(E_EincZT_l_name(v_mode),                           &
                                                combine_delays_l_Enode,vref_node,              &
                                                Vv_end2_node(v_mode),Vv_end1_ref_node(v_mode), &
                                                Hp_integrate_ZT_filter_with_cancellation,source_scale,  &
                                                vref_node,next_free_node)         ! note: gain set to 1.0

! deallocate the temporary filter data
  CALL deallocate_Sfilter(ZT_filter_minus_Rdc)
  CALL deallocate_Sfilter(integrator_filter)
  CALL deallocate_Sfilter(integrate_ZT_filter)
  CALL deallocate_Sfilter(Hp_integrate_ZT_filter)
  CALL deallocate_Sfilter(Hp_integrate_ZT_filter_with_cancellation)
 
end do ! next victim mode

if(verbose) write(*,*)'FINISHED write_ZT_incient_field_excitation_circuit'

RETURN

END SUBROUTINE write_ZT_incident_field_circuit
