write_incident_field_excitation_circuit.F90 22.4 KB
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554
!
! 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_incident_field_circuit
!
! NAME
!     write_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.8 
!     and the Spice circuit is seen in Theory_Manual_Figures 3.17, and 3.18
!
!     INPUTS REQUIRED
!     1.    nodes for connection of the incident field sources
!     2.    time delays for the viictim domain modes
!     3.    bundle length
!     4.    modal decomposition matrices
!     5.    victim domain propagation correction filter functions
!     6.    pre-calculated coefficients related to the incident plane wave field excitation 
!     7.    nodes for the specification of the incident field function
!     8.    next free node and reference node for the spice sub-circuit
!
!
!     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.
!
! HISTORY
!
!     STAGE 6 developments started 16/06/2016 CJS.
!     Include ground plane 28/6/2016 CJS
!     Replace Tz transmission lines with a very small delay by a very small resistance - causes long runtimes in Pspice. 6/12/2016.
!     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_incident_field_circuit(n_victim_domain_modes, &
                             Vv_end1_node,Vv_end2_node, &
                             T_victim, &
                             length, &
                             TVI_victim,TI_victim, &
                             Hpv_filter, &
                             Tz,alpha0,alphaL, &
                             next_free_node,Einc_node1,Einc_node2,vref_node,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_victim_domain_modes   ! number of modes in the victim domain

! victim domain node lists for connection of the incident field source terms at each end.
integer,intent(IN) :: Vv_end1_node(1:n_victim_domain_modes)
integer,intent(IN) :: Vv_end2_node(1:n_victim_domain_modes)

! mode delays for the victim domain transmission line modes
real(dp),intent(IN) :: T_victim(1:n_victim_domain_modes)

! modal decomposition matrices for the victim domain
real(dp),intent(IN) :: TVI_victim(1:n_victim_domain_modes,1:n_victim_domain_modes)
real(dp),intent(IN) :: TI_victim(1:n_victim_domain_modes,1:n_victim_domain_modes)
 
real(dp),intent(IN) :: length                                   ! bundle length
TYPE(Sfilter),intent(IN) :: Hpv_filter(1:n_victim_domain_modes) ! victim mode proopagation correction terms

real(dp),intent(IN) :: Tz                                    ! incident field excitation delay time
real(dp),intent(IN) :: alpha0(1:n_victim_domain_modes)       ! pre-calculated model constants 
real(dp),intent(IN) :: alphaL(1:n_victim_domain_modes)       ! pre-calculated model constants 

integer,intent(INOUT) :: next_free_node      ! the next free node number for the spice sub-circuit 
integer,intent(IN) :: Einc_node1             ! excitation function sub-circuit terminal node 1
integer,intent(IN) :: Einc_node2             ! excitation function sub-circuit terminal node 2 
integer,intent(IN) :: vref_node              ! spice sub-circuit reference node number
integer,intent(IN) :: victim_domain          ! victim domain number

! local variables

integer :: v_mode

real(dp) :: Z0_delay=50d0   ! transmission line impedance for pure delay lines

! Names for delay lines

character(len=spice_name_length) :: delay_line_Tz_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_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_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) :: E_Einc_l_name(1:n_victim_domain_modes)
character(len=spice_name_length) :: E_Einc_s_name(1:n_victim_domain_modes)

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

! names for special case Tz=Tvictim components
character(len=spice_name_length) :: G_Einc_derivative_name(1:n_victim_domain_modes)
character(len=spice_name_length) :: L_Einc_derivative_name(1:n_victim_domain_modes)

! working strings

character(len=spice_name_length) :: name1,name2,name3
character(len=spice_name_length) :: Einc_string

! String used if a transmission line delay is very small and the T element is replaced by a resistance
character(len=spice_name_length) :: Rtempstring

! circuit to combine incident field terms

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

real(dp)       :: Evalue

! 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)

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            ! gain value for s-domain transfer functions

real(dp) :: Tz_minus_Tv     ! difference between source and victim mode velocities
logical  :: Tz_equal_Tv     ! if incident field delay and victim mode delay are the same then we have a different circuit topology

! loop variables
integer :: row,i

! nodes for special case Tz=Tvictim components
integer :: Einc_derivative_node(1:n_victim_domain_modes)

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

! START

if (verbose) write(*,*)'CALLED write_incient_field_excitation_circuit'

source_scale=c0  

! check for the case where Tz<0 and stop with an error for now. 
! Eventually we will cope with this situation by reversing the model.
if (Tz.LT.0d0) then
  run_status='ERROR in write_incient_field_excitation_circuit. Tz is less than zero'
  CALL write_program_status()
  STOP 1
end if

CALL write_spice_comment('START OF INCIDENT FIELD EXCITATION MODELS')

! 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
      
! Set Tz delay line component names

delay_line_Tz_name='T_Tz_Einc'
  
! mode impedance
  
delay_line_ZC_Tz_name='RZC_Tz_Einc'

! source terms
  
delay_line_E1_Tz_name='E1_Tz_Einc'
      
! Write Tz delay line components
! Theory_Manual_Fig 3.17, Theory_Manual_Eqn 3.168b

if (Tz.GT.Tz_min_delay) then
! significant delay so use a T element
  
  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
! the delay is very small so just use a very small series resistance instead of a delay line
  
  CALL write_spice_comment('Incident field delay lines, Tz delay=0. Delay line replaced by small resistance') 

  Rtempstring='R'//trim(delay_line_Tz_name)
  write(spice_model_file_unit,'(A30,2I6,E16.6)')Rtempstring,  &
                                           delay_line_Tz_s_nodes,delay_line_Tz_l_nodes,Rsmall
end if
 
! modal impedances on modal delay lines, Tz
 
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 Einc_string which labels the transfer impedance model number plus the source mode and victim mode numbers
  name1='EINC_vm_'
  CALL add_integer_to_string(name1,v_mode,Einc_string)

! we always need the Tz+T_victim(v_mode) delay lines so write these components now
! Theory_Manual_Fig 3.17, Theory_Manual_Eqn 3.168a

! 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(Einc_string)
  
! mode impedance
  
  delay_line_ZC_TzPTv_name(v_mode)='RZC_TzPTv_'//trim(Einc_string)

! source terms
  
  delay_line_E1_TzPTv_name(v_mode)='E1_TzPTv_'//trim(Einc_string)  
    
! Write TzPTv delay lines
! Theory_Manual_Fig 3.17, Theory_Manual_Eqn 3.168a

  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      

! 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
! Theory_Manual_Section 3.8.1, Theory_Manual_Eqn 3.169

  Tz_minus_Tv=Tz-T_victim(v_mode)

  if (abs(Tz_minus_Tv).GT.Einc_min_delay) then
  
    Tz_equal_Tv=.FALSE.
    
! The whole incident field excitation circuit is implemented using delay lines
! Theory_Manual_Fig 3.17

! 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(Einc_string) 
  
! mode impedance
  
    delay_line_ZC_Tv_name(v_mode)='RZC_Tv_'//trim(Einc_string)

! source terms
  
    delay_line_E1_Tv_name(v_mode)='E1_Tv_'//trim(Einc_string)
  
! Write Tv delay lines

    CALL write_spice_comment('Incident field 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, T_victim(v_mode)

  
    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 for positive z propagation, T_victim(v_mode)
  
    CALL write_spice_comment('Controlled source for 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 
       
  else
! Special case Tz=T_victim(v_mode)as see in Theory_Manual_Eqn 3.169  Theory_Manual_Fig 3.18 

    Tz_equal_Tv=.TRUE.
      
! **** The special case required time derivative circuits operating on the delayed incident field excitation****
! new node for time derivative delayed incident field function
    next_free_node=next_free_node+1
    Einc_derivative_node(v_mode)=next_free_node
  
! ****** Names for the special case circuit for Tsource=Tvictim  
    G_Einc_derivative_name(v_mode)='G_Vp_ddt_'//trim(Einc_string)
  
    L_Einc_derivative_name(v_mode)='L_Vp_ddt_'//trim(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-Tz)') 

    write(spice_model_file_unit,'(A30,4I6,E16.6)')G_Einc_derivative_name(v_mode), &
                             Einc_derivative_node(v_mode),vref_node, &
                             delay_line_Tz_l_nodes,vref_node,1.0 
                       
    CALL write_spice_comment('Inductor for derivative  of delayed Einc i.e. Einc(0,t-Tz)') 
  
    write(spice_model_file_unit,'(A30,2I6,E16.6)')L_Einc_derivative_name(v_mode), &
                                 Einc_derivative_node(v_mode),vref_node,1d0/source_scale

  end if ! Special case Tz-Tv =0
          
! The remaining part of the circuit combines all the contributions to the victim mode voltage source
! So as to evaluate the incident field terms in Theory_Manual_Eqn 3.168a,b
! We create the nodes for the summation circuit as we go
    
! create Einc_string which labels the transfer impedance model number plus victim mode number
  name1='EINC_vm_'
  CALL add_integer_to_string(name1,v_mode,Einc_string)

! first add the contributions which are common to both forms of circuit  

! Incident field excitation voltage source names      

  combine_delays_s_E_name(v_mode,1)='E_zt_dsum_s_'//trim(Einc_string)//'_E1'
  combine_delays_s_E_name(v_mode,2)='E_zt_dsum_s_'//trim(Einc_string)//'_E2'
  
 ! START OF CIRCUIT TO COMBINE INCIDENT FIELD EXCITATION TERMS

!  calculation of V_victim at z=0. Theory_Manual_Eqn 3.168a
 
  CALL write_spice_comment('Circuit to combine incident field excitation terms') 

! Einc, no delay:   Theory_Manual_Eqn 3.168a, line 2, term 1.
  next_free_node=next_free_node+1
  Evalue=+alpha0(v_mode)/(Tz+T_victim(v_mode))
  write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_s_E_name(v_mode,1) &
                                ,next_free_node,vref_node &
                                ,Einc_node1,Einc_node2 &
                                ,Evalue
  combine_delays_s_Enode=next_free_node
    
! Einc, delay=T_victim+Tz:  Theory_Manual_Eqn 3.168a, line 2, term 2
  next_free_node=next_free_node+1
  Evalue=-alpha0(v_mode)/(Tz+T_victim(v_mode))
  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_TzPTv_l_nodes(v_mode),vref_node &
                                 ,Evalue 
  combine_delays_s_Enode=next_free_node
    
!  calculation of V_victim at z=L. Theory_Manual_Eqn 3.168 a
                             			      			      
  if (.NOT.Tz_equal_Tv) then

! normal form based on delay lines,     Theory_Manual_Eqn 3.168b
    combine_delays_l_E_name(v_mode,1)='E_zt_dsum_l_'//trim(Einc_string)//'_E1'
    combine_delays_l_E_name(v_mode,2)='E_zt_dsum_l_'//trim(Einc_string)//'_E2'

! Einc, delay=T_victim: Theory_Manual_Eqn 3.168b, line 2, term 1.
    next_free_node=next_free_node+1
    Evalue=+alphaL(v_mode)/(T_victim(v_mode)-Tz)
    write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_l_E_name(v_mode,1) &
                                ,next_free_node,vref_node &
                                ,delay_line_Tv_l_nodes(v_mode),vref_node &
                                ,Evalue
    combine_delays_l_Enode=next_free_node
    
! Einc, delay=Tz: Theory_Manual_Eqn 3.168b, line 2, term 2.
    next_free_node=next_free_node+1
    Evalue=-alphaL(v_mode)/(T_victim(v_mode)-Tz)
    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_Tz_l_nodes,vref_node &
                                 ,Evalue 
    combine_delays_l_Enode=next_free_node

  else
! we need to calcuate the contribution to z=L using the time derivative of the delayed incident field  
! Theory_Manual_Eqn 3.169, Theory_Manual_Fig 3.18

    combine_delays_l_E_name(v_mode,1)='E_zt_dsum_l_'//trim(Einc_string)//'_E1'
     
    next_free_node=next_free_node+1
    Evalue=-alphaL(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,combine_delays_l_Enode  &
                                 ,Einc_derivative_node(v_mode),vref_node &
                                 ,Evalue 
    combine_delays_l_Enode=next_free_node
    
  end if  ! special case T_victim=Tz
    
  CALL write_spice_comment('Large resistance to complete the circuit for the series voltage sources') 
  
! Resistance to complete the circuit for the series voltage sources
  R_combine_delays_l_name(v_mode)='R_Einc_dsum_l_'//trim(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
  
! Resistance to complete the circuit for the series voltage sources
  R_combine_delays_s_name(v_mode)='R_Einc_dsum_s_'//trim(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 
  
! Write filter function for propagation correction filter 
! Theory_Manual_Equation 3.170

  CALL write_spice_comment('Incident field excitation sources, end 1')
  
  E_Einc_s_name(v_mode)='Einc_s_'//trim(Einc_string)
  CALL write_s_domain_controlled_voltage_source(E_Einc_s_name(v_mode),                          &
                                                combine_delays_s_Enode,vref_node,               &
                                                Vv_end1_node(v_mode),Vref_node,                 &
                                                Hpv_filter(v_mode),1d0,vref_node,next_free_node) ! note: gain set to 1.0


  CALL write_spice_comment('Incident field excitation sources, end 2')
  
  E_Einc_l_name(v_mode)='Einc_l_'//trim(Einc_string)
  CALL write_s_domain_controlled_voltage_source(E_Einc_l_name(v_mode),                          &
                                                combine_delays_l_Enode,vref_node,               &
                                                Vv_end2_node(v_mode),Vref_node,                 &
                                                Hpv_filter(v_mode),1d0,vref_node,next_free_node) ! note: gain set to 1.0

end do ! next victim mode

RETURN

END SUBROUTINE write_incident_field_circuit