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