!
! 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_FH2
!
! NAME
!     frequency_domain_analysis_FH2
!
! 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
!
!     26/10/2023      CJS: use FastHenry impedance matrix in analytic solution
!
SUBROUTINE frequency_domain_analysis_FH2(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

! FastHenry2 stuff...

character(LEN=256) :: line
character(LEN=256) :: freq_and_dim_string

integer :: local_line_length

integer :: n_freq_in

real(dp) :: freq
integer :: nrows,ncols,r,c

complex(dp),allocatable :: Zc_in(:,:,:)
real(dp),allocatable    :: freq_in(:)
real(dp),allocatable    :: values(:),re,im

integer :: loop,floop


! 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)
  
   
  if (verbose) then
    write(*,*)
    write(*,*)' Processing frequency dependent impedance file from FastHenry2'
    write(*,*)
  end if
  
  OPEN(unit=fh2_output_file_unit,file='Zc.mat')  
  
  do loop=1,2

    n_freq_in=0

10  CONTINUE
    read(fh2_output_file_unit,'(A)',END=1000,ERR=1000)line
    if (line(1:32).NE.'Impedance matrix for frequency =') GOTO 10

!     if we get here then we have some impedance matrix data to read
!      write(*,*)'Found impedance matrix:',trim(line)

! increase the number of frequencies for which we have an impedance matrix    
      n_freq_in=n_freq_in+1

      local_line_length=len(trim(line))
      freq_and_dim_string=line(33:local_line_length)
  
!      write(*,*)'reading matrix for:',trim(freq_and_dim_string)

! replace the 'x' by a space then read the frequency, nrows and ncols
      local_line_length=len(trim(freq_and_dim_string)) 
      do i=1,local_line_length
        if (freq_and_dim_string(i:i).EQ.'x') freq_and_dim_string(i:i)=' '
      end do
  
!      write(*,*)'data string:',trim(freq_and_dim_string)
  
      read(freq_and_dim_string,*)freq,nrows,ncols
      if (nrows.NE.ncols) then
        write(*,*)'****** ERROR: nrows.NE.ncols ******'
      end if
!      write(*,*)'frequency=',freq,' n=',nrows
    
      if (loop.EQ.2) then
        freq_in(n_freq_in)=freq
      end if
        
      do r=1,nrows
  
        read(fh2_output_file_unit,'(A)',END=1000,ERR=1000)line

        if (loop.EQ.2) then
      
!   replace the 'j's by a space then read the complex data
        local_line_length=len(trim(line)) 
          
        do i=1,local_line_length
          if (line(i:i).EQ.'j') line(i:i)=' '
        end do
  
        read(line,*)(values(i),i=1,2*ncols)
    
        do c=1,ncols
    
          re=values(c*2-1)
          im=values(c*2)
                
          Zc_in(n_freq_in,r,c)=cmplx(re,im)
	
        end do ! next column
	 
      end if
  
    end do ! next row
     
    GOTO 10 ! continue to read the file

! Jump here when the file has been read

1000 CONTINUE
 
    if (loop.Eq.1) then
      if (verbose) write(*,*)'Number of frequencies=',n_freq_in
      ALLOCATE( freq_in(n_freq_in) )
      ALLOCATE( Zc_in(n_freq_in,nrows,ncols) )
      ALLOCATE( values(2*ncols) )
      rewind(unit=fh2_output_file_unit)
    end if
  
  end do ! next file read loop
  
  CLOSE(unit=fh2_output_file_unit)
  
! Check dimensions
  
  if (nrows.NE.dim) then
    write(*,*)'ERROR in frequency_domain_analysis_FH2: nrows.NE.dim'
    write(*,*)'nrows from FH2 =',nrows,' dim=',dim
    write(run_status,*)'ERROR in frequency_domain_analysis_FH2: dimension error:',dim,nrows
    CALL write_program_status()
    STOP 1  
  end if
 
! ************* LOOP OVER FREQUENCY   ***************

! Loop over the specified frequencies
  do frequency_loop=1,n_freq_in
    
! get the frequency and angular frequency values
    f=freq_in(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)=Zc_in(frequency_loop,row,col)

        
!! TEMP TEST FOR STP CABLES: CALCULATE DM AND CM IMPEDANCE MATRIX
!        Z_domain(1,1)=2d0*(Zc_in(frequency_loop,1,1)-Zc_in(frequency_loop,1,2))
!        Z_domain(1,2)=0d0
!        Z_domain(2,1)=0d0
!        Z_domain(2,2)=0.5d0*(Zc_in(frequency_loop,1,1)+Zc_in(frequency_loop,1,2))

! 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

! Contribution to the impedance matrices from the conductor based impedance models is zero here

    Z_domain_conductor_impedance_correction(1:dim,1:dim)=(0d0,0d0)
    
    if (verbose) write(*,*)'Add transfer impedance contributions'
    
! add transfer impedance contributions   *********** NOT SURE WHETHER WE NEED TO DO ANYTHING HERE *********

! 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 )
 
  DEALLOCATE( freq_in )
  DEALLOCATE( Zc_in )
  DEALLOCATE( values )

  RETURN

END SUBROUTINE frequency_domain_analysis_FH2