!
! This file is part of SACAMOS, State of the Art CAble MOdels for Spice. 
! It was developed by the University of Nottingham and the Netherlands Aerospace 
! Centre (NLR) for ESA under contract number 4000112765/14/NL/HK.
! 
! Copyright (C) 2016-2018 University of Nottingham
! 
! SACAMOS is free software: you can redistribute it and/or modify it under the 
! terms of the GNU General Public License as published by the Free Software 
! Foundation, either version 3 of the License, or (at your option) any later 
! version.
! 
! SACAMOS is distributed in the hope that it will be useful, but 
! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 
! or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
! for more details.
! 
! A copy of the GNU General Public License version 3 can be found in the 
! file GNU_GPL_v3 in the root or at <http://www.gnu.org/licenses/>.
! 
! SACAMOS uses the EISPACK library (in /SRC/EISPACK). EISPACK is subject to 
! the GNU Lesser General Public License. A copy of the GNU Lesser General Public 
! License version can be found in the file GNU_LGPL in the root of EISPACK 
! (/SRC/EISPACK ) or at <http://www.gnu.org/licenses/>.
! 
! The University of Nottingham can be contacted at: ggiemr@nottingham.ac.uk
!
!
! File Contents:
! SUBROUTINE frequency_domain_analysis
!
! NAME
!     frequency_domain_analysis
!
! AUTHORS
!     Chris Smartt
!
! DESCRIPTION
!     This subroutine controls the analytic solution for the frequency domain analysis of
!     multi-conductor transmission lines. 
!     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 frequency domain termination voltage for the specified validation test case written to file
!     
! 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
!     Include general conductor impedance model 12/05/2016 CJS
!     Fix bug with conductor impedance contributions 12/05/2016 CJS
!     25/08/2016 CJS Include revised transfer impedance/ condcutor impedance model for shields
!     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 frequency_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

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(:,:)

! temporary variables
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    :: row,col
integer    :: i

integer :: ierr    ! error code for matrix inverse calls

! START

! Open output file 
  open(unit=analytic_soln_file_unit,file=trim(analytic_soln_filename))
  
! write the file header line
  if (spice_validation_test%analysis_freq_spec%freq_range_type.EQ.'log') then
    write(analytic_soln_file_unit,'(A)')log_freq_header
  else if (spice_validation_test%analysis_freq_spec%freq_range_type.EQ.'lin') then
    write(analytic_soln_file_unit,'(A)')lin_freq_header
  end if
  
  dim=spice_bundle_model%bundle%system_dimension
  
! allocate memory
  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 of 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(:,:)=cmplx(spice_bundle_model%bundle%global_MI%mat(:,:))
  MV(:,:)=cmplx(spice_bundle_model%bundle%global_MV%mat(:,:))
  
  if (verbose) write(*,*)'Invert MI'
  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)
  
  if (verbose) then
    write(*,*)'Transpose[MII]'
    do row=1,dim
      write(*,8000)(real(MII(col,row)),col=1,dim)
    end do
  
    write(*,*)'[MV]'
    do row=1,dim
      write(*,8000)(real(MV(row,col)),col=1,dim)
    end do
8000 format(20F4.1)
  
  end if ! verbose
  
  if (verbose) write(*,*)'Invert MV'
  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 the specified frequencies
  do frequency_loop=1,spice_validation_test%analysis_freq_spec%n_frequencies
    
! get the frequency and angular frequency values
    f=spice_validation_test%analysis_freq_spec%freq_list(frequency_loop)
    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. Initailly set to zero
! See Theory_Manual_Section 2.2.3 

    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
      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,'Rdc=',Rdc_c
        write(*,*)'Zc reference:',Zint_c_ref,'Rdc=',Rdc_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) then
        
      write(*,*)'[R_domain]=Re[Z_domain]'
      do row=1,dim
        write(*,8020)(real(Z_domain_conductor_impedance_correction(row,col)),col=1,dim)
      end do
    
    end if ! verbose 
    
    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(:,:)
    
    if (verbose) then
        
      write(*,*)'[R_domain]=Re[Z_domain]'
      do row=1,dim
        write(*,8020)(real(Z_domain(row,col)),col=1,dim)
      end do
    
    end if ! verbose 
    
    if (verbose) then
    
      write(*,*)'Im[Z_domain]'
      do row=1,dim
        write(*,8010)(aimag(Z_domain(row,col)),col=1,dim)
      end do
    
      write(*,*)'[R_domain]=Re[Z_domain]'
      do row=1,dim
        write(*,8020)(real(Z_domain(row,col)),col=1,dim)
      end do
  
      write(*,*)'Im[Y_domain]'
      do row=1,dim
        write(*,8010)(aimag(Y_domain(row,col)),col=1,dim)
      end do
8010 format(20ES10.2)
8020 format(20F12.4)
  
    end if ! verbose 

! Get the incident field amplitude    
    Einc=cmplx(spice_bundle_model%Eamplitude)

! Solve the frequency domain multi-conductor transmission line equations with the specified termination circuit and 
! incident field excitation, return the required conductor voltage in Vout.        
    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
! Output the result to file  

    if (spice_validation_test%output_type.EQ.'li') then
      
      if (plot_real) then
        write(analytic_soln_file_unit,*)f,real(Vout),aimag(Vout)    
      else
        write(analytic_soln_file_unit,*)f,abs(Vout),atan2(aimag(Vout),real(Vout))
      end if
       
    else if (spice_validation_test%output_type.EQ.'dB') then
    
      write(analytic_soln_file_unit,*)f,20d0*log10(abs(Vout))
      
    end if   ! output format (linear or dB)

  end do ! next frequency in frequency loop
      
! Close output file 
  Close(unit=analytic_soln_file_unit) 

! deallocate memory
  DEALLOCATE( Z_domain )
  DEALLOCATE( Y_domain )
  
  DEALLOCATE( Z_domain_conductor_impedance_correction )
  
  DEALLOCATE( Vs1 )
  DEALLOCATE( Z1 )
  DEALLOCATE( Vs2 )
  DEALLOCATE( Z2 )
  
! 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 )

  RETURN

END SUBROUTINE frequency_domain_analysis