diff --git a/SRC/MTL_ANALYTIC_SOLUTION/MTL_analytic_solution.F90 b/SRC/MTL_ANALYTIC_SOLUTION/MTL_analytic_solution.F90
index d4dea31..83c12d3 100644
--- a/SRC/MTL_ANALYTIC_SOLUTION/MTL_analytic_solution.F90
+++ b/SRC/MTL_ANALYTIC_SOLUTION/MTL_analytic_solution.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 'modal_decomposition.F90'
include 'frequency_domain_analysis.F90'
+include 'frequency_domain_analysis_FH2.F90'
+
include 'time_domain_analysis.F90'
include 'frequency_domain_MTL_solution.F90'
diff --git a/SRC/MTL_ANALYTIC_SOLUTION/Makefile b/SRC/MTL_ANALYTIC_SOLUTION/Makefile
index 7d68cab..02c7b5d 100644
--- a/SRC/MTL_ANALYTIC_SOLUTION/Makefile
+++ b/SRC/MTL_ANALYTIC_SOLUTION/Makefile
@@ -1,7 +1,7 @@
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
diff --git a/SRC/MTL_ANALYTIC_SOLUTION/frequency_domain_analysis_FH2.F90 b/SRC/MTL_ANALYTIC_SOLUTION/frequency_domain_analysis_FH2.F90
new file mode 100644
index 0000000..b64fd14
--- /dev/null
+++ b/SRC/MTL_ANALYTIC_SOLUTION/frequency_domain_analysis_FH2.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 .
+!
+! 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 .
+!
+! 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
diff --git a/SRC/compare_results.F90 b/SRC/compare_results.F90
index f37d737..87c8f15 100644
--- a/SRC/compare_results.F90
+++ b/SRC/compare_results.F90
@@ -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
diff --git a/SRC/spice_cable_bundle_model_builder.F90 b/SRC/spice_cable_bundle_model_builder.F90
index 46f58dd..7ef84c5 100644
--- a/SRC/spice_cable_bundle_model_builder.F90
+++ b/SRC/spice_cable_bundle_model_builder.F90
@@ -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'
diff --git a/SRC/write_FH_input_file.F90 b/SRC/write_FH_input_file.F90
index 1ec161d..d646f0d 100644
--- a/SRC/write_FH_input_file.F90
+++ b/SRC/write_FH_input_file.F90
@@ -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 "WRITE_FH2_IPFILE/get_grid_type.F90"
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 "WRITE_FH2_IPFILE/create_grid.F90"
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
--
libgit2 0.21.2