!
! 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:
! MODULE cable_bundle_module
! CONTAINS
! SUBROUTINE read_cable_bundle
! SUBROUTINE write_cable_bundle
! SUBROUTINE deallocate_cable_bundle
! SUBROUTINE check_cables_wrt_ground_plane
! SUBROUTINE check_cable_intersections
!
! NAME
! MODULE cable_bundle_module
!
! Data strunctures and subroutines relating to cable bundles
!
! COMMENTS
! Need to improve error handling in file reading
! Not sure that some arrays in the type_specification need to be in this structure as they are only used in one subroutine...
!
! HISTORY
! started 25/11/2015 CJS
! stage 3 developments started 21/3/2015 CJS
! generalise conductor impedance for transfer impedance model (stage 5) CJS 12/5/2016
! CJS 25/08/2016 updates for revised transfer impedance and conductor impedance model for shields
! 16/11/2017 CJS Include network synthesis process to replace s-domain transfer functions
!
MODULE cable_bundle_module
USE type_specifications
USE cable_module
IMPLICIT NONE
! This is a general structure used to hold conductor based lists,
TYPE::conductor_list_type
integer :: n_elements
integer,allocatable :: element(:)
END TYPE conductor_list_type
! This is the main structure used to hold all the bundle information
TYPE::bundle_specification_type
character(LEN=line_length) :: version
character(LEN=line_length) :: bundle_name
integer :: tot_n_conductors
integer :: tot_n_shielded_conductors
integer :: tot_n_external_conductors
! dimension of the global matrix system (total number of conductors-1)
integer :: system_dimension
! cable information, total number of cables including and excluding the ground plane (if it exists)
integer :: n_cables
integer :: n_cables_without_ground_plane
! Array of cable specification data for all the cables in this bundle
type(cable_specification_type),allocatable :: cable(:)
! list of cable positions in the bundle cross section and rotation angles
real(dp),allocatable :: cable_x_offset(:)
real(dp),allocatable :: cable_y_offset(:)
real(dp),allocatable :: cable_angle(:)
! ground plane information: note we have two notations at the moment.
! This will need to be sorted out and one only adopted
logical :: ground_plane_present
real(dp) :: ground_plane_angle
real(dp) :: ground_plane_offset
real(dp) :: ground_plane_sigma ! ground plane conductivity for FastHenry2
real(dp) :: ground_plane_w ! ground plane width for FastHenry2
real(dp) :: ground_plane_h ! ground plane height for FastHenry2
real(dp) :: ground_plane_Rdc ! ground plane dc resistance for FastHenry2
integer :: ground_plane_nsegx ! ground plane number of segments in x for FastHenry2
integer :: ground_plane_nsegz ! ground plane number of segments in z for FastHenry2
integer :: ground_plane_nh ! ground plane number of layers in h for FastHenry2
real(dp) :: ground_plane_x,ground_plane_y,ground_plane_theta
real(dp) :: ground_plane_nx,ground_plane_ny
integer :: ground_plane_cable_side ! +1 if the cables are on the normal direction side of the ground plane, -1 otherwise
! domain based information
integer :: tot_n_domains
integer,allocatable :: n_conductors(:) ! number of conductors in each domain
type(conductor_list_type),allocatable :: terminal_conductor_list(:) ! external conductor number for coonductors in each domain
type(matrix),allocatable :: L(:) ! 'high frequency' inductance matrix for each domain
type(matrix),allocatable :: C(:) ! 'high frequency' capacitance matrix for each domain
type(Sfilter_matrix),allocatable :: Z(:) ! frequency dependent impedance matrix for each domain
type(Sfilter_matrix),allocatable :: Y(:) ! frequency dependent admittance matrix for each domain
! global conductor based information
type(matrix) :: global_MI ! global domain current transformation matrix
type(matrix) :: global_MV ! global domain voltage transformation matrix
type(matrix) :: global_L ! global 'high frequency' inductance matrix
type(matrix) :: global_C ! global 'high frequency' capacitance matrix
type(Sfilter_matrix) :: global_Z ! global frequency dependent impedance matrix
type(Sfilter_matrix) :: global_Y ! global frequency dependent admittance matrix
! global conductor impedance model information for each conductor
type(conductor_impedance_model),allocatable :: conductor_impedance(:)
! global conductor position in the bundle cross-section
real(dp),allocatable :: conductor_x_offset(:)
real(dp),allocatable :: conductor_y_offset(:)
! arrays for cross referecing numbering information from the global (external) conductor numbering and the domain based numbering systems
logical,allocatable :: terminal_conductor_is_shield_flag(:)
integer,allocatable :: terminal_conductor_to_inner_domain(:)
integer,allocatable :: terminal_conductor_to_outer_domain(:)
integer,allocatable :: terminal_conductor_to_global_domain_conductor(:)
integer,allocatable :: terminal_conductor_to_local_domain_conductor(:)
integer,allocatable :: terminal_conductor_to_reference_terminal_conductor(:)
! Y matrix element function fitting information for frequency dependent dielectrics
! model order for filter fitting
integer :: Y_fit_model_order
! frequency range specification
type(frequency_specification) :: Y_fit_freq_spec
character(LEN=line_length),allocatable :: conductor_label(:)
END TYPE bundle_specification_type
CONTAINS
! The following subroutines apply to all bundles and cable types
! NAME
! SUBROUTINE read_cable_bundle(unit)
!
! read cable bundle structure from a specified unit
!
! COMMENTS
!
!
! HISTORY
! started 2/12/2015 CJS
! stage 3 developments started 21/3/2015 CJS
! generalise conductor impedance for transfer impedance model (stage 5) CJS 12/5/2016
!
SUBROUTINE read_cable_bundle(bundle,file_unit)
USE type_specifications
USE constants
USE general_module
USE cable_module
USE maths
IMPLICIT NONE
! variables passed to subroutine
type(bundle_specification_type),intent(INOUT) :: bundle
integer ,intent(IN) :: file_unit
! local variables
character(len=filename_length) :: filename
logical :: file_exists
character(len=line_length) :: line
integer :: cable
integer :: n_domains
integer :: domain
integer :: matrix_dimension
integer :: n_conductors
integer :: conductor
integer :: ierr
! START
! Open the (.bundle) file to read
filename=trim(MOD_bundle_lib_dir)//trim(bundle%bundle_name)//bundle_file_extn
inquire(file=trim(filename),exist=file_exists)
if (.NOT.file_exists) then
run_status='ERROR, Cannot find the required bundle file:'//trim(filename)
CALL write_program_status()
STOP 1
end if
open(unit=file_unit,file=trim(filename))
if (verbose) write(*,*)'Opened file:',trim(filename)
! read .bundle file
read(file_unit,'(A)',ERR=9000)bundle%version
read(file_unit,'(A)',ERR=9000)bundle%bundle_name
! read cable information
read(file_unit,*,ERR=9000)bundle%n_cables_without_ground_plane
read(file_unit,*,ERR=9000)bundle%n_cables
! allocate and read cable information and cable positions in bundle
ALLOCATE( bundle%cable(1:bundle%n_cables) )
ALLOCATE( bundle%cable_x_offset(1:bundle%n_cables) )
ALLOCATE( bundle%cable_y_offset(1:bundle%n_cables) )
ALLOCATE( bundle%cable_angle(1:bundle%n_cables) )
do cable=1,bundle%n_cables_without_ground_plane
read(file_unit,'(A)')bundle%cable(cable)%cable_name
CALL read_cable(bundle%cable(cable),cable_file_unit)
read(file_unit,*,ERR=9000)bundle%cable_x_offset(cable),bundle%cable_y_offset(cable),bundle%cable_angle(cable)
! convert angle to radians
bundle%cable_angle(cable)=bundle%cable_angle(cable)*pi/180d0
end do ! next cable
! read ground plane specification
read(file_unit,'(A)',IOSTAT=ierr)line
if (ierr.NE.0) then
run_status='ERROR reading ground plane present/ absent information'
CALL write_program_status()
STOP 1
end if
CALL convert_to_lower_case(line,line_length)
bundle%ground_plane_present=(line.eq.'ground_plane')
if (bundle%ground_plane_present) then
read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_x,bundle%ground_plane_y,bundle%ground_plane_angle
read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_nx,bundle%ground_plane_ny
read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_cable_side
if (ierr.NE.0) then
run_status='ERROR reading ground plane position and angle'
CALL write_program_status()
STOP 1
end if
! convert ground plane angle to radians
bundle%ground_plane_angle=bundle%ground_plane_angle*pi/180d0
! calculate offset
bundle%ground_plane_offset=-bundle%ground_plane_x*sin(bundle%ground_plane_angle)+ &
bundle%ground_plane_y*cos(bundle%ground_plane_angle)
read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_sigma
read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_w
read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_h
read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_Rdc
read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_nsegx
read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_nsegz
read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_nh
if (ierr.NE.0) then
run_status='ERROR reading ground plane sigma, w h and Rdc'
CALL write_program_status()
STOP 1
end if
! set the geometric data for the last cable in the list i.e. the ground plane
cable=bundle%n_cables
bundle%cable_x_offset(cable)=bundle%ground_plane_x
bundle%cable_y_offset(cable)=bundle%ground_plane_y
bundle%cable_angle(cable)=bundle%ground_plane_angle-pi/2d0 ! note there are two ground plane spec formats at the moment...
else
! set x, y, angle and offset to 0
bundle%ground_plane_x=0d0
bundle%ground_plane_y=0d0
bundle%ground_plane_nx=0d0
bundle%ground_plane_ny=0d0
bundle%ground_plane_cable_side=0
bundle%ground_plane_angle=0d0
bundle%ground_plane_offset=0d0
end if
! read system dimension information
read(file_unit,*)bundle%tot_n_conductors
read(file_unit,*)bundle%tot_n_external_conductors
read(file_unit,*)bundle%system_dimension
! read domain information
read(file_unit,*,ERR=9000)bundle%tot_n_domains
ALLOCATE( bundle%n_conductors(1:bundle%tot_n_domains) )
ALLOCATE( bundle%L(1:bundle%tot_n_domains) )
ALLOCATE( bundle%C(1:bundle%tot_n_domains) )
ALLOCATE( bundle%Z(1:bundle%tot_n_domains))
ALLOCATE( bundle%Y(1:bundle%tot_n_domains))
ALLOCATE( bundle%terminal_conductor_list(1:bundle%tot_n_domains))
n_domains=bundle%tot_n_domains
do domain=1,bundle%tot_n_domains
read(file_unit,*,ERR=9000) ! comment line
read(file_unit,*)bundle%n_conductors(domain)
read(file_unit,*,ERR=9000) ! comment line
read(file_unit,*,ERR=9000)bundle%terminal_conductor_list(domain)%n_elements
ALLOCATE( bundle%terminal_conductor_list(domain)%element(1:bundle%terminal_conductor_list(domain)%n_elements) )
do conductor=1,bundle%terminal_conductor_list(domain)%n_elements
read(file_unit,*,ERR=9000)bundle%terminal_conductor_list(domain)%element(conductor)
end do
read(file_unit,*,ERR=9000) ! comment line
read(file_unit,*)matrix_dimension
bundle%L(domain)%dim=matrix_dimension
ALLOCATE( bundle%L(domain)%mat(1:matrix_dimension,1:matrix_dimension) )
CALL dread_matrix(bundle%L(domain)%mat,matrix_dimension,matrix_dimension,matrix_dimension,file_unit)
write(*,*)'L='
CALL dwrite_matrix(bundle%L(domain)%mat,matrix_dimension,matrix_dimension,matrix_dimension,0) ! ******
read(file_unit,*,ERR=9000) ! comment line
read(file_unit,*)matrix_dimension
bundle%C(domain)%dim=matrix_dimension
ALLOCATE( bundle%C(domain)%mat(1:matrix_dimension,1:matrix_dimension) )
CALL dread_matrix(bundle%C(domain)%mat,matrix_dimension,matrix_dimension,matrix_dimension,file_unit)
write(*,*)'C='
CALL dwrite_matrix(bundle%C(domain)%mat,matrix_dimension,matrix_dimension,matrix_dimension,0) ! ******
read(file_unit,*,ERR=9000) ! comment line
CALL read_Sfilter_matrix( bundle%Z(domain),file_unit )
read(file_unit,*,ERR=9000) ! comment line
CALL read_Sfilter_matrix( bundle%Y(domain),file_unit )
end do
read(file_unit,*,ERR=9000) ! comment line
read(file_unit,*)matrix_dimension
bundle%global_MI%dim=matrix_dimension
ALLOCATE( bundle%global_MI%mat(1:matrix_dimension,1:matrix_dimension) )
CALL dread_matrix(bundle%global_MI%mat,matrix_dimension,matrix_dimension,matrix_dimension,file_unit)
read(file_unit,*,ERR=9000) ! comment line
read(file_unit,*)matrix_dimension
bundle%global_MV%dim=matrix_dimension
ALLOCATE( bundle%global_MV%mat(1:matrix_dimension,1:matrix_dimension) )
CALL dread_matrix(bundle%global_MV%mat,matrix_dimension,matrix_dimension,matrix_dimension,file_unit)
read(file_unit,*,ERR=9000) ! comment line
read(file_unit,*)matrix_dimension
bundle%global_L%dim=matrix_dimension
ALLOCATE( bundle%global_L%mat(1:matrix_dimension,1:matrix_dimension) )
CALL dread_matrix(bundle%global_L%mat,matrix_dimension,matrix_dimension,matrix_dimension,file_unit)
read(file_unit,*,ERR=9000) ! comment line
read(file_unit,*)matrix_dimension
bundle%global_C%dim=matrix_dimension
ALLOCATE( bundle%global_C%mat(1:matrix_dimension,1:matrix_dimension) )
CALL dread_matrix(bundle%global_C%mat,matrix_dimension,matrix_dimension,matrix_dimension,file_unit)
read(file_unit,*,ERR=9000) ! comment line
CALL read_Sfilter_matrix( bundle%global_Z,file_unit )
read(file_unit,*,ERR=9000) ! comment line
CALL read_Sfilter_matrix( bundle%global_Y,file_unit )
! read the loss model information for each conductor
ALLOCATE( bundle%conductor_impedance(1:bundle%tot_n_conductors) )
read(file_unit,*) ! comment line
do conductor=1,bundle%tot_n_conductors
CALL read_conductor_impedance_model(bundle%conductor_impedance(conductor),file_unit)
end do
! read the x and y position for each conductor: this is needed for incident field excitation
ALLOCATE( bundle%conductor_x_offset(1:bundle%tot_n_conductors) )
ALLOCATE( bundle%conductor_y_offset(1:bundle%tot_n_conductors) )
read(file_unit,*) ! comment line
do conductor=1,bundle%tot_n_conductors
read(file_unit,*)bundle%conductor_x_offset(conductor),bundle%conductor_y_offset(conductor)
end do
! numbering information
read(file_unit,*) ! comment line
ALLOCATE( bundle%terminal_conductor_is_shield_flag(1:bundle%tot_n_conductors) )
ALLOCATE( bundle%terminal_conductor_to_inner_domain(1:bundle%tot_n_conductors) )
ALLOCATE( bundle%terminal_conductor_to_outer_domain(1:bundle%tot_n_conductors) )
ALLOCATE( bundle%terminal_conductor_to_global_domain_conductor(1:bundle%tot_n_conductors) )
ALLOCATE( bundle%terminal_conductor_to_local_domain_conductor(1:bundle%tot_n_conductors) )
ALLOCATE( bundle%terminal_conductor_to_reference_terminal_conductor(1:bundle%tot_n_conductors) )
do conductor=1,bundle%tot_n_conductors
read(file_unit,*)bundle%terminal_conductor_is_shield_flag(conductor), &
bundle%terminal_conductor_to_inner_domain(conductor), &
bundle%terminal_conductor_to_outer_domain(conductor), &
bundle%terminal_conductor_to_global_domain_conductor(conductor), &
bundle%terminal_conductor_to_local_domain_conductor(conductor), &
bundle%terminal_conductor_to_reference_terminal_conductor(conductor)
end do
! read the conductor labels
ALLOCATE( bundle%conductor_label(1:bundle%tot_n_conductors) )
read(file_unit,*) ! comment line
do conductor=1,bundle%tot_n_conductors
read(file_unit,'(A)')bundle%conductor_label(conductor)
end do
! close .bundle file
CLOSE(unit=file_unit)
if (verbose) write(*,*)'Closed file:',trim(filename)
RETURN
9000 run_status='ERROR reading the bundle file:'//trim(filename)
CALL write_program_status()
STOP 1
END SUBROUTINE read_cable_bundle
! NAME
! SUBROUTINE write_cable_bundle
!
! write the cable bundle structure to a specified unit
!
! COMMENTS
!
!
! HISTORY
! started 2/12/2015 CJS
! stage 3 developments started 21/3/2015 CJS
! generalise conductor impedance for transfer impedance model (stage 5) CJS 12/5/2016
!
SUBROUTINE write_cable_bundle(bundle,file_unit)
USE type_specifications
USE constants
USE general_module
USE cable_module
USE maths
IMPLICIT NONE
! variables passed to subroutine
type(bundle_specification_type),intent(INOUT) :: bundle
integer,intent(IN) :: file_unit
! local variables
character(len=filename_length) :: filename
character(len=line_length) :: line
integer :: cable
integer :: domain
integer :: matrix_dimension
integer :: conductor
integer :: n_cables_to_write
! START
! Open the output (.bundle) file
filename=trim(MOD_bundle_lib_dir)//trim(bundle%bundle_name)//bundle_file_extn
open(unit=file_unit,file=trim(filename))
! write .bundle file
write(file_unit,'(A)')trim(bundle%version)
write(file_unit,'(A)')trim(bundle%bundle_name)
write(file_unit,*)bundle%n_cables_without_ground_plane,' ! number of cables not including ground plane'
write(file_unit,*)bundle%n_cables,' ! number of cables, cable name and x y coordinates follow...'
! write cable names and coordinates
do cable=1,bundle%n_cables_without_ground_plane
write(file_unit,'(A)')trim(bundle%cable(cable)%cable_name)
write(file_unit,*)bundle%cable_x_offset(cable),bundle%cable_y_offset(cable), &
bundle%cable_angle(cable),' x y coordinates and angle of cable '
end do ! next cable
! write ground plane specification
if (bundle%ground_plane_present) then
write(file_unit,'(A)')'ground_plane'
! note: convert angle back to degrees
write(file_unit,*)bundle%ground_plane_x,bundle%ground_plane_y,bundle%ground_plane_angle*180d0/pi, &
' x y coordinates and angle of ground plane '
write(file_unit,*)bundle%ground_plane_nx,bundle%ground_plane_ny,' ground plane normal direction'
write(file_unit,*)bundle%ground_plane_cable_side,' orientation of cables wrt ground plane'
write(file_unit,*)bundle%ground_plane_sigma,' ground plane conductivity'
write(file_unit,*)bundle%ground_plane_w,' ground plane width '
write(file_unit,*)bundle%ground_plane_h,' ground plane height (thickness) '
write(file_unit,*)bundle%ground_plane_Rdc,' ground plane Rdc '
write(file_unit,*)bundle%ground_plane_nsegx,' ground plane nsegx '
write(file_unit,*)bundle%ground_plane_nsegz,' ground plane nsegz '
write(file_unit,*)bundle%ground_plane_nh,' ground plane nh '
else
write(file_unit,'(A)')'no_ground_plane'
end if
write(file_unit,*)bundle%tot_n_conductors,' # total number of conductors'
write(file_unit,*)bundle%tot_n_external_conductors,' # total number of external conductors'
write(file_unit,*)bundle%system_dimension,' # dimension of the matrix system characterising the MTL propagation'
write(file_unit,*)bundle%tot_n_domains,' # number of domains'
do domain=1,bundle%tot_n_domains
write(file_unit,*)'Domain number',domain
write(file_unit,*)bundle%n_conductors(domain),' number of conductors in this domain'
write(file_unit,*)'terminal_conductor_list'
write(file_unit,*)bundle%terminal_conductor_list(domain)%n_elements,' ! number of elements '
do conductor=1,bundle%terminal_conductor_list(domain)%n_elements
write(file_unit,*)bundle%terminal_conductor_list(domain)%element(conductor)
end do
matrix_dimension=bundle%L(domain)%dim
write(file_unit,*)'Per-Unit-length Inductance Matrix, [L]'
write(file_unit,*)matrix_dimension,' Dimension of [L]'
CALL dwrite_matrix(bundle%L(domain)%mat,matrix_dimension,matrix_dimension,matrix_dimension,file_unit)
write(file_unit,*)'Per-Unit-length Capacitance Matrix, [C]'
write(file_unit,*)matrix_dimension,' Dimension of [C]'
CALL dwrite_matrix(bundle%C(domain)%mat,matrix_dimension,matrix_dimension,matrix_dimension,file_unit)
write(file_unit,*)'Per-Unit-length Impedance Matrix, Z'
CALL write_Sfilter_matrix( bundle%Z(domain),file_unit )
write(file_unit,*)'Per-Unit-length Admittance Matrix, Y'
CALL write_Sfilter_matrix( bundle%Y(domain),file_unit )
end do
write(file_unit,*)'# Global current domain transformation matrix, [MI]'
matrix_dimension=bundle%global_MI%dim-1
write(file_unit,*)matrix_dimension,' Dimension of [MI]'
CALL dwrite_matrix(bundle%global_MI%mat,matrix_dimension,matrix_dimension,bundle%global_MI%dim,file_unit)
write(file_unit,*)'# Global voltage domain transformation matrix, [MV]'
matrix_dimension=bundle%global_MV%dim-1
write(file_unit,*)matrix_dimension,' Dimension of [MV]'
CALL dwrite_matrix(bundle%global_MV%mat,matrix_dimension,matrix_dimension,bundle%global_MI%dim,file_unit)
write(file_unit,*)'# Global domain based inductance matrix, [L]'
matrix_dimension=bundle%global_L%dim
write(file_unit,*)matrix_dimension,' Dimension of [L]'
CALL dwrite_matrix(bundle%global_L%mat,matrix_dimension,matrix_dimension,matrix_dimension,file_unit)
write(file_unit,*)'# Global domain based capacitance matrix, [C]'
matrix_dimension=bundle%global_C%dim
write(file_unit,*)matrix_dimension,' Dimension of [C]'
CALL dwrite_matrix(bundle%global_C%mat,matrix_dimension,matrix_dimension,matrix_dimension,file_unit)
write(file_unit,*)'Global domain based Per-Unit-length Impedance Matrix, Z'
CALL write_Sfilter_matrix( bundle%global_Z,file_unit )
write(file_unit,*)'Global domain based Per-Unit-length Admittance Matrix, Y'
CALL write_Sfilter_matrix( bundle%global_Y,file_unit )
! write the loss model information for each conductor
write(file_unit,*)' # conductor impedance models'
do conductor=1,bundle%tot_n_conductors
CALL write_conductor_impedance_model(bundle%conductor_impedance(conductor),file_unit)
end do
! write the x and y position for each conductor: this is needed for incident field excitation
write(file_unit,*)' # conductor x y positions'
do conductor=1,bundle%tot_n_conductors
write(file_unit,*)bundle%conductor_x_offset(conductor),bundle%conductor_y_offset(conductor)
end do
! numbering information
write(file_unit,'(A)')'is_shield tc_to_in_domain tc_to_out_domain tc_to_gdc tc_to_ldc tc_to_ref_tc '
do conductor=1,bundle%tot_n_conductors
write(file_unit,8000)bundle%terminal_conductor_is_shield_flag(conductor), &
bundle%terminal_conductor_to_inner_domain(conductor), &
bundle%terminal_conductor_to_outer_domain(conductor), &
bundle%terminal_conductor_to_global_domain_conductor(conductor), &
bundle%terminal_conductor_to_local_domain_conductor(conductor), &
bundle%terminal_conductor_to_reference_terminal_conductor(conductor)
8000 format(L5,5I15)
end do
! write the conductor labels
write(file_unit,*)' # Conductor labels'
do conductor=1,bundle%tot_n_conductors
write(file_unit,*)trim(bundle%conductor_label(conductor))
end do
! close .bundle file
CLOSE(unit=file_unit)
RETURN
END SUBROUTINE write_cable_bundle
! NAME
! SUBROUTINE deallocate_cable_bundle
!
! deallocate cable_bundle structure data
!
! COMMENTS
!
!
! HISTORY
! started 2/12/2015 CJS
! generalise conductor impedance for transfer impedance model (stage 5) CJS 12/5/2016
!
SUBROUTINE deallocate_cable_bundle(bundle)
USE type_specifications
USE general_module
USE cable_module
IMPLICIT NONE
! variables passed to subroutine
type(bundle_specification_type),intent(INOUT) :: bundle
! local variables
integer :: cable
integer :: domain
integer :: conductor
! START
do cable=1,bundle%n_cables
CALL deallocate_cable(bundle%cable(cable))
end do
if (allocated(bundle%cable)) DEALLOCATE( bundle%cable )
if (allocated(bundle%cable_x_offset )) DEALLOCATE( bundle%cable_x_offset )
if (allocated(bundle%cable_y_offset )) DEALLOCATE( bundle%cable_y_offset )
if (allocated(bundle%cable_angle )) DEALLOCATE( bundle%cable_angle )
if (allocated(bundle%n_conductors )) DEALLOCATE( bundle%n_conductors )
if (allocated(bundle%terminal_conductor_list )) then
do domain=1,bundle%tot_n_domains
if (allocated(bundle%terminal_conductor_list(domain)%element )) &
DEALLOCATE( bundle%terminal_conductor_list(domain)%element )
end do
DEALLOCATE( bundle%terminal_conductor_list )
end if
if (allocated(bundle%L )) then
do domain=1,bundle%tot_n_domains
if (allocated(bundle%L(domain)%mat )) DEALLOCATE( bundle%L(domain)%mat )
end do
DEALLOCATE( bundle%L )
end if
if (allocated(bundle%C )) then
do domain=1,bundle%tot_n_domains
if (allocated(bundle%C(domain)%mat )) DEALLOCATE( bundle%C(domain)%mat )
end do
DEALLOCATE( bundle%C )
end if
if (ALLOCATED(bundle%Z)) then
do domain=1,bundle%tot_n_domains
CALL deallocate_Sfilter_matrix( bundle%Z(domain) )
end do
DEALLOCATE(bundle%Z)
end if
if (ALLOCATED(bundle%Y)) then
do domain=1,bundle%tot_n_domains
CALL deallocate_Sfilter_matrix( bundle%Y(domain) )
end do
DEALLOCATE(bundle%Y)
end if
if (allocated(bundle%global_MI%mat)) DEALLOCATE( bundle%global_MI%mat )
if (allocated(bundle%global_MV%mat)) DEALLOCATE( bundle%global_MV%mat )
if (allocated(bundle%global_L%mat)) DEALLOCATE( bundle%global_L%mat )
if (allocated(bundle%global_C%mat)) DEALLOCATE( bundle%global_C%mat )
CALL deallocate_Sfilter_matrix( bundle%global_Z )
CALL deallocate_Sfilter_matrix( bundle%global_Y )
if (ALLOCATED(bundle%conductor_impedance)) then
do conductor=1,bundle%tot_n_conductors
CALL deallocate_conductor_impedance_model(bundle%conductor_impedance(conductor))
end do
DEALLOCATE(bundle%conductor_impedance)
end if
if (ALLOCATED( bundle%conductor_x_offset )) DEALLOCATE( bundle%conductor_x_offset )
if (ALLOCATED( bundle%conductor_y_offset )) DEALLOCATE( bundle%conductor_y_offset )
! Numbering information required for transfer impedance calculation
if(ALLOCATED( bundle%terminal_conductor_is_shield_flag )) DEALLOCATE( bundle%terminal_conductor_is_shield_flag )
if(ALLOCATED( bundle%terminal_conductor_to_inner_domain )) DEALLOCATE( bundle%terminal_conductor_to_inner_domain )
if(ALLOCATED( bundle%terminal_conductor_to_outer_domain )) DEALLOCATE( bundle%terminal_conductor_to_outer_domain )
if(ALLOCATED( bundle%terminal_conductor_to_global_domain_conductor )) &
DEALLOCATE( bundle%terminal_conductor_to_global_domain_conductor )
if(ALLOCATED( bundle%terminal_conductor_to_local_domain_conductor )) &
DEALLOCATE( bundle%terminal_conductor_to_local_domain_conductor )
if(ALLOCATED( bundle%terminal_conductor_to_reference_terminal_conductor )) &
DEALLOCATE( bundle%terminal_conductor_to_reference_terminal_conductor )
if (ALLOCATED( bundle%conductor_label )) then
DEALLOCATE(bundle%conductor_label)
end if
RETURN
END SUBROUTINE deallocate_cable_bundle
!
! NAME
! SUBROUTINE check_cables_wrt_ground_plane
!
! check that the cables are all on one side of the ground plane to ensure consistency
! and indicate which side the cables are on by setting bundle%ground_plane_cable_side
!
! COMMENTS
!
!
! HISTORY
! started 29/06/2016 CJS
!
SUBROUTINE check_cables_wrt_ground_plane(bundle)
USE type_specifications
USE general_module
USE cable_module
IMPLICIT NONE
! variables passed to subroutine
type(bundle_specification_type),intent(INOUT) :: bundle
! local variables
integer :: cable
real(dp) :: norm_x,norm_y
real(dp) :: p
real(dp) :: gp_offset
real(dp) :: gp_offset_cable_1
! START
! calculate the ground plane normal direction
norm_x=bundle%ground_plane_nx
norm_y=bundle%ground_plane_ny
! The equation of the ground plane is r.norm-P=0, calculate P here
P=bundle%ground_plane_x*norm_x+bundle%ground_plane_y*norm_y
! loop over all of the cables, calculating the offset from the ground plane in the normal direction
! offset=r.norm-P, and check that all offsets have the same sign
if (verbose) then
write(*,*)'Checking the orientation of cables wrt the ground plane'
write(*,*)'GP normal:',norm_x,norm_y
write(*,*)'GP offset:',bundle%ground_plane_x+bundle%ground_plane_y
write(*,*)'P=',P
end if
do cable=1,bundle%n_cables_without_ground_plane
gp_offset=bundle%cable_x_offset(cable)*norm_x+bundle%cable_y_offset(cable)*norm_y-P
if (verbose) then
write(*,*)'Cable',cable,' offset',gp_offset
end if
if (cable.EQ.1) then
gp_offset_cable_1=gp_offset
else
! check consistency of ground plane offsets
if (gp_offset_cable_1*gp_offset.LE.0d0) then
run_status='ERROR there are cables both sides of the ground plane'
CALL write_program_status()
STOP 1
end if
end if ! cable.NE.1
end do ! next cable to check
if (gp_offset_cable_1.GT.0d0) then
bundle%ground_plane_cable_side=1
else
bundle%ground_plane_cable_side=-1
end if
RETURN
END SUBROUTINE check_cables_wrt_ground_plane
!
! NAME
! SUBROUTINE check_cable_intersections
!
! check that the cables do not intersect
! the intersection test checks intersection of the outer shape of each cable
! which is usually the outer dielectric layer apart from overshields or D connectors which
! are assumed not to have dielectric
!
! COMMENTS
! Revised to work with all cable types
!
! HISTORY
! started 8/10/2016 CJS
!
!
SUBROUTINE check_cable_intersection(bundle)
USE type_specifications
USE general_module
USE cable_module
IMPLICIT NONE
! variables passed to subroutine
type(bundle_specification_type),intent(IN) :: bundle
! local variables
integer :: cable1,cable2
integer :: nec1,nec2,ec1,ec2
integer :: shape1,shape2
integer :: type1,type2
real(dp) :: ox,oy,theta
real(dp) :: r
real(dp) :: wd,hd
integer :: nc,ncrow(2)
real(dp) :: rw,p,s,o,W(2)
integer :: npts1
real(dp),allocatable :: shape1_x(:)
real(dp),allocatable :: shape1_y(:)
integer :: npts2
real(dp),allocatable :: shape2_x(:)
real(dp),allocatable :: shape2_y(:)
logical :: intersect
logical :: intersect2
logical :: nested_1_in_2
logical :: nested_2_in_1
logical :: gp_intersect
logical :: intersection_found
! START
intersection_found=.FALSE.
! loop over cable 1
do cable1=1,bundle%n_cables_without_ground_plane
ox=bundle%cable_x_offset(cable1)
oy=bundle%cable_y_offset(cable1)
theta=bundle%cable_angle(cable1)
type1=bundle%cable(cable1)%cable_type
if ((bundle%cable(cable1)%cable_type.NE.cable_geometry_type_flex_cable).AND. &
(bundle%cable(cable1)%cable_type.NE.cable_geometry_type_ML_flex_cable) ) then
nec1=bundle%cable(cable1)%n_external_conductors
else
nec1=1
end if
! loop over the external conductors of cable 1
do ec1=1,nec1
shape1=bundle%cable(cable1)%external_model(ec1)%conductor_type
! generate a list of points on cable 1 outer surface
if (shape1.EQ.rectangle) then
wd=bundle%cable(cable1)%external_model(1)%dielectric_width
hd=bundle%cable(cable1)%external_model(1)%dielectric_height
CALL generate_rectangle_points(npts1,shape1_x,shape1_y,ox,oy,theta,wd,hd)
else if (shape1.EQ.circle) then
r=bundle%cable(cable1)%external_model(ec1)%dielectric_radius
CALL generate_circle_points(npts1,shape1_x,shape1_y,ox,oy,r)
else if (shape1.EQ.Dshape) then
nc=bundle%cable(cable1)%tot_n_conductors
rw=bundle%cable(cable1)%parameters(1)
p=bundle%cable(cable1)%parameters(2)
s=bundle%cable(cable1)%parameters(3)
o=bundle%cable(cable1)%parameters(4)
ncrow(1)=nc/2
ncrow(2)=(nc-1)-ncrow(1)
w(1)=(ncrow(1)-1)*p
w(2)=(ncrow(2)-1)*p
CALL generate_Dshape_points(npts1,shape1_x,shape1_y,ox,oy,w(1)/2d0,w(2)/2d0,s/2d0,rw+o,theta)
end if
if (bundle%ground_plane_present) then
! check for intersection of the cables with the ground plane
gp_intersect=.FALSE.
if (verbose) write(*,*)'TESTING: intersection of cable',cable1,' with the ground plane'
CALL gptest(npts1,shape1_x,shape1_y,gp_intersect)
if (gp_intersect) then
if (verbose) write(*,*)'Cable ',cable1,' intersects the ground plane'
intersection_found=.TRUE.
end if
end if
! loop over cable 2
do cable2=cable1+1,bundle%n_cables_without_ground_plane
ox=bundle%cable_x_offset(cable2)
oy=bundle%cable_y_offset(cable2)
theta=bundle%cable_angle(cable2)
type2=bundle%cable(cable2)%cable_type
if ((bundle%cable(cable1)%cable_type.NE.cable_geometry_type_flex_cable).AND. &
(bundle%cable(cable1)%cable_type.NE.cable_geometry_type_ML_flex_cable) ) then
nec2=bundle%cable(cable2)%n_external_conductors
else
nec2=1
end if
! loop over the external conductors of cable 2
do ec2=1,nec2
shape2=bundle%cable(cable2)%external_model(ec2)%conductor_type
! generate a list of points on cable 1 outer surface
if (shape2.EQ.rectangle) then
wd=bundle%cable(cable2)%external_model(1)%dielectric_width
hd=bundle%cable(cable2)%external_model(1)%dielectric_height
CALL generate_rectangle_points(npts2,shape2_x,shape2_y,ox,oy,theta,wd,hd)
else if (shape2.EQ.circle) then
r=bundle%cable(cable2)%external_model(ec2)%dielectric_radius
CALL generate_circle_points(npts2,shape2_x,shape2_y,ox,oy,r)
else if (shape2.EQ.Dshape) then
nc=bundle%cable(cable2)%tot_n_conductors
rw=bundle%cable(cable2)%parameters(1)
p=bundle%cable(cable2)%parameters(2)
s=bundle%cable(cable2)%parameters(3)
o=bundle%cable(cable2)%parameters(4)
ncrow(1)=nc/2
ncrow(2)=(nc-1)-ncrow(1)
w(1)=(ncrow(1)-1)*p
w(2)=(ncrow(2)-1)*p
CALL generate_Dshape_points(npts2,shape2_x,shape2_y,ox,oy,w(1)/2d0,w(2)/2d0,s/2d0,rw+o,theta)
end if
! assume that there is no intersection and the shapes are not nested
intersect=.FALSE.
nested_1_in_2=.FALSE.
nested_2_in_1=.FALSE.
! check for the intersection of cable 1 with cable 2 and flag if cable 1 is nested within cable 2.
if (verbose) write(*,*)'TESTING: intersection of cable',cable1,' with cable ',cable2
CALL shape_test(npts1,shape1_x,shape1_y,npts2,shape2_x,shape2_y,intersect,nested_1_in_2)
if (intersect) then
if (verbose) write(*,*)'Cable ',cable1,' intersects Cable',cable2
intersection_found=.TRUE.
end if
! Cables can only be nested if the outer cable is an overshield
if ( nested_1_in_2 ) then
if (type1.NE.cable_geometry_type_overshield) then
if (verbose) write(*,*)'Cable ',cable2,' is inside Cable',cable1
intersection_found=.TRUE.
else
if (verbose) write(*,*)'Cable ',cable2,' is OK inside Overshield Cable',cable1
end if
end if
! check for the intersection of cable 1 with cable 2 and flag if cable 1 is nested within cable 2.
if (verbose) write(*,*)'TESTING: intersection of cable',cable2,' with cable ',cable1
CALL shape_test(npts2,shape2_x,shape2_y,npts1,shape1_x,shape1_y,intersect2,nested_2_in_1)
if (intersect) then
if (verbose) write(*,*)'Cable ',cable2,' intersects Cable',cable1
intersection_found=.TRUE.
end if
! Cables can only be nested if the outer cable is an overshield
if ( nested_2_in_1 ) then
if (type2.NE.cable_geometry_type_overshield) then
if (verbose) write(*,*)'Cable ',cable1,' is inside Cable',cable2
intersection_found=.TRUE.
else
if (verbose) write(*,*)'Cable ',cable1,' is OK inside Overshield Cable',cable2
end if
end if
DEALLOCATE( shape2_x )
DEALLOCATE( shape2_y )
end do ! next external conductor of cable 2
end do ! next cable 2
DEALLOCATE( shape1_x )
DEALLOCATE( shape1_y )
end do ! next external conductor of cable1
end do ! next cable1
if (intersection_found) then
run_status='ERROR there are cables in the bundle which intersect'
CALL write_program_status()
STOP 1
end if
RETURN
END SUBROUTINE check_cable_intersection
!
! NAME
! SUBROUTINE gptest
!
! check that a cable does not intersect the ground plane
!
! COMMENTS
!
!
! HISTORY
! started 20/4/2017 CJS
!
!
SUBROUTINE gptest(npts,shape_x,shape_y,gp_intersect)
USE type_specifications
IMPLICIT NONE
integer,intent(IN) :: npts
real(dp),intent(IN) :: shape_x(npts)
real(dp),intent(IN) :: shape_y(npts)
logical,intent(OUT) :: gp_intersect
! local variables
integer :: i
! START
do i=1,npts
if (shape_y(i).LE.0D0) then
gp_intersect=.TRUE.
RETURN
end if
end do ! next point
RETURN
END SUBROUTINE gptest
!
! NAME
! SUBROUTINE shape_test
!
! check that cables do not intersect, also flag whether the first cable is nested
! i.e. lies entirely within the second
!
! COMMENTS
! We assume that the shape of the cable is convex and that the shape is closed
!
! HISTORY
! started 20/4/2017 CJS
!
!
SUBROUTINE shape_test(npts1,shape1_x,shape1_y,npts2,shape2_x,shape2_y,intersect,nested)
USE type_specifications
IMPLICIT NONE
integer,intent(IN) :: npts1
real(dp),intent(IN) :: shape1_x(npts1)
real(dp),intent(IN) :: shape1_y(npts1)
integer,intent(IN) :: npts2
real(dp),intent(IN) :: shape2_x(npts2)
real(dp),intent(IN) :: shape2_y(npts2)
logical,intent(OUT) :: intersect ! flag to indicate that shapes 1 and 2 intersect
logical,intent(OUT) :: nested ! flag to indicate shape 2 lies within shape 1
! local variables
integer :: i,j
real(dp) :: vx1,vy1
real(dp) :: vx2,vy2
real(dp) :: vp
real(dp) :: dirn
logical :: inside_point_found
logical :: outside_point_found
logical :: local_outside_point_found
! START
intersect=.FALSE.
nested=.FALSE.
inside_point_found=.FALSE.
outside_point_found=.FALSE.
! loop around test points in shape 2
do j=1,npts2
local_outside_point_found=.FALSE.
dirn=0.0 ! set to zero initially to indicate that it is not known
! loop around line segments in shape 1. Note there is one less line segment than the number of points
do i=1,npts1-1
! vector along line segment in shape 1
vx1=shape1_x(i+1)-shape1_x(i)
vy1=shape1_y(i+1)-shape1_y(i)
! vector from first point in line segment to test point
vx2=shape2_x(j)-shape1_x(i)
vy2=shape2_y(j)-shape1_y(i)
! calculate the z component of the vector product of the two vectors
vp=vx1*vy2-vx2*vy1
if (dirn.EQ.0D0) then ! maybe we need to use abs(dirn).LT.small here...
! this is the first value so this sets the direction to test for all other sets of points
dirn=vp
else
! test whether the sign of the direction has changed
if (dirn*vp.LT.0.0) then
! the signs are different then this indicates a test point outside shape 1
local_outside_point_found=.TRUE.
end if ! change in sign of direction
end if ! first test line segment.
end do ! next line segment in shape 1
if (local_outside_point_found) then
outside_point_found=.TRUE.
else
inside_point_found=.TRUE.
end if
end do ! next test point in shape 2
! work out whether we have an intersection
if (inside_point_found.AND.outside_point_found) intersect=.TRUE.
! work out whether we have a nested shape
if (inside_point_found.AND.(.NOT.outside_point_found)) nested=.TRUE.
RETURN
END SUBROUTINE
END MODULE cable_bundle_module