!
! 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 PUL_RL_FastHenry2
!
! NAME
! SUBROUTINE PUL_RL_FastHenry2
!
! DESCRIPTION
! Wrapping subroutine to control the calculation of frequency dependent impedance
! using FastHenry2
!
! The process is divided into the following stages:
! STAGE 1: work out the configuration for the calculation i.e. is there a ground plane, is the outer boundary free space or a conductor (overshield)
! STAGE 2: Create the input file for FastHenry2
! STAGE 4: Call FastHenry2
! STAGE 5: Read the frequency dependent impedance matrices and fit rational functions to these
!
! COMMENTS
!
!
! HISTORY
! started 23/10/2023
!
SUBROUTINE PUL_RL_FastHenry2(PUL,name,fit_order,freq_spec,domain)
!
USE type_specifications
USE constants
USE general_module
USE maths
USE filter_module
USE filter_module
USE Sfilter_fit_module
USE frequency_spec
!
IMPLICIT NONE
! variables passed to subroutine
type(PUL_type), intent(INOUT) :: PUL ! per-unit-length parameter calculation structure
character(LEN=line_length),intent(IN) :: name ! string used as the base for filenames
integer, intent(IN) :: fit_order ! filter fit_order for frequency dependent dielectrics
type(frequency_specification),intent(IN) :: freq_spec ! filter frequency range specification for frequency dependent dielectrics
integer, intent(IN) :: domain ! domain number used to label the mesh files
! parameters
integer,parameter :: inside =1
integer,parameter :: outside=2
! local variables
integer :: first_conductor,last_conductor,conductor
character(LEN=filename_length) :: command ! string used for running external commands
integer :: exit_stat
! Process_Zc variables:
integer :: matrix_dimension
complex(dp) :: Zc_temp
complex(dp),allocatable :: function_to_fit(:) ! complex function to be fitted using Sfilter_fit
real(dp) :: ferr
real(dp) :: gp_half_width
character(LEN=256) :: line
character(LEN=256) :: freq_and_dim_string
integer :: local_line_length
integer :: i
integer :: n_freq_in
real(dp) :: freq,w
integer :: nrows,ncols,r,c,row,col
complex(dp),allocatable :: Zc_in(:,:,:)
real(dp),allocatable :: freq_in(:)
real(dp),allocatable :: Rdc_domain(:,:)
real(dp),allocatable :: values(:),re,im
logical :: subtract_dc_resistance
integer :: loop,floop
integer :: ierr
! START
if (verbose) write(*,*)'CALLED: PUL_RL_FastHenry2'
! work out the local frequency sampling: FastHenry uses a specific form based on log frequency
! and number of samples per decade
write(*,*)'Filter fit order=',fit_order
write(*,*)'fmin=',freq_spec%fmin
write(*,*)'fmax=',freq_spec%fmax
write(*,*)'n_frequencies=',freq_spec%n_frequencies
write(*,*)'ndec=',freq_spec%ndec
if (freq_spec%freq_range_type.EQ.'lin') then
write(run_status,*)'ERROR in PUL_RL_FastHenry2: We must have a logarithmic frequency range specified'
CALL write_program_status()
STOP 1
end if
! STAGE 1: work out the configuration for the calculation i.e. is there a ground plane, is the outer boundary free space or a conductor (overshield)
! write the file which is used to generate to FastHenry2 input file
OPEN(unit=fh2_input_file_unit,file='fh2.txt')
write(fh2_input_file_unit,'(A)')'fh2.inp'
if ((.NOT.PUL%overshield_present).AND.(PUL%ground_plane_present)) then
! SOLUTION TYPE 1: NO OVERSHIELD, WITH GROUND PLANE
first_conductor=1
last_conductor=PUL%n_conductors-1
write(fh2_input_file_unit,'(A)')'ground_plane'
! ********* STILL SOME HARD WIRED GROUND PLANE SPEC: TO BE SORTED OUT ************
! ALSO NEED TO CHECK THAT EVERYTHING WE NEED IS SET...
write(*,*)'Ground plane parameters:'
write(*,*)'PUL%ground_plane_w=',PUL%ground_plane_w
write(*,*)'PUL%ground_plane_h=',PUL%ground_plane_h
write(*,*)'PUL%ground_plane_sigma=',PUL%ground_plane_sigma
write(*,*)'PUL%ground_plane_nsegx=',PUL%ground_plane_nsegx
write(*,*)'PUL%ground_plane_nsegz=',PUL%ground_plane_nsegz
write(*,*)'PUL%ground_plane_nh=',PUL%ground_plane_nh
if ( (PUL%ground_plane_w.EQ.0d0).OR.(PUL%ground_plane_h.EQ.0d0).OR.(PUL%ground_plane_sigma.EQ.0d0) &
.OR.(PUL%ground_plane_nsegx.EQ.0).OR.(PUL%ground_plane_nsegz.EQ.0).OR.(PUL%ground_plane_nh.EQ.0) )then
write(run_status,*)'ERROR in PUL_RL_FastHenry2: ground plane parameters not defined'
CALL write_program_status()
STOP 1
end if
gp_half_width=PUL%ground_plane_w/2.0
write(fh2_input_file_unit,'(2ES12.4,A)')-gp_half_width,-PUL%ground_plane_h/2.0,' 0.0e-3 ! x y z of ground plane point 1'
write(fh2_input_file_unit,'(2ES12.4,A)') gp_half_width,-PUL%ground_plane_h/2.0,' 0.0e-3 ! x y z of ground plane point 2'
write(fh2_input_file_unit,'(2ES12.4,A)') gp_half_width,-PUL%ground_plane_h/2.0,' 1000.0e-3 ! x y z of ground plane point 3'
write(fh2_input_file_unit,'(ES12.4,A)')PUL%ground_plane_h,' ! thickness '
write(fh2_input_file_unit,'(2I6,A)')PUL%ground_plane_nsegx,PUL%ground_plane_nsegz, &
' ! discretisation in p1-p2 and p2=p3 directions'
write(fh2_input_file_unit,'(ES12.4,A)')PUL%ground_plane_sigma,' ! conductivity, sigma '
write(fh2_input_file_unit,'(I4,A)')PUL%ground_plane_nh,' ! gp discretisation in thickness, nhinc '
write(fh2_input_file_unit,'(A)')'2.0 ! gp discretisation in thickness ratio, rh '
write(fh2_input_file_unit,'(A)')'0e-3 0.0e-3 0.0e-3 ! x y z of ground plane node at end 1 '
write(fh2_input_file_unit,'(A)')'0e-3 0.0e-3 1000.0e-3 ! x y z of ground plane node at end 2 '
! ********* STILL SOME HARD WIRED GROUND PLANE SPEC: TO BE SORTED OUT ************
else if ((PUL%overshield_present).AND.(.NOT.PUL%ground_plane_present)) then
! SOLUTION TYPE 2: OVERSHIELD, NO GROUND PLANE
first_conductor=1
last_conductor=PUL%n_conductors
else if ((.NOT.PUL%overshield_present).AND.(.NOT.PUL%ground_plane_present)) then
! SOLUTION TYPE 3: NO OVERSHIELD, NO GROUND PLANE
first_conductor=1
last_conductor=PUL%n_conductors
else
write(run_status,*)'ERROR in PUL_RL_FastHenry2 cannot handle this configuration'
CALL write_program_status()
STOP 1
end if
write(fh2_input_file_unit,'(I2,A)')last_conductor-first_conductor+1,' ! number of conductors'
do conductor=first_conductor,last_conductor
write(fh2_input_file_unit,'(I2,A)')conductor,' ! conductor number'
if ( (conductor.EQ.last_conductor).AND.PUL%overshield_present) then
write(fh2_input_file_unit,'(A)')'annulus ! conductor 1 shape (cylinder, rectangle, annulus)'
write(fh2_input_file_unit,'(2ES12.4,A)')PUL%overshield_x,PUL%overshield_y,' ! centre xc yc'
write(fh2_input_file_unit,'(ES12.4,A)')PUL%overshield_r,' ! inner radius rc'
write(fh2_input_file_unit,'(ES12.4,A)')PUL%overshield_r*1.01,' ! ***** outer radius rc: 1.01*rinner *****'
write(fh2_input_file_unit,'(ES12.4,A)')PUL%overshield_r/2,' ! ***** conductor discretisation size, dl=r/2 *****'
write(fh2_input_file_unit,'(ES12.4,A)')1e20,' ! **** conductivity, sigma: high conductivity ****'
! write(fh2_input_file_unit,'(2I4,A)')FH2_nw,FH2_nh,' ! segment discretisation in x and y, nwinc nhinc '
! write(fh2_input_file_unit,'(2ES12.4,A)')FH2_rw,FH2_rh,' ! segment discretisation ratio in x and y, rw rh '
write(fh2_input_file_unit,'(A)')' 3 1 ! segment discretisation in x and y, nwinc nhinc '
write(fh2_input_file_unit,'(A)')'1.0 1.0 ! segment discretisation ratio in x and y, rw rh '
else ! not an overshield
if (PUL%shape(conductor).EQ.circle) then
if (PUL%conductor_is_shield(conductor)) then
! ********* USED FOR SHIELDS IN THE EXTERNAL DOMAIN: STILL SOME HARD WIRED PARAMETERS ************
write(fh2_input_file_unit,'(A)')'annulus ! conductor 1 shape (cylinder, rectangle, annulus)'
write(fh2_input_file_unit,'(2ES12.4,A)')PUL%x(conductor)+PUL%ox(conductor),&
PUL%y(conductor)+PUL%oy(conductor),' ! centre xc yc'
write(fh2_input_file_unit,'(ES12.4,A)')PUL%r(conductor)/1.01d0,' ! inner radius rs/1.01'
write(fh2_input_file_unit,'(ES12.4,A)')PUL%r(conductor),' ! outer radius rs'
write(fh2_input_file_unit,'(ES12.4,A)')PUL%overshield_r/2,' ! ***** conductor discretisation size, dl=r/2 *****'
write(fh2_input_file_unit,'(ES12.4,A)')1e20,' ! **** conductivity, sigma: high conductivity ****'
write(fh2_input_file_unit,'(A)')' 3 1 ! segment discretisation in x and y, nwinc nhinc '
write(fh2_input_file_unit,'(A)')'1.0 1.0 ! segment discretisation ratio in x and y, rw rh '
else
if(auto_cyl_grid) then
write(fh2_input_file_unit,'(A)')'cylindrical grid_block auto ! conductor 1 shape (cylinder, rectangle, annulus)'
else
write(fh2_input_file_unit,'(A)')'cylindrical ! conductor 1 shape (cylinder, rectangle, annulus)'
end if
write(fh2_input_file_unit,'(2ES12.4,A)')PUL%x(conductor)+PUL%ox(conductor),&
PUL%y(conductor)+PUL%oy(conductor),' ! centre xc yc'
write(fh2_input_file_unit,'(ES12.4,A)')PUL%r(conductor),' ! radius rc'
write(fh2_input_file_unit,'(ES12.4,A)')PUL%r(conductor)/FH2_nlayers_radius, &
' ! cylindrical conductor discretisation size, dl'
write(fh2_input_file_unit,'(ES12.4,A)')PUL%sigma(conductor),' ! conductivity, sigma '
write(fh2_input_file_unit,'(2I4,A)')FH2_nw,FH2_nh,' ! segment discretisation in x and y, nwinc nhinc '
write(fh2_input_file_unit,'(2ES12.4,A)')FH2_rw,FH2_rh,' ! segment discretisation ratio in x and y, rw rh '
end if ! conductor_is_shield
else if (PUL%shape(conductor).EQ.rectangle) then
write(fh2_input_file_unit,'(A)')'rectangle ! conductor 1 shape (cylinder, rectangle, annulus)'
write(fh2_input_file_unit,'(2ES12.4,A)')PUL%x(conductor)+PUL%ox(conductor),&
PUL%y(conductor)+PUL%oy(conductor),' ! centre xc yc'
write(fh2_input_file_unit,'(2ES12.4,A)')PUL%rw(conductor),PUL%rh(conductor),' ! width and height'
write(fh2_input_file_unit,'(ES12.4,A)')PUL%sigma(conductor),' ! conductivity, sigma '
write(fh2_input_file_unit,'(2I4,A)')FH2_nw,FH2_nh,' ! segment discretisation in x and y, nwinc nhinc '
write(fh2_input_file_unit,'(2ES12.4,A)')FH2_rw,FH2_rh,' ! segment discretisation ratio in x and y, rw rh '
else
write(run_status,*)'ERROR in PUL_RL_FastHenry2 cannot handle shape:',PUL%shape(conductor)
CALL write_program_status()
STOP 1
end if ! shape
end if ! overshield
end do ! next conductor
write(fh2_input_file_unit,'(ES12.4,A)')freq_spec%fmin,' ! Minimum frequency, fmin'
write(fh2_input_file_unit,'(ES12.4,A)')freq_spec%fmax,' ! Maximum frequency, fmax'
write(fh2_input_file_unit,'(ES12.4,A)')freq_spec%ndec,' ! Number of points per decade, ndec'
CLOSE(unit=fh2_input_file_unit)
write(*,*)'CALLING write_FH_input_file'
command='write_FH_input_file < fh2.txt'
write(*,*)'COMMAND:'
write(*,'(A)')command
CALL execute_command_line(command,EXITSTAT=exit_stat)
! Check that the FastHenry2 input file generated correctly
if (exit_stat.NE.0) then
! write_FH_input_file returned with a non zdero exit code indicating an error
write(run_status,*)'ERROR in PUL_RL_FastHenry2 running write_FH_input_file: exit_stat=',exit_stat
CALL write_program_status()
STOP 1
end if
write(*,*)'CALLING FastHenry2'
if(fh2_no_refinement) then
command='fasthenry -a off fh2.inp'
else
command='fasthenry fh2.inp'
end if
write(*,*)'COMMAND:'
write(*,'(A)')command
CALL execute_command_line(command,EXITSTAT=exit_stat)
! Check that the fasthenry ran correctly
if (exit_stat.NE.0) then
! fasthenry returned with a non zdero exit code indicating an error
write(run_status,*)'ERROR in PUL_RL_FastHenry2 running fasthenry: exit_stat=',exit_stat
CALL write_program_status()
STOP 1
end if
command='cp fh2.txt fh2.txt_'//trim(name)
CALL execute_command_line(command,EXITSTAT=exit_stat)
command='cp fh2.inp fh2.inp_'//trim(name)
CALL execute_command_line(command,EXITSTAT=exit_stat)
command='cp Zc.mat Zc.mat_'//trim(name)
CALL execute_command_line(command,EXITSTAT=exit_stat)
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))
! Changed to avoid compilation warning
! freq_and_dim_string=line(33:local_line_length)
freq_and_dim_string(1:local_line_length-32)=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( Rdc_domain(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)
matrix_dimension=PUL%n_conductors-1
if (matrix_dimension.NE.nrows) then
write(run_status,*)'ERROR in PUL_RL_FastHenry2 matrix dimension',matrix_dimension,nrows
CALL write_program_status()
STOP 1
end if
! Make the impedance matrix explicitly symmetric
do floop=1,n_freq_in
do r=1,nrows
do c=r+1,ncols
Zc_temp=(Zc_in(floop,r,c)+Zc_in(floop,c,r))/2d0
Zc_in(floop,r,c)=Zc_temp
Zc_in(floop,c,r)=Zc_temp
end do
end do
end do
if (L_from_fasthenry) then
! For the modal decomposition, use the inductance from the highest frequency calculated
! ******************* TO BE USED WITH CARE DUE TO OBSERVED FASTER THAN LIGHT MODES AND *************
! ********* MODE BEATING IN SOME PLANE PROBLEMS E.G. TWO WIRES OVER GND IN FREE SPACE *********
do r=1,nrows
do c=1,ncols
freq=freq_in(n_freq_in)
w=2.0*pi*freq
PUL%L%mat(r,c)=real(Zc_in(n_freq_in,r,c)/(j*w))
end do
end do
end if
! Save the resistance matrix from the lowest frequency calculated
do r=1,nrows
do c=1,ncols
Rdc_domain(r,c)=real(Zc_in(1,r,c))
end do
end do
! We calculate the inductance and resistance matrix at a number of frequencies before fitting filter functions
! to each of the frequency dependent impedance matrix entries.
write(*,*)'**************************'
write(*,*)'Process FastHenry2 results'
write(*,*)'**************************'
write(*,*)'Number of frequencies (freq_spec) =',freq_spec%n_frequencies
write(*,*)'Number of frequencies (FastHenry2)=',n_freq_in
! Frequency check
if (n_freq_in.NE.freq_spec%n_frequencies) then
write(run_status,*)'ERROR in PUL_RL_FastHenry2 n_frequencies specification',n_freq_in,freq_spec%n_frequencies
CALL write_program_status()
STOP 1
end if
do floop=1,n_freq_in
ferr=abs(freq_spec%freq_list(floop)-freq_in(floop))/(freq_spec%freq_list(floop)+freq_in(floop))
if (ferr.GT.0.001) then
write(run_status,*)'ERROR in PUL_RL_FastHenry2 frequency samples, floop=',floop, &
' freq_spec=',freq_spec%freq_list(floop),' freq_FH2=',freq_in(floop)
CALL write_program_status()
STOP 1
end if
end do
! loop over frequency
if (verbose) then
do floop=1,n_freq_in
write(*,*)freq_spec%freq_list(floop),freq_in(floop)
end do
end if
! we have frequency domain impedance matrix, Zc
! loop over the elements of Zc and the Zfilter matrix by filter fitting
ALLOCATE( function_to_fit(1:freq_spec%n_frequencies) )
do row=1,matrix_dimension
do col=row,matrix_dimension
! get the function values for this matrix element function_to_fit=R-Rdc+jwL
do floop=1,freq_spec%n_frequencies
function_to_fit(floop)=Zc_in(floop,row,col)-Rdc_domain(row,col)
end do
! calculate the Zfilter matrix using the filter fitting process
! note aorder=border and no restrictions are put on the function
CALL Calculate_Sfilter(function_to_fit,freq_spec%freq_list,freq_spec%n_frequencies, &
PUL%Zfilter%sfilter_mat(row,col),fit_order+1,-1,2) ! note type 2 filter fit=> a0=0.0
if (col.NE.row) then
! make the matrix symmetrical
PUL%Zfilter%sfilter_mat(col,row)=PUL%Zfilter%sfilter_mat(row,col)
end if
end do ! next col
end do ! next row
DEALLOCATE( function_to_fit )
DEALLOCATE( freq_in )
DEALLOCATE( Zc_in )
DEALLOCATE( Rdc_domain )
DEALLOCATE( values )
if (verbose) then
write(*,*)'FINISHED: PUL_RL_FastHenry2'
write(*,*)'Inductance matrix, L'
CALL dwrite_matrix(PUL%L%mat,nrows,ncols,nrows,0)
end if
if (verbose) write(*,*)'FINISHED PUL_RL_FastHenry2'
! STOP 1
END SUBROUTINE PUL_RL_FastHenry2