Include software updates for infinite ground plane model, new flex cable and ite…

…rative solver in Laplace'

Include software updates for infinite ground plane model, new flex cable and ite…
…rative solver in Laplace'
Chris Smartt
1 parent 5f44abed
Showing 32 changed files with 2229 additions and 187 deletions Show diff stats
SRC/BUNDLE_DOMAIN_CREATION/create_global_domain_structure.F90
SRC/BUNDLE_DOMAIN_CREATION/plot_bundle_cross_section.F90
SRC/BUNDLE_DOMAIN_CREATION/set_conductor_positions_for_Einc.F90
SRC/CABLE_BUNDLE_MODULES/cable_bundle_module.F90
SRC/CABLE_MODULES/ML_flex_cable.F90
SRC/CABLE_MODULES/Makefile
SRC/CABLE_MODULES/cable_checks.F90
SRC/CABLE_MODULES/cable_module.F90
SRC/CABLE_MODULES/flex_cable.F90
SRC/GENERAL_MODULES/constants.F90
SRC/GENERAL_MODULES/general_module.F90
SRC/GENERAL_MODULES/generate_shapes.F90
SRC/GENERAL_MODULES/plot_geometry.F90
SRC/MATHS_MODULES/cmatrix.F90
SRC/MATHS_MODULES/dmatrix.F90
SRC/PUL_PARAMETER_CALCULATION/Aprod.F90
SRC/PUL_PARAMETER_CALCULATION/CG_solve.F90
SRC/PUL_PARAMETER_CALCULATION/Laplace.F90
SRC/PUL_PARAMETER_CALCULATION/Makefile
SRC/PUL_PARAMETER_CALCULATION/PUL_LC_Laplace.F90
@@ -72,6 +72,7 @@
 !                    so that we can work out the is_shield flag properly in all circumstances.
 !     18/10/2017 CJS: include 8b. Copy the cable based conductor labels to the bundle structure
 !     16/11/2017 CJS Include network synthesis process to replace s-domain transfer functions
+!    16/3/2018 CJS add y offset for ML_flex_cable
 !
 SUBROUTINE create_global_domain_structure(bundle)
@@ -994,7 +995,8 @@ USE maths
             PUL%y(conductor)=bundle%cable_y_offset(cable)   
             PUL%rtheta(conductor)=bundle%cable_angle(cable)
-            PUL%o(conductor)=bundle%cable(cable)%external_model(local_conductor)%conductor_ox
+            PUL%ox(conductor)=bundle%cable(cable)%external_model(local_conductor)%conductor_ox
+            PUL%oy(conductor)=bundle%cable(cable)%external_model(local_conductor)%conductor_oy
             PUL%r(conductor)=bundle%cable(cable)%external_model(local_conductor)%conductor_radius
             PUL%rw(conductor)=bundle%cable(cable)%external_model(local_conductor)%conductor_width
             PUL%rw2(conductor)=bundle%cable(cable)%external_model(local_conductor)%conductor_width2
@@ -1002,7 +1004,8 @@ USE maths
             PUL%rd(conductor)=bundle%cable(cable)%external_model(local_conductor)%dielectric_radius
             PUL%rdw(conductor)=bundle%cable(cable)%external_model(local_conductor)%dielectric_width
             PUL%rdh(conductor)=bundle%cable(cable)%external_model(local_conductor)%dielectric_height
-            PUL%rdo(conductor)=bundle%cable(cable)%external_model(local_conductor)%dielectric_ox
+            PUL%rdox(conductor)=bundle%cable(cable)%external_model(local_conductor)%dielectric_ox
+            PUL%rdoy(conductor)=bundle%cable(cable)%external_model(local_conductor)%dielectric_oy
             PUL%epsr(conductor)=bundle%cable(cable)%external_model(local_conductor)%dielectric_epsr 
           end if ! not a ground plane    
@@ -45,7 +45,7 @@
 !    14/6/2016 CJS  Create the arrays of x and y coordinates of each conductor which are used in the incident field excitation
 !    5/9/2016 CJS   include additional shielded cable types with surface impedance loss models
 !    20/4/2017 CJS  Separate the set_conductor_positions_for_Einc process from plot_bundle_cross_section and make use of generate_shapes.F90
-!    
+!    16/3/2018 CJS plot ML_flex_cable 
 !
 SUBROUTINE plot_bundle_cross_section(bundle)
@@ -151,6 +151,13 @@ USE maths
                       bundle%cable_x_offset(cable), &
                       bundle%cable_y_offset(cable), &
                       bundle%cable_angle(cable),xmin,xmax,ymin,ymax )
+                                                                                           
+    else if (bundle%cable(cable)%cable_type.EQ.cable_geometry_type_ML_flex_cable) then
+    
+      CALL ML_flex_cable_plot( bundle%cable(cable),                 &
+                      bundle%cable_x_offset(cable), &
+                      bundle%cable_y_offset(cable), &
+                      bundle%cable_angle(cable),xmin,xmax,ymin,ymax )
     else if (bundle%cable(cable)%cable_type.EQ.cable_geometry_type_dconnector) then
@@ -46,7 +46,7 @@
 !    14/6/2016 CJS  Create the arrays of x and y coordinates of each conductor which are used in the incident field excitation
 !    5/9/2016 CJS   include additional shielded cable types with surface impedance loss models
 !    20/4/2017 CJS Separate the set_conductor_positions_for_Einc process from plot_bundle_cross_section
-!
+!    16/3/2018 CJS add y offset for ML_flex_cable
 !
 SUBROUTINE set_conductor_positions_for_Einc(bundle)
@@ -98,7 +98,8 @@ USE maths
   do cable=1,bundle%n_cables
 ! loop over the conductors in the cable
-    if (bundle%cable(cable)%cable_type.EQ.cable_geometry_type_flex_cable) then
+    if ((bundle%cable(cable)%cable_type.EQ.cable_geometry_type_flex_cable).OR.    &
+        (bundle%cable(cable)%cable_type.EQ.cable_geometry_type_ML_flex_cable) ) then
 ! We include a fix here for flex cables which have to have their conductor offsets calculated here
 ! work out the centre of this conductor before rotation
@@ -107,7 +108,7 @@ USE maths
       do i=1,bundle%cable(cable)%tot_n_conductors
         FC_x0=bundle%cable(cable)%external_model(i)%conductor_ox
-        FC_y0=0.0    
+        FC_y0=bundle%cable(cable)%external_model(i)%conductor_oy
 ! work out the centre of this conductor when the flex cable is rotated and offset
         conductor=conductor+1
         bundle%conductor_x_offset(conductor)=FC_x0*cos(theta)-FC_y0*sin(theta)+bundle%cable_x_offset(cable)
@@ -892,7 +892,8 @@ CONTAINS
     theta=bundle%cable_angle(cable1)
     type1=bundle%cable(cable1)%cable_type
-    if (bundle%cable(cable1)%cable_type.NE.cable_geometry_type_flex_cable) then
+    if ((bundle%cable(cable1)%cable_type.NE.cable_geometry_type_flex_cable).AND. &
+        (bundle%cable(cable1)%cable_type.NE.cable_geometry_type_ML_flex_cable) )  then
       nec1=bundle%cable(cable1)%n_external_conductors
     else
       nec1=1
@@ -953,7 +954,8 @@ CONTAINS
         theta=bundle%cable_angle(cable2)
         type2=bundle%cable(cable2)%cable_type
-        if (bundle%cable(cable2)%cable_type.NE.cable_geometry_type_flex_cable) then
+        if ((bundle%cable(cable1)%cable_type.NE.cable_geometry_type_flex_cable).AND. &
+            (bundle%cable(cable1)%cable_type.NE.cable_geometry_type_ML_flex_cable) )  then
           nec2=bundle%cable(cable2)%n_external_conductors
         else
           nec2=1
@@ -0,0 +1,471 @@
+!
+! 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 ML_flex_cable_set_parameters
+! SUBROUTINE ML_flex_cable_set_internal_domain_information
+! SUBROUTINE ML_flex_cable_plot
+!
+! NAME
+!     ML_flex_cable_set_parameters
+!
+! AUTHORS
+!     Chris Smartt
+!
+! DESCRIPTION
+!     Set the overall parameters for a ML_flex_cable cable
+!
+! COMMENTS
+!      
+!
+! HISTORY
+!
+!     started 26/9/2016 CJS 
+!     16/11/2017 CJS Include network synthesis process to replace s-domain transfer functions
+!
+!
+SUBROUTINE ML_flex_cable_set_parameters(cable)
+
+USE type_specifications
+
+IMPLICIT NONE
+
+! variables passed to subroutine
+
+  type(cable_specification_type),intent(INOUT)    :: cable
+
+! local variables
+
+! START
+
+  cable%cable_type=cable_geometry_type_ML_flex_cable
+  cable%tot_n_conductors=0        ! this is set in the cable specification file
+  cable%tot_n_domains=1
+  cable%n_external_conductors=0   ! this is set in the cable specification file
+  cable%n_internal_conductors=0  
+  cable%n_internal_domains=0
+  cable%n_parameters=0
+  cable%n_dielectric_filters=1
+  cable%n_transfer_impedance_models=0
+  
+END SUBROUTINE ML_flex_cable_set_parameters
+!
+! NAME
+!     ML_flex_cable_set_internal_domain_information
+!
+! AUTHORS
+!     Chris Smartt
+!
+! DESCRIPTION
+!     Set the overall parameters for a ML_flex_cable cable
+!
+! COMMENTS
+!     Set the dimension of the domain transformation matrices to include an external reference conductor for the cable 
+!      
+!
+! HISTORY
+!
+!     started 13/4/2016 CJS 
+!     8/5/2017         CJS: Include references to Theory_Manual
+!
+!
+SUBROUTINE ML_flex_cable_set_internal_domain_information(cable)
+
+USE type_specifications
+USE general_module
+USE constants
+
+IMPLICIT NONE
+
+! variables passed to subroutine
+
+  type(cable_specification_type),intent(INOUT)    :: cable
+
+! local variables
+
+  integer :: nc
+  integer :: dim
+  integer :: i,prow,prow2,row2,row,col,conductivity_row
+  
+  real(dp) :: full_width,ML_flex_cable_xmin
+  
+! variables for cable parameter checks 
+  logical :: cable_spec_error
+  real(dp) :: wd      ! dielectric width
+  real(dp) :: hd      ! dielectric height
+  integer  :: nrows   ! number of rows of conductors
+  
+  real(dp) :: w      ! conductor width
+  real(dp) :: h      ! conductor height
+  real(dp) :: s      ! conductor separation
+  real(dp) :: ox     ! conductor offset x
+  real(dp) :: oy     ! conductor offset y
+  integer  :: ncrow  ! number of conductors in this row
+  integer  :: conductor  ! counter for conductors
+  
+  real(dp) :: w1,h1,w2,h2,ox1,ox2,oy1,oy2
+  
+  real(dp) :: sigma
+  
+  type(Sfilter) :: epsr
+  
+  integer :: nclocal
+  
+  character(LEN=error_message_length) :: message 
+  
+  character(LEN=2) :: conductor_string
+  
+  integer :: check_type
+
+! START
+
+  nc=cable%tot_n_conductors 
+
+! set the information related ot the number of conductors
+
+  cable%n_external_conductors=nc
+    
+! Set the domain decomposition matrices ! Theory_Manual_Eqn 6.17, 6.18
+
+! The dimension of the domain transformation matrices is the number of conductors+1
+  dim=nc+1
+
+  cable%MI%dim=dim
+  ALLOCATE(cable%MI%mat(dim,dim))
+  cable%MI%mat(:,:)=0d0
+  
+  cable%MV%dim=dim
+  ALLOCATE(cable%MV%mat(dim,dim))
+  cable%MV%mat(:,:)=0d0
+
+  do i=1,nc
+    row=i
+    col=row
+    cable%MI%mat(row,col)=1d0
+  end do
+  do i=1,dim
+    row=dim
+    col=i
+    cable%MI%mat(row,col)=1d0
+  end do
+
+  do i=1,nc
+    row=i
+    col=row
+    cable%MV%mat(row,col)=1d0
+    cable%MV%mat(row,dim)=-1d0
+  end do
+  cable%MV%mat(dim,dim)=1d0
+  
+! Set the local reference conductor numbering  
+  ALLOCATE( cable%local_reference_conductor(nc) )
+  cable%local_reference_conductor(1:nc)=0              ! external domain conductor, reference not known
+
+! Set the local domain information: include a reference conductor in the count
+  ALLOCATE( cable%local_domain_n_conductors(1:cable%tot_n_domains) )
+  cable%local_domain_n_conductors(1)=nc+1     ! external domain 
+
+! Read the parameter list and set the conductor information
+
+  nrows=NINT(cable%parameters(3))
+  
+! check that the number of conductors is consistent
+  prow=3
+  nclocal=0  
+  
+  do row=1,nrows   ! loop over rows of conductors
+  
+    ncrow=NINT(cable%parameters(prow+6))    
+    nclocal=nclocal+ncrow    
+    prow=prow+6
+  
+  end do ! next row of conductors
+  
+!  check that the number of conductors is consistent
+  if (nclocal.NE.nc) then 
+    write(message,*)' Nc in .cable file=',nc,' conductor count=',nclocal
+    run_status='ERROR in cable_model_builder, inconsistent conductor count for flex_cable:'   &
+                //trim(cable%cable_name)//'. '//trim(message)
+    CALL write_program_status()
+    STOP 1
+  end if
+
+! conductivity for the conductor loss models
+  conductivity_row=prow+1
+  sigma=cable%parameters(conductivity_row)   ! conductivity
+  
+  ALLOCATE( cable%external_model(cable%n_external_conductors) )
+  
+  prow=3
+  conductor=0      ! reset the conductor count
+   
+  do row=1,nrows   ! loop over rows of conductors
+  
+! get the parameters for this row of conductors
+    ox=cable%parameters(prow+1)   ! conductor offset in x
+    oy=cable%parameters(prow+2)   ! conductor offset in y
+    w=cable%parameters(prow+3)   ! conductor width
+    h=cable%parameters(prow+4)   ! conductor height
+    s=cable%parameters(prow+5)   ! conductor separation
+    ncrow=NINT(cable%parameters(prow+6))
+    
+    full_width=w*ncrow+s*(ncrow-1)
+    ML_flex_cable_xmin=-full_width/2d0
+    
+    do i=1,ncrow    ! loop over the conductors in this row
+    
+      conductor=conductor+1
+      
+! set the conductor impedance model for this conductor
+      cable%conductor_impedance(conductor)%impedance_model_type=impedance_model_type_rectangular_with_conductivity
+      cable%conductor_impedance(conductor)%width=w
+      cable%conductor_impedance(conductor)%height=h
+      cable%conductor_impedance(conductor)%conductivity=sigma
+       
+! Set the external domain conductor information
+      
+      CALL reset_external_conductor_model(cable%external_model(conductor))
+      cable%external_model(conductor)%conductor_type=rectangle
+      cable%external_model(conductor)%conductor_width=w
+      cable%external_model(conductor)%conductor_width2=w
+      cable%external_model(conductor)%conductor_height=h
+
+! work out the offset of the ith conductor
+      cable%external_model(conductor)%conductor_ox=ox+ML_flex_cable_xmin+(w+s)*(i-1)+w/2d0
+      cable%external_model(conductor)%conductor_oy=oy
+   
+    end do ! next conductor in this row
+        
+    prow=prow+6
+  
+  end do ! next row of conductors
+
+! dielectric data
+
+  epsr=cable%dielectric_filter(1)
+  wd=cable%parameters(1)   ! dielectric width
+  hd=cable%parameters(2)   ! dielectric height
+
+! do some consistency checks
+  cable_spec_error=.FALSE.    ! assume no errors initially
+  message=''
+    
+! Do some intersection checks on the flex cable 
+
+  prow=3
+  do row=1,nrows   ! loop over rows of conductors
+  
+! get the parameters for this row of conductors
+    ox1=cable%parameters(prow+1)   ! conductor offset in x
+    oy1=cable%parameters(prow+2)   ! conductor offset in y
+    w=cable%parameters(prow+3)   ! conductor width
+    h1=cable%parameters(prow+4)   ! conductor height
+    s=cable%parameters(prow+5)   ! conductor separation
+    ncrow=NINT(cable%parameters(prow+6))
+    
+    w1=w*ncrow+s*(ncrow-1)  ! w1 is the full width of the row of conductors
+    
+! check whether this row of conductors is well specified - conductor width, height and spearation >0
+    CALL flex_cable_check(ncrow,w,h1,s,0d0,0d0,cable_spec_error,cable%cable_name,message)
+
+    if (cable_spec_error) then       
+      run_status='ERROR in cable_model_builder, error on parameters for cable:'//trim(cable%cable_name)//'. '//trim(message)
+      CALL write_program_status()
+      STOP 1
+    end if
+    
+! check whether this row of conductors intersects the dielectric boundary of the flex cable
+    check_type=2
+    CALL ML_flex_cable_check(w1,h1,ox1,oy1,wd,hd,0d0,0d0,check_type,cable_spec_error,cable%cable_name,message)
+   
+    if (cable_spec_error) then       
+      run_status='ERROR in cable_model_builder, error on parameters for cable:'//trim(cable%cable_name)//'. '//trim(message)
+      CALL write_program_status()
+      STOP 1
+    end if
+    
+! loop over all the other rows to check for intersections between the rows of conductors
+        
+    prow=prow+6
+    
+    prow2=prow
+    
+    do row2=row+1,nrows
+    
+! get the parameters for this row of conductors
+      ox2=cable%parameters(prow2+1)   ! conductor offset in x
+      oy2=cable%parameters(prow2+2)   ! conductor offset in y
+      w=cable%parameters(prow2+3)   ! conductor width
+      h2=cable%parameters(prow2+4)   ! conductor height
+      s=cable%parameters(prow2+5)   ! conductor separation
+      ncrow=NINT(cable%parameters(prow2+6))
+    
+      w2=w*ncrow+s*(ncrow-1)  ! w2 is the full width of the second row of conductors
+    
+    
+! check whether the two rows of conductors intersect
+      check_type=1
+      CALL ML_flex_cable_check(w1,h1,ox1,oy1,w2,h2,ox2,oy2,check_type,cable_spec_error,cable%cable_name,message)
+   
+      if (cable_spec_error) then       
+        run_status='ERROR in cable_model_builder, error on parameters for cable:'//trim(cable%cable_name)//'. '//trim(message)
+        CALL write_program_status()
+        STOP 1
+      end if 
+    
+      prow2=prow2+6
+    
+    end do ! next row2 for checking intersection between rows of conductors
+  
+  end do ! next row of conductors to check
+  
+  CALL dielectric_check(epsr,cable_spec_error,cable%cable_name,message)
+   
+  if (cable_spec_error) then       
+    run_status='ERROR in cable_model_builder, error on parameters for cable:'//trim(cable%cable_name)//'. '//trim(message)
+    CALL write_program_status()
+    STOP 1
+  end if
+  
+! add a dielectric region to the first conductor which encloses the whole cable
+
+! write the dielectric which is offset from the conductors      
+! The dielectric is centred at the cable centre 
+  
+  cable%external_model(1)%dielectric_width=wd
+  cable%external_model(1)%dielectric_height=hd
+  cable%external_model(1)%dielectric_ox=0d0
+  cable%external_model(1)%dielectric_oy=0d0
+  cable%external_model(1)%dielectric_epsr=epsr
+
+  CALL deallocate_Sfilter(epsr)
+  
+  ALLOCATE( cable%conductor_label(1:cable%tot_n_conductors) )
+  do i=1,cable%tot_n_conductors
+    write(conductor_string,'(I2)')i
+    cable%conductor_label(i)='Cable name: '//trim(cable%cable_name)//   &
+    '. type: '//trim(cable%cable_type_string)//'. conductor '//conductor_string//' : ML_flex_cable conductor'
+  end do
+    
+END SUBROUTINE ML_flex_cable_set_internal_domain_information
+!
+! NAME
+!     ML_flex_cable_plot
+!
+! AUTHORS
+!     Chris Smartt
+!
+! DESCRIPTION
+!     plot ML_flex_cable cable
+!
+! COMMENTS
+!      
+!
+! HISTORY
+!
+!     started 23/9/2016 CJS 
+!
+!
+SUBROUTINE ML_flex_cable_plot(cable,x_offset,y_offset,theta,xmin,xmax,ymin,ymax)
+
+USE type_specifications
+USE general_module
+
+IMPLICIT NONE
+
+! variables passed to subroutine
+
+  type(cable_specification_type),intent(IN)    :: cable
+  
+  real(dp),intent(IN) :: x_offset,y_offset,theta
+  real(dp),intent(INOUT) ::  xmin,xmax,ymin,ymax
+  
+! local variables
+
+  integer nc
+
+  real(dp) :: full_width,ML_flex_cable_xmin
+  real(dp) :: xoff,yoff,x,y,w,h,s,wd,hd,ox,oy
+  
+  integer nrows,ncrow,row,prow
+  integer i
+
+! START
+
+! plot ML_flex_cable conductor
+
+  nc=cable%tot_n_conductors
+  
+  nrows=NINT(cable%parameters(3))
+  
+  prow=3
+   
+  do row=1,nrows   ! loop over rows of conductors
+  
+! get the parameters for this row of conductors
+    ox=cable%parameters(prow+1)   ! conductor offset in x
+    oy=cable%parameters(prow+2)   ! conductor offset in y
+    w=cable%parameters(prow+3)   ! conductor width
+    h=cable%parameters(prow+4)   ! conductor height
+    s=cable%parameters(prow+5)   ! conductor separation
+    ncrow=NINT(cable%parameters(prow+6))
+    
+    full_width=w*ncrow+s*(ncrow-1)
+    ML_flex_cable_xmin=-full_width/2d0
+    
+    do i=1,ncrow    ! loop over the conductors in this row
+    
+! work out the centre of this conductor before rotation
+      xoff=ox+ML_flex_cable_xmin+(w+s)*(i-1)+w/2d0
+      yoff=oy
+    
+! work out the centre of this conductor when the flex cable is rotated and offset
+      x=x_offset+xoff*cos(theta)-yoff*sin(theta)
+      y=y_offset+xoff*sin(theta)+yoff*cos(theta)
+    
+! write the conductor    
+      CALL write_rectangle(x,y,w,h,theta,conductor_geometry_file_unit,xmin,xmax,ymin,ymax)
+        
+    end do ! next conductor in this row
+        
+    prow=prow+6
+  
+  end do ! next row of conductors
+
+! write the dielectric which is offset from the conductors    
+  
+  wd=cable%external_model(1)%dielectric_width
+  hd=cable%external_model(1)%dielectric_height
+  CALL write_rectangle(x_offset,y_offset,wd,hd,theta,dielectric_geometry_file_unit,xmin,xmax,ymin,ymax)
+   
+  RETURN
+  
+END SUBROUTINE ML_flex_cable_plot
+
+
+
+
@@ -8,6 +8,7 @@ $(OBJ_MOD_DIR)/cable_module.o: cable_module.F90  $(TYPE_SPEC_MODULE) $(GENERAL_M
                                twinax.F90 \
                                spacewire.F90 \
                                flex_cable.F90 \
+                               ML_flex_cable.F90 \
                                Dconnector.F90 \
                                overshield.F90 \
                                ground_plane.F90 \
@@ -627,6 +627,142 @@ RETURN
 END SUBROUTINE flex_cable_check
 !
 ! NAME
+!     ML_flex_cable_check
+!
+! AUTHORS
+!     Chris Smartt
+!
+! DESCRIPTION
+!     check that the parameters defining a multi-layer flex cable are consistent
+!     This check takes the width, height and offset of two rectangles and checks whether they intersect
+!     It is used for checking intersection of rows of conductors and conducotrs and dielectric in flex cable models
+!     If the check fails the cable_spec_error flag is set on return, otherwise it is left unchanged
+!
+! COMMENTS
+!      
+!
+! HISTORY
+!
+!  
+!     started 12/6/2018 CJS, based on flex_cable_check
+!
+!
+SUBROUTINE ML_flex_cable_check(w1,h1,ox1,oy1,w2,h2,ox2,oy2,check_type,cable_spec_error,cable_name,message)
+
+  real(dp),intent(IN) :: w1,h1
+  real(dp),intent(IN) :: ox1,oy1
+  real(dp),intent(IN) :: w2,h2
+  real(dp),intent(IN) :: ox2,oy2
+  
+  integer,intent(IN) :: check_type
+
+  logical,intent(INOUT)  :: cable_spec_error
+  
+  character(LEN=line_length),intent(IN) :: cable_name  
+  character(LEN=error_message_length),intent(INOUT) :: message 
+
+! local variables
+
+logical :: intersect,internal
+
+real(dp) :: p1,p2,p3,p4
+logical :: internalp3,internalp4
+
+logical :: intersectx,internalx,externalx
+logical :: intersecty,internaly,externaly
+
+! START
+
+if (cable_spec_error) RETURN ! return if an error has already been flagged
+
+intersect=.FALSE.
+internal=.FALSE.
+intersectx=.FALSE.
+internalx=.FALSE.
+externalx=.FALSE.
+intersecty=.FALSE.
+internaly=.FALSE.
+externaly=.FALSE.
+
+p1=ox1-w1/2.0
+p2=ox1+w1/2.0
+p3=ox2-w2/2.0
+p4=ox2+w2/2.0
+
+internalp3=.FALSE.
+internalp4=.FALSE.
+
+if ((p3.GE.p1).AND.(p3.LE.p2)) internalp3=.TRUE.
+if ((p4.GE.p1).AND.(p4.LE.p2)) internalp4=.TRUE.
+
+if (internalp3.AND.internalp4) then
+  intersectx=.TRUE.
+  internalx=.TRUE.                        ! rectangle 2 is inside rectangel1
+else if (internalp3.OR.internalp4) then
+  intersectx=.TRUE.
+  internalx=.FALSE.
+else
+! p3 and p4 are outside the p1-p2 range, we now nead to decide if p1-p2 is inside the range p3-p4
+  if ((p3.LT.p1).AND.(p4.GT.p2)) then
+    intersectx=.TRUE.
+    externalx=.TRUE.                      ! rectangle 2 is outside rectangel1
+  end if
+end if
+
+p1=oy1-h1/2.0
+p2=oy1+h1/2.0
+p3=oy2-h2/2.0
+p4=oy2+h2/2.0
+
+internalp3=.FALSE.
+internalp4=.FALSE.
+
+if ((p3.GE.p1).AND.(p3.LE.p2)) internalp3=.TRUE.
+if ((p4.GE.p1).AND.(p4.LE.p2)) internalp4=.TRUE.
+
+if (internalp3.AND.internalp4) then
+  intersecty=.TRUE.
+  internaly=.TRUE.                        ! rectangle 2 is inside rectangel1
+else if (internalp3.OR.internalp4) then
+  intersecty=.TRUE.
+  internaly=.FALSE.
+else
+! p3 and p4 are outside the p1-p2 range, we now nead to decide if p1-p2 is inside the range p3-p4
+  if ((p3.LT.p1).AND.(p4.GT.p2)) then
+    intersecty=.TRUE.
+    externaly=.TRUE.                      ! rectangle 2 is outside rectangel1
+  end if
+end if
+
+if (intersectx.AND.intersecty) intersect=.TRUE.
+if (internalx.AND.internaly) internal=.TRUE.
+if (externalx.AND.externaly) internal=.TRUE. ! rectangle 2 is outside rectangel1
+
+if ( (check_type.EQ.1).AND.(intersect) ) then
+! we are checking conductor intersections
+  write(*,*)'Error in cable:',trim(cable_name)
+  message='conductors intersect'
+  write(*,'(A)')trim(message)
+  cable_spec_error=.TRUE.
+  RETURN
+
+end if
+
+! this is the chek for dielectrc being outside the condustor row, the dielectric is shape 2 when 
+! the subroutine is called
+if ( (check_type.EQ.2).AND.( .NOT.(externalx.AND.externaly) ) ) then
+  write(*,*)'Error in cable:',trim(cable_name)
+  message='Dielectric intersects conductors'
+  write(*,'(A)')trim(message)
+  cable_spec_error=.TRUE.
+  RETURN
+end if 
+
+RETURN
+
+END SUBROUTINE ML_flex_cable_check    
+!
+! NAME
 !     dielectric_check
 !
 ! AUTHORS
@@ -92,6 +92,10 @@
 ! flex_cable.F90:  SUBROUTINE flex_cable_set_internal_domain_information
 ! flex_cable.F90:  SUBROUTINE flex_cable_plot
 ! 
+! ML_flex_cable.F90:  SUBROUTINE ML_flex_cable_set_parameters
+! ML_flex_cable.F90:  SUBROUTINE ML_flex_cable_set_internal_domain_information
+! ML_flex_cable.F90:  SUBROUTINE ML_flex_cable_plot
+! 
 ! Dconnector.F90:  SUBROUTINE Dconnector_set_parameters
 ! Dconnector.F90:  SUBROUTINE Dconnector_set_internal_domain_information
 ! Dconnector.F90:  SUBROUTINE Dconnector_plot
@@ -120,6 +124,7 @@
 !    Include general conductor impedance model 12/05/2016 CJS
 !    put the external condcutor model information into its own structure 6/10/2016 CJS
 !     16/11/2017 CJS Include network synthesis process to replace s-domain transfer functions
+!     16/3/2018 CJS Add ML_flex_cable
 !
 MODULE cable_module
@@ -171,10 +176,12 @@ TYPE:: external_conductor_model
   real(dp) :: conductor_width2
   real(dp) :: conductor_height
   real(dp) :: conductor_ox
+  real(dp) :: conductor_oy
   real(dp) :: dielectric_radius
   real(dp) :: dielectric_width
   real(dp) :: dielectric_height
   real(dp) :: dielectric_ox
+  real(dp) :: dielectric_oy
   type(Sfilter) :: dielectric_epsr
 END TYPE external_conductor_model
@@ -261,6 +268,7 @@ integer,parameter    :: cable_geometry_type_twinax                   =7
 integer,parameter    :: cable_geometry_type_flex_cable               =8
 integer,parameter    :: cable_geometry_type_Dconnector               =9
 integer,parameter    :: cable_geometry_type_ground_plane             =10 
+integer,parameter    :: cable_geometry_type_ML_flex_cable            =11
 integer,parameter    :: impedance_model_type_PEC                                =0
 integer,parameter    :: impedance_model_type_cylindrical_with_conductivity      =1
@@ -280,6 +288,7 @@ include &#39;shielded_twisted_pair.F90&#39;
 include 'spacewire.F90'
 include 'twinax.F90'
 include 'flex_cable.F90'
+include 'ML_flex_cable.F90'
 include 'Dconnector.F90'
 include 'overshield.F90'
 include 'ground_plane.F90'
@@ -329,10 +338,12 @@ include &#39;conductor_impedance_model.F90&#39;
   ext_model%conductor_width2=0d0
   ext_model%conductor_height=0d0
   ext_model%conductor_ox=0d0
+  ext_model%conductor_oy=0d0
   ext_model%dielectric_radius=0d0
   ext_model%dielectric_width=0d0
   ext_model%dielectric_height=0d0
   ext_model%dielectric_ox=0d0
+  ext_model%dielectric_oy=0d0
   ext_model%dielectric_epsr=1d0
 END SUBROUTINE reset_external_conductor_model
@@ -372,6 +383,7 @@ END SUBROUTINE reset_external_conductor_model
   integer        :: matrix_dimension
   logical        :: file_exists
+  character(len=line_length)    :: line
 ! START
@@ -502,11 +514,27 @@ END SUBROUTINE reset_external_conductor_model
     read(file_unit,*) cable%external_model(i)%conductor_width
     read(file_unit,*) cable%external_model(i)%conductor_width2
     read(file_unit,*) cable%external_model(i)%conductor_height
-    read(file_unit,*) cable%external_model(i)%conductor_ox
+
+! format change. This way of reading allows backward compatibility    
+    read(file_unit,'(A)')line    
+    read(line,*,err=100) cable%external_model(i)%conductor_ox,cable%external_model(i)%conductor_oy
+    GOTO 110
+100 read(line,*) cable%external_model(i)%conductor_ox
+    cable%external_model(i)%conductor_oy=0d0
+110 CONTINUE
+    
     read(file_unit,*) cable%external_model(i)%dielectric_radius
     read(file_unit,*) cable%external_model(i)%dielectric_width
     read(file_unit,*) cable%external_model(i)%dielectric_height
-    read(file_unit,*) cable%external_model(i)%dielectric_ox
+
+! format change. This way of reading allows backward compatibility    
+    read(file_unit,'(A)')line    
+    read(line,*,err=200) cable%external_model(i)%dielectric_ox,cable%external_model(i)%dielectric_oy
+    GOTO 210
+200 read(line,*) cable%external_model(i)%dielectric_ox
+    cable%external_model(i)%dielectric_oy=0d0
+210 CONTINUE
+   
     CALL read_Sfilter(cable%external_model(i)%dielectric_epsr,cable_file_unit)
   end do
@@ -677,11 +705,13 @@ END SUBROUTINE reset_external_conductor_model
     write(file_unit,*) cable%external_model(i)%conductor_width   ,' conductor_width  '
     write(file_unit,*) cable%external_model(i)%conductor_width2   ,' conductor_width2  '
     write(file_unit,*) cable%external_model(i)%conductor_height  ,' conductor_height '
-    write(file_unit,*) cable%external_model(i)%conductor_ox      ,' conductor_ox     '
+    write(file_unit,*) cable%external_model(i)%conductor_ox,cable%external_model(i)%conductor_oy     &
+                      ,' conductor_ox,  conductor_oy    '
     write(file_unit,*) cable%external_model(i)%dielectric_radius ,' dielectric_radius'
     write(file_unit,*) cable%external_model(i)%dielectric_width  ,' dielectric_width '
     write(file_unit,*) cable%external_model(i)%dielectric_height ,' dielectric_height'
-    write(file_unit,*) cable%external_model(i)%dielectric_ox     ,' dielectric_ox    '
+    write(file_unit,*) cable%external_model(i)%dielectric_ox,cable%external_model(i)%dielectric_oy   &
+                      ,' dielectric_ox, dielectric_oy    '
     CALL write_Sfilter(cable%external_model(i)%dielectric_epsr,cable_file_unit)
   end do
@@ -249,7 +249,7 @@ IMPLICIT NONE
   do i=1,cable%tot_n_conductors
     write(conductor_string,'(I2)')i
     cable%conductor_label(i)='Cable name: '//trim(cable%cable_name)//   &
-    '. type: '//trim(cable%cable_type_string)//'. conductor '//conductor_string//' : flex_cable conductor'
+    '. type: '//trim(cable%cable_type_string)//'. conductor '//conductor_string//' : original_flex_cable conductor'
   end do
 END SUBROUTINE flex_cable_set_internal_domain_information
@@ -104,6 +104,8 @@ IMPLICIT NONE
   real(dp) ::  Laplace_surface_mesh_constant=3d0 ! the mesh edge length on boundaries is calculated as radius/Laplace_surface_mesh_constant 
                                                  ! if Laplace_surface_mesh_constant=5 we have just over 30 elements on the circumference
+                                           
+  real(dp) ::  max_mesh_edge_length=1d30          ! the maximum mesh edge length on internal boundaries
   real(dp) ::  Twisted_pair_equivalent_radius=1.5d0 ! The twisted pair commmon mode interaction is calculated by treating the
                                                     ! common mode as being carried on an 'equivalent cylindrical conductor'
@@ -132,6 +132,8 @@ integer,parameter :: Pspice =3
 integer,parameter :: circle =1 
 integer,parameter :: rectangle=2
 integer,parameter :: Dshape=3
+integer,parameter :: semicircle=4
+integer,parameter :: semicircle_diameter=5
 integer :: spice_version 
@@ -168,6 +170,13 @@ logical :: verbose=.FALSE.
 logical :: plot_propagation_correction_filter_fit_data=.FALSE.
+logical :: direct_solver=.FALSE.
+
+logical :: inf_gnd=.TRUE.
+
+!logical :: use_ABC=.TRUE.
+logical :: use_ABC=.FALSE.
+
 ! There are some differences between windows and unix relating to
 ! file separators and making directories. The different formats for
 ! the two ooperating systems supported are included below 
@@ -30,6 +30,7 @@
 !SUBROUTINE generate_rectangle_points
 !SUBROUTINE generate_Dshape_points
 !SUBROUTINE generate_arc_points
+!SUBROUTINE generate_semicircle_points
 !
 !SUBROUTINE generate_circle_points
 !
@@ -545,3 +546,76 @@ IMPLICIT NONE
   RETURN
 END SUBROUTINE generate_line_points
+!
+!SUBROUTINE generate_semicircle_points
+!
+! NAME
+!     SUBROUTINE generate_semicircle_points
+!
+! DESCRIPTION
+!     write a semicircle with specified x,y centre and radius to file for plotting with gnuplot
+!     the semicircle is in the region y>=ycentre
+!
+!     
+! COMMENTS
+!     return the extent of the plotting area...
+!
+! HISTORY
+!
+!     started 10/05/2013 CJS
+!
+!
+SUBROUTINE generate_semicircle_points(npts,shape_x,shape_y,x,y,r)
+   
+USE type_specifications
+USE constants
+
+IMPLICIT NONE
+
+  integer,intent(OUT)                :: npts
+  real(dp),allocatable,intent(INOUT)  :: shape_x(:)
+  real(dp),allocatable,intent(INOUT)  :: shape_y(:)
+  real(dp),intent(IN)                :: x,y,r                 ! centre x and y coordinates and radius
+  
+! local variables  
+
+  real(dp) x1,y1
+  real(dp) x2,y2
+  real(dp) x3,y3
+  
+  integer :: loop
+  
+! START
+
+! write the semicircle as two arcs and a straight line to close the shape
+! the straght line is formed by the first and last points of the arc
+
+  x1=x+r
+  y1=y
+  
+  x2=x
+  y2=y+r
+  
+  x3=x-r
+  y3=y
+
+  do loop=1,2  ! first pass to count the points, second pass to set the point coordinates
+    
+    npts=0
+    
+    CALL generate_arc_points(npts,shape_x,shape_y,loop,x,y,x1,y1,x2,y2)
+  
+    CALL generate_arc_points(npts,shape_x,shape_y,loop,x,y,x2,y2,x3,y3)
+        
+    CALL generate_line_points(npts,shape_x,shape_y,loop,x3,y3,x1,y1)
+
+    if (loop.EQ.1) then
+      ALLOCATE( shape_x(1:npts) )
+      ALLOCATE( shape_y(1:npts) )
+    end if
+  
+  end do ! next loop
+   
+  RETURN
+  
+END SUBROUTINE generate_semicircle_points
@@ -29,6 +29,7 @@
 !SUBROUTINE write_circle
 !SUBROUTINE write_rectangle
 !SUBROUTINE write_Dshape
+!SUBROUTINE write_semicircle
 !
 !SUBROUTINE write_circle
 !
@@ -220,3 +221,67 @@ IMPLICIT NONE
   RETURN
 END SUBROUTINE write_Dshape
+!
+!SUBROUTINE write_semicircle
+!
+! NAME
+!     SUBROUTINE write_semicircle
+!
+! DESCRIPTION
+!     write a semicircle with specified x,y centre and radius to file for plotting with gnuplot
+!     the semicircle is in the region y>=0
+!
+!     
+! COMMENTS
+!     return the extent of the plotting area...
+!
+! HISTORY
+!
+!     started 10/05/2013 CJS
+!     using generate_shapes.F90 20/4/2017 CJS
+!
+!
+SUBROUTINE write_semicircle(x,y,r,unit,xmin,xmax,ymin,ymax)
+   
+USE type_specifications
+USE constants
+
+IMPLICIT NONE
+
+  real(dp),intent(IN)    :: x,y,r                 ! centre x and y coordinates and radius
+  real(dp),intent(INOUT) :: xmin,xmax,ymin,ymax   ! extent of the plotting area
+  integer, intent(IN)    :: unit                  ! unit to write to  
+  
+! local variables  
+  
+  integer :: npts
+  real(dp),allocatable :: shape_x(:)
+  real(dp),allocatable :: shape_y(:)
+  
+  integer :: i
+  
+! START
+
+  CALL generate_semicircle_points(npts,shape_x,shape_y,x,y,r)
+  
+  do i=1,npts
+    
+    write(unit,8000)shape_x(i),shape_y(i)
+8000 format (4E14.6)
+
+    xmin=min(xmin,shape_x(i))
+    xmax=max(xmax,shape_x(i))
+    ymin=min(ymin,shape_y(i))
+    ymax=max(ymax,shape_y(i))
+  
+  end do
+  
+  write(unit,*)
+  write(unit,*)
+
+  DEALLOCATE( shape_x )
+  DEALLOCATE( shape_y )
+   
+  RETURN
+  
+END SUBROUTINE write_semicircle
@@ -32,6 +32,7 @@
 !       SUBROUTINE write_cmatrix_abs(Mat,dim,unit)
 !       SUBROUTINE cinvert_Gauss_Jordan(A,n,AI,dim) 
 !       SUBROUTINE c_condition_number (A,n,condition_number,dim)
+!
 !     16/11/2017 CJS Include network synthesis process to replace s-domain transfer functions
 !
 ! NAME
@@ -454,3 +455,5 @@ IMPLICIT NONE
   RETURN
   END SUBROUTINE c_condition_number
+
+
@@ -28,6 +28,7 @@
 ! SUBROUTINE dwrite_matrix(a,ar,ac,dim,unit)
 ! SUBROUTINE dread_matrix(a,ar,ac,dim,unit)
 ! SUBROUTINE dinvert_Gauss_Jordan(A,ar,AI,dim) 
+!
 !_____________________________________________________________________
 !
 ! NAME
@@ -283,3 +284,4 @@ IMPLICIT NONE
   RETURN
   END SUBROUTINE dinvert_Gauss_Jordan
+
@@ -0,0 +1,207 @@
+! 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 Aprod
+!     SUBROUTINE ATprod
+!     SUBROUTINE zAprod
+!     SUBROUTINE zATprod
+!
+! NAME
+!     SUBROUTINE Aprod
+!
+! DESCRIPTION
+!     Matrix vector multiplication for sparse, real matrices
+!
+! COMMENTS
+!     
+!
+! HISTORY
+!    19/3/2018 CJS 
+!
+!
+SUBROUTINE Aprod ( n, x, y )
+
+! form the vector y=K*x
+
+USE type_specifications
+
+IMPLICIT NONE
+
+integer n
+real(dp) :: x(n), y(n)
+
+integer :: i,row,col
+
+! START
+
+y(:)=0d0
+
+do i=1,n_entry
+  row=i_K(i)
+  col=j_K(i)
+  
+!  if (row.GT.n) then
+!    write(*,*)'Error in Aprod, row=',row,' n=',n
+!  end if
+!  if (col.GT.n) then
+!    write(*,*)'Error in Aprod, col=',col,' n=',n
+!  end if
+  
+  y(row)=y(row)+s_K_re(i)*x(col)
+  
+end do
+
+RETURN
+
+END SUBROUTINE Aprod 
+!
+! NAME
+!     SUBROUTINE ATprod
+!
+! DESCRIPTION
+!     Transpose Matrix vector multiplication for sparse, real matrices
+!
+! COMMENTS
+!     
+!
+! HISTORY
+!    19/3/2018 CJS 
+!
+!
+SUBROUTINE ATprod ( n, x, y )
+
+! form the vector y=AT*x
+
+USE type_specifications
+
+IMPLICIT NONE
+
+integer n
+real(dp) :: x(n), y(n)
+
+integer :: i,row,col
+
+! START
+
+y(:)=0d0
+
+do i=1,n_entry
+  row=j_K(i)
+  col=i_K(i)
+    
+  y(row)=y(row)+s_K_re(i)*x(col)
+  
+end do
+
+RETURN
+
+END SUBROUTINE ATprod 
+!
+! NAME
+!     SUBROUTINE zAprod
+!
+! DESCRIPTION
+!     Matrix vector multiplication for sparse, real matrices
+!
+! COMMENTS
+!     
+!
+! HISTORY
+!    19/3/2018 CJS 
+!
+!
+SUBROUTINE zAprod ( n, x, y )
+
+! form the vector y=A*x
+
+USE type_specifications
+
+IMPLICIT NONE
+
+integer n
+complex(dp) :: x(n), y(n)
+
+integer :: i,row,col
+
+! START
+
+y(:)=0d0
+
+do i=1,n_entry
+  row=i_K(i)
+  col=j_K(i)
+  
+  y(row)=y(row)+s_K(i)*x(col)
+  
+end do
+
+RETURN
+
+END SUBROUTINE zAprod 
+!     
+! NAME
+!     SUBROUTINE zATprod
+!
+! DESCRIPTION
+!     Transpose Matrix vector multiplication for sparse, real matrices
+!
+! COMMENTS
+!     
+!
+! HISTORY
+!    19/3/2018 CJS 
+!
+!
+SUBROUTINE zATprod ( n, x, y )
+
+! form the vector y=AT*x
+
+USE type_specifications
+
+IMPLICIT NONE
+
+integer n
+complex(dp) :: x(n), y(n)
+
+integer :: i,row,col
+
+! START
+
+y(:)=0d0
+
+do i=1,n_entry
+  row=j_K(i)
+  col=i_K(i)
+    
+  y(row)=y(row)+(s_K(i))*x(col)
+  
+end do
+
+RETURN
+
+END SUBROUTINE zATprod 
@@ -0,0 +1,382 @@
+!
+! 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
+!
+! SUBROUTINE solve_real_symm(n, b, x,tol)
+! SUBROUTINE dot ( n, x, y ,res)
+! SUBROUTINE solve_complex_symm(n, b, x,tol)
+! SUBROUTINE zdot ( n, x, y ,res)
+!
+! NAME
+!    solve_real_symm
+!
+! DESCRIPTION
+!
+! Solve the system Ax=b using a conjugate gradient method. 
+! The iterative solution halts when the norm of the residual is less than tol
+!
+! HISTORY
+!
+!     started 19/6/2018 CJS
+!
+! COMMENTS
+! 
+SUBROUTINE solve_real_symm(n, b, x,tol)
+
+USE type_specifications
+
+IMPLICIT NONE
+
+! variables passed to subroutine
+
+integer,intent(IN) :: n
+
+real(dp),intent(IN) :: b(n)
+real(dp),intent(INOUT) :: x(n)   ! x contains the initial guess for x
+real(dp),intent(IN) :: tol
+
+! local variables
+
+real(dp) :: ak
+real(dp) :: rk(n),rkb(n)
+real(dp) :: pk(n),pkb(n)
+real(dp) :: bk
+
+real(dp) :: t1(n),t2(n)
+real(dp) :: res(n),res_norm
+real(dp) :: num,den
+
+integer :: i,ii,itloop
+
+logical :: minres
+
+! START
+
+write(*,*)'CALLED solve_real_symm with n=',n
+
+! choice between normal conjugate gradient and minimum residual forms.
+
+minres=.TRUE.
+
+! set starting values
+    
+! calculate the residual res=b-A.xk
+  CALL Aprod( n, x, t1 )    
+  res(:)=b(:)-t1(:)
+
+! calculate the size of the residual
+  res_norm=0d0
+  do i=1,n
+    res_norm=res_norm+res(i)*res(i)
+  end do
+  res_norm=sqrt(res_norm)
+    
+  write(*,*)'iteration ',0,' norm=',res_norm
+  
+! initial r1=b-A.x
+  CALL Aprod( n, x, t1 )    
+  rk(:)=b(:)-t1(:) 
+  
+  if (minres) then
+! For the minimum residual algorithm r1b=A.r1
+    CALL Aprod( n, rk, rkb )    
+  else
+! For the normal algorithm  r1b=r1
+    rkb(:)=rk(:)    
+  end if
+
+! p1=r1
+  pk=rk
+  
+! p1b=r1b
+  pkb=rkb
+
+! iteration loop
+  do itloop=1,2*n
+  
+! calculate ak=(rkb.rk)/(pkb.A.pk)
+    
+    CALL dot(n,rkb,rk,num)
+
+    CALL Aprod( n, pk, t1 )    
+    CALL dot(n,pkb,t1,den)
+    
+    ak=num/den    
+       
+! calculate the next value of the solution vector, xk+1=xk+ak*pk
+    x(:)=x(:)+ak*pk(:)
+    
+! calculate the residual res=b-A.xk
+    CALL Aprod( n, x, t2 )    
+    res(:)=b(:)-t2(:)
+
+! calculate the size of the residual
+    res_norm=0d0
+    do i=1,n
+      res_norm=res_norm+res(i)*res(i)
+    end do
+    res_norm=sqrt(res_norm)
+    
+    if (res_norm.LT.tol) then
+      write(*,*)'iteration ',itloop,' norm=',res_norm,' ak=',ak
+      RETURN
+    end if
+    
+! calculate rk+1=rk-ak*A.pk
+    rk(:) =rk(:) -ak*t1(:)
+   
+! calculate rk+1b=rkb-ak*AT.pkb, note for the symmetric matrices here AT=A
+    CALL ATprod( n, pkb, t2 )    
+    rkb(:)=rkb(:)-ak*t2(:)
+    
+! calculate bk=(rk+1b.rk+1)/(rkb.rk)
+    den=num
+    CALL dot(n,rkb,rk,num)
+    bk=num/den
+    
+    write(*,*)'iteration ',itloop,' norm=',res_norm,' ak=',ak,' bk=',bk
+        
+! calculate pk+1=rk+1-bk*pk
+    pk(:)=rk(:)+bk*pk(:)
+    
+! calculate pk+1b=rk+1b-bk*pkb
+    pkb(:)=rkb(:)+bk*pkb(:)
+
+  end do ! next iteration
+  
+! finish
+
+  RETURN
+
+
+END SUBROUTINE solve_real_symm
+!
+! _________________________________________________________________________
+!
+!
+
+SUBROUTINE dot ( n, x, y ,res)
+
+! form the dot product of x and y
+
+USE type_specifications
+
+IMPLICIT NONE
+
+integer n
+real(dp) :: x(n), y(n), res
+
+integer :: i
+
+! START
+
+res=0d0
+
+do i=1,n
+
+  res=res+x(i)*y(i)
+  
+end do
+
+RETURN
+
+END SUBROUTINE dot
+!
+! NAME
+!    solve_complex_symm
+!
+! DESCRIPTION
+!
+! Solve the system Ax=b using a conjugate gradient method. 
+! The iterative solution halts when the norm of the residual is less than tol
+!
+! HISTORY
+!
+!     started 19/6/2018 CJS
+!
+! COMMENTS
+! 
+ SUBROUTINE solve_complex_symm(n, b, x,tol)
+
+USE type_specifications
+
+IMPLICIT NONE
+
+! variables passed to subroutine
+
+integer,intent(IN) :: n
+
+complex(dp),intent(IN) :: b(n)
+complex(dp),intent(INOUT) :: x(n)   ! x contains the initial guess for x
+real(dp),intent(IN) :: tol
+
+! local variables
+
+complex(dp) :: ak
+complex(dp) :: rk(n),rkb(n)
+complex(dp) :: pk(n),pkb(n)
+complex(dp) :: bk
+
+complex(dp) :: t1(n),t2(n)
+complex(dp) :: res(n)
+real(dp)    :: res_norm
+complex(dp) :: num,den
+
+integer :: i,ii,itloop
+
+logical :: minres
+
+! START
+
+write(*,*)'CALLED solve_real_symm with n=',n
+
+! choice between normal conjugate gradient and minimum residual forms.
+
+minres=.TRUE.
+
+! set starting values
+    
+! calculate the residual res=b-A.xk
+  CALL zAprod( n, x, t1 )    
+  res(:)=b(:)-t1(:)
+
+! calculate the size of the residual
+  res_norm=0d0
+  do i=1,n
+    res_norm=res_norm+res(i)*res(i)
+  end do
+  res_norm=sqrt(res_norm)
+    
+  write(*,*)'iteration ',0,' norm=',res_norm
+  
+! initial r1=b-A.x
+  CALL zAprod( n, x, t1 )    
+  rk(:)=b(:)-t1(:) 
+  
+  if (minres) then
+! For the minimum residual algorithm r1b=A.r1
+    CALL zAprod( n, rk, rkb )    
+  else
+! For the normal algorithm  r1b=r1
+    rkb(:)=rk(:)    
+  end if
+
+! p1=r1
+  pk=rk
+  
+! p1b=r1b
+  pkb=rkb
+
+! iteration loop
+  do itloop=1,2*n
+  
+! calculate ak=(rkb.rk)/(pkb.A.pk)
+    
+    CALL zdot(n,rkb,rk,num)
+
+    CALL zAprod( n, pk, t1 )    
+    CALL zdot(n,pkb,t1,den)
+    
+    ak=num/den    
+       
+! calculate the next value of the solution vector, xk+1=xk+ak*pk
+    x(:)=x(:)+ak*pk(:)
+    
+! calculate the residual res=b-A.xk
+    CALL zAprod( n, x, t2 )    
+    res(:)=b(:)-t2(:)
+
+! calculate the size of the residual
+    res_norm=0d0
+    do i=1,n
+      res_norm=res_norm+abs(res(i))**2
+    end do
+    res_norm=sqrt(res_norm)
+    
+    if (res_norm.LT.tol) then
+      write(*,*)'iteration ',itloop,' norm=',res_norm,' ak=',ak
+      RETURN
+    end if
+    
+! calculate rk+1=rk-ak*A.pk
+    rk(:) =rk(:) -ak*t1(:)
+   
+! calculate rk+1b=rkb-ak*AT.pkb, note for the symmetric matrices here AT=A
+    CALL zATprod( n, pkb, t2 )    
+    rkb(:)=rkb(:)-ak*t2(:)
+    
+! calculate bk=(rk+1b.rk+1)/(rkb.rk)
+    den=num
+    CALL zdot(n,rkb,rk,num)
+    bk=num/den
+    
+    write(*,*)'iteration ',itloop,' norm=',res_norm,' ak=',ak,' bk=',bk
+        
+! calculate pk+1=rk+1-bk*pk
+    pk(:)=rk(:)+bk*pk(:)
+    
+! calculate pk+1b=rk+1b-bk*pkb
+    pkb(:)=rkb(:)+bk*pkb(:)
+
+  end do ! next iteration
+  
+! finish
+
+  RETURN
+
+
+END SUBROUTINE solve_complex_symm
+!
+! _________________________________________________________________________
+!
+!
+
+SUBROUTINE zdot ( n, x, y ,res)
+
+! form the dot product of x and y
+
+USE type_specifications
+
+IMPLICIT NONE
+
+integer n
+complex(dp) :: x(n), y(n), res
+
+integer :: i
+
+! START
+
+res=0d0
+
+do i=1,n
+
+  res=res+x(i)*y(i)     ! note no complex conjugate in the inner product here. This gives the biconjugate gradient method.
+  
+end do
+
+RETURN
+
+END SUBROUTINE zdot
@@ -106,6 +106,7 @@
 !    started 5/7/2016 CJS. This subroutineis based on software from NLR and has been 
 !                          translated from Matlab to Fortran.
 !     8/5/2017         CJS: Include references to Theory_Manual
+!     19/6/2018         CJS: Include iterative solver based on the biconjugate gradient method
 !
 SUBROUTINE Laplace(mesh_filename,dim,first_surface_is_free_space_boundary,   &
                    n_dielectric_regions,dielectric_region_epsr,frequency,Lmat,Cmat,Gmat,ox,oy)
@@ -169,7 +170,7 @@ integer :: N_nodes_in                                     ! number of nodes in t
 integer :: N_elements_in                                  ! number of elements in the gmsh file (this is not necessarily the number of triangular elements in the FE mesh)
 integer :: n_boundaries                                   ! number of boundaries, not including dielectric (internal) boundaries
-integer :: N_boundaries_max                               ! maximum boundary number generated by PUL_LC_Laplace and found in the mesh file
+integer :: N_boundaries_max                               ! maximum boundary number generated by PUL_LC_Laplace and found insolve_real_symm(n, b, x,tol) the mesh file
 integer :: N_boundary                                     ! number of viable boundaries in the mesh i.e. boundaries with two or more points on
 integer :: boundary_number
 integer :: N_elements_boundary_temp
@@ -202,7 +203,7 @@ integer     :: N_unknown                     ! number of known node voltages i.e
 integer     :: N_known                       ! number of unknown node voltages
 integer     :: jmax                          ! maximum boundary number (i.e. number of conductors)
-                                             ! this is used to determine the number of finite element solutions (right hand sides to solve for) to
+                                             ! this is used to determine the number of finite element solutions (right hand solve_real_symm(n, b, x,tol)sides to solve for) to
                                              ! fill the capacitance/conductance matrix
 integer     :: total_n_boundary_nodes        ! total number of boundary nodes i.e. the number of known node values
@@ -216,10 +217,12 @@ complex(dp),allocatable     :: c(:,:)      ! element based array of constants re
 complex(dp),allocatable     :: delta(:)    ! element based array with a value related to the element geometry
 complex(dp),allocatable     :: eps_r(:)    ! element based array with the complex relative permittivity of the element
-! 1D arrays used in the construction of the K matrix ( K(i_K(:),j_K(:))=K(i_K(:),j_K(:))+s_K(:) )
-integer,allocatable         :: i_K(:)
-integer,allocatable         :: j_K(:)
-complex(dp),allocatable     :: s_K(:)
+!
+!integer     :: n_entrysolve_real_symm(n, b, x,tol)
+!! 1D arrays used in the construction of the K matrix ( K(i_K(:),j_K(:))=K(i_K(:),j_K(:))+s_K(:) )
+!integer,allocatable         :: i_K(:)
+!integer,allocatable         :: j_K(:)
+!complex(dp),allocatable     :: s_K(:)
 ! 1D arrays used in the construction of the right hand side vectors ( K_rhs(i_K_rhs(:),j_K_rhs(:))=K_rhs(i_K_rhs(:),j_K_rhs(:))+s_K_rhs(:) )
 integer,allocatable         :: i_K_rhs(:)        ! row number
@@ -248,7 +251,6 @@ real(dp),allocatable    :: Conductance_energy(:,:)
 integer :: node,new_node
 integer :: element
-integer     :: n_entry
 integer     :: n_entry_K_rhs
 integer     :: n1,n2,n3
 complex(dp) :: x1,x2,x3,y1,y2,y3
@@ -280,6 +282,37 @@ integer :: nr
 integer      :: ierr   ! error code for matrix inversion 
+! variables for iterative solver
+
+logical  checkA 
+integer  itnlim 
+real(dp)   rtol   
+integer  nout
+logical  goodb
+logical  precon 
+real(dp)   shift
+
+integer istop, itn
+real(dp)  anorm, acond, rnorm, ynorm
+      
+real(dp),allocatable :: r1(:)
+real(dp),allocatable :: r2(:)
+real(dp),allocatable :: vt(:)
+real(dp),allocatable :: wt(:)
+real(dp),allocatable :: yt(:)
+real(dp),allocatable :: bt(:)
+real(dp),allocatable :: xt(:)
+
+logical :: lossy_dielectric                  ! flag to indicate lossy dielectric i.e. we must solve a complex problem
+
+logical :: wantse
+integer :: n,m
+real(dp) :: atol, btol, conlim, damp
+real(dp) :: Arnorm, xnorm
+real(dp),allocatable :: se(:)
+
+real(dp) :: Vout
+
 ! START
 if (verbose) write(*,*)'CALLED: Laplace'
@@ -542,12 +575,16 @@ if (verbose) write(*,*)&#39; Number of materials=&#39;,N_materials
 ALLOCATE ( Mat_prop(1:N_materials,1:4) )
+lossy_dielectric=.FALSE.
+
 do i=1,N_materials
   Mat_prop(i,1)=i    ! material number
   Mat_prop(i,2)=dble(dielectric_region_epsr(i-1))   ! Re{epsr}
   Mat_prop(i,3)=aimag(dielectric_region_epsr(i-1))  ! Im{epsr}
   Mat_prop(i,4)=frequency
+  if ( abs(Mat_prop(i,3)).GT.small ) lossy_dielectric=.TRUE.
+    
   if (verbose) then
     count=0
     do element=1,N_elements
@@ -942,7 +979,16 @@ do loop=1,2
               ls=sqrt((x2-x1)**2+(y2-y1)**2)      ! l2= boundary element edge length
-              gamma=-eps_r(el)/(log(rho)*rho)     ! gamma (see equation 4.3 with equation 4.93)
+              if (use_ABC) then
+! This implements an asymptotic boundary condition. This apears to cause convergence
+! issues for the Capacitance matrix calculation so should be used with caution.
+! The default is the Neumann boundary condition (gamma=0)
+                gamma=-eps_r(el)/(log(rho)*rho)     ! gamma (see equation 4.3 with equation 4.93)
+              else
+! Neumann boundary condition on the outer boundary. This effectively sets the charge to zero on the outer
+! boundary and the convergence of the capacitance matrix is improved compared with the ABC
+                gamma=0d0   
+              end if
               ! K11 contributions according to equation 4.51 (note delta_ij=1 if i=j, 0 otherwise)
               n_entry=n_entry+1
@@ -1002,8 +1048,6 @@ ALLOCATE ( v_tmp(1:N_known) )
 ! solution based on a full matrix inverse 
-ALLOCATE ( K(1:N_unknown,1:N_unknown) )
-ALLOCATE ( KI(1:N_unknown,1:N_unknown) )
 ALLOCATE ( K_rhs(1:N_unknown,1:N_known) )
 if(verbose) then
@@ -1011,12 +1055,6 @@ if(verbose) then
   write(*,*)'Number of entries in K_rhs',n_entry_K_rhs
 end if
-! STAGE 11a: fill the K matrix
-K(1:N_unknown,1:N_unknown)=0d0
-do i=1,n_entry
-  K(i_K(i),j_K(i))=K(i_K(i),j_K(i))+s_K(i)
-end do
-
 ! STAGE 11b: fill the K_rhs matrix
 K_rhs(1:N_unknown,1:N_known)=0d0
 do i=1,n_entry_K_rhs
@@ -1030,20 +1068,108 @@ else
   write(*,*)'Dimension of K in Laplace is',N_unknown,N_unknown
 end if
+if (direct_solver) then
+  ALLOCATE ( K(1:N_unknown,1:N_unknown) )
+  ALLOCATE ( KI(1:N_unknown,1:N_unknown) )
+
+! STAGE 11a: fill the K matrix
+  K(1:N_unknown,1:N_unknown)=0d0
+  do i=1,n_entry
+    K(i_K(i),j_K(i))=K(i_K(i),j_K(i))+s_K(i)
+  end do
+
 ! STAGE 11c: Invert the K matrix
-if(verbose) write(*,*)'Invert the K matrix'
-ierr=0   ! set ierr=0 to cause an error within cinvert_Gauss_Jordan if there is a problem calculating the inverse
-CALL cinvert_Gauss_Jordan(K,N_unknown,KI,N_unknown,ierr) 
+  if(verbose) write(*,*)'Invert the K matrix'
+  ierr=0   ! set ierr=0 to cause an error within cinvert_Gauss_Jordan if there is a problem calculating the inverse
+  CALL cinvert_Gauss_Jordan(K,N_unknown,KI,N_unknown,ierr) 
 ! STAGE 11d loop over all the RHS vectors solving the matrix equation 
-do j1=1,jmax-1
-    do j2=j1,jmax-1
+  do j1=1,jmax-1
+      do j2=j1,jmax-1
+          v_tmp(1:n_known)=V(j1,j2,1:n_known)
+          b_tmp(1:N_unknown)=-matmul(K_rhs,v_tmp)
+          x_tmp=matmul(KI,b_tmp)
+          x(j1,j2,1:N_unknown)=x_tmp(1:N_unknown)
+      end do
+  end do
+  
+else if(.NOT.lossy_dielectric) then
+    
+  checkA = .true.
+  itnlim = N_unknown * 2
+  rtol   = 1.0D-12
+  nout=6
+  goodb=.FALSE.
+  precon = .FALSE.
+  shift=0d0
+     
+  allocate(r1(1:N_unknown))
+  allocate(r2(1:N_unknown))
+  allocate(vt(1:N_unknown))
+  allocate(wt(1:N_unknown))
+  allocate(yt(1:N_unknown))
+  allocate(xt(1:N_unknown))
+  allocate(bt(1:N_unknown))
+  
+  ALLOCATE( s_K_re(1:n_entry) )
+  s_K_re(1:n_entry)=dble(s_K(1:n_entry))
+  
+! Iterative solution
+   do j1=1,jmax-1
+      do j2=j1,jmax-1
+      
         v_tmp(1:n_known)=V(j1,j2,1:n_known)
         b_tmp(1:N_unknown)=-matmul(K_rhs,v_tmp)
-        x_tmp=matmul(KI,b_tmp)
+        bt(1:N_unknown)=b_tmp(1:N_unknown)
+      
+! UoN conjugate gradient solution
+        CALL solve_real_symm(N_unknown, bt, xt, rtol)
+
+        x(j1,j2,1:N_unknown)=xt(1:N_unknown)
+     
+      end do
+  end do
+      
+  deallocate(r1)
+  deallocate(r2)
+  deallocate(vt)
+  deallocate(wt)
+  deallocate(yt)
+  deallocate(xt)
+  deallocate(bt)
+  
+else if (lossy_dielectric) then
+
+  itnlim=4*N_unknown
+  nout=6
+  wantse=.FALSE.
+  atol=1D-8
+  btol=1d-8
+  conlim=1d10
+  damp=0d0
+  allocate(se(1:N_unknown))
+        
+! Iterative solution
+   do j1=1,jmax-1
+      do j2=j1,jmax-1
+      
+        v_tmp(1:n_known)=V(j1,j2,1:n_known)
+        b_tmp(1:N_unknown)=-matmul(K_rhs,v_tmp)
+        n=N_unknown
+        m=N_unknown
+        
+! UoN conjugate gradient solution
+        rtol   = 1.0D-12
+        CALL solve_complex_symm(N_unknown, b_tmp, x_tmp, rtol)
+ 
         x(j1,j2,1:N_unknown)=x_tmp(1:N_unknown)
-    end do
-end do
+          
+      end do
+  end do
+  
+  deallocate(se)
+       
+end if
 ! STAGE 12 Determine the voltage phi in each node of the mesh
@@ -1321,6 +1447,7 @@ DEALLOCATE( eps_r )
 DEALLOCATE( i_K )
 DEALLOCATE( j_K )
 DEALLOCATE( s_K )
+if ( ALLOCATED(s_K_re) ) DEALLOCATE( s_K_re )
 DEALLOCATE( i_K_rhs )
 DEALLOCATE( j_K_rhs )
@@ -1330,8 +1457,8 @@ DEALLOCATE ( x )
 DEALLOCATE ( x_tmp )
 DEALLOCATE ( b_tmp )
 DEALLOCATE ( v_tmp )
-DEALLOCATE ( K )
-DEALLOCATE ( KI )
+if ( ALLOCATED(K) ) DEALLOCATE ( K )
+if ( ALLOCATED(KI) ) DEALLOCATE ( KI )
 DEALLOCATE ( K_rhs )
 DEALLOCATE( phi )
 default:  $(PUL_PARAMETER_CALCULATION_OBJS)
 # 
 $(OBJ_MOD_DIR)/%.o: %.F90 $(TYPE_SPEC_MODULE) $(CONSTANTS_MODULE) $(GENERAL_MODULE)  \
-                          PUL_analytic.F90 PUL_LC_Laplace.F90 Laplace.F90
+                          PUL_analytic.F90 PUL_LC_Laplace.F90 Laplace.F90 Aprod.F90 CG_solve.F90
 	$(FC) $(FLAGS) -c -o $@ $<
@@ -57,8 +57,9 @@
 !    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)
 !
@@ -89,6 +90,8 @@ IMPLICIT NONE
 ! flags to indicate the configuration
   logical :: ground_plane
+  logical :: finite_ground_plane
+  logical :: infinite_ground_plane
   logical :: overshield
   logical :: open_boundary
@@ -115,16 +118,20 @@ IMPLICIT NONE
 ! 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 :: ro(:)          ! offset in the x direction of this shape from the centre (x():),y(:) above)
+  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
@@ -151,6 +158,7 @@ IMPLICIT NONE
   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.
@@ -171,6 +179,8 @@ IMPLICIT NONE
   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
@@ -181,9 +191,11 @@ IMPLICIT NONE
 ! 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
@@ -200,10 +212,16 @@ IMPLICIT NONE
 ! 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
-    tot_n_boundaries_without_dielectric=PUL%n_conductors+1   ! add a free space boundary
+    end_conductor_to_dielectric_loop=PUL%n_conductors    
     first_conductor_shape=2                    ! shape 1 is for the outer boundary
   else
@@ -221,6 +239,8 @@ IMPLICIT NONE
   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
@@ -304,13 +324,16 @@ IMPLICIT NONE
   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( ro(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) )
@@ -332,17 +355,17 @@ IMPLICIT NONE
     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%o(i))  
-        xmax=max(xmax,PUL%x(i)+PUL%r(i)+PUL%o(i))  
-        ymin=min(ymin,PUL%y(i)-PUL%r(i))
-        ymax=max(ymax,PUL%y(i)+PUL%r(i))
+        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%o(i)) 
-        xmax=max(xmax,PUL%x(i)+PUL%rw(i)+PUL%o(i)) 
-        ymin=min(ymin,PUL%y(i)-PUL%rh(i))
-        ymax=max(ymax,PUL%y(i)+PUL%rh(i))
+        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
@@ -352,13 +375,17 @@ IMPLICIT NONE
         ymin=min(ymin,PUL%y(i)-PUL%rh(i))
         ymax=max(ymax,PUL%y(i)+PUL%rh(i))
-      else
+! 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  
+! 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)
@@ -382,21 +409,48 @@ IMPLICIT NONE
       write(*,*)'boundary radius=',rboundary
     end if
-! set the first shape information to relate to the outer boundary, this is a circle, centred on the bundle centre
+! 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.
-    x(1)=xc
-    y(1)=yc
-    r(1)=rboundary
-    rw(1)=0d0
-    rw2(1)=0d0
-    rh(1)=0d0
-    ro(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
+    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
@@ -418,7 +472,8 @@ IMPLICIT NONE
     rw(shape)=PUL%rw(i)
     rw2(shape)=PUL%rw2(i)
     rh(shape)=PUL%rh(i)
-    ro(shape)=PUL%o(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
@@ -426,6 +481,8 @@ IMPLICIT NONE
       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
@@ -434,22 +491,46 @@ IMPLICIT NONE
   if (ground_plane) then ! Add the ground plane information
     shape=PUL%n_conductors+1   ! ground plane
-    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
-    ro(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   
+    
+    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 
@@ -461,7 +542,8 @@ IMPLICIT NONE
     rw(shape)=PUL%overshield_w
     rw2(shape)=PUL%overshield_w2
     rh(shape)=PUL%overshield_h
-    ro(shape)=0d0   
+    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
@@ -489,7 +571,8 @@ IMPLICIT NONE
       rw(shape)=PUL%rdw(i)          
       rw2(shape)=PUL%rw2(i)         
       rh(shape)=PUL%rdh(i)          
-      ro(shape)=PUL%rdo(i)          
+      rox(shape)=PUL%rdox(i)          
+      roy(shape)=PUL%rdoy(i)          
       rtheta(shape)=PUL%rtheta(i)   
       if (shape_type(shape).EQ.circle) then                 
@@ -525,8 +608,8 @@ IMPLICIT NONE
 ! 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=ro(conductor_shape)
-    yp=0d0
+    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))
@@ -537,7 +620,8 @@ IMPLICIT NONE
       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),ro(shape),rtheta(shape))) then
+      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
@@ -581,56 +665,121 @@ IMPLICIT NONE
       write(gmsh_geometry_file_unit,*)' // circle ',i
       point=point+1
-      write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',x(i)-xc+ro(i),','     ,y(i)-yc,',',0.0,','     ,dl(i),' };'
+      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+ro(i),','     ,y(i)-yc-r(i),',',0.0,',',dl(i),' };'
+      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+ro(i)+r(i),',',y(i)-yc,',',0.0,','     ,dl(i),' };'
+      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+ro(i),','     ,y(i)-yc+r(i),',',0.0,',',dl(i),' };'
+      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+ro(i)-r(i),',',y(i)-yc,',',0.0,','     ,dl(i),' };'
+      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
-      point=point+1
-      
-      xt=rw(i)/2d0+ro(i)
-      yt=rh(i)/2d0
+
+! 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))
-      write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl(i),' };'
-      point=point+1
-      
-      xt=-rw(i)/2d0+ro(i)
-      yt=rh(i)/2d0
+      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))
-      write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl(i),' };'
-      point=point+1
-      
-      xt=-rw(i)/2d0+ro(i)
-      yt=-rh(i)/2d0
+      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))
-      write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl(i),' };'
-      point=point+1
-      
-      xt=rw(i)/2d0+ro(i)
-      yt=-rh(i)/2d0
+      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))
-      write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl(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),ro(i),rtheta(i),dl(i),point,i,gmsh_geometry_file_unit)
+      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
-    end if
+      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,*)' '
@@ -666,6 +815,10 @@ IMPLICIT NONE
       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
@@ -675,9 +828,17 @@ IMPLICIT NONE
       line=line+1
       write(gmsh_geometry_file_unit,*)'Line(',line,' ) = {',p3,',',p4,' };'
       line=line+1
-      write(gmsh_geometry_file_unit,*)'Line(',line,' ) = {',p4,',',p1,' };'
+      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+4
+      point=point+8
     else if (shape_type(i).EQ.Dshape) then
@@ -713,7 +874,36 @@ IMPLICIT NONE
       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
@@ -728,10 +918,11 @@ IMPLICIT NONE
   line_loop=0      ! initialise line_loop counter 
   line=0           ! reset the line counter to the first line
   do i=1,n_shapes
-    
-    line_loop=line_loop+1  
+          
+    if ( (shape_type(i).EQ.circle) ) then
-    if ( (shape_type(i).EQ.circle).OR.(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,' };'  
@@ -739,11 +930,46 @@ IMPLICIT NONE
       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,' };'  
@@ -753,6 +979,18 @@ IMPLICIT NONE
       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
@@ -766,6 +1004,7 @@ IMPLICIT NONE
   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
@@ -785,32 +1024,39 @@ IMPLICIT NONE
 ! 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,i,string2)
+          CALL add_integer_to_string(string1,line_loop,string2)
           string1=trim(string2)//', '        
         end if    
       end if ! overshield
-
-      do i=1,n_shapes
+      
+      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         
-            CALL add_integer_to_string(string1,i,string2)
-            string1=trim(string2)//', '        
+              (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
+      end do  ! next shape      
-    else
+    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         
-          CALL add_integer_to_string(string1,i,string2)
-          string1=trim(string2)//', '        
+            (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
@@ -834,13 +1080,15 @@ IMPLICIT NONE
 ! 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( ro ))  DEALLOCATE( ro )
+  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 )
@@ -40,6 +40,14 @@
 !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
@@ -47,13 +55,18 @@
 !     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
@@ -74,14 +87,16 @@ TYPE::PUL_type
   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        :: o(:)            ! offset in the x direction of this conductor from the cable centre (x():),y(:) above) 
+  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        :: rdo(:)          ! offset of dielectric in the x direction of this conductor from the cable centre
+  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
@@ -112,6 +127,15 @@ TYPE::PUL_type
 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"
@@ -120,6 +144,10 @@ include &quot;PUL_LC_Laplace.F90&quot;
 include "Laplace.F90"
+include "Aprod.F90"
+
+include "CG_solve.F90"
+
 ! NAME
 !     SUBROUTINE allocate_and_reset_PUL_data
 !
@@ -156,11 +184,13 @@ include &quot;Laplace.F90&quot;
   ALLOCATE( PUL%x( PUL%n_conductors) ) 
   ALLOCATE( PUL%y( PUL%n_conductors) ) 
   ALLOCATE( PUL%r( PUL%n_conductors) ) 
-  ALLOCATE( PUL%o( 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%rdo( 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) ) 
@@ -174,14 +204,16 @@ include &quot;Laplace.F90&quot;
     PUL%x(i)=0d0 
     PUL%y(i)=0d0 
     PUL%r(i)=0d0 
-    PUL%o(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%rdo(i)=0d0 
+    PUL%rdox(i)=0d0 
+    PUL%rdoy(i)=0d0 
     PUL%rtheta(i)=0d0
     PUL%epsr(i)=1d0
@@ -240,14 +272,16 @@ include &quot;Laplace.F90&quot;
   if (allocated(PUL%x )    ) DEALLOCATE( PUL%x )
   if (allocated(PUL%y )    ) DEALLOCATE( PUL%y )
   if (allocated(PUL%r )    ) DEALLOCATE( PUL%r )
-  if (allocated(PUL%o )    ) DEALLOCATE( PUL%o )
+  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%rdo )   ) DEALLOCATE( PUL%rdo )
+  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
@@ -281,9 +315,10 @@ include &quot;Laplace.F90&quot;
 !
 ! 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,ro,rtheta)
+  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
@@ -293,7 +328,8 @@ include &quot;Laplace.F90&quot;
   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) ::   ro           ! offset in the x direction of this test object from the centre
+  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
@@ -335,8 +371,8 @@ include &quot;Laplace.F90&quot;
     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-ro
-    yt_ro=yt_r  
+    xt_ro=xt_r-rox
+    yt_ro=yt_r-roy
     write(*,*)'Offset Inverse rotation on test point',xt_ro,yt_ro
@@ -345,9 +381,9 @@ include &quot;Laplace.F90&quot;
     dx=xt_ro
     dy=yt_ro
-    write(*,*)'dielectric centre point',x,y,' offset',ro,' angle',rtheta
-    write(*,*)'dx',dx,' rw/2',rw/2d0
-    write(*,*)'dy',dy,' rh/2',rh/2d0
+!    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.
@@ -442,7 +478,7 @@ include &quot;Laplace.F90&quot;
 !    started 15/11/2016 CJS
 !
-  SUBROUTINE write_Dshape_gmsh(x,y,w_in,w2_in,h_in,r,ox,theta,dl,point,number,unit)
+  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
@@ -458,6 +494,7 @@ include &quot;Laplace.F90&quot;
   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
@@ -498,7 +535,7 @@ include &quot;Laplace.F90&quot;
 ! POINT 1      ! top right
       xt=w1+ox
-      yt=h+r
+      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,' };'
@@ -506,7 +543,7 @@ include &quot;Laplace.F90&quot;
 ! POINT 2      ! top left
       xt=-w1+ox
-      yt=h+r
+      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,' };'
@@ -514,7 +551,7 @@ include &quot;Laplace.F90&quot;
 ! POINT 3      ! top left centre
       xt=-w1+ox
-      yt=h
+      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,' };'
@@ -522,7 +559,7 @@ include &quot;Laplace.F90&quot;
 ! POINT 4      ! top left edge
       xt=-w1+ox+voxl
-      yt=h+voyl
+      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,' };'
@@ -530,7 +567,7 @@ include &quot;Laplace.F90&quot;
 ! POINT 5      ! bottom left edge
       xt=-w2+ox+voxl
-      yt=-h+voyl
+      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,' };'
@@ -538,7 +575,7 @@ include &quot;Laplace.F90&quot;
 ! POINT 6      ! bottom left centre
       xt=-w2+ox
-      yt=-h
+      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,' };'
@@ -546,7 +583,7 @@ include &quot;Laplace.F90&quot;
 ! POINT 7      ! bottom left
       xt=-w2+ox
-      yt=-(h+r)
+      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,' };'
@@ -554,7 +591,7 @@ include &quot;Laplace.F90&quot;
 ! POINT 8      ! bottom right
       xt=w2+ox
-      yt=-(h+r)
+      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,' };'
@@ -562,7 +599,7 @@ include &quot;Laplace.F90&quot;
 ! POINT 9      ! bottom right centre
       xt=w2+ox
-      yt=-h
+      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,' };'
@@ -570,7 +607,7 @@ include &quot;Laplace.F90&quot;
 ! POINT 10      ! bottom right edge
       xt=w2+ox-voxl
-      yt=-h+voyl
+      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,' };'
@@ -578,7 +615,7 @@ include &quot;Laplace.F90&quot;
 ! POINT 11      ! top right edge
       xt=w1+ox-voxl
-      yt=h+voyl
+      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,' };'
@@ -586,7 +623,7 @@ include &quot;Laplace.F90&quot;
 ! POINT 12      ! top right centre
       xt=w1+ox
-      yt=h
+      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,' };'
@@ -74,6 +74,8 @@
 !     29/6/2016 CJS: Check that the cables are all on one side of the ground plane if it is present
 !     December 2016 CJS Version 2: Rationalise cable types so that there is only a single version of each type of cable
 !     24/2/2017 CJS Allow the input name to include a path i.e. the _spec file does not need to be local.
+!     13/3/2018 CJS Add flag for direct/ iterative matrix solver in Laplace solution and inf/finite ground plane
+!     19/6/2018 CJS Add flag for Neumann/ Asymptotic boundary condition in Laplace solver. Default is Neumann
 !
 PROGRAM cable_bundle_model_builder
@@ -333,6 +335,15 @@ integer        :: dim
     if (INDEX(line,'no_plot_potential').NE.0) plot_potential=.FALSE.
     if (INDEX(line,'plot_mesh').NE.0) plot_mesh=.TRUE.
     if (INDEX(line,'no_plot_mesh').NE.0) plot_mesh=.FALSE.
+    
+    if (INDEX(line,'direct_solver').NE.0) direct_solver=.TRUE.
+    if (INDEX(line,'iterative_solver').NE.0) direct_solver=.FALSE.
+    
+    if (INDEX(line,'inf_gnd').NE.0) inf_gnd=.TRUE.
+    if (INDEX(line,'finite_gnd').NE.0) inf_gnd=.FALSE.
+    
+    if (INDEX(line,'abc').NE.0) use_ABC=.TRUE.
+    if (INDEX(line,'neumann').NE.0) use_ABC=.FALSE.
 ! redefine mesh generation parameters if required
     if (INDEX(line,'laplace_boundary_constant').NE.0) then
@@ -346,6 +357,10 @@ integer        :: dim
     if (INDEX(line,'twisted_pair_equivalent_radius').NE.0) then
       read(bundle_spec_file_unit,*,END=100,ERR=100)Twisted_pair_equivalent_radius
     end if
+    
+    if (INDEX(line,'max_mesh_edge_length').NE.0) then
+      read(bundle_spec_file_unit,*,END=100,ERR=100)max_mesh_edge_length
+    end if
   end do  ! continue until all flags are read - indicated by an end of file.
@@ -84,6 +84,7 @@
 !     24/2/2017 CJS Allow the input name to include a path i.e. the _spec file does not need to be local.
 !     23/10/2017 CJS Check the pole stability of the dielectric and transfer impedance filters.
 !     16/11/2017 CJS Include network synthesis process to replace s-domain transfer functions
+!     13/3/2018 CJS Add flag for direct/ iterative matrix solver in Laplace solution
 !
 !
 PROGRAM cable_model_builder
@@ -214,10 +215,14 @@ integer        :: i
     CALL spacewire_set_parameters(cable_spec)
-  else if (cable_spec%cable_type_string.eq.'flex_cable') then
+  else if (cable_spec%cable_type_string.eq.'original_flex_cable') then
     CALL flex_cable_set_parameters(cable_spec)
+  else if (cable_spec%cable_type_string.eq.'flex_cable') then
+
+    CALL ML_flex_cable_set_parameters(cable_spec)
+
   else if (cable_spec%cable_type_string.eq.'dconnector') then
     CALL Dconnector_set_parameters(cable_spec)
@@ -274,11 +279,20 @@ integer        :: i
     STOP 1
   end if
+  if (cable_spec%n_parameters.NE.0) then
+  
 ! Check that the number of parameters is consistent with this cable type
-  if (np_in.NE.cable_spec%n_parameters) then
-    run_status='ERROR, Wrong number of parameters'
-    CALL write_program_status()
-    STOP 1
+    if (np_in.NE.cable_spec%n_parameters) then
+      run_status='ERROR, Wrong number of parameters'
+      CALL write_program_status()
+      STOP 1
+    end if
+  
+  else
+
+! for a ML_flex cable the number of parameters depends on the cable specification  
+    cable_spec%n_parameters=np_in
+    
   end if
 ! allocate and read parameters
@@ -410,6 +424,9 @@ integer        :: i
     if (INDEX(line,'no_plot_potential').NE.0) plot_potential=.FALSE.
     if (INDEX(line,'plot_mesh').NE.0) plot_mesh=.TRUE.
     if (INDEX(line,'no_plot_mesh').NE.0) plot_mesh=.FALSE.
+    
+    if (INDEX(line,'direct_solver').NE.0) direct_solver=.TRUE.
+    if (INDEX(line,'iterative_solver').NE.0) direct_solver=.FALSE.
 ! redefine mesh generation parameters if required
     if (INDEX(line,'laplace_boundary_constant').NE.0) then
@@ -423,6 +440,10 @@ integer        :: i
     if (INDEX(line,'twisted_pair_equivalent_radius').NE.0) then
       read(cable_spec_file_unit,*,END=100,ERR=100)Twisted_pair_equivalent_radius
     end if
+    
+    if (INDEX(line,'max_mesh_edge_length').NE.0) then
+      read(cable_spec_file_unit,*,END=100,ERR=100)max_mesh_edge_length
+    end if
   end do  ! continue until all flags are read - indicated by an end of file.
@@ -482,6 +503,10 @@ integer        :: i
     CALL flex_cable_set_internal_domain_information(cable_spec)
+  else if (cable_spec%cable_type.EQ.cable_geometry_type_ML_flex_cable) then
+
+    CALL ML_flex_cable_set_internal_domain_information(cable_spec)
+
   else if (cable_spec%cable_type.EQ.cable_geometry_type_dconnector) then
     CALL dconnector_set_internal_domain_information(cable_spec)
@@ -2,13 +2,17 @@
 .
 flex_cable
 3		# number of conductors
-6		# number of parameters
-1.0e-3  	# parameter 1: conductor width (x dimension)
-0.25e-3  	# parameter 2: conductor height (y dimension)
-0.5e-3  	# parameter 3: conductor separation (x dimension)
-0.5e-3  	# parameter 4: dielectric offset x
-0.5e-3  	# parameter 5: dielectric offset y
-0.0  	        # parameter 6: conductivity
+10		# number of parameters
+5.0e-3          # parameter 1: dielectric width (x dimension)
+1.25e-3  	# parameter 2: dielectric height (y dimension)
+1               # parameter 3: number of rows of conductors
+0.0e-3  	# parameter 4: row 1 centre offset x
+0.0e-3  	# parameter 5: row 1 centre offset y
+1.0e-3          # parameter 6: row 1 conductor width (x dimension)
+0.25e-3         # parameter 7: row 1 conductor height (y dimension)
+0.5e-3          # parameter 8: row 1 conductor separation
+3               # parameter 9: row 1 number of conductors
+0.0  	        # parameter 10: conductivity
 1		# number of frequency dependent parameters
 # dielectric relative permittivity model follows
    1E9        # w normalisation constant
@@ -2,6 +2,33 @@
 .
 flex_cable
 3		# number of conductors
+10		# number of parameters
+4.4e-3          # parameter 1: dielectric width (x dimension)
+0.45e-3  	# parameter 2: dielectric height (y dimension)
+1               # parameter 3: number of rows of conductors
+0.0e-3  	# parameter 4: row 1 centre offset x
+0.0e-3  	# parameter 5: row 1 centre offset y
+1.0e-3          # parameter 6: row 1 conductor width (x dimension)
+0.25e-3         # parameter 7: row 1 conductor height (y dimension)
+0.5e-3          # parameter 8: row 1 conductor separation
+3               # parameter 9: row 1 number of conductors
+0.0  	        # parameter 10: conductivity
+1		# number of frequency dependent parameters
+# dielectric relative permittivity model follows
+   1E9        # w normalisation constant
+           0   # a order, a coefficients follow below:
+  1.0         
+           0   # b order, b coefficients follow below:
+  1.0      
+
+
+
+# OLD model
+
+#MOD_cable_lib_dir
+.
+flex_cable
+3		# number of conductors
 6		# number of parameters
 1.0e-3  	# parameter 1: conductor width (x dimension)
 0.25e-3  	# parameter 2: conductor height (y dimension)
@@ -16,11 +43,3 @@ flex_cable
   1.0         
            0   # b order, b coefficients follow below:
   1.0      
-
-
-# dielectric relative permittivity model follows
-   1E9        # w normalisation constant
-           1   # a order, a coefficients follow below:
-  4.0      2.0     
-           1   # b order, b coefficients follow below:
-  1.0      1.0
@@ -2,6 +2,33 @@
 .
 flex_cable
 3		# number of conductors
+10		# number of parameters
+5.0e-3          # parameter 1: dielectric width (x dimension)
+1.25e-3  	# parameter 2: dielectric height (y dimension)
+1               # parameter 3: number of rows of conductors
+0.0e-3  	# parameter 4: row 1 centre offset x
+0.0e-3  	# parameter 5: row 1 centre offset y
+1.0e-3          # parameter 6: row 1 conductor width (x dimension)
+0.25e-3         # parameter 7: row 1 conductor height (y dimension)
+0.5e-3          # parameter 8: row 1 conductor separation
+3               # parameter 9: row 1 number of conductors
+0.0  	        # parameter 10: conductivity
+1		# number of frequency dependent parameters
+# dielectric relative permittivity model follows
+   1E9        # w normalisation constant
+           0   # a order, a coefficients follow below:
+  2.2         
+           0   # b order, b coefficients follow below:
+  1.0      
+
+
+
+
+
+#MOD_cable_lib_dir
+.
+flex_cable
+3		# number of conductors
 6		# number of parameters
 1.0e-3  	# parameter 1: conductor width (x dimension)
 0.25e-3  	# parameter 2: conductor height (y dimension)
@@ -2,6 +2,32 @@
 .
 flex_cable
 3		# number of conductors
+10		# number of parameters
+5.0e-3          # parameter 1: dielectric width (x dimension)
+1.25e-3  	# parameter 2: dielectric height (y dimension)
+1               # parameter 3: number of rows of conductors
+0.0e-3  	# parameter 4: row 1 centre offset x
+0.0e-3  	# parameter 5: row 1 centre offset y
+1.0e-3          # parameter 6: row 1 conductor width (x dimension)
+0.25e-3         # parameter 7: row 1 conductor height (y dimension)
+0.5e-3          # parameter 8: row 1 conductor separation
+3               # parameter 9: row 1 number of conductors
+0.0  	        # parameter 10: conductivity
+1		# number of frequency dependent parameters
+# dielectric relative permittivity model follows
+   1E9        # w normalisation constant
+           0   # a order, a coefficients follow below:
+  2.2         
+           0   # b order, b coefficients follow below:
+  1.0      
+
+
+
+
+#MOD_cable_lib_dir
+.
+flex_cable
+3		# number of conductors
 6		# number of parameters
 1.0e-3  	# parameter 1: conductor width (x dimension)
 0.25e-3  	# parameter 2: conductor height (y dimension)
@@ -2,6 +2,32 @@
 .
 flex_cable
 3		# number of conductors
+10		# number of parameters
+5.0e-3          # parameter 1: dielectric width (x dimension)
+1.25e-3  	# parameter 2: dielectric height (y dimension)
+1               # parameter 3: number of rows of conductors
+0.0e-3  	# parameter 4: row 1 centre offset x
+0.0e-3  	# parameter 5: row 1 centre offset y
+1.0e-3          # parameter 6: row 1 conductor width (x dimension)
+0.25e-3         # parameter 7: row 1 conductor height (y dimension)
+0.5e-3          # parameter 8: row 1 conductor separation
+3               # parameter 9: row 1 number of conductors
+0.0  	        # parameter 10: conductivity
+1		# number of frequency dependent parameters
+# dielectric relative permittivity model follows
+   1E9        # w normalisation constant
+           0   # a order, a coefficients follow below:
+  2.2         
+           0   # b order, b coefficients follow below:
+  1.0      
+
+
+
+
+#MOD_cable_lib_dir
+.
+flex_cable
+3		# number of conductors
 6		# number of parameters
 1.0e-3  	# parameter 1: conductor width (x dimension)
 0.25e-3  	# parameter 2: conductor height (y dimension)
@@ -2,6 +2,32 @@
 .
 flex_cable
 3		# number of conductors
+10		# number of parameters
+5.0e-3          # parameter 1: dielectric width (x dimension)
+1.25e-3  	# parameter 2: dielectric height (y dimension)
+1               # parameter 3: number of rows of conductors
+0.0e-3  	# parameter 4: row 1 centre offset x
+0.0e-3  	# parameter 5: row 1 centre offset y
+1.0e-3          # parameter 6: row 1 conductor width (x dimension)
+0.25e-3         # parameter 7: row 1 conductor height (y dimension)
+0.5e-3          # parameter 8: row 1 conductor separation
+3               # parameter 9: row 1 number of conductors
+0.0  	        # parameter 10: conductivity
+1		# number of frequency dependent parameters
+# dielectric relative permittivity model follows
+   1E9        # w normalisation constant
+           0   # a order, a coefficients follow below:
+  2.2         
+           0   # b order, b coefficients follow below:
+  1.0      
+
+
+
+
+#MOD_cable_lib_dir
+.
+flex_cable
+3		# number of conductors
 6		# number of parameters
 1.0e-3  	# parameter 1: conductor width (x dimension)
 0.25e-3  	# parameter 2: conductor height (y dimension)
@@ -16,3 +42,29 @@ flex_cable
   2.2         
            0   # b order, b coefficients follow below:
   1.0      
+#MOD_cable_lib_dir
+.
+ML_flex_cable
+3		# number of conductors
+10		# number of parameters
+5.0e-3          # parameter 1: dielectric width (x dimension)
+1.25e-3  	# parameter 2: dielectric height (y dimension)
+1               # parameter 3: number of rows of conductors
+0.0e-3  	# parameter 4: row 1 centre offset x
+0.0e-3  	# parameter 5: row 1 centre offset y
+1.0e-3          # parameter 6: row 1 conductor width (x dimension)
+0.25e-3         # parameter 7: row 1 conductor height (y dimension)
+0.5e-3          # parameter 8: row 1 conductor separation
+3               # parameter 9: row 1 number of conductors
+0.0  	        # parameter 10: conductivity
+1		# number of frequency dependent parameters
+# dielectric relative permittivity model follows
+   1E9        # w normalisation constant
+           0   # a order, a coefficients follow below:
+  2.2         
+           0   # b order, b coefficients follow below:
+  1.0      
+
+
+
+
@@ -2,6 +2,32 @@
 .
 flex_cable
 2		# number of conductors
+10		# number of parameters
+2.9e-3          # parameter 1: dielectric width (x dimension)
+0.45e-3  	# parameter 2: dielectric height (y dimension)
+1               # parameter 3: number of rows of conductors
+0.0e-3  	# parameter 4: row 1 centre offset x
+0.0e-3  	# parameter 5: row 1 centre offset y
+1.0e-3          # parameter 6: row 1 conductor width (x dimension)
+0.25e-3         # parameter 7: row 1 conductor height (y dimension)
+0.5e-3          # parameter 8: row 1 conductor separation
+2               # parameter 9: row 1 number of conductors
+5e7  	        # parameter 10: conductivity
+1		# number of frequency dependent parameters
+# dielectric relative permittivity model follows
+   1E9        # w normalisation constant
+           0   # a order, a coefficients follow below:
+  2.25         
+           0   # b order, b coefficients follow below:
+  1.0      
+
+
+
+
+#MOD_cable_lib_dir
+.
+flex_cable
+2		# number of conductors
 6		# number of parameters
 1.0e-3  	# parameter 1: conductor width (x dimension)
 0.25e-3  	# parameter 2: conductor height (y dimension)
@@ -163,7 +163,6 @@ CM_TO_DM_TWISTED_PAIR_OVER_GROUND_PLANE
 LOW_MASS_SPACEWIRE
 "
-
 ERROR_STRING="\
 Run using the following command:
@@ -0,0 +1,16 @@
+# Script to clean object files, mod files and 
+# all temporary test case files before updating 
+# the git repository.
+
+cd SRC 
+make clean
+rm -f compilation_date.inc 
+cd ..
+#
+cd TEST_CASES
+./generate_spice_cable_bundle_model clean_all
+cd ..
+#
+cd DOCUMENTATION
+make clean
+cd ..