Blame view

SRC/CREATE_SPICE_CIRCUIT_MODEL/write_transfer_impedance_circuit.F90 46.5 KB
886c558b   Steve Greedy   SACAMOS Public Re...
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

! 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 
189467e4   Steve Greedy   First Public Release
63
64
65
!     14/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
886c558b   Steve Greedy   SACAMOS Public Re...
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
!
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

189467e4   Steve Greedy   First Public Release
130
131
real(dp) :: source_scale     ! scaling factor for controlled sources to ensure that 
                             ! voltage/ currents stay in a sensible range
886c558b   Steve Greedy   SACAMOS Public Re...
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

! 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

189467e4   Steve Greedy   First Public Release
263
264
265
! scaling factor for controlled sources to ensure that voltage/ currents stay in a sensible range
source_scale=c0   

886c558b   Steve Greedy   SACAMOS Public Re...
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
! 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.
  
189467e4   Steve Greedy   First Public Release
372
373
374
375
376
377
378
379
!    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_delay_line(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, &
                          Z_source(s_mode),T_source(s_mode)+T_victim(v_mode),length)
886c558b   Steve Greedy   SACAMOS Public Re...
380
381
382

    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.
  
189467e4   Steve Greedy   First Public Release
383
384
385
386
387
388
389
390
391
!    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)
    CALL write_delay_line(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, &
                          Z_source(s_mode),T_source(s_mode)+T_victim(v_mode),length)
    
886c558b   Steve Greedy   SACAMOS Public Re...
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
! 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.

189467e4   Steve Greedy   First Public Release
455
456
457
458
459
460
461
462
!      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_delay_line(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, &
                            Z_source(s_mode),T_source(s_mode),length)
886c558b   Steve Greedy   SACAMOS Public Re...
463
464
465

      CALL write_spice_comment('Delay lines for negative z propagation, T_source(s_mode) delay') ! Theory_Manaul_Eqn 3.132, line 3, term2.

189467e4   Steve Greedy   First Public Release
466
467
468
469
470
471
472
473
!      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)
      CALL write_delay_line(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, &
                            Z_source(s_mode),T_source(s_mode),length)
886c558b   Steve Greedy   SACAMOS Public Re...
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
 
! 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
   
189467e4   Steve Greedy   First Public Release
549
550
551
552
553
554
555
556
!      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_delay_line(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, &
                            Z_source(s_mode),T_victim(v_mode),length)
886c558b   Steve Greedy   SACAMOS Public Re...
557
558
559

      CALL write_spice_comment('Delay lines for negative z propagation, T_victim(v_mode)') ! Theory_Manaul_Eqn 3.132, line 3, term1
  
189467e4   Steve Greedy   First Public Release
560
561
562
563
564
565
566
567
!      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)
      CALL write_delay_line(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, &
                            Z_source(s_mode),T_victim(v_mode),length)
886c558b   Steve Greedy   SACAMOS Public Re...
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620

! 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, &
189467e4   Steve Greedy   First Public Release
621
                               delay_line_pz_Ts_l_nodes(s_mode,v_mode),vref_node,1d0
886c558b   Steve Greedy   SACAMOS Public Re...
622
                        
189467e4   Steve Greedy   First Public Release
623
      CALL write_spice_comment('Inductor for derivative of positive z propagating voltage wave, Vs+(0,t-Ts)')
886c558b   Steve Greedy   SACAMOS Public Re...
624
625
  
      write(spice_model_file_unit,'(A30,2I6,E16.6)')L_Vplus_derivative_name(s_mode,v_mode), &
189467e4   Steve Greedy   First Public Release
626
                                   Vplus_derivative_node(s_mode,v_mode),vref_node,1d0/source_scale
886c558b   Steve Greedy   SACAMOS Public Re...
627
628
629
630
631
632
633
634
635

! ***** 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, &
189467e4   Steve Greedy   First Public Release
636
                                  delay_line_mz_Ts_l_nodes(s_mode,v_mode),vref_node,1D0
886c558b   Steve Greedy   SACAMOS Public Re...
637
                               
189467e4   Steve Greedy   First Public Release
638
      CALL write_spice_comment('Inductor for derivative of negative z propagating voltage wave, Vs-(0,t-Ts)')
886c558b   Steve Greedy   SACAMOS Public Re...
639
640
      
      write(spice_model_file_unit,'(A30,2I6,E16.6)')L_Vminus_derivative_name(s_mode,v_mode), &
189467e4   Steve Greedy   First Public Release
641
                                  Vminus_derivative_node(s_mode,v_mode),vref_node,1D0 /source_scale
886c558b   Steve Greedy   SACAMOS Public Re...
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691

    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)))
189467e4   Steve Greedy   First Public Release
692
    Evalue=Evalue/source_scale
886c558b   Steve Greedy   SACAMOS Public Re...
693
694
695
    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 &
189467e4   Steve Greedy   First Public Release
696
                                ,Evalue
886c558b   Steve Greedy   SACAMOS Public Re...
697
698
699
700
701
    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)))
189467e4   Steve Greedy   First Public Release
702
    Evalue=Evalue/source_scale
886c558b   Steve Greedy   SACAMOS Public Re...
703
704
705
706
707
708
709
710
711
712
713
    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)))    
189467e4   Steve Greedy   First Public Release
714
    Evalue=Evalue/source_scale
886c558b   Steve Greedy   SACAMOS Public Re...
715
716
717
718
719
720
721
722
723
    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)))     
189467e4   Steve Greedy   First Public Release
724
    Evalue=Evalue/source_scale
886c558b   Steve Greedy   SACAMOS Public Re...
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
    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)))     
189467e4   Steve Greedy   First Public Release
747
      Evalue=Evalue/source_scale
886c558b   Steve Greedy   SACAMOS Public Re...
748
749
750
751
752
753
754
755
756
      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)))     
189467e4   Steve Greedy   First Public Release
757
      Evalue=Evalue/source_scale
886c558b   Steve Greedy   SACAMOS Public Re...
758
759
760
761
762
763
764
765
766
      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)))
189467e4   Steve Greedy   First Public Release
767
      Evalue=Evalue/source_scale
886c558b   Steve Greedy   SACAMOS Public Re...
768
769
770
771
772
773
774
775
776
777
      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)))
189467e4   Steve Greedy   First Public Release
778
      Evalue=Evalue/source_scale
886c558b   Steve Greedy   SACAMOS Public Re...
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
      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
  
189467e4   Steve Greedy   First Public Release
817
! Resistance to complete the circuit for the series voltage sources
886c558b   Steve Greedy   SACAMOS Public Re...
818
  R_combine_delays_l_name(v_mode)='R_zt_dsum_l_'//trim(ZT_string)
189467e4   Steve Greedy   First Public Release
819
  write(spice_model_file_unit,'(A30,2I6,E16.6)')R_combine_delays_l_name(v_mode),combine_delays_l_Enode,vref_node,Rcombine_sources
886c558b   Steve Greedy   SACAMOS Public Re...
820
  
189467e4   Steve Greedy   First Public Release
821
! Resistance to complete the circuit for the series voltage sources
886c558b   Steve Greedy   SACAMOS Public Re...
822
  R_combine_delays_s_name(v_mode)='R_zt_dsum_s_'//trim(ZT_string)
189467e4   Steve Greedy   First Public Release
823
  write(spice_model_file_unit,'(A30,2I6,E16.6)')R_combine_delays_s_name(v_mode),combine_delays_s_Enode,vref_node,Rcombine_sources
886c558b   Steve Greedy   SACAMOS Public Re...
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
  
! 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')
189467e4   Steve Greedy   First Public Release
897
898
  E_ZT_s_name(v_mode)='ZT_s_'//trim(ZT_string)    

886c558b   Steve Greedy   SACAMOS Public Re...
899
900
901
  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), &
189467e4   Steve Greedy   First Public Release
902
903
904
905
                                                Hp_integrate_ZT_filter_with_cancellation,source_scale,  &
                                                vref_node,next_free_node)         
                                                
                                                ! note: gain set to 1.0 but with source scaling applied
886c558b   Steve Greedy   SACAMOS Public Re...
906
907
908
909
910
911
912
913


! 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), &
189467e4   Steve Greedy   First Public Release
914
915
916
917
                                                Hp_integrate_ZT_filter_with_cancellation,source_scale,  &
                                                vref_node,next_free_node) 
                                                        
                                                ! note: gain set to 1.0 but with source scaling applied
886c558b   Steve Greedy   SACAMOS Public Re...
918
919
920
921
922
923
924
925
926
927
928
929

! 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