!
! 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_LC_Laplace
!
! NAME
! SUBROUTINE PUL_LC_Laplace
!
! DESCRIPTION
! Wrapping subroutine to control the calculation of PUL_parameters by the Finite Element method
! which is implemented in Laplace.F90
!
! The finite element mesh is greated by the open source software gmsh: see http://gmsh.info/
!
! 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: Work out the numbers of conductors, points, lines, line loops and surfaces in the cross section geometry
! STAGE 3: Create the input file for the mesh generator, gmsh.
! STAGE 4: Call the mesh generator, gmsh
! STAGE 5: Call Laplace to calculate the L,C and G matrices, using the mesh that we have just generated
!
! COMMENTS
! This subroutine is never called in the case of a homogeneous frequency dependent dielectric region
! If this were to change then the logic of how Laplace is called in stage 5 will have to change to
! take this into account.
!
! HISTORY
! started 5/7/2016 CJS.
! add dielectric regions 13/7/2016 CJS.
! Consolidated ground_plane, overhsield and open_boundary subroutines into a single subroutine here 13/12/2016 CJS.
! 14/10/2017 CJS make the filter fitting minimum aorder=1, border=0 and
! ensure that border=aorder-1 to make the choice of model order more sensible
! 16/11/2017 CJS Include network synthesis process to replace s-domain transfer functions
! 4/3/2018 CJS Start to include the infinite ground plane option
!
! 16/3/2018 CJS add offset to x and y for shapes. Relates to ML_flex cables
!
SUBROUTINE PUL_LC_Laplace(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
! flags to indicate the configuration
logical :: ground_plane
logical :: finite_ground_plane
logical :: infinite_ground_plane
logical :: overshield
logical :: open_boundary
! counters and loop limits which depend on the configuration required
integer :: end_conductor_loop
integer :: end_conductor_to_dielectric_loop
integer :: tot_n_boundaries_without_dielectric
integer :: first_conductor_shape
! local variables
character(LEN=filename_length) :: geom_filename ! filename for the geometry - i.e. the input file for gmsh
character(LEN=filename_length) :: mesh_filename ! filename for the mesh - i.e. the output file from gmsh
! temposrary strings used to construct filenames
character(LEN=filename_length) :: string1,string2
character(LEN=filename_length) :: temp_string
character(LEN=filename_length) :: mesh_name
real(dp) :: xmin,xmax,ymin,ymax ! extent of bundle conductors
real(dp) :: xc,yc ! centre of bundle conductors
real(dp) :: rboundary ! free space boundary radius
! variables required for generating the gmsh geometry input file: we may be able to simplify
! things here as there is some overlap with the PUL structure
integer,allocatable :: shape_type(:) ! holds a nuber which indicates the shape type
integer,allocatable :: loop_number(:) ! holds the loop number for this shape
real(dp),allocatable :: x(:) ! x coordinate of the centre of the cable to which this shape belongs
real(dp),allocatable :: y(:) ! y coordinate of the centre of the cable to which this shape belongs
real(dp),allocatable :: r(:) ! radius of a circular shape or curve radius for a Dshape
real(dp),allocatable :: rw(:) ! width1 (x dimension) of rectangular/ Dshape
real(dp),allocatable :: rw2(:) ! width2 (x dimension) of rectangular/ Dshape
real(dp),allocatable :: rh(:) ! height (y dimension) of rectangular shape
real(dp),allocatable :: rox(:) ! offset in the x direction of this shape from the centre (x():),y(:) above)
real(dp),allocatable :: roy(:) ! offset in the y direction of this shape from the centre (x():),y(:) above)
real(dp),allocatable :: rtheta(:) ! rotation angle of conductor
real(dp),allocatable :: dl(:) ! edge length for the mesh on this shape
logical,allocatable :: conductor_has_dielectric(:) ! flag to indicate that a conductor has a dielectric around it
real(dp) :: dl_gp ! edge length at the centre of the ground plane
real(dp) :: dl_local ! edge length for the mesh at rectangle edge centre
complex(dp) :: epsr ! relative permittivity value
integer,allocatable :: shape_to_region(:,:) ! array which stores the inner and outer region numbers for each shape
! undefined regions (inside PEC or outside the outer boundary) are set to -1
logical :: first_surface_is_free_space_boundary ! flag to indicate properties on the first surface defined
! counters for conductors, points, lines etc. Should be self explanatory
integer :: n_conductors_without_ground_plane
integer :: tot_n_conductors,conductor
integer :: tot_n_boundaries
integer :: n_points,point
integer :: n_lines,line
integer :: n_line_loops,line_loop
integer :: n_surfaces,surface
integer :: n_dielectric_regions,tot_n_dielectric_regions,dielectric_region
integer :: n_shapes,shape,conductor_shape
integer :: matrix_dimension
complex(dp),allocatable :: dielectric_region_epsr(:) ! list of the relative permittivity of each region
! temporary variables used in writing the geometry input file for gmsh
real(dp) :: rl,rt,rdl,rx1,rx2,ry1,ry2
real(dp) :: xp,yp,xt,yt
integer :: p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12
real(dp) :: edge_length
character(LEN=filename_length) :: command ! string used for the command which runs gmsh from the shell
integer :: gmsh_exit_stat ! exit status from the shell command which runs gmsh.
type(matrix) :: mat1 ! temporary matrix
type(matrix) :: mat2 ! temporary matrix
! frequency dependent dielectric variables
integer :: floop ! frequency loop variable
real(dp) :: f ! frequency
type(matrix),allocatable :: C_freq(:) ! frequency dependent Capacitance matrices for each frequency of evaluation
type(matrix),allocatable :: G_freq(:) ! frequency dependent Conductance matrices for each frequency of evaluation
complex(dp),allocatable :: function_to_fit(:) ! complex function to be fitted using Sfilter_fit
logical :: dielectric_defined ! flag to indicate whether there is a dielectric
logical :: frequency_dependent_dielectric ! flag to indicate whether there is a frequency dependent dielectric
integer :: gpline1,gpline2 ! line segments for infinite ground plane
integer :: row,col ! loop variables for matrix elements
integer :: i,ii ! temporary loop variables
! START
if (verbose) write(*,*)'CALLED: PUL_LC_Laplace'
! 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)
ground_plane=.FALSE.
infinite_ground_plane=.FALSE.
finite_ground_plane=.FALSE.
overshield=.FALSE.
open_boundary=.FALSE.
if ((.NOT.PUL%overshield_present).AND.(.NOT.PUL%ground_plane_present)) then
! SOLUTION TYPE 1: NO OVERSHIELD, NO GROUND PLANE
open_boundary=.TRUE.
first_surface_is_free_space_boundary=.TRUE.
end_conductor_loop=PUL%n_conductors
end_conductor_to_dielectric_loop=PUL%n_conductors
tot_n_boundaries_without_dielectric=PUL%n_conductors+1 ! add a free space boundary
first_conductor_shape=2 ! shape 1 is for the outer boundary
else if ((.NOT.PUL%overshield_present).AND.(PUL%ground_plane_present)) then
! SOLUTION TYPE 2: NO OVERSHIELD, WITH GROUND PLANE
ground_plane=.TRUE.
if (inf_gnd) then
infinite_ground_plane=.TRUE.
else
finite_ground_plane=.TRUE.
end if
tot_n_boundaries_without_dielectric=PUL%n_conductors+1 ! add a free space boundary
first_surface_is_free_space_boundary=.TRUE.
end_conductor_loop=PUL%n_conductors-1 ! the ground plane is not included here - added separately
end_conductor_to_dielectric_loop=PUL%n_conductors
first_conductor_shape=2 ! shape 1 is for the outer boundary
else
! SOLUTION TYPE 3: THERE IS AN OVERSHIELD WHICH FORMS THE OUTER BOUNDARY OF THE DOMAIN
overshield=.TRUE.
first_surface_is_free_space_boundary=.FALSE.
end_conductor_loop=PUL%n_conductors-1
end_conductor_to_dielectric_loop=PUL%n_conductors-1 ! overshield has no dielectric in the internal domain
tot_n_boundaries_without_dielectric=PUL%n_conductors
first_conductor_shape=1 ! shape 1 is a conductor
end if ! no overshield
if (verbose) then
write(*,*)'ground_plane:',ground_plane
write(*,*)'infinite_ground_plane:',infinite_ground_plane
write(*,*)'finite_ground_plane:',finite_ground_plane
write(*,*)'open_boundary:',open_boundary
write(*,*)'overshield:',overshield
write(*,*)'first_surface_is_free_space_boundary:',first_surface_is_free_space_boundary
write(*,*)'N_conductors:',PUL%n_conductors
write(*,*)'end_conductor_loop:',end_conductor_loop
write(*,*)'end_conductor_to_dielectric_loop:',end_conductor_to_dielectric_loop
write(*,*)'tot_n_boundaries_without_dielectric:',tot_n_boundaries_without_dielectric
write(*,*)'first_conductor_shape:',first_conductor_shape
end if ! verbose
! STAGE 2: Work out the numbers of conductors, points, lines, line loops and surfaces in the cross section geometry
tot_n_conductors = PUL%n_conductors ! the conductor count includes the ground plane if present
! Count the dielectric regions:
! Each dielectric region adds a boundary and another surface to mesh so work out whether a dielectric
! is present for each conductor.
! A dielectric region is assumed to exist if the dielectric radius is larger than the conductor radius and
! the relative permittivity at d.c. is not equal to 1 or the permittivity model order is not equal to zero
ALLOCATE ( conductor_has_dielectric(1:PUL%n_conductors) )
n_dielectric_regions=0
do conductor=1,end_conductor_loop
! assume no dielectric initially
conductor_has_dielectric(conductor)=.FALSE.
! If the low freq permittivity is not 1.0 or the order of the filter model is not 0 then assume we have a dielectric
dielectric_defined=.FALSE.
! check the low frequency value of epsr
epsr=PUL%epsr(conductor)%a%coeff(0)/PUL%epsr(conductor)%b%coeff(0)
if ( abs(epsr-1d0).GT.small ) dielectric_defined=.TRUE.
! check the model order
if ( (PUL%epsr(conductor)%a%order.GT.0).OR. (PUL%epsr(conductor)%b%order.GT.0) ) dielectric_defined=.TRUE.
! for now use the high frequency limit of epsr
epsr=evaluate_Sfilter_high_frequency_limit(PUL%epsr(conductor))
if (PUL%shape(conductor).EQ.circle) THEN
if ( ( (PUL%rd(conductor)-PUL%r(conductor)).GT.small ).AND.( dielectric_defined ) ) then
n_dielectric_regions= n_dielectric_regions+1
conductor_has_dielectric(conductor)=.TRUE.
end if
else if (PUL%shape(conductor).EQ.rectangle) THEN
! include dielectric region if the dielectric is not air and there is a dielectric region surrounding the conductor
if ( ( dielectric_defined ).AND. ( ( PUL%rdw(conductor).GT.small ).AND.( PUL%rdw(conductor).GT.small ) ) ) then
n_dielectric_regions= n_dielectric_regions+1
conductor_has_dielectric(conductor)=.TRUE.
end if
end if
end do ! next conductor
tot_n_boundaries = tot_n_boundaries_without_dielectric+n_dielectric_regions
tot_n_dielectric_regions = n_dielectric_regions+1 ! add a dielectric region for the 'free space' region
n_points = 5*tot_n_boundaries ! Each region is enclosed by a shape, defined by 5 points
n_lines = 4*tot_n_boundaries ! Each region is enclosed by a shape, defined by 4 lines
n_line_loops= tot_n_boundaries ! Each region is enclosed by a shape, defined by a single line loop
n_surfaces = n_dielectric_regions+1 ! Each dielectric region plus the free space region are a separate surface
if (verbose) then
write(*,*)'Number of dielectric regions =',n_dielectric_regions
write(*,*)'Total number of regions =',tot_n_dielectric_regions
write(*,*)'Total number of conductors =',tot_n_conductors
write(*,*)'Total number of boundaries =',tot_n_boundaries
write(*,*)'first_surface_is_free_space_boundary? ',first_surface_is_free_space_boundary
end if
! STAGE 3: Create the input file for the mesh generator, gmsh.
! Allocate the memory required for all the boundaries in the cross section geometry
n_shapes=tot_n_boundaries
ALLOCATE( shape_type(1:n_shapes) )
shape_type(1:n_shapes)=circle ! set all boudaries to type circle initially
ALLOCATE( loop_number(1:n_shapes) )
loop_number(1:n_shapes)=-1 ! set all boudaries to -1 i.e. not set
ALLOCATE( x(1:n_shapes) )
ALLOCATE( y(1:n_shapes) )
ALLOCATE( r(1:n_shapes) )
ALLOCATE( rw(1:n_shapes) )
ALLOCATE( rw2(1:n_shapes) )
ALLOCATE( rh(1:n_shapes) )
ALLOCATE( rox(1:n_shapes) )
ALLOCATE( roy(1:n_shapes) )
ALLOCATE( rtheta(1:n_shapes) )
ALLOCATE( dl(1:n_shapes) )
ALLOCATE( shape_to_region(1:n_shapes,1:2) )
shape_to_region(1:n_shapes,1:2)=-1 ! set all regions to undefined initially
ALLOCATE( dielectric_region_epsr(0:n_dielectric_regions) )
if (.NOT.overshield) then
! Set the position of the outer free space boundary. This is related to the size of the cable bundle enclosed
! by the Laplace_boundary_constant
! First get the extent of the bundle
xmin=large
xmax=-large
ymin=large
ymax=-large
do i=1,end_conductor_loop ! if the last conductor is the ground plane it is included in x/y,max/min calculation later
if (PUL%shape(i).EQ.circle) then
xmin=min(xmin,PUL%x(i)-PUL%r(i)+PUL%ox(i))
xmax=max(xmax,PUL%x(i)+PUL%r(i)+PUL%ox(i))
ymin=min(ymin,PUL%y(i)-PUL%r(i)+PUL%oy(i))
ymax=max(ymax,PUL%y(i)+PUL%r(i)+PUL%oy(i))
else if (PUL%shape(i).EQ.rectangle) then
xmin=min(xmin,PUL%x(i)-PUL%rw(i)+PUL%ox(i))
xmax=max(xmax,PUL%x(i)+PUL%rw(i)+PUL%ox(i))
ymin=min(ymin,PUL%y(i)-PUL%rh(i)+PUL%oy(i))
ymax=max(ymax,PUL%y(i)+PUL%rh(i)+PUL%oy(i))
else if (PUL%shape(i).EQ.Dshape) then
!*** Does this have to take rotation into account ? *****
xmin=min(xmin,PUL%x(i)-PUL%rw(i))
xmax=max(xmax,PUL%x(i)+PUL%rw(i))
ymin=min(ymin,PUL%y(i)-PUL%rh(i))
ymax=max(ymax,PUL%y(i)+PUL%rh(i))
! note don't check semicircle as this is only for infinite ground only
else if (PUL%shape(i).NE.semicircle) then
write(*,*)'Unknown shape type',PUL%shape(i)
end if
end do ! next conductor
if (ground_plane) then
! include the ground plane in the bundle sizing by adding the origin point
! here we assume that the ground plane is along the x axis.
xmin=min(xmin,0d0)
xmax=max(xmax,0d0)
ymin=min(ymin,0d0)
ymax=max(ymax,0d0)
end if ! ground_plane
! Centre of conductor bundle
xc=(xmax+xmin)/2d0
yc=(ymax+ymin)/2d0
rboundary=max((xmax-xmin),(ymax-ymin))*Laplace_boundary_constant
if (verbose) then
write(*,*)'bundle xmin=',xmin,' bundle xmax=',xmax
write(*,*)'bundle ymin=',ymin,' bundle ymax=',ymax
write(*,*)'bundle centre=',xc,yc
write(*,*)'boundary radius=',rboundary
end if
! set the first shape information to relate to the outer boundary
!
! If we have a finite gound plane of free space boundary then this is a circle, centred on the bundle centre
! note that for the Laplace calculation the geometry will be shifted so that the origin is at xc,yc so the
! free space outer boundary is centred on the origin. This is required for the free space boundary condition in Laplace to work correctly.
if (.NOT.infinite_ground_plane) then
x(1)=xc
y(1)=yc
r(1)=rboundary
rw(1)=0d0
rw2(1)=0d0
rh(1)=0d0
rox(1)=0d0
roy(1)=0d0
rtheta(1)=0d0
dl(1)=r(1)/(2*Laplace_surface_mesh_constant)
shape_to_region(1,inside) =0 ! inside is the dielectric region
shape_to_region(1,outside)=-1 ! no outside region
else
!
! If we have an infinite ground plane then we have a semicircle centred on the origin
shape_type(1)=semicircle
xc=0d0
yc=0d0
x(1)=xc
y(1)=yc
r(1)=rboundary
rw(1)=0d0
rw2(1)=0d0
rh(1)=0d0
rox(1)=0d0
roy(1)=0d0
rtheta(1)=0d0
dl(1)=r(1)/(2*Laplace_surface_mesh_constant)
shape_to_region(1,inside) =0 ! inside is the dielectric region
shape_to_region(1,outside)=-1 ! no outside region
dl_gp=min(gp_edge_length,dl(1))
end if ! infinite ground plane
else ! there is an overshield specified so we do not need to offset the bundle to define a free space outer boundary
xc=0d0
yc=0d0
end if ! overshield
dielectric_region_epsr(0)= evaluate_Sfilter_high_frequency_limit(PUL%epsr_background) ! set material properties in region 0, the background region
! copy the conductor information from the PUL structure to the shape based structure
do i=1,end_conductor_loop ! all conductors except the ground plane if it exists
shape=first_conductor_shape-1+i
shape_type(shape)=PUL%shape(i)
x(shape)=PUL%x(i)
y(shape)=PUL%y(i)
r(shape)=PUL%r(i)
rw(shape)=PUL%rw(i)
rw2(shape)=PUL%rw2(i)
rh(shape)=PUL%rh(i)
rox(shape)=PUL%ox(i)
roy(shape)=PUL%oy(i)
rtheta(shape)=PUL%rtheta(i)
if (shape_type(shape).EQ.circle) then
dl(shape)=r(shape)/Laplace_surface_mesh_constant
else if (shape_type(shape).EQ.rectangle) then
dl(shape)=min(rw(shape),rh(shape))/Laplace_surface_mesh_constant
else if (shape_type(shape).EQ.Dshape) then
dl(shape)=r(shape)/(2D0*Laplace_surface_mesh_constant)
else if (shape_type(shape).EQ.semicircle) then
dl(shape)=r(shape)/Laplace_surface_mesh_constant
else
write(*,*)'Unknown shape type',shape_type(shape),' i=',i,' shape=',shape
end if
end do
if (ground_plane) then ! Add the ground plane information
shape=PUL%n_conductors+1 ! ground plane
if (finite_ground_plane) then
shape_type(shape)=rectangle
rl=rboundary*Laplace_ground_plane_ratio ! width of the ground plane (x dimension) see constants.F90 for parameter
rt=rl*Laplace_ground_plane_thickness_ratio ! height of the ground plane (y dimension) see constants.F90 for parameter
rdl=rt/2d0
x(shape)=xc ! x centre of ground plane rectangle
y(shape)=-rt ! y centre of ground plane rectangle
r(shape)=0d0 ! not used
rw(shape)=rl*2d0 ! x extent of the ground plane
rw2(shape)=rw(shape) ! x extent of the ground plane
rh(shape)=rt*2d0 ! y extent of the ground plane
rox(shape)=0d0 ! no offset from centre
roy(shape)=0d0 ! no offset from centre
rtheta(shape)=0d0 ! always zero for ground plane
dl(shape)=rdl ! this is the thickness of the ground plane
shape_to_region(shape,inside) =-1 ! no internal region
shape_to_region(shape,outside) = 0 ! no dielectric so a free space region
else
shape_type(shape)=semicircle_diameter
rl=0d0 ! not used
rt=0d0 ! not used
rdl=0d0 ! not used
x(shape)=0d0 ! x centre of ground plane
y(shape)=0d0 ! y centre of ground plane
r(shape)=rboundary ! semicircle radius
rw(shape)=0d0 ! not used
rw2(shape)=0d0 ! not used
rh(shape)=0d0 ! not used
rox(shape)=0d0 ! no offset from centre
roy(shape)=0d0 ! no offset from centre
rtheta(shape)=0d0 ! always zero for infinite ground plane
dl(shape)=r(1)/(2*Laplace_surface_mesh_constant)
shape_to_region(shape,inside) =-1 ! no internal region
shape_to_region(shape,outside) = 0 ! no dielectric so a free space region
end if ! infinite ground plane
else if (overshield) then ! Add the overshield conductor information
shape=PUL%n_conductors
shape_type(shape)=PUL%overshield_shape
x(shape)=PUL%overshield_x
y(shape)=PUL%overshield_y
r(shape)=PUL%overshield_r
rw(shape)=PUL%overshield_w
rw2(shape)=PUL%overshield_w2
rh(shape)=PUL%overshield_h
rox(shape)=0d0
roy(shape)=0d0
rtheta(shape)=0d0
dl(shape)=r(shape)/(2d0*Laplace_surface_mesh_constant)
shape_to_region(shape,inside) = 0 ! inside the overshield is the free space region
end if ! overshield
! loop over conductors copying the associated dielectric information from the PUL structure,
! also associate shapes with regions
dielectric_region=0
do i=1,end_conductor_loop ! exclude the ground plane if it exists as it has no dielectric region
conductor_shape=first_conductor_shape-1+i ! shape number for the conductor
shape_to_region(conductor_shape,inside) =-1 ! conductors have no inner surface
if (conductor_has_dielectric(i)) then
dielectric_region=dielectric_region+1
shape=shape+1 ! create a new shape number for the dielectric
shape_type(shape)=PUL%shape(i) ! the dielectric is the same shape as the conductor
x(shape)=PUL%x(i)
y(shape)=PUL%y(i)
r(shape)=PUL%rd(i)
rw(shape)=PUL%rdw(i)
rw2(shape)=PUL%rw2(i)
rh(shape)=PUL%rdh(i)
rox(shape)=PUL%rdox(i)
roy(shape)=PUL%rdoy(i)
rtheta(shape)=PUL%rtheta(i)
if (shape_type(shape).EQ.circle) then
dl(shape)=r(shape)/Laplace_surface_mesh_constant
else if (shape_type(shape).EQ.rectangle) then
dl(shape)=min(rw(shape),rh(shape))/Laplace_surface_mesh_constant
else if (shape_type(shape).EQ.Dshape) then
dl(shape)=r(shape)/(2D0*Laplace_surface_mesh_constant)
else
write(*,*)'Unknown shape type',shape_type(shape)
end if
shape_to_region(shape,inside) =dielectric_region
shape_to_region(shape,outside) =0 ! free space region
dielectric_region_epsr(dielectric_region)=evaluate_Sfilter_high_frequency_limit(PUL%epsr(i))
else
! no dielectric
shape_to_region(conductor_shape,outside) =0 ! no dielectric so a free space region
end if ! this conductor has a dielectric or not
end do ! next conductor
! Work out which conductors are associated with which dielectrics
do i=1,end_conductor_loop
conductor_shape=first_conductor_shape-1+i ! shape number for the conductor
dielectric_region=0
! calculate the centre of the conductor (the test point)
! xp,yp is the coordinate of the centre of the conductor relative to the centre of the cable
xp=rox(conductor_shape)
yp=roy(conductor_shape)
! rotate the cable about it's origin then translate to the cable position to give the final conductor position
xt=x(conductor_shape)+xp*cos(rtheta(conductor_shape))-yp*sin(rtheta(conductor_shape))
yt=y(conductor_shape)+xp*sin(rtheta(conductor_shape))+yp*cos(rtheta(conductor_shape))
! loop over dielectrc shapes
do shape=tot_n_boundaries_without_dielectric+1,tot_n_boundaries
dielectric_region=dielectric_region+1
! check whether the conductor is inside the dielectric
if (point_is_inside(xt,yt,shape_type(shape),x(shape),y(shape),r(shape),rh(shape),rw(shape), &
rox(shape),roy(shape),rtheta(shape))) then
shape_to_region(conductor_shape,outside) =dielectric_region
write(*,*)'Conductor ',i,' is inside dielectric',dielectric_region,'(shape',shape,'):'
endif
end do ! next dielectric shape
end do ! next conductor
if (verbose) then
write(*,*)' shape shape_type x(shape) y(shape) r(shape) rh(shape) rw(shape) rtheta(shape) dl(shape)'
do shape=1,n_shapes
write(*,'(2I10,7F12.6)')shape,shape_type(shape),x(shape),y(shape),r(shape),rw(shape),rh(shape),rtheta(shape),dl(shape)
end do
end if
! Open a file for the gmsh geometry input data
temp_string=trim(name)//'_geometry_domain_'
CALL add_integer_to_string(temp_string,domain,geom_filename)
geom_filename=trim(geom_filename)//gmsh_geometry_file_extn
open (unit=gmsh_geometry_file_unit,file=trim(geom_filename))
if (verbose) write(*,*)'Opened file:',trim(geom_filename)
! open a file for the boundary information
open(unit=boundary_file_unit,file=boundary_file_name)
if (verbose) write(*,*)'Opened file:',boundary_file_name
! write the points required to define the each shape
! Translate all coordinates by (-xc,-yc) such that a free space outer boundary is centred at the origin
point=0
do i=1,n_shapes
if (shape_type(i).EQ.circle) then
write(gmsh_geometry_file_unit,*)' // circle ',i
point=point+1
write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',x(i)-xc+rox(i),',' &
,y(i)-yc+roy(i),',',0.0,',' ,dl(i),' };'
point=point+1
write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',x(i)-xc+rox(i),',' &
,y(i)-yc-r(i)+roy(i),',',0.0,',',dl(i),' };'
point=point+1
write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',x(i)-xc+rox(i)+r(i),',' &
,y(i)-yc+roy(i),',',0.0,',' ,dl(i),' };'
point=point+1
write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',x(i)-xc+rox(i),',' &
,y(i)-yc+r(i)+roy(i),',',0.0,',',dl(i),' };'
point=point+1
write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',x(i)-xc+rox(i)-r(i),',' &
,y(i)-yc+roy(i),',',0.0,',' ,dl(i),' };'
else if (shape_type(i).EQ.rectangle) then
write(gmsh_geometry_file_unit,*)' // rectangle ',i
! corner point
point=point+1
xt=rw(i)/2d0+rox(i)
yt=rh(i)/2d0+roy(i)
xp=x(i)-xc+xt*cos(rtheta(i))-yt*sin(rtheta(i))
yp=y(i)-yc+xt*sin(rtheta(i))+yt*cos(rtheta(i))
edge_length=min(max_mesh_edge_length,dl(i))
write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',edge_length,' };'
! edge centre point
point=point+1
xt=rox(i)
yt=rh(i)/2d0+roy(i)
xp=x(i)-xc+xt*cos(rtheta(i))-yt*sin(rtheta(i))
yp=y(i)-yc+xt*sin(rtheta(i))+yt*cos(rtheta(i))
dl_local=min(max_mesh_edge_length,rw(i)/Laplace_surface_mesh_constant)
write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl_local,' };'
! corner point
point=point+1
xt=-rw(i)/2d0+rox(i)
yt=rh(i)/2d0+roy(i)
xp=x(i)-xc+xt*cos(rtheta(i))-yt*sin(rtheta(i))
yp=y(i)-yc+xt*sin(rtheta(i))+yt*cos(rtheta(i))
edge_length=min(max_mesh_edge_length,dl(i))
write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',edge_length,' };'
! edge centre point
point=point+1
xt=-rw(i)/2d0+rox(i)
yt=roy(i)
xp=x(i)-xc+xt*cos(rtheta(i))-yt*sin(rtheta(i))
yp=y(i)-yc+xt*sin(rtheta(i))+yt*cos(rtheta(i))
dl_local=min(max_mesh_edge_length,rh(i)/Laplace_surface_mesh_constant)
write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl_local,' };'
! corner point
point=point+1
xt=-rw(i)/2d0+rox(i)
yt=-rh(i)/2d0+roy(i)
xp=x(i)-xc+xt*cos(rtheta(i))-yt*sin(rtheta(i))
yp=y(i)-yc+xt*sin(rtheta(i))+yt*cos(rtheta(i))
edge_length=min(max_mesh_edge_length,dl(i))
write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',edge_length,' };'
! edge centre point
point=point+1
xt=rox(i)
yt=-rh(i)/2d0+roy(i)
xp=x(i)-xc+xt*cos(rtheta(i))-yt*sin(rtheta(i))
yp=y(i)-yc+xt*sin(rtheta(i))+yt*cos(rtheta(i))
dl_local=min(max_mesh_edge_length,rw(i)/Laplace_surface_mesh_constant)
write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl_local,' };'
! corner point
point=point+1
xt=rw(i)/2d0+rox(i)
yt=-rh(i)/2d0+roy(i)
xp=x(i)-xc+xt*cos(rtheta(i))-yt*sin(rtheta(i))
yp=y(i)-yc+xt*sin(rtheta(i))+yt*cos(rtheta(i))
edge_length=min(max_mesh_edge_length,dl(i))
write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',edge_length,' };'
! edge centre point
point=point+1
xt=rw(i)/2d0+rox(i)
yt=roy(i)
xp=x(i)-xc+xt*cos(rtheta(i))-yt*sin(rtheta(i))
yp=y(i)-yc+xt*sin(rtheta(i))+yt*cos(rtheta(i))
dl_local=min(max_mesh_edge_length,rh(i)/Laplace_surface_mesh_constant)
write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl_local,' };'
else if (shape_type(i).EQ.Dshape) then
CALL write_Dshape_gmsh(x(i),y(i),rw(i),rw2(i),rh(i),r(i),rox(i),roy(i),rtheta(i),dl(i),point,i,gmsh_geometry_file_unit)
else if (shape_type(i).EQ.semicircle) then
write(gmsh_geometry_file_unit,*)' // semicircle ',i
point=point+1
write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',x(i)-xc+rox(i),',' &
,y(i)-yc+roy(i),',',0.0,',' ,dl_gp,' };' ! middle of ground plane point
point=point+1
write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',x(i)-xc+rox(i)+r(i),',' &
,y(i)-yc+roy(i),',',0.0,',',dl(i),' };'
point=point+1
write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',x(i)-xc+rox(i) ,',' &
,y(i)-yc+r(i)+roy(i),',',0.0,',' ,dl(i),' };'
point=point+1
write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',x(i)-xc+rox(i)-r(i),',' &
,y(i)-yc+roy(i),',',0.0,',',dl(i),' };'
end if
end do ! next shape
! write the boudary line segment definitions for each shape, these are constructed from the previously defined points
write(gmsh_geometry_file_unit,*)' '
point=0 ! initialise point counter
line=0 ! initialise line counter
do i=1,n_shapes
if (shape_type(i).EQ.circle) then
p1=point+1
p2=point+2
p3=point+3
p4=point+4
p5=point+5
write(gmsh_geometry_file_unit,*)' // circle ',i
line=line+1
write(gmsh_geometry_file_unit,*)'Circle(',line,' ) = {',p2,',',p1,',',p3,' };'
line=line+1
write(gmsh_geometry_file_unit,*)'Circle(',line,' ) = {',p3,',',p1,',',p4,' };'
line=line+1
write(gmsh_geometry_file_unit,*)'Circle(',line,' ) = {',p4,',',p1,',',p5,' };'
line=line+1
write(gmsh_geometry_file_unit,*)'Circle(',line,' ) = {',p5,',',p1,',',p2,' };'
point=point+5
else if (shape_type(i).EQ.rectangle) then
p1=point+1
p2=point+2
p3=point+3
p4=point+4
p5=point+5
p6=point+6
p7=point+7
p8=point+8
write(gmsh_geometry_file_unit,*)' // rectangle line segments ',i
line=line+1
write(gmsh_geometry_file_unit,*)'Line(',line,' ) = {',p1,',',p2,' };'
line=line+1
write(gmsh_geometry_file_unit,*)'Line(',line,' ) = {',p2,',',p3,' };'
line=line+1
write(gmsh_geometry_file_unit,*)'Line(',line,' ) = {',p3,',',p4,' };'
line=line+1
write(gmsh_geometry_file_unit,*)'Line(',line,' ) = {',p4,',',p5,' };'
line=line+1
write(gmsh_geometry_file_unit,*)'Line(',line,' ) = {',p5,',',p6,' };'
line=line+1
write(gmsh_geometry_file_unit,*)'Line(',line,' ) = {',p6,',',p7,' };'
line=line+1
write(gmsh_geometry_file_unit,*)'Line(',line,' ) = {',p7,',',p8,' };'
line=line+1
write(gmsh_geometry_file_unit,*)'Line(',line,' ) = {',p8,',',p1,' };'
point=point+8
else if (shape_type(i).EQ.Dshape) then
p1=point+1
p2=point+2
p3=point+3
p4=point+4
p5=point+5
p6=point+6
p7=point+7
p8=point+8
p9=point+9
p10=point+10
p11=point+11
p12=point+12
write(gmsh_geometry_file_unit,*)' // Dshape line segments',i
line=line+1
write(gmsh_geometry_file_unit,*)'Line(',line,' ) = {',p1,',',p2,' };'
line=line+1
write(gmsh_geometry_file_unit,*)'Circle(',line,' ) = {',p2,',',p3,',',p4,' };'
line=line+1
write(gmsh_geometry_file_unit,*)'Line(',line,' ) = {',p4,',',p5,' };'
line=line+1
write(gmsh_geometry_file_unit,*)'Circle(',line,' ) = {',p5,',',p6,',',p7,' };'
line=line+1
write(gmsh_geometry_file_unit,*)'Line(',line,' ) = {',p7,',',p8,' };'
line=line+1
write(gmsh_geometry_file_unit,*)'Circle(',line,' ) = {',p8,',',p9,',',p10,' };'
line=line+1
write(gmsh_geometry_file_unit,*)'Line(',line,' ) = {',p10,',',p11,' };'
line=line+1
write(gmsh_geometry_file_unit,*)'Circle(',line,' ) = {',p11,',',p12,',',p1,' };'
point=point+12
else if (shape_type(i).EQ.semicircle) then
p1=point+1
p2=point+2
p3=point+3
p4=point+4
write(gmsh_geometry_file_unit,*)' // semicircle ',i
line=line+1
write(gmsh_geometry_file_unit,*)'Circle(',line,' ) = {',p2,',',p1,',',p3,' };'
line=line+1
write(gmsh_geometry_file_unit,*)'Circle(',line,' ) = {',p3,',',p1,',',p4,' };'
point=point+4
else if (shape_type(i).EQ.semicircle_diameter) then
! this is the infinite ground plane line so pick points from the outer boundary circle
p1=4 ! semicircle end point 2
p2=1 ! centre
p3=2 ! semicircle end point 1
write(gmsh_geometry_file_unit,*)' // semicircle_diameter ',i
line=line+1
write(gmsh_geometry_file_unit,*)'Line(',line,' ) = {',p1,',',p2,' };'
gpline1=line
line=line+1
write(gmsh_geometry_file_unit,*)'Line(',line,' ) = {',p2,',',p3,' };'
gpline2=line
end if
end do ! next shape
! write the number of boundary segments to the boundary file
write(boundary_file_unit,*)line,' # number of boundary segments'
! Create the closed boundary line loops for each shape, these are constructed from the preciously defined line segments.
write(gmsh_geometry_file_unit,*)' '
line_loop=0 ! initialise line_loop counter
line=0 ! reset the line counter to the first line
do i=1,n_shapes
if ( (shape_type(i).EQ.circle) ) then
line_loop=line_loop+1
loop_number(i)=line_loop
write(gmsh_geometry_file_unit,*)'Line Loop(',line_loop,' ) = {',line+1,',',line+2,',',line+3,',',line+4,' };'
! write the boundary segment and the boundary number to the boundary file
do ii=1,4
write(boundary_file_unit,*)line+ii,i,' # boundary segment number and boundary number'
end do
line=line+4
else if ( shape_type(i).EQ.rectangle ) then
line_loop=line_loop+1
loop_number(i)=line_loop
write(gmsh_geometry_file_unit,*)'Line Loop(',line_loop,' ) = {',line+1,',',line+2,',',line+3,',',line+4,',', &
line+5,',',line+6,',',line+7,',',line+8,' };'
! write the boundary segment and the boundary number to the boundary file
do ii=1,8
write(boundary_file_unit,*)line+ii,i,' # boundary segment number and boundary number'
end do
line=line+8
else if (shape_type(i).EQ.semicircle) then
! this is the outer boundary with infinite ground plane.
! the line loop consists of 4 lines, 2 for the semicircle and 2 for the ground plane
line_loop=line_loop+1
loop_number(i)=line_loop
write(gmsh_geometry_file_unit,*)'Line Loop(',line_loop,' ) = {',line+1,',',line+2,',',gpline1,',',gpline2,' };'
! write the boundary segment and the boundary number to the boundary file
do ii=1,2
write(boundary_file_unit,*)line+ii,i,' # boundary segment number and boundary number'
end do
line=line+2
else if (shape_type(i).EQ.Dshape) then
line_loop=line_loop+1
loop_number(i)=line_loop
write(gmsh_geometry_file_unit,*)'Line Loop(',line_loop,' ) = {',line+1,',',line+2,',',line+3,',',line+4,',', &
line+5,',',line+6,',',line+7,',',line+8,' };'
! write the boundary segment and the boundary number to the boundary file
do ii=1,8
write(boundary_file_unit,*)line+ii,i,' # boundary segment number and boundary number'
end do
line=line+8
else if (shape_type(i).EQ.semicircle_diameter) then
! note there is no line loop for the semicircle_diameter shape as this is dealt with along with the semicircle
! just skip these lines
! write the boundary segment and the boundary number to the boundary file
do ii=1,2
write(boundary_file_unit,*)line+ii,i,' # boundary segment number and boundary number'
end do
line=line+2
end if
end do ! next shape
! Create the plane Surfaces, one for each dielectric region. Surfaces are defined by their boundling line loops
! Note that the order of the line loop specification is important for the orientation of the elements in the mesh
write(gmsh_geometry_file_unit,*)' '
surface=0 ! initialise surface counter
! loop over the dielectric regions (including free space) as each one contributes a surface to be meshed.
do dielectric_region=0,n_dielectric_regions ! 0 is the free space region which always exists
surface=surface+1
string1=' Plane Surface('
CALL add_integer_to_string(string1,surface,string2)
string1=trim(string2)//') = {'
! for this surface, write a list of all the line loops bounding the surface
! in order to construct this list we loop over all shapes (line loops)
! and check whether the inside or outside region of the line loop is associated with the current region
if (surface.eq.1) then
! order the line loops for the free space region outer boundary first
if (overshield) then
! The overshield reference conductor should be the last conductor
i=PUL%n_conductors
line_loop=loop_number(i)
if( (shape_to_region(i,inside).EQ.dielectric_region).OR. &
(shape_to_region(i,outside).EQ.dielectric_region) ) then
CALL add_integer_to_string(string1,line_loop,string2)
string1=trim(string2)//', '
end if
end if ! overshield
do i=1,n_shapes ! don't include ground plane in check
line_loop=loop_number(i)
if ((.NOT.overshield).OR.(i.NE.PUL%n_conductors)) then
if( (shape_to_region(i,inside).EQ.dielectric_region).OR. &
(shape_to_region(i,outside).EQ.dielectric_region) ) then
if (line_loop.GT.0) then
CALL add_integer_to_string(string1,line_loop,string2)
string1=trim(string2)//', '
end if
end if
end if
end do ! next shape
else ! not the fisrt surface
! dielectric boundaries are defined second so reverse the order of the loop for dielectric regions
do i=n_shapes,1,-1
line_loop=loop_number(i)
if( (shape_to_region(i,inside).EQ.dielectric_region).OR. &
(shape_to_region(i,outside).EQ.dielectric_region) ) then
if (line_loop.GT.0) then
CALL add_integer_to_string(string1,line_loop,string2)
string1=trim(string2)//', '
end if
end if
end do
end if
string1=trim(string2)//'} '
write(gmsh_geometry_file_unit,'(A,A)')trim(string1),' ;'
end do ! next region (surface) to mesh
write(gmsh_geometry_file_unit,*)' '
! Close the file for the gmsh input data
close(gmsh_geometry_file_unit)
if (verbose) write(*,*)'Closed file:',trim(geom_filename)
! close the boundary information file
close(unit=boundary_file_unit)
if (verbose) write(*,*)'Closed file:',boundary_file_name
! Dealllocate the local shape geometry data.
if (allocated( shape_type )) DEALLOCATE( shape_type )
if (allocated( loop_number )) DEALLOCATE( loop_number )
if (allocated( x )) DEALLOCATE( x )
if (allocated( y )) DEALLOCATE( y )
if (allocated( r )) DEALLOCATE( r )
if (allocated( rw )) DEALLOCATE( rw )
if (allocated( rw2 )) DEALLOCATE( rw2 )
if (allocated( rh )) DEALLOCATE( rh )
if (allocated( rox )) DEALLOCATE( rox )
if (allocated( roy )) DEALLOCATE( roy )
if (allocated( rtheta )) DEALLOCATE( rtheta )
if (allocated( dl )) DEALLOCATE( dl )
if (allocated( shape_to_region )) DEALLOCATE( shape_to_region )
! STAGE 4: Call the mesh generator, gmsh
temp_string=trim(name)//'_mesh_domain_'
CALL add_integer_to_string(temp_string,domain,mesh_filename)
mesh_filename=trim(mesh_filename)//mesh_file_extn
write(*,*)'CALLING gmsh for domain',domain
command='gmsh -2 -o '//trim(mesh_filename)//' '//trim(geom_filename)
CALL execute_command_line(command,EXITSTAT=gmsh_exit_stat)
! Check that the mesh generator finished correctly, if not exit and write the gmsh error code to the screen
if (gmsh_exit_stat.NE.0) then
! gmsh returned with a non zdero exit code indicating an error
write(run_status,*)'ERROR in PUL_LC_Laplace. gmsh exit status=',gmsh_exit_stat
CALL write_program_status()
STOP 1
end if
! STAGE 5: Call Laplace to calculate the L,C and G matrices, using the mesh that we have just generated
! Work out the dimension of the per-unit-length matrices
matrix_dimension=PUL%n_conductors-1
if (verbose) then
write(*,*)'Number of conductors=',PUL%n_conductors
write(*,*)'matrix dimension=',matrix_dimension
end if
! check we have a valid system
if (matrix_dimension.LT.1) then
run_status='ERROR in PUL_LC_Laplace, Matrix dimension is less than 1 '
CALL write_program_status()
STOP 1
end if
! work out if we have any frequency dependent dielectric
frequency_dependent_dielectric=.FALSE.
dielectric_region=0
if ( (PUL%epsr_background%a%order.GT.0).OR. (PUL%epsr_background%b%order.GT.0) ) then
frequency_dependent_dielectric=.TRUE.
end if
! check the external dielectric first
do i=1,end_conductor_loop ! don't include the ground plane in the conductor loop for the dielectric region count
write(*,*)'Check dielectric loop',i,' conductor has dielectric:',conductor_has_dielectric(i)
if (conductor_has_dielectric(i)) then
write(*,*)'aorder,border',PUL%epsr(i)%a%order,PUL%epsr(i)%b%order
if ( (PUL%epsr(i)%a%order.GT.0).OR. (PUL%epsr(i)%b%order.GT.0) ) then
frequency_dependent_dielectric=.TRUE.
end if
end if
end do ! next conductor
if (verbose) write(*,*)'Frequency dependent dielectric:',frequency_dependent_dielectric
! Check for special case which is not currently used but could happen in future developments maybe...
if ( (frequency_dependent_dielectric).AND.(n_dielectric_regions.EQ.0) ) then
run_status='Error in PUL_LC_Laplace. Conductors are situated in a homogeneous frequency dependent dielectric'
CALL write_program_status()
STOP 1
end if
! Allocate the per-unit-length matrices
PUL%L%dim=matrix_dimension
ALLOCATE( PUL%L%mat(1:PUL%L%dim,1:PUL%L%dim) )
PUL%C%dim=matrix_dimension
ALLOCATE( PUL%C%mat(1:PUL%C%dim,1:PUL%C%dim) )
PUL%G%dim=matrix_dimension
ALLOCATE( PUL%G%mat(1:PUL%G%dim,1:PUL%G%dim) )
! Allocate the impedance and admittance filter matrices
PUL%Zfilter%dim=matrix_dimension
ALLOCATE(PUL%Zfilter%sfilter_mat(1:PUL%Zfilter%dim,1:PUL%Zfilter%dim))
PUL%Yfilter%dim=matrix_dimension
ALLOCATE(PUL%Yfilter%sfilter_mat(1:PUL%Yfilter%dim,1:PUL%Yfilter%dim))
if (verbose) write(*,*)'CALLING: Laplace'
if (n_dielectric_regions.EQ.0) then
! no dielectrics so we can calculate L, C and G in the same solution
f=1D0 ! arbitrary frequency here for lossless dielectric but a frequency must be set. The value will have no effect in this case.
CALL Laplace(mesh_filename,matrix_dimension,first_surface_is_free_space_boundary, &
n_dielectric_regions,dielectric_region_epsr,f,PUL%L%mat,PUL%C%mat,PUL%G%mat,xc,yc)
! calculate the impedance and admittance filters for a frequency independent model
CALL Z_Y_from_L_C(PUL%L,PUL%C,PUL%Zfilter,PUL%Yfilter)
else
! there are dielectrics present
ALLOCATE( mat1%mat(1:matrix_dimension,1:matrix_dimension) )
mat1%dim=matrix_dimension
ALLOCATE( mat2%mat(1:matrix_dimension,1:matrix_dimension) )
mat2%dim=matrix_dimension
! Firstly solve using the 'high frequency' dielectric material properties
dielectric_region=0
do i=1,end_conductor_loop ! don't include the ground plane in the dielectric region count
if (conductor_has_dielectric(i)) then
dielectric_region=dielectric_region+1
dielectric_region_epsr(dielectric_region)=evaluate_Sfilter_high_frequency_limit(PUL%epsr(i))
end if
end do ! next conductor
! call first with the dielectric materials and solve for C and G first
f=1D0 ! arbitrary frequency here for lossless dielectric
CALL Laplace(mesh_filename,matrix_dimension,first_surface_is_free_space_boundary, &
n_dielectric_regions,dielectric_region_epsr,f,mat1%mat,PUL%C%mat,PUL%G%mat,xc,yc)
! call second with all materials set to free space properties and solve for L
dielectric_region_epsr(:)=(1d0,0d0)
f=1D0 ! arbitrary frequency here for lossless dielectric
CALL Laplace(mesh_filename,matrix_dimension,first_surface_is_free_space_boundary, &
n_dielectric_regions,dielectric_region_epsr,f,PUL%L%mat,mat1%mat,mat2%mat,xc,yc)
if ((.NOT.frequency_dependent_dielectric).OR.(freq_spec%n_frequencies.EQ.1)) then
! calculate the impedance and admittance filters
CALL Z_Y_from_L_C(PUL%L,PUL%C,PUL%Zfilter,PUL%Yfilter)
else
! There are frequency dependent dielectrics present
! Secondly we calculate the capacitance matrix at a number of frequencies before fitting filter functions
! to each of the frequency dependent admittance matrix entries.
! allocate the frequency dependent matrices
ALLOCATE( C_freq(1:freq_spec%n_frequencies) )
ALLOCATE( G_freq(1:freq_spec%n_frequencies) )
! loop over frequency
do floop=1,freq_spec%n_frequencies
! allocate the matrices for this frequency
ALLOCATE( C_freq(floop)%mat(1:matrix_dimension,1:matrix_dimension) )
C_freq(floop)%dim=matrix_dimension
ALLOCATE( G_freq(floop)%mat(1:matrix_dimension,1:matrix_dimension) )
G_freq(floop)%dim=matrix_dimension
! evaluate the complex dielectric constant at the current frequency
f=freq_spec%freq_list(floop)
! set material properties in region 0, the background region
dielectric_region=0
dielectric_region_epsr(dielectric_region)=evaluate_Sfilter_frequency_response(PUL%epsr_background,f)
do i=1,end_conductor_loop ! don't include the ground plane in the dielectric region count
if (conductor_has_dielectric(i)) then
dielectric_region=dielectric_region+1
dielectric_region_epsr(dielectric_region)=evaluate_Sfilter_frequency_response(PUL%epsr(i),f)
end if
end do ! next conductor
! solve for C and G matrices at this frequency
! call Laplace to solve for C and G at this frequency
CALL Laplace(mesh_filename,matrix_dimension,first_surface_is_free_space_boundary, &
n_dielectric_regions,dielectric_region_epsr,f,mat1%mat,C_freq(floop)%mat,G_freq(floop)%mat,xc,yc)
end do ! next frequency
! we now have frequency domain C and G matrices.
! loop over the elements of C and G and calculate the Y matrix by filter fitting to Y=G(f)+jwC(f), and the Z matrix as jwL
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=G+jwC
do floop=1,freq_spec%n_frequencies
function_to_fit(floop)=G_freq(floop)%mat(row,col)+ &
j*2d0*pi*freq_spec%freq_list(floop)*C_freq(floop)%mat(row,col)
end do
! calculate the Yfilter 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%Yfilter%sfilter_mat(row,col),fit_order+1,-1,0)
! calcualte the Zfilter matrix as jwL
PUL%Zfilter%sfilter_mat(row,col)=jwA_filter( PUL%L%mat(row,col) )
if (col.NE.row) then
! make the matrix symmetrical
PUL%Zfilter%sfilter_mat(col,row)=PUL%Zfilter%sfilter_mat(row,col)
PUL%Yfilter%sfilter_mat(col,row)=PUL%Yfilter%sfilter_mat(row,col)
end if
end do ! next col
end do ! next row
DEALLOCATE( function_to_fit )
! deallocate the temporary matrices associated with the frequency dependent C and G matrices
do floop=1,freq_spec%n_frequencies
! allocate the matrices for this frequency
DEALLOCATE( C_freq(floop)%mat )
DEALLOCATE( G_freq(floop)%mat )
end do ! next frequency
DEALLOCATE( C_freq )
DEALLOCATE( G_freq )
end if ! frequency dependent dielectric
DEALLOCATE( mat1%mat )
DEALLOCATE( mat2%mat )
end if ! any dielectrics
if (verbose) then
write(*,*)'FINISHED: Laplace'
write(*,*)'Inductance matrix, L'
CALL dwrite_matrix(PUL%L%mat,matrix_dimension,matrix_dimension,matrix_dimension,0)
write(*,*)'Capacitance matrix, C'
CALL dwrite_matrix(PUL%C%mat,matrix_dimension,matrix_dimension,matrix_dimension,0)
write(*,*)'Conductance matrix, G'
CALL dwrite_matrix(PUL%G%mat,matrix_dimension,matrix_dimension,matrix_dimension,0)
end if
if (allocated( dielectric_region_epsr )) DEALLOCATE( dielectric_region_epsr )
if (verbose) write(*,*)'FINISHED PUL_LC_Laplace'
END SUBROUTINE PUL_LC_Laplace