Blame view

SRC/MTL_ANALYTIC_SOLUTION/time_domain_analysis.F90 29.2 KB
886c558b   Steve Greedy   SACAMOS Public Re...
1
!
fe64b32b   Chris Smartt   Update file heade...
2
! This file is part of SACAMOS, State of the Art CAble MOdels for Spice. 
886c558b   Steve Greedy   SACAMOS Public Re...
3
4
5
! It was developed by the University of Nottingham and the Netherlands Aerospace 
! Centre (NLR) for ESA under contract number 4000112765/14/NL/HK.
! 
fe64b32b   Chris Smartt   Update file heade...
6
! Copyright (C) 2016-2018 University of Nottingham
886c558b   Steve Greedy   SACAMOS Public Re...
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
! 
! 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
!
fe64b32b   Chris Smartt   Update file heade...
28
!
886c558b   Steve Greedy   SACAMOS Public Re...
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
555
556
557
558
559
560
561
562
563
564
565
566
567
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
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
! File Contents:
! SUBROUTINE time_domain_analysis
!
! NAME
!     time_domain_analysis
!
! AUTHORS
!     Chris Smartt
!
! DESCRIPTION
!     This subroutine controls the analytic solution for the transient analysis of
!     multi-conductor transmission lines. The solution is obtained using the Inverse 
!     Fourier Transform (IFT) method.
!     The solution is obtained using the full dimension transmission line equations 
!     i.e. we are NOT using the weak form of transfer impedance coupling
!     Note also that frequency dependent quantities are evaluated separately at
!     each frequency of analysis, i.e. the frequency dependence of the solution is rigorous 
!     given only the frequency dependence of the dielectrics is modelled using impedance/ admittance 
!     matrices whose elements are rational frequency dependent filter functions.
!
!     INPUTS:
!       spice_bundle_model structure
!       spice_validation_test structure 
!    
!     OUTPUT
!        analytic time domain termination voltage for the specified validation test case written to file
!
!    Note that the analysis assumes that the excitation (and hence the response) are periodic
!    with period equal to the maximimum time of the simulation. If the response of the circuit 
!    is longer than this then aliasing will occur and the comparison between Spice model results
!    and analytic results will be in error so please ensure that all transients have reduced
!    to insignificant levels within the simulation time.
!     
! COMMENTS
!     STAGE_1: frequency independent parameter solution
!     STAGE_2: multi-conductor solution
!     STAGE_3: shielded cable solution
!     STAGE_4: frequency dependent model 
!     STAGE_5: transfer impedance model 
!
! HISTORY
!
!     started 7/12/2015 CJS: STAGE_1 developments
!     24/03/2016 CJS: STAGE_3 developments -shielded cables
!     22/04/2016 CJS: STAGE_4 developments -frequency dependent model
!     04/05/2016 CJS: Only write time domain output over the time period requested,not the full FFT time period.
!    Include general conductor impedance model 12/05/2016 CJS
!     Fix bug with conductor impedance contributions 12/05/2016 CJS
!     8/09/2016 CJS Correct the common mode/ differential mode loss terms for twisted pairs
!     13/10/2016 CJS Correct transfer impedance for multiple modes in external domain
!     7/3/2017         CJS: Add resistance and voltage source onto the reference coonductor 
!     8/5/2017         CJS: Include references to Theory_Manual
!
!
SUBROUTINE time_domain_analysis(spice_bundle_model,spice_validation_test)

USE type_specifications
USE general_module
USE constants
USE cable_module
USE cable_bundle_module
USE spice_cable_bundle_module
USE maths
USE frequency_spec

IMPLICIT NONE

! variables passed to subroutine

TYPE(spice_model_specification_type),intent(IN):: spice_bundle_model    ! Spice cable bundle model structure

TYPE(spice_validation_test_type),intent(IN)    :: spice_validation_test ! Spice validation circuit structure

! local variables

! variables for time domain analytic solution using IFT

integer            :: n_timesteps             ! number of timesteps requested
integer            :: n_FFT                   ! nummber of samples in the FFT
complex(dp),allocatable    :: Vs_time(:)      ! Source voltage as a function of time
complex(dp),allocatable    :: Vs_freq(:)      ! Source voltage as a function of frequency
complex(dp),allocatable    :: VTL_time(:)     ! Output voltage as a function of time
complex(dp),allocatable    :: VTL_freq(:)     ! Output voltage as a function of frequency
real(dp),allocatable       :: time(:)         ! time values
real(dp),allocatable       :: freq(:)         ! frequency values

! variables for the frequency domain MTL solution

real(dp)    :: f,w            ! frequency and angular frequency
integer     :: frequency_loop ! frequency loop variable

integer        :: dim         ! dimension of matrix system to solve

! domain based impedance and admittance matrices
complex(dp),allocatable     :: Z_domain(:,:)
complex(dp),allocatable     :: Y_domain(:,:)

! domain based conductor impedance terms
complex(dp),allocatable     ::Z_domain_conductor_impedance_correction(:,:)

! Vectors and matrices used in the frequency domain solution of the transmission line equations with termination conditions 
complex(dp),allocatable     :: Vs1(:)
complex(dp),allocatable     :: Z1(:,:)
complex(dp),allocatable     :: Vs2(:)
complex(dp),allocatable     :: Z2(:,:)

complex(dp)    :: Vout ! complex output voltage value

! domain transformation matrices
complex(dp),allocatable     :: MI(:,:)
complex(dp),allocatable     :: MII(:,:)
complex(dp),allocatable     :: MV(:,:)
complex(dp),allocatable     :: MVI(:,:)

! temporary working matrices
complex(dp),allocatable     :: TM1(:,:)

integer    :: conductor,inner_domain,outer_domain

integer    :: domain1,inner_domain1,outer_domain1
integer    :: conductor1,reference_conductor1
integer    :: domain_conductor1,domain_reference_conductor1
logical    :: is_shield1

integer    :: domain2,inner_domain2,outer_domain2
integer    :: conductor2,reference_conductor2
integer    :: domain_conductor2,domain_reference_conductor2
logical    :: is_shield2

! conductor based impedance (loss) and transfer impedance model data
complex(dp) :: Zint_c          ! conductor surface impedance
complex(dp) :: Zint_c_ref      ! reference conductor surface impedance
real(dp)    :: Rdc_c           ! d.c. resistance of conductor
real(dp)    :: Rdc_c_ref       ! d.c. resistance of reference conductor
complex(dp) :: Zint_t          ! conductor transfer impedance
complex(dp) :: Zint_t_ref      ! reference conductor transfer impedance
real(dp)    :: Rdc_t           ! d.c. resistance of conductor (from transfer impedance)
real(dp)    :: Rdc_t_ref       ! d.c. resistance of reference conductor (from transfer impedance)

! complex amplitude of incident field
complex(dp) :: Einc

logical,allocatable  :: is_shielded_flag(:)            ! flags conductors which are not exposed to the incident field
integer              :: shield_conductor               ! temporary variable, shield conductor number for shielded conductors
real(dp),allocatable :: local_conductor_x_offset(:)    ! x coordinate in bundle cross section of conductors
real(dp),allocatable :: local_conductor_y_offset(:)    ! y coordinate in bundle cross section of conductors

integer             :: n_conductors_outer_domain ! for shield conductors, the number of conductors in the domain outside the shield
integer             :: shield_conductor_number_in_outer_domain ! for shield conductors, the conductor number in the domain outside the shield

! loop variables
integer    :: i
integer    :: row,col

integer :: ierr

! START

! initialise FFT stuff

! set the number of frequencies to be equal to the next power of 2 above the number of timesteps requested
! Theory_Manual_Section 2.3.3

  n_timesteps=NINT(spice_validation_test%runtime/spice_validation_test%timestep)+1
  
  if(verbose) write(*,*)'Number of timesteps in transient analysis:',n_timesteps
  
  i=1
10  CONTINUE
    i=i*2
    if (i.ge.n_timesteps) then
      n_FFT=i
    else
      GOTO 10
    end if
    
  if(verbose) write(*,*)'Number of frequency in transient analysis:',n_FFT

! allocate the arrays used in the FFT
  ALLOCATE( time(n_FFT) ) 
  ALLOCATE( Vs_time(n_FFT) ) 
  ALLOCATE( VTL_time(n_FFT) ) 
  ALLOCATE( freq(n_FFT) ) 
  ALLOCATE( Vs_freq(n_FFT) ) 
  ALLOCATE( VTL_freq(n_FFT) )   

! Generate the time domain voltage source function array
! The time domain waveform is the Exponential Pulse, the same as that used in Spice
! See Ngspice manual, section 4.1.3

  do i=1,n_FFT
  
    time(i)=(i-1)*spice_validation_test%timestep
        
    if (time(i).LT.spice_validation_test%width) then
      Vs_time(i)=(1d0-exp(-time(i)/spice_validation_test%risetime))
    else if (time(i).GE.spice_validation_test%width) then
      Vs_time(i)=(exp(-(time(i)-spice_validation_test%width)/spice_validation_test%risetime))
    end if
      
  end do

! FFT the time domain voltage source array to give the frequency domain excitation function
! Theory_Manual_Equation 2.48

  CALL FFT_TIME_TO_FREQ(n_FFT,time,Vs_time,freq,Vs_freq)
  
! allocate memory

  dim=spice_bundle_model%bundle%system_dimension
  ALLOCATE( Z_domain(dim,dim) )
  ALLOCATE( Y_domain(dim,dim) )
  
  ALLOCATE( Z_domain_conductor_impedance_correction(dim,dim) )
  
  ALLOCATE( Vs1(dim) )
  ALLOCATE( Z1(dim,dim) )
  ALLOCATE( Vs2(dim) )
  ALLOCATE( Z2(dim,dim) )
  
! domain transformation matrices
  ALLOCATE( MI(dim,dim) )
  ALLOCATE( MII(dim,dim) )
  ALLOCATE( MV(dim,dim) )
  ALLOCATE( MVI(dim,dim) )
  
! temporary working matrices
  ALLOCATE( TM1(dim,dim) )
   
   
  ALLOCATE( is_shielded_flag(1:dim+1) )
  ALLOCATE( local_conductor_x_offset(1:dim+1) )
  ALLOCATE( local_conductor_y_offset(1:dim+1) )
  
! loop over conductors to work out which are in shielded domains and which are in the external domain
! also get the position of the conductor in the bundle cross section for incident field excitation

  do i=1,dim+1
    if (spice_bundle_model%bundle%terminal_conductor_to_outer_domain(i).EQ.spice_bundle_model%bundle%tot_n_domains) then
      is_shielded_flag(i)=.FALSE.
      local_conductor_x_offset(i)=spice_bundle_model%bundle%conductor_x_offset(i)
      local_conductor_y_offset(i)=spice_bundle_model%bundle%conductor_y_offset(i)
    else
      is_shielded_flag(i)=.TRUE.
! work out the conductor number of the shield
      shield_conductor=spice_bundle_model%bundle%terminal_conductor_to_reference_terminal_conductor(i)
! shielded conductors pick up the coordinate of the shield for the purposes on incident field excitation
      local_conductor_x_offset(i)=spice_bundle_model%bundle%conductor_x_offset(shield_conductor)
      local_conductor_y_offset(i)=spice_bundle_model%bundle%conductor_y_offset(shield_conductor)
    end if
  end do
  
! build the termination specifications and convert to complex
  Vs1(1:dim)=cmplx( spice_validation_test%Vs_end1(1:dim)-spice_validation_test%Vs_end1(dim+1) )
  Vs2(1:dim)=cmplx( spice_validation_test%Vs_end2(1:dim)-spice_validation_test%Vs_end2(dim+1) )
  
  Z1(:,:) =cmplx( spice_validation_test%R_end1(dim+1)  )
  Z2(:,:) =cmplx( spice_validation_test%R_end2(dim+1)  )
  do i=1,dim
    Z1(i,i) =Z1(i,i)+cmplx( spice_validation_test%R_end1(i)  )
    Z2(i,i) =Z2(i,i)+cmplx( spice_validation_test%R_end2(i)  )
  end do

! Copy the domain transformation matrices and calculate the inverses
  MI(:,:)=spice_bundle_model%bundle%global_MI%mat(:,:)
  MV(:,:)=spice_bundle_model%bundle%global_MV%mat(:,:)
  
  ierr=0   ! set ierr=0 on input to matrix inverse to cause the program to stop if we have a singular matrix
  CALL cinvert_Gauss_Jordan(MI,dim,MII,dim,ierr)
  ierr=0   ! set ierr=0 on input to matrix inverse to cause the program to stop if we have a singular matrix
  CALL cinvert_Gauss_Jordan(MV,dim,MVI,dim,ierr)

! loop over frequency to calculate the frequency domain transmission line transfer function

  write(6,8100,advance='no')'Frequency ',0,' of ',n_FFT/2+1
  flush(6)
  
  do frequency_loop=1,n_FFT/2+1
  
    write(6,'(A)',advance='no')char(13)
    write(6,8100,advance='no')'Frequency ',frequency_loop,' of ',n_FFT/2+1
    flush(6)
    
8100 format(A10,I10,A4,I10)    

    f=freq(frequency_loop)

! shift frequency from d.c. as solution fails for lossless systems at f=0 Hz
    if (f.eq.0d0) then
      f=1d0
    end if 

    w=2d0*pi*f

! Use the global domain based L and C matrices and the domain voltage and current 
! domain transformation matrices to calculate the impedance [Z] and admittance [Y] matrices
    do row=1,dim
      do col=1,dim
      
! Evaluate the cable impedance filter function
        Z_domain(row,col)=evaluate_Sfilter_frequency_response(spice_bundle_model%bundle%global_Z%sfilter_mat(row,col),f)

! Evaluate the cable admittance filter function        
        Y_domain(row,col)=evaluate_Sfilter_frequency_response(spice_bundle_model%bundle%global_Y%sfilter_mat(row,col),f)

      end do  ! next col
       
    end do ! next row

! calculate the contribution to the matrices from the conductor based impedance models. 
! See Theory_Manual_Section 2.2.3 

!Initailly set to zero

    Z_domain_conductor_impedance_correction(1:dim,1:dim)=(0d0,0d0)

!   new domain based model

    do conductor1=1,dim
   
      domain1=spice_bundle_model%bundle%terminal_conductor_to_outer_domain(conductor1)
      reference_conductor1=spice_bundle_model%bundle%terminal_conductor_to_reference_terminal_conductor(conductor1)
      domain_conductor1=spice_bundle_model%bundle%terminal_conductor_to_global_domain_conductor(conductor1)
      domain_reference_conductor1=spice_bundle_model%bundle%terminal_conductor_to_global_domain_conductor(reference_conductor1)
      is_shield1=spice_bundle_model%bundle%terminal_conductor_is_shield_flag(conductor1)
      
! evaluate the surface impedance for this conductor conductor 
      CALL evaluate_conductor_impedance_model(spice_bundle_model%bundle%conductor_impedance(conductor1),    &
                                              f,Zint_c,Rdc_c,Zint_t,Rdc_t)
                                              
! Apply multiplication factor to the conductor impedance to correct for common mode/ differential modes in twisted pairs
! See note at the end of Theory_Manual_Section 3.5.4

      Zint_c=Zint_c*spice_bundle_model%bundle%conductor_impedance(conductor1)%Resistance_multiplication_factor
                                              
! evaluate the surface impedance for the reference conductor 
      CALL evaluate_conductor_impedance_model(spice_bundle_model%bundle%conductor_impedance(reference_conductor1),  &
                                              f,Zint_c_ref,Rdc_c_ref,Zint_t_ref,Rdc_t_ref)

! Apply multiplication factor to the conductor impedance to correct for common mode/ differential modes in twisted pairs
      Zint_c_ref=Zint_c_ref*spice_bundle_model%bundle%conductor_impedance(reference_conductor1)%Resistance_multiplication_factor
              
! The surface impedance of the conductor and the reference conductor contribute to the diagonal element
      Z_domain_conductor_impedance_correction(domain_conductor1,domain_conductor1)=Zint_c+Zint_c_ref
      
      if (verbose) then
        write(*,*)'conductor:',conductor1,' reference conductor',reference_conductor1
        write(*,*)'conductor loss model type',spice_bundle_model%bundle%conductor_impedance(conductor1)%impedance_model_type
        write(*,*)'refconductor loss model type',  &
                  spice_bundle_model%bundle%conductor_impedance(reference_conductor1)%impedance_model_type
        write(*,*)'radius',  &
                  spice_bundle_model%bundle%conductor_impedance(reference_conductor1)%radius
        write(*,*)'conductivity',  &
                  spice_bundle_model%bundle%conductor_impedance(reference_conductor1)%conductivity
        write(*,*)'Resistance_multiplication_factor',  &
                  spice_bundle_model%bundle%conductor_impedance(reference_conductor1)%Resistance_multiplication_factor
        write(*,*)'thickness',  &
                  spice_bundle_model%bundle%conductor_impedance(reference_conductor1)%thickness

        write(*,*)'domain conductor:',domain_conductor1,' domain reference conductor',domain_reference_conductor1
        write(*,*)'Contribution to Zc(',domain_conductor1,domain_conductor1,')'
        write(*,*)'Zc conductor:',Zint_c
        write(*,*)'Zc reference:',Zint_c_ref
      end if ! verbose
      
! conductor always gets the contribution from its own surface impedance
    
      do conductor2=conductor1+1,dim
   
        domain2=spice_bundle_model%bundle%terminal_conductor_to_outer_domain(conductor2)
        reference_conductor2=spice_bundle_model%bundle%terminal_conductor_to_reference_terminal_conductor(conductor2)
        domain_conductor2=spice_bundle_model%bundle%terminal_conductor_to_global_domain_conductor(conductor2)
        domain_reference_conductor2=spice_bundle_model%bundle%terminal_conductor_to_global_domain_conductor(reference_conductor2)
        is_shield2=spice_bundle_model%bundle%terminal_conductor_is_shield_flag(conductor2)

! if the two conductors belong to the same domain then add the reference conductor impedance
        if (domain1.EQ.domain2) then
          Z_domain_conductor_impedance_correction(domain_conductor1,domain_conductor2)=Zint_c_ref
          Z_domain_conductor_impedance_correction(domain_conductor2,domain_conductor1)=Zint_c_ref
          
          if (verbose) then
            write(*,*)'Contribution to Zc(',domain_conductor1,domain_conductor2,')'
            write(*,*)'Contribution to Zc(',domain_conductor2,domain_conductor1,')'
            write(*,*)'Zc conductor:',Zint_c_ref
          end if ! verbose
          
       end if
      
      end do ! next conductor2
      
    end do ! next conductor1
    
    if (verbose) write(*,*)'Add transfer impedance contributions'
! add transfer impedance contributions

! loop over conductors looking for shields. Note include all conductors including the reference here
    do conductor=1,dim+1
   
      is_shield1=spice_bundle_model%bundle%terminal_conductor_is_shield_flag(conductor)
      
      if (is_shield1) then
! add transfer impedance contributions to inner and outer domain conductors
 
        inner_domain=spice_bundle_model%bundle%terminal_conductor_to_inner_domain(conductor)
        outer_domain=spice_bundle_model%bundle%terminal_conductor_to_outer_domain(conductor)
        
        CALL evaluate_conductor_impedance_model(spice_bundle_model%bundle%conductor_impedance(conductor),  &
                                                f,Zint_c,Rdc_c,Zint_t,Rdc_t)
                                                
! Check whether the shield is the reference conductor in the outer domain - the contributions
! are different if this is the case.

        n_conductors_outer_domain=spice_bundle_model%bundle%n_conductors(outer_domain)
        shield_conductor_number_in_outer_domain=spice_bundle_model%bundle%terminal_conductor_to_local_domain_conductor(conductor)
         
! number of conductors in a domain is spice_bundle_model%bundle%n_conductors(domain)
        if (shield_conductor_number_in_outer_domain.NE.n_conductors_outer_domain) then
               
! loop over all conductors
          do row=1,dim
         
! get the domain of row conductor
            domain1=spice_bundle_model%bundle%terminal_conductor_to_outer_domain(row)

            if (domain1.EQ.inner_domain) then
! The row conductor is in the inner shield domain and so gets a transfer impedance contribution from the shield conductor
                 
! the shield couples these two domains so add the transfer impedance term - also include term to make the matrix symmetric
              domain_conductor1=spice_bundle_model%bundle%terminal_conductor_to_global_domain_conductor(row)
                 
              domain_conductor2=spice_bundle_model%bundle%terminal_conductor_to_global_domain_conductor(conductor)
              Z_domain_conductor_impedance_correction(domain_conductor1,domain_conductor2)=        &
              Z_domain_conductor_impedance_correction(domain_conductor1,domain_conductor2) -Zint_t
              Z_domain_conductor_impedance_correction(domain_conductor2,domain_conductor1)=        &
              Z_domain_conductor_impedance_correction(domain_conductor2,domain_conductor1) -Zint_t
                 
              if (verbose) then
                write(*,*)'Shield conductor',conductor,' inner domain',inner_domain,' outer domain',outer_domain
                write(*,*)'row',row,' col',col,' row domain',domain1,' col domain',domain2
                write(*,*)'Contribution to Zt(',domain_conductor1,domain_conductor2,')'
                write(*,*)'Contribution to Zt(',domain_conductor2,domain_conductor1,')'
                write(*,*)'Zt conductor:',-Zint_t
              end if ! verbose
                                                 
            end if  ! transfer impedance term required
            
          end do ! next row conductor
          
        else ! shield IS reference conductor in outer domain
               
! loop over all conductors
          do row=1,dim
         
! get the domain of row conductor
            domain1=spice_bundle_model%bundle%terminal_conductor_to_outer_domain(row)

            if (domain1.EQ.inner_domain) then
! The row conductor is in the inner shield domain and so gets a transfer impedance contribution from the shield conductor
                 
! the shield couples these two domains so add the transfer impedance term - also include term to make the matrix symmetric
              domain_conductor1=spice_bundle_model%bundle%terminal_conductor_to_global_domain_conductor(row)
                 
! As the shield conductor is the reference we need to find all the conductors contributing to the shield current
! note that the contribution is then -ve of the normal transfer impedance contribution as the currents are in the
! opposite direction

              do col=1,dim
              
                domain2=spice_bundle_model%bundle%terminal_conductor_to_outer_domain(col)
! Check the domain of the col conductor. If it is an outer domain conductor of the shield then it contributes

                if (domain2.EQ.outer_domain) then

                  domain_conductor2=spice_bundle_model%bundle%terminal_conductor_to_global_domain_conductor(col)
                  
                  Z_domain_conductor_impedance_correction(domain_conductor1,domain_conductor2)=        &
                  Z_domain_conductor_impedance_correction(domain_conductor1,domain_conductor2) +Zint_t
                  Z_domain_conductor_impedance_correction(domain_conductor2,domain_conductor1)=        &
                  Z_domain_conductor_impedance_correction(domain_conductor2,domain_conductor1) +Zint_t
                 
                  if (verbose) then
                    write(*,*)'Shield conductor',conductor,' inner domain',inner_domain,' outer domain',outer_domain
                    write(*,*)'row',row,' col',col,' row domain',domain1,' col domain',domain2
                    write(*,*)'Contribution to Zt(',domain_conductor1,domain_conductor2,')'
                    write(*,*)'Contribution to Zt(',domain_conductor2,domain_conductor1,')'
                    write(*,*)'Zt conductor:',-Zint_t
                  end if ! verbose
                  
                end if  ! transfer impedance term required for this col conductor
                 
              end do ! next condutor to check
                                                
            end if  ! transfer impedance term required for this row conductor
            
          end do ! next row conductor
         
        end if ! shield is/ is not reference conductor in outer domain
      
      end if ! conductor is a shield

    end do ! next conductor
    
! Add the conductor impedance contributions to the domain based impedance matrix
    Z_domain(:,:)=Z_domain(:,:)+Z_domain_conductor_impedance_correction(:,:)
    
    Einc=cmplx(spice_bundle_model%Eamplitude)
    
! Solve the frequency domain multi-conductor transmission line equations with the specified termination circuit and 
! return the required conductor voltage in Vout.   
! For the purposes of the frequency domain solution, the excitations are assumed to be constant in frequency
! The frequency response of the excitation is included afterwards         
    if (.NOT.run_validation_test_Vbased) then
                                          
      CALL frequency_domain_MTL_solution(dim,Z_domain,Y_domain,MV,MVI,MI,MII, &
                                       Einc,spice_bundle_model%Ex,spice_bundle_model%Ey,spice_bundle_model%Ez,    &
                                       spice_bundle_model%Hx,spice_bundle_model%Hy,spice_bundle_model%Hz,         &
                                       spice_bundle_model%kx,spice_bundle_model%ky,spice_bundle_model%kz,         &
                                       local_conductor_x_offset,                                                  &
                                       local_conductor_y_offset,                                                  &
                                       spice_bundle_model%bundle%ground_plane_present,                            &
                                       spice_bundle_model%bundle%ground_plane_x,                                  &
                                       spice_bundle_model%bundle%ground_plane_y,                                  &
                                       spice_bundle_model%bundle%ground_plane_theta,                              &
                                       spice_bundle_model%length,Vs1,Z1,Vs2,Z2,                                   &
                                       is_shielded_flag,                                                          &
                                       f,spice_validation_test%output_end,spice_validation_test%output_conductor, &
                                       spice_validation_test%output_conductor_ref,Vout)
    else
                                        
      CALL frequency_domain_MTL_solution_V(dim,Z_domain,Y_domain,MV,MVI,MI,MII, &
                                       Einc,spice_bundle_model%Ex,spice_bundle_model%Ey,spice_bundle_model%Ez,    &
                                       spice_bundle_model%Hx,spice_bundle_model%Hy,spice_bundle_model%Hz,         &
                                       spice_bundle_model%kx,spice_bundle_model%ky,spice_bundle_model%kz,         &
                                       local_conductor_x_offset,                                                  &
                                       local_conductor_y_offset,                                                  &
                                       spice_bundle_model%bundle%ground_plane_present,                            &
                                       spice_bundle_model%bundle%ground_plane_x,                                  &
                                       spice_bundle_model%bundle%ground_plane_y,                                  &
                                       spice_bundle_model%bundle%ground_plane_theta,                              &
                                       spice_bundle_model%length,Vs1,Z1,Vs2,Z2,                                   &
                                       is_shielded_flag,                                                          &
                                       f,spice_validation_test%output_end,spice_validation_test%output_conductor, &
                                       spice_validation_test%output_conductor_ref,Vout)
    end if
    
! multiply the Frequency domain source and transmission line transfer function 
! Save transfer function multiplied by frequency domain source voltage function i.e. include the time response of
! the excitation here. Theory_Manual_Section 2.3.3

    VTL_freq(frequency_loop)=Vout*Vs_freq(frequency_loop)
    
    if ( (frequency_loop.NE.1).AND.(frequency_loop.NE.N_FFT/2+1) ) then
      VTL_freq(N_FFT-frequency_loop+2)=conjg(VTL_freq(frequency_loop))
    end if
    
  end do ! next frequency in frequency loop

! Inverse FFT to give the time domain transmission line response
! Theory_Manual_Equation 2.50

  CALL FFT_FREQ_TO_TIME(n_FFT,time,VTL_time,freq,VTL_freq)

! write the time domain response to file

! Open output file 
  open(unit=analytic_soln_file_unit,file=trim(analytic_soln_filename))
  
! write the file header line
  write(analytic_soln_file_unit,'(A)')time_header

  do i=1,n_timesteps                                              ! changed from n_FFT 04/5/2016 CJS
    write(analytic_soln_file_unit,*)time(i),real(VTL_time(i))
  end do
  
! Close output file 
  Close(unit=analytic_soln_file_unit) 

! deallocate memory
  DEALLOCATE( Z_domain )
  DEALLOCATE( Y_domain )
  
  DEALLOCATE( Z_domain_conductor_impedance_correction )
  
  DEALLOCATE( Vs_time ) 
  DEALLOCATE( Vs_freq ) 
  DEALLOCATE( VTL_time ) 
  DEALLOCATE( VTL_freq ) 
  DEALLOCATE( time ) 
  DEALLOCATE( freq ) 
  
! domain transformation matrices
  DEALLOCATE( MI )
  DEALLOCATE( MII )
  DEALLOCATE( MV )
  DEALLOCATE( MVI )
    
  DEALLOCATE( is_shielded_flag )
  DEALLOCATE( local_conductor_x_offset )
  DEALLOCATE( local_conductor_y_offset )

! temporary working matrices
  DEALLOCATE( TM1 )
  
  write(6,*)

  RETURN

END SUBROUTINE time_domain_analysis