!
! This file is part of SACAMOS, State of the Art CAble MOdels in 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-2017 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 PUL_parameter_module
! CONTAINS
! SUBROUTINE allocate_and_reset_PUL_data
! SUBROUTINE deallocate_PUL_data
! FUNCTION point_is_inside
! SUBROUTINE Z_Y_from_L_C
! SUBROUTINE write_Dshape_gmsh
!
! plus the following in include files:
!Laplace.F90: SUBROUTINE Laplace
!PUL_analytic.F90: SUBROUTINE PUL_LC_calc_wide_separation_approximation
!PUL_analytic.F90: SUBROUTINE PUL_LC_calc_overshield_wide_separation_approximation
!PUL_analytic.F90: SUBROUTINE calculate_height_over_ground_plane
!PUL_LC_Laplace.F90: SUBROUTINE PUL_LC_Laplace
!Aprod.F90: SUBROUTINE Aprod
!Aprod.F90: SUBROUTINE ATprod
!Aprod.F90: SUBROUTINE zAprod
!Aprod.F90: SUBROUTINE zATprod
!CG_solve.F90: SUBROUTINE solve_real_symm(n, b, x,tol)
!CG_solve.F90: SUBROUTINE solve_complex_symm(n, b, x,tol)
!
!
!
! NAME
! MODULE PUL_parameters
!
! Data relating to the calculation of PUL_parameters
!
! COMMENTS
! The conjugate gradient solutions are included in this module rather than the maths module
! as they are intimately linked with the matrix-vector product subroutine which in turn relies
! on the sparse matrix storage structure in this module
!
! HISTORY
! started 25/11/2015 CJS
! 30/3/2016 CJS fix factor of 2 in mutual inductance of wires over ground plane
! 18_4_2016 CJS include calculation for conductors within a cylindrical shield for oversheild domains
! 5_7_2016 CJS include numerical Laplace solver for L,C,G in inhomogeneous regions
! 13/3/2018CJS include iterative solver stuff based on Stanford University symmlq.f
! 16/3/2018 CJS add y offset for ML_flex_cable (PUL%ox, and PUL%oy)(PUL%rdox, and PUL%rdoy)
! 19/6/2018 CJS include iterative sparse matrix solver based on the biconjugate gradient method
!
MODULE PUL_parameter_module
USE type_specifications
USE filter_module
USE frequency_spec
USE general_module
IMPLICIT NONE
! Main structure used to hold data for the per-unit-length parameter calculation
TYPE::PUL_type
integer :: n_conductors ! number of conductors
! conductor based arrays:
integer,allocatable :: shape(:) ! holds a nuber which indicates the shape of the conductor
real(dp),allocatable :: r(:) ! radius of a circular conductor
real(dp),allocatable :: x(:) ! x coordinate of the centre of the cable to which this conductor belongs
real(dp),allocatable :: y(:) ! y coordinate of the centre of the cable to which this conductor belongs
real(dp),allocatable :: ox(:) ! offset in the x direction of this conductor from the cable centre (x():),y(:) above)
real(dp),allocatable :: oy(:) ! offset in the y direction of this conductor from the cable centre (x():),y(:) above)
real(dp),allocatable :: rh(:) ! height (y dimension) of rectangular conductor
real(dp),allocatable :: rw(:) ! width1 (x dimension) of rectangular conductor/ Dshape
real(dp),allocatable :: rw2(:) ! width2 (x dimension) of rectangular conductor/ Dshape
real(dp),allocatable :: rd(:) ! radius of dielectric surrounding a circular conductor
real(dp),allocatable :: rdh(:) ! height (y dimension) of rectangular dielectric around conductor
real(dp),allocatable :: rdw(:) ! width (x dimension) of rectangular dielectric around conductor
real(dp),allocatable :: rdox(:) ! offset of dielectric in the x direction of this conductor from the cable centre
real(dp),allocatable :: rdoy(:) ! offset of dielectric in the y direction of this conductor from the cable centre
real(dp),allocatable :: rtheta(:) ! rotation angle of conductor
type(Sfilter),allocatable :: epsr(:) ! relative permittivity of dielecrtric surrounding the conductor
logical :: ground_plane_present ! flag to indicate the presence of a ground plane
real(dp) :: ground_plane_angle ! angle of ground plane normal from the x axis
real(dp) :: ground_plane_offset ! ground plane offset
logical :: overshield_present ! flag to indicate the presence of an overshield
integer :: overshield_shape ! holds a nuber which indicates the shape of the opvershield
real(dp) :: overshield_x ! x coordinate of the centre of the overshield
real(dp) :: overshield_y ! y coordinate of the centre of the overshield
real(dp) :: overshield_r ! radius of a circular overshield
real(dp) :: overshield_h ! height (y dimension) of Dshape overshield
real(dp) :: overshield_w ! width1 (x) parameter of Dshape overshield
real(dp) :: overshield_w2 ! width2 (x) parameter of Dshape overshield
type(Sfilter) :: epsr_background ! relative permittivity of background dielecrtric (unually air)
type(matrix) :: L ! inductance matrix
type(matrix) :: C ! capacitance matrix
type(matrix) :: G ! conductance matrix
type(Sfilter_matrix):: Zfilter ! frequency dependent impedance matrix
type(Sfilter_matrix):: Yfilter ! frequency dependent impedance matrix
integer :: filter_fit_order ! order for filter fitting for frequency dependent models
type(frequency_specification) :: filter_fit_freq ! frequency range for filter fitting for frequency dependent models
END TYPE PUL_type
integer,public :: n_entry
! 1D arrays used in the construction of the K matrix ( K(i_K(:),j_K(:))=K(i_K(:),j_K(:))+s_K(:) )
integer,allocatable,public :: i_K(:)
integer,allocatable,public :: j_K(:)
real(dp),allocatable,public :: s_K_re(:)
complex(dp),allocatable,public :: s_K(:)
CONTAINS
include "PUL_analytic.F90"
include "PUL_LC_Laplace.F90"
include "Laplace.F90"
include "Aprod.F90"
include "CG_solve.F90"
! NAME
! SUBROUTINE allocate_and_reset_PUL_data
!
! allocate and reset per-unit-length calculation data
!
! COMMENTS
!
!
! HISTORY
! started 10/10/2015 CJS
!
SUBROUTINE allocate_and_reset_PUL_data(PUL,n_conductors)
USE type_specifications
USE general_module
IMPLICIT NONE
! variables passed to subroutine
type(PUL_type),intent(OUT) :: PUL
integer,intent(IN) :: n_conductors
! local variables
integer :: i
! START
! allocate memory for PUL structure
PUL%n_conductors=n_conductors
ALLOCATE( PUL%shape( PUL%n_conductors) )
ALLOCATE( PUL%x( PUL%n_conductors) )
ALLOCATE( PUL%y( PUL%n_conductors) )
ALLOCATE( PUL%r( PUL%n_conductors) )
ALLOCATE( PUL%ox( PUL%n_conductors) )
ALLOCATE( PUL%oy( PUL%n_conductors) )
ALLOCATE( PUL%rd( PUL%n_conductors) )
ALLOCATE( PUL%rdw( PUL%n_conductors) )
ALLOCATE( PUL%rdh( PUL%n_conductors) )
ALLOCATE( PUL%rdox( PUL%n_conductors) )
ALLOCATE( PUL%rdoy( PUL%n_conductors) )
ALLOCATE( PUL%epsr( PUL%n_conductors) )
ALLOCATE( PUL%rh( PUL%n_conductors) )
ALLOCATE( PUL%rw( PUL%n_conductors) )
ALLOCATE( PUL%rw2( PUL%n_conductors) )
ALLOCATE( PUL%rtheta( PUL%n_conductors) )
! reset all the elements of the PUL structure
do i=1,PUL%n_conductors
PUL%shape(i)=0
PUL%x(i)=0d0
PUL%y(i)=0d0
PUL%r(i)=0d0
PUL%ox(i)=0d0
PUL%oy(i)=0d0
PUL%rh(i)=0d0
PUL%rw(i)=0d0
PUL%rw2(i)=0d0
PUL%rd(i)=0d0
PUL%rdw(i)=0d0
PUL%rdh(i)=0d0
PUL%rdox(i)=0d0
PUL%rdoy(i)=0d0
PUL%rtheta(i)=0d0
PUL%epsr(i)=1d0
end do
PUL%ground_plane_present=.FALSE.
PUL%ground_plane_angle=0d0
PUL%ground_plane_offset=0d0
PUL%overshield_present=.FALSE.
PUL%overshield_shape=circle ! circle by default
PUL%overshield_x=0d0
PUL%overshield_y=0d0
PUL%overshield_r=0d0
PUL%overshield_w=0d0
PUL%overshield_w2=0d0
PUL%overshield_h=0d0
PUL%epsr_background=1d0
RETURN
END SUBROUTINE allocate_and_reset_PUL_data
! NAME
! SUBROUTINE deallocate_PUL_data
!
! deallocate per-unit-length calculation data
!
! COMMENTS
!
!
! HISTORY
! started 2/12/2015 CJS
! 5/6/2016 CJS include additional variables for numerical Laplace solution.
!
SUBROUTINE deallocate_PUL_data(PUL)
USE type_specifications
USE general_module
IMPLICIT NONE
! variables passed to subroutine
type(PUL_type),intent(INOUT) :: PUL
! local variables
integer :: i
! START
if (allocated(PUL%shape ) ) DEALLOCATE( PUL%shape )
if (allocated(PUL%x ) ) DEALLOCATE( PUL%x )
if (allocated(PUL%y ) ) DEALLOCATE( PUL%y )
if (allocated(PUL%r ) ) DEALLOCATE( PUL%r )
if (allocated(PUL%ox ) ) DEALLOCATE( PUL%ox )
if (allocated(PUL%oy ) ) DEALLOCATE( PUL%oy )
if (allocated(PUL%rh ) ) DEALLOCATE( PUL%rh )
if (allocated(PUL%rw ) ) DEALLOCATE( PUL%rw )
if (allocated(PUL%rw2 ) ) DEALLOCATE( PUL%rw2 )
if (allocated(PUL%rd ) ) DEALLOCATE( PUL%rd )
if (allocated(PUL%rdw ) ) DEALLOCATE( PUL%rdw )
if (allocated(PUL%rdh ) ) DEALLOCATE( PUL%rdh )
if (allocated(PUL%rdox ) ) DEALLOCATE( PUL%rdox )
if (allocated(PUL%rdoy ) ) DEALLOCATE( PUL%rdoy )
if (allocated(PUL%rtheta ) ) DEALLOCATE( PUL%rtheta )
if (allocated(PUL%epsr ) ) then
do i=1,PUL%n_conductors
CALL deallocate_Sfilter(PUL%epsr(i))
end do
DEALLOCATE( PUL%epsr )
end if
CALL deallocate_Sfilter(PUL%epsr_background)
if (allocated(PUL%L%mat) ) DEALLOCATE( PUL%L%mat )
if (allocated(PUL%C%mat) ) DEALLOCATE( PUL%C%mat )
if (allocated(PUL%G%mat) ) DEALLOCATE( PUL%G%mat )
CALL deallocate_Sfilter_matrix( PUL%Zfilter )
CALL deallocate_Sfilter_matrix( PUL%Yfilter )
RETURN
END SUBROUTINE deallocate_PUL_data
! NAME
! FUNCTION point_is_inside(xt,yt,shape_type,x,y,r,rh,rw,rtheta)
!
! function to determine whether conductors are inside dielectrics - used to create region boundary
! lists for mesh generation
!
! COMMENTS
!
!
! HISTORY
! started 6/10/2016 CJS
! add rox, roy instead of only the x offset, ro 16/3/2018 CJS
!
logical FUNCTION point_is_inside(xt,yt,shape_type,x,y,r,rh,rw,rox,roy,rtheta)
real(dp),intent(IN) :: xt ! test point x
real(dp),intent(IN) :: yt ! test point y
integer,intent(IN) :: shape_type ! holds a nuber which indicates the shape of the test object
real(dp),intent(IN) :: x ! x coordinate of the centre of the test object
real(dp),intent(IN) :: y ! y coordinate of the centre of the test object
real(dp),intent(IN) :: r ! radius of cylindrical test object
real(dp),intent(IN) :: rh ! height (y dimension) of rectangular test object
real(dp),intent(IN) :: rw ! width (x dimension) of test object
real(dp),intent(IN) :: rox ! offset in the x direction of this test object from the centre
real(dp),intent(IN) :: roy ! offset in the y direction of this test object from the centre
real(dp),intent(IN) :: rtheta ! rotation angle of test object
! local variables
real(dp):: d
real(dp):: xt_r,yt_r
real(dp):: xt_ro,yt_ro
real(dp):: xt_o,yt_o
real(dp):: dx,dy
! START
if (shape_type.EQ.circle) then
! calculate the distance between the test point and the centre of the dielectric cylinder
d=sqrt((x-xt)**2+(y-yt)**2)
if (d.LT.r) then
point_is_inside=.TRUE.
else
point_is_inside=.FALSE.
end if
else
! rectangular conductor
write(*,*)'Test point',xt,yt
! apply the inverse of the offset to the shape i.e. move it into the shifted, un-rotated coordinate system of the rectangle
xt_o=xt-x
yt_o=yt-y
write(*,*)'Offset on test point',xt_o,yt_o
! apply the inverse of the rotation to the test point i.e. move it into the un-rotated coordinate system of the rectangle
xt_r= xt_o*cos(rtheta)+yt_o*sin(rtheta)
yt_r=-xt_o*sin(rtheta)+yt_o*cos(rtheta)
write(*,*)'Inverse rotation on test point',xt_r,yt_r
! apply the inverse of the offset to the shape i.e. move it into the shifted, un-rotated coordinate system of the rectangle
xt_ro=xt_r-rox
yt_ro=yt_r-roy
write(*,*)'Offset Inverse rotation on test point',xt_ro,yt_ro
! xt_ro,yt_ro is the the distance from the transformed test point to the centre of the rectangle in x and y
dx=xt_ro
dy=yt_ro
! write(*,*)'dielectric centre point',x,y,' offset',rox,roy,' angle',rtheta
! write(*,*)'dx',dx,' rw/2',rw/2d0
! write(*,*)'dy',dy,' rh/2',rh/2d0
if ( (abs(dx).LT.rw/2d0).AND.(abs(dy).LT.rh/2d0) ) then
point_is_inside=.TRUE.
write(*,*)'POINT IS INSIDE'
else
point_is_inside=.FALSE.
write(*,*)'POINT IS OUTSIDE'
end if
end if
END FUNCTION point_is_inside
!
! NAME
! SUBROUTINE Z_Y_from_L_C(unit)
!
! calculate frequency dependent impedance and admittance (filter function) matrices
! for inductance and capacitance matrices with frequency independent properties
! Z=jwL Y=jwC
!
! COMMENTS
!
!
! HISTORY
! started 2/12/2015 CJS
!
SUBROUTINE Z_Y_from_L_C(L,C,Z,Y)
USE type_specifications
USE general_module
USE maths
USE filter_module
IMPLICIT NONE
! variables passed to subroutine
type(matrix),intent(IN) :: L
type(matrix),intent(IN) :: C
type(Sfilter_matrix),intent(OUT) :: Z
type(Sfilter_matrix),intent(OUT) :: Y
! local variables
integer :: dim
integer :: row,col
type(Sfilter) :: jw
! START
jw=jwA_filter(1d0)
! Impedance matrix
dim=L%dim
Z%dim=dim
if (.NOT.ALLOCATED(Z%sfilter_mat)) ALLOCATE(Z%sfilter_mat(dim,dim))
do row=1,dim
do col=1,dim
Z%sfilter_mat(row,col)=L%mat(row,col)*jw
end do
end do
! Capacitance matrix
dim=C%dim
Y%dim=dim
if (.NOT.ALLOCATED(Y%sfilter_mat)) ALLOCATE(Y%sfilter_mat(dim,dim))
do row=1,dim
do col=1,dim
Y%sfilter_mat(row,col)=C%mat(row,col)*jw
end do
end do
! Deallocate work filter
CALL deallocate_Sfilter(jw)
RETURN
END SUBROUTINE Z_Y_from_L_C
!
! NAME
! SUBROUTINE write_Dshape_gmsh
!
! write the points required for a Dshape connector shell model in gmsh
!
! COMMENTS
!
!
! HISTORY
! started 15/11/2016 CJS
!
SUBROUTINE write_Dshape_gmsh(x,y,w_in,w2_in,h_in,r,ox,oy,theta,dl,point,number,unit)
USE type_specifications
USE general_module
IMPLICIT NONE
! variables passed to subroutine
real(dp),intent(IN) :: x ! x coordinate of the centre of the Dshape
real(dp),intent(IN) :: y ! y coordinate of the centre of the Dshape
real(dp),intent(IN) :: w_in ! width1 (x dimension) of the Dshape
real(dp),intent(IN) :: w2_in ! width2 (x dimension) of the Dshape
real(dp),intent(IN) :: h_in ! height (x dimension) of the Dshape
real(dp),intent(IN) :: r ! radius of curves in Dshape
real(dp),intent(IN) :: ox ! x offset
real(dp),intent(IN) :: oy ! y offset
real(dp),intent(IN) :: theta ! rotation angle of Dshape
real(dp),intent(IN) :: dl ! edge length for the line segments making up the Dshape
integer,intent(INOUT) :: point ! point counter
integer,intent(IN) :: number ! number of the Dshape - used as a label in the gmsh file
integer,intent(IN) :: unit ! unit to write to
! local variables
real(dp) :: w1,w2,h
real(dp) :: vx,vy
real(dp) :: voxl,voyl
real(dp) :: norm
real(dp) :: xt,yt
real(dp) :: xp,yp
integer :: i
! START
w1=w_in/2d0-r
w2=w2_in/2d0-r
h=h_in/2d0-r
! vector from top left conductor to bottem left conductor
vx=w1-w2
vy=-2d0*h
! perpendicular vector to -x edge
norm=sqrt(vx*vx+vy*vy)
voxl=vy*r/norm
voyl=-vx*r/norm
write(unit,*)' // Dshape ',number
point=point+1
! POINT 1 ! top right
xt=w1+ox
yt=h+r+oy
xp=x+xt*cos(theta)-yt*sin(theta)
yp=y+xt*sin(theta)+yt*cos(theta)
write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };'
point=point+1
! POINT 2 ! top left
xt=-w1+ox
yt=h+r+oy
xp=x+xt*cos(theta)-yt*sin(theta)
yp=y+xt*sin(theta)+yt*cos(theta)
write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };'
point=point+1
! POINT 3 ! top left centre
xt=-w1+ox
yt=h+oy
xp=x+xt*cos(theta)-yt*sin(theta)
yp=y+xt*sin(theta)+yt*cos(theta)
write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };'
point=point+1
! POINT 4 ! top left edge
xt=-w1+ox+voxl
yt=h+voyl+oy
xp=x+xt*cos(theta)-yt*sin(theta)
yp=y+xt*sin(theta)+yt*cos(theta)
write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };'
point=point+1
! POINT 5 ! bottom left edge
xt=-w2+ox+voxl
yt=-h+voyl+oy
xp=x+xt*cos(theta)-yt*sin(theta)
yp=y+xt*sin(theta)+yt*cos(theta)
write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };'
point=point+1
! POINT 6 ! bottom left centre
xt=-w2+ox
yt=-h+oy
xp=x+xt*cos(theta)-yt*sin(theta)
yp=y+xt*sin(theta)+yt*cos(theta)
write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };'
point=point+1
! POINT 7 ! bottom left
xt=-w2+ox
yt=-(h+r)+oy
xp=x+xt*cos(theta)-yt*sin(theta)
yp=y+xt*sin(theta)+yt*cos(theta)
write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };'
point=point+1
! POINT 8 ! bottom right
xt=w2+ox
yt=-(h+r)+oy
xp=x+xt*cos(theta)-yt*sin(theta)
yp=y+xt*sin(theta)+yt*cos(theta)
write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };'
point=point+1
! POINT 9 ! bottom right centre
xt=w2+ox
yt=-h+oy
xp=x+xt*cos(theta)-yt*sin(theta)
yp=y+xt*sin(theta)+yt*cos(theta)
write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };'
point=point+1
! POINT 10 ! bottom right edge
xt=w2+ox-voxl
yt=-h+voyl+oy
xp=x+xt*cos(theta)-yt*sin(theta)
yp=y+xt*sin(theta)+yt*cos(theta)
write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };'
point=point+1
! POINT 11 ! top right edge
xt=w1+ox-voxl
yt=h+voyl+oy
xp=x+xt*cos(theta)-yt*sin(theta)
yp=y+xt*sin(theta)+yt*cos(theta)
write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };'
point=point+1
! POINT 12 ! top right centre
xt=w1+ox
yt=h+oy
xp=x+xt*cos(theta)-yt*sin(theta)
yp=y+xt*sin(theta)+yt*cos(theta)
write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };'
RETURN
END SUBROUTINE write_Dshape_gmsh
END MODULE PUL_parameter_module