!
! 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 <http://www.gnu.org/licenses/>.
! 
! 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 <http://www.gnu.org/licenses/>.
! 
! 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_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

    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(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),' };'
    
    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