Use FastHenry2 impedance matrices in analytic MTL solution

Chris Smartt
1 parent ca8ab9e9
Showing 6 changed files with 828 additions and 7 deletions Show diff stats
SRC/MTL_ANALYTIC_SOLUTION/MTL_analytic_solution.F90
SRC/MTL_ANALYTIC_SOLUTION/Makefile
SRC/MTL_ANALYTIC_SOLUTION/frequency_domain_analysis_FH2.F90
SRC/compare_results.F90
SRC/spice_cable_bundle_model_builder.F90
SRC/write_FH_input_file.F90
@@ -29,6 +29,7 @@
 ! FILE CONTENTS (within include files)
 !
 !frequency_domain_analysis.F90:     SUBROUTINE frequency_domain_analysis
+!frequency_domain_analysis_FH2.F90: SUBROUTINE frequency_domain_analysis_FH2
 !time_domain_analysis.F90:          SUBROUTINE time_domain_analysis
 !frequency_domain_MTL_solution.F90: SUBROUTINE frequency_domain_MTL_solution
 !incident_field_excitation.F90:     SUBROUTINE calc_incident_field_components
@@ -73,6 +74,8 @@ include &#39;modal_decomposition.F90&#39;
 include 'frequency_domain_analysis.F90'
+include 'frequency_domain_analysis_FH2.F90'
+
 include 'time_domain_analysis.F90'
 include 'frequency_domain_MTL_solution.F90'
 default: $(MTL_ANALYTIC_SOLUTION_OBJS)
 #
 $(OBJ_MOD_DIR)/%.o: %.F90 $(TYPE_SPEC_MODULE)  $(CONSTANTS_MODULE)  $(EISPACK_MODULE)  $(MATHS_MODULE)\
-                    frequency_domain_analysis.F90 frequency_domain_MTL_solution.F90 \
+                    frequency_domain_analysis.F90 frequency_domain_analysis_FH2.F90 frequency_domain_MTL_solution.F90 \
                     modal_decomposition_LC.F90 modal_decomposition.F90 \
                     time_domain_analysis.F90 propagation_correction_filters.F90 incident_field_excitation.F90
@@ -0,0 +1,644 @@
+!
+! 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)
+
+! 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
@@ -51,6 +51,8 @@
 ! HISTORY
 !
 !     started 24/11/2015 CJS 
+!     26/10/2023 CJS    don't send exit code 1 if there is an error. This allows results to be plotted using
+!                       the generate_spice_cable_bundle_model process
 !
 ! _________________________________________________________________
 !
@@ -432,7 +434,8 @@ IMPLICIT NONE
     run_status='ERROR: Sample not found in file 2'
     CALL write_program_status()
-    STOP 1
+!    STOP 1
+    GOTO 2000
 1000  CONTINUE
@@ -447,6 +450,8 @@ IMPLICIT NONE
     local_compare_results_value=local_compare_results_value+    &
       error*( ((dataset1%x(sample+1))-(dataset1%x(sample-1)) )/2d0)     &
             /( (dataset1%x(sample_max+1))-(dataset1%x(sample_min-1)) )
+            
+2000 CONTINUE
   end do ! next sample in the integration over x
@@ -571,6 +571,9 @@ integer    :: i
     if (INDEX(line,'no_high_freq_zt_model').NE.0) high_freq_Zt_model=.FALSE.
     if (INDEX(line,'plot_real').NE.0) plot_real=.TRUE.
     if (INDEX(line,'plot_mag').NE.0) plot_real=.FALSE.
+
+    if (INDEX(line,'no_fasthenry').NE.0) use_FastHenry=.FALSE.   
+    if (INDEX(line,'use_fasthenry').NE.0) use_FastHenry=.TRUE.
     if (INDEX(line,'use_analytic_i').NE.0) run_validation_test_Vbased=.FALSE.
     if (INDEX(line,'use_analytic_v').NE.0) run_validation_test_Vbased=.TRUE.
@@ -667,10 +670,20 @@ integer    :: i
     if (spice_validation_test%analysis_type.EQ.analysis_type_AC) then
-      if (verbose) write(*,*)'CALLING: frequency_domain_analysis'
-      CALL frequency_domain_analysis(spice_bundle_model,spice_validation_test)
-      if (verbose) write(*,*)'FINISHED: frequency_domain_analysis'
-    
+      if (use_FastHenry) then
+  
+        if (verbose) write(*,*)'CALLING: frequency_domain_analysis using FastHenry2 impedance matrices'
+        CALL frequency_domain_analysis_FH2(spice_bundle_model,spice_validation_test)
+        if (verbose) write(*,*)'FINISHED: frequency_domain_analysis using FastHenry2 impedance matrices'
+      
+      else
+  
+        if (verbose) write(*,*)'CALLING: frequency_domain_analysis'
+        CALL frequency_domain_analysis(spice_bundle_model,spice_validation_test)
+        if (verbose) write(*,*)'FINISHED: frequency_domain_analysis'
+        
+      end if
+      
     else
       if (verbose) write(*,*)'CALLING: time_domain_analysis'
@@ -16,6 +16,9 @@ real(dp)    :: gp_t,gp_sigma,gp_rh,gp_ex1,gp_ey1,gp_ez1,gp_ex2,gp_ey2,gp_ez2
 real(dp)    :: gp_skin_depth
 integer :: gp_nhinc,gp_seg1,gp_seg2
+real(dp) :: gp_xmin,gp_xmax
+real(dp) :: gp_xc,gp_yc,gp_w,gp_h
+
 character(len=80),allocatable    :: gp_node1
 character(len=80),allocatable    :: gp_node2
@@ -35,6 +38,14 @@ TYPE::conductor_type
   real(dp)    :: yc
   real(dp)    :: width
   real(dp)    :: height
+  
+  real(dp)    :: rss
+  real(dp)    :: rss_outer
+  real(dp)    :: xcss(7)
+  real(dp)    :: ycss(7)
+  real(dp)    :: rot_angle
+  integer     :: n_layers2ss
+  
   integer :: n_layers2
   integer :: tot_n_layers
@@ -74,6 +85,7 @@ integer,parameter :: type_cyl=1
 integer,parameter :: type_rect=2
 integer,parameter :: type_annulus=3
 integer,parameter :: type_gnd=4
+integer,parameter :: type_seven_strand=5
 integer,parameter :: mesh_type_layer=1
 integer,parameter :: mesh_type_grid=2
@@ -121,6 +133,10 @@ character :: type_ch
 integer :: tot_n_segments,tot_n_filaments
 integer :: gp_n_segments,gp_n_filaments
+integer :: i
+
+real(dp) :: xpt,ypt
+
 real(dp),parameter :: pi=3.1415926535
 ! START
@@ -235,7 +251,33 @@ INCLUDE &quot;WRITE_FH2_IPFILE/get_grid_type.F90&quot;
     write(*,*)'Enter the cylindrical conductor discretisation, dl in metres'
     line=line+1
     read(*,*,ERR=9000)conductor_data(conductor)%dl
-        
+    
+  else if ( (type_ch.EQ.'s').OR.(type_ch.EQ.'S') ) then
+    conductor_data(conductor)%type=type_seven_strand   
+
+! work out the grid type    
+
+INCLUDE "WRITE_FH2_IPFILE/get_grid_type.F90"
+
+    write(*,*)Conductor,conductor,' mesh type=',conductor_data(conductor)%mesh_type
+
+    write(*,*)'Enter the seven strand conductor centre coordinates, xc yc in metres'
+    line=line+1
+    read(*,*,ERR=9000)conductor_data(conductor)%xc,conductor_data(conductor)%yc
+    
+    write(*,*)'Enter the seven strand conductor equivalent radius, rc in metres'
+    line=line+1
+    read(*,*,ERR=9000)conductor_data(conductor)%rc
+    
+    write(*,*)'Enter the seven strand conductor rotation angle in degrees'
+    line=line+1
+    read(*,*,ERR=9000)conductor_data(conductor)%rot_angle
+    conductor_data(conductor)%rot_angle=conductor_data(conductor)%rot_angle*pi/180d0
+    
+    write(*,*)'Enter the seven strand conductor discretisation, dl in metres'
+    line=line+1
+    read(*,*,ERR=9000)conductor_data(conductor)%dl
+       
   else if ( (type_ch.EQ.'r').OR.(type_ch.EQ.'R') ) then
     conductor_data(conductor)%type=type_rect
@@ -376,6 +418,71 @@ INCLUDE &quot;WRITE_FH2_IPFILE/create_grid.F90&quot;
       write(*,*)'Unknown grid type'
       STOP 1
     end if
+  
+  else if (conductor_data(conductor)%type.EQ.type_seven_strand) then
+
+! radius of each strand for the same conductor area  
+    conductor_data(conductor)%rss=conductor_data(conductor)%rc/sqrt(7.0)
+! maximum radius of the combined conductors
+    conductor_data(conductor)%rss_outer=conductor_data(conductor)%rss*3.0
+
+! seven conductor centre coordinates
+! central conductor
+    conductor_data(conductor)%xcss(1)=conductor_data(conductor)%xc
+    conductor_data(conductor)%ycss(1)=conductor_data(conductor)%yc
+    do i=1,6
+      angle=real(i-1)*2.0*pi/6.0+conductor_data(conductor)%rot_angle
+      conductor_data(conductor)%xcss(i+1)=conductor_data(conductor)%xc+2d0*conductor_data(conductor)%rss*cos(angle)
+      conductor_data(conductor)%ycss(i+1)=conductor_data(conductor)%yc+2d0*conductor_data(conductor)%rss*sin(angle)    
+    end do
+    
+    if (conductor_data(conductor)%mesh_type.EQ.mesh_type_layer) then
+    
+      conductor_data(conductor)%n_layers2ss=NINT(conductor_data(conductor)%rss/conductor_data(conductor)%dl)
+      conductor_data(conductor)%tot_n_layers=7*2*conductor_data(conductor)%n_layers2ss
+    
+    else if (conductor_data(conductor)%mesh_type.EQ.mesh_type_grid) then
+    
+! allocate grid for the circular geometry
+      conductor_data(conductor)%grid_dim=conductor_data(conductor)%rss_outer
+     
+INCLUDE "WRITE_FH2_IPFILE/create_grid.F90"
+         
+! Loop over the grid and set segments within the conductor
+    
+      conductor_data(conductor)%tot_n_layers=0
+
+! loop over 7 strands      
+      do i=1,7
+
+! loop over grid    
+        do ix=nxmin,nxmax
+          do iy=nymin,nymax
+          
+! calculate dpt, the distance to the centre of strand i
+            xpt=conductor_data(conductor)%dl*real(ix)- &
+                (conductor_data(conductor)%xcss(i)-conductor_data(conductor)%xcss(1))
+            ypt=conductor_data(conductor)%dl*real(iy)- &
+                (conductor_data(conductor)%ycss(i)-conductor_data(conductor)%ycss(1))
+            
+            dpt=sqrt(xpt**2+ypt*2)
+            
+            if ( dpt.LE.conductor_data(conductor)%rss ) then
+              conductor_data(conductor)%grid(ix,iy)=1
+              conductor_data(conductor)%depth(ix,iy)=conductor_data(conductor)%rss-dpt
+              conductor_data(conductor)%tot_n_layers=conductor_data(conductor)%tot_n_layers+1
+            end if
+          end do
+        end do
+    
+        conductor_data(conductor)%n_layers2=0
+
+      end do ! next strand
+
+    else
+      write(*,*)'Unknown grid type'
+      STOP 1
+    end if
   else if (conductor_data(conductor)%type.EQ.type_annulus) then
@@ -457,6 +564,17 @@ if (ground_plane) then
   write(20,'(9A)')'+ ',trim(gp_node2),' (',trim(adjustl(x_string)),  &
                                        ',',trim(adjustl(y_string)),  &
                        ',',trim(adjustl(z_string)),')'
+  
+  gp_xmin=min(gp_x(1),gp_x(2),gp_x(3))
+  gp_xmax=max(gp_x(1),gp_x(2),gp_x(3))
+  gp_xc=(gp_xmin+gp_xmax)/2d0
+  gp_yc=gp_y(1)
+  gp_w=(gp_xmax-gp_xmin)
+  gp_h=gp_t
+  
+  CALL plot_layer(gp_xc,gp_yc,gp_w,gp_h,1d0,0d0,10)
+    
+  CALL plot_grid(gp_xc,gp_yc,gp_w,gp_h,1d0,0d0,1d0,gp_rh,gp_seg1,gp_nhinc,12)
 end if
@@ -520,6 +638,44 @@ do conductor=1,n_conductors
 INCLUDE "WRITE_FH2_IPFILE/set_segments_from_grid.F90"
     end if   ! mesh_type_grid
+  
+  else if (conductor_data(conductor)%type.EQ.type_seven_strand) then 
+    
+    if (conductor_data(conductor)%mesh_type.EQ.mesh_type_layer) then
+    
+      layer_number=0
+
+! loop over 7 strands      
+      do i=1,7
+
+        do layer=-conductor_data(conductor)%n_layers2ss+1,conductor_data(conductor)%n_layers2ss
+  
+          layer_number=layer_number+1
+    
+          ymin=conductor_data(conductor)%dl*real(layer)
+          ymax=conductor_data(conductor)%dl*real(layer-1)
+          y=(ymin+ymax)/2.0
+          x=sqrt(conductor_data(conductor)%rss**2-y**2)
+    
+          conductor_data(conductor)%x(layer_number)=conductor_data(conductor)%xcss(i)
+          conductor_data(conductor)%y(layer_number)=conductor_data(conductor)%ycss(i)+y
+          conductor_data(conductor)%w(layer_number)=2.0*x
+          conductor_data(conductor)%h(layer_number)=conductor_data(conductor)%dl
+          conductor_data(conductor)%d(layer_number)=0.0
+          conductor_data(conductor)%anwinc(layer_number)=conductor_data(conductor)%nwinc
+          conductor_data(conductor)%anhinc(layer_number)=conductor_data(conductor)%nwinc
+   
+        end do  ! next layer
+   
+      end do  ! next strand
+        
+    else if (conductor_data(conductor)%mesh_type.EQ.mesh_type_grid) then
+    
+! Loop over the grid and set segments within the circular conductor
+
+INCLUDE "WRITE_FH2_IPFILE/set_segments_from_grid.F90"
+     
+    end if   ! mesh_type_grid
   else if (conductor_data(conductor)%type.EQ.type_rect) then