! ! 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