!
! This file is part of SACAMOS, State of the Art CAble MOdels for Spice. 
! It was developed by the University of Nottingham and the Netherlands Aerospace 
! Centre (NLR) for ESA under contract number 4000112765/14/NL/HK.
! 
! Copyright (C) 2016-2018 University of Nottingham
! 
! SACAMOS is free software: you can redistribute it and/or modify it under the 
! terms of the GNU General Public License as published by the Free Software 
! Foundation, either version 3 of the License, or (at your option) any later 
! version.
! 
! SACAMOS is distributed in the hope that it will be useful, but 
! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 
! or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
! for more details.
! 
! A copy of the GNU General Public License version 3 can be found in the 
! file GNU_GPL_v3 in the root or at <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:
! MODULE cable_bundle_module
!   CONTAINS
!     SUBROUTINE read_cable_bundle
!     SUBROUTINE write_cable_bundle
!     SUBROUTINE deallocate_cable_bundle
!     SUBROUTINE check_cables_wrt_ground_plane
!     SUBROUTINE check_cable_intersections
!
! NAME
!     MODULE cable_bundle_module
!
!     Data strunctures and subroutines relating to cable bundles
!     
! COMMENTS
!     Need to improve error handling in file reading 
!  Not sure that some arrays in the type_specification need to be in this structure as they are only used in one subroutine...    
!
! HISTORY
!    started 25/11/2015 CJS
!    stage 3 developments started 21/3/2015 CJS
!    generalise conductor impedance for transfer impedance model (stage 5) CJS 12/5/2016
!    CJS 25/08/2016   updates for revised transfer impedance and conductor impedance model for shields
!     16/11/2017 CJS Include network synthesis process to replace s-domain transfer functions
!
MODULE cable_bundle_module

USE type_specifications
USE cable_module

IMPLICIT NONE

! This is a general structure used to hold conductor based lists, 

TYPE::conductor_list_type

  integer             :: n_elements
  integer,allocatable :: element(:)

END TYPE conductor_list_type

! This is the main structure used to hold all the bundle information

TYPE::bundle_specification_type 

  character(LEN=line_length)    :: version 
  character(LEN=line_length)    :: bundle_name
  
  integer            :: tot_n_conductors
  integer            :: tot_n_shielded_conductors
  integer            :: tot_n_external_conductors
  
! dimension of the global matrix system (total number of conductors-1)
  integer            :: system_dimension

! cable information, total number of cables including and excluding the ground plane (if it exists)
  integer            :: n_cables
  integer            :: n_cables_without_ground_plane
  
! Array of cable specification data for all the cables in this bundle
  type(cable_specification_type),allocatable    :: cable(:)
  
! list of cable positions in the bundle cross section and rotation angles
  real(dp),allocatable        :: cable_x_offset(:)
  real(dp),allocatable        :: cable_y_offset(:)
  real(dp),allocatable        :: cable_angle(:)
  
! ground plane information: note we have two notations at the moment. 
! This will need to be sorted out and one only adopted

  logical             :: ground_plane_present
  
  real(dp)            :: ground_plane_angle
  real(dp)            :: ground_plane_offset
  
  real(dp)            :: ground_plane_sigma      ! ground plane conductivity for FastHenry2
  real(dp)            :: ground_plane_w          ! ground plane width for FastHenry2
  real(dp)            :: ground_plane_h          ! ground plane height for FastHenry2
  real(dp)            :: ground_plane_Rdc        ! ground plane dc resistance for FastHenry2
  integer             :: ground_plane_nsegx      ! ground plane number of segments in x for FastHenry2
  integer             :: ground_plane_nsegz      ! ground plane number of segments in z for FastHenry2
  integer             :: ground_plane_nh         ! ground plane number of layers in h for FastHenry2
  
  real(dp)            :: ground_plane_x,ground_plane_y,ground_plane_theta
  real(dp)            :: ground_plane_nx,ground_plane_ny
  
  integer             :: ground_plane_cable_side    ! +1 if the cables are on the normal direction side of the ground plane, -1 otherwise

! domain based information  
  integer             :: tot_n_domains

  integer,allocatable :: n_conductors(:)  ! number of conductors in each domain
  
  type(conductor_list_type),allocatable  :: terminal_conductor_list(:) ! external conductor number for coonductors in each domain
  
  type(matrix),allocatable        :: L(:)  ! 'high frequency' inductance matrix for each domain
  type(matrix),allocatable        :: C(:)  ! 'high frequency' capacitance matrix for each domain
  
  type(Sfilter_matrix),allocatable    :: Z(:)  ! frequency dependent impedance matrix for each domain
  type(Sfilter_matrix),allocatable    :: Y(:)  ! frequency dependent admittance matrix for each domain
  
! global conductor based information  
  type(matrix)        :: global_MI     ! global domain current transformation matrix
  type(matrix)        :: global_MV     ! global domain voltage transformation matrix
  
  type(matrix)        :: global_L       ! global 'high frequency' inductance matrix 
  type(matrix)        :: global_C       ! global 'high frequency' capacitance matrix
  
  type(Sfilter_matrix)    :: global_Z    ! global frequency dependent impedance matrix 
  type(Sfilter_matrix)    :: global_Y    ! global frequency dependent admittance matrix

! global conductor impedance model information for each conductor
  type(conductor_impedance_model),allocatable    :: conductor_impedance(:)
  
! global conductor position in the bundle cross-section
  real(dp),allocatable        :: conductor_x_offset(:)
  real(dp),allocatable        :: conductor_y_offset(:)

! arrays for cross referecing numbering information from the global (external) conductor numbering and the domain based numbering systems
  logical,allocatable :: terminal_conductor_is_shield_flag(:)
  integer,allocatable :: terminal_conductor_to_inner_domain(:)
  integer,allocatable :: terminal_conductor_to_outer_domain(:)
  integer,allocatable :: terminal_conductor_to_global_domain_conductor(:)
  integer,allocatable :: terminal_conductor_to_local_domain_conductor(:)
  integer,allocatable :: terminal_conductor_to_reference_terminal_conductor(:)
  
! Y matrix element function fitting information for frequency dependent dielectrics

! model order for filter fitting
  integer    :: Y_fit_model_order
      
! frequency range specification 
  type(frequency_specification) :: Y_fit_freq_spec
  
  character(LEN=line_length),allocatable   :: conductor_label(:)
    
END TYPE bundle_specification_type

CONTAINS

! The following subroutines apply to all bundles and cable types

! NAME
!     SUBROUTINE read_cable_bundle(unit)
!
!     read cable bundle structure from a specified unit
!     
! COMMENTS
!     
!
! HISTORY
!    started 2/12/2015 CJS
!    stage 3 developments started 21/3/2015 CJS
!    generalise conductor impedance for transfer impedance model (stage 5) CJS 12/5/2016
!

  SUBROUTINE read_cable_bundle(bundle,file_unit)

  USE type_specifications
  USE constants
  USE general_module
  USE cable_module
  USE maths

  IMPLICIT NONE

! variables passed to subroutine

  type(bundle_specification_type),intent(INOUT) :: bundle
  integer ,intent(IN)                           :: file_unit
  
! local variables
    
  character(len=filename_length)    :: filename
  logical                :: file_exists
  character(len=line_length)        :: line

  integer                :: cable
  integer                :: n_domains
  integer                :: domain
  integer                :: matrix_dimension
  integer                :: n_conductors
  integer                :: conductor
  
  integer                :: ierr

! START

! Open the (.bundle) file to read

  filename=trim(MOD_bundle_lib_dir)//trim(bundle%bundle_name)//bundle_file_extn

  inquire(file=trim(filename),exist=file_exists)
  if (.NOT.file_exists) then
    run_status='ERROR, Cannot find the required bundle file:'//trim(filename)
    CALL write_program_status()
    STOP 1
  end if 
  open(unit=file_unit,file=trim(filename))

  if (verbose) write(*,*)'Opened file:',trim(filename)

! read .bundle file
  read(file_unit,'(A)',ERR=9000)bundle%version
  read(file_unit,'(A)',ERR=9000)bundle%bundle_name
  
! read cable information
  read(file_unit,*,ERR=9000)bundle%n_cables_without_ground_plane
  read(file_unit,*,ERR=9000)bundle%n_cables
  
! allocate and read cable information and cable positions in bundle 
  ALLOCATE( bundle%cable(1:bundle%n_cables) )
  ALLOCATE( bundle%cable_x_offset(1:bundle%n_cables) )
  ALLOCATE( bundle%cable_y_offset(1:bundle%n_cables) )
  ALLOCATE( bundle%cable_angle(1:bundle%n_cables) )

  do cable=1,bundle%n_cables_without_ground_plane
    read(file_unit,'(A)')bundle%cable(cable)%cable_name
    CALL read_cable(bundle%cable(cable),cable_file_unit)
    read(file_unit,*,ERR=9000)bundle%cable_x_offset(cable),bundle%cable_y_offset(cable),bundle%cable_angle(cable)
! convert angle to radians
    bundle%cable_angle(cable)=bundle%cable_angle(cable)*pi/180d0
  end do ! next cable
  
! read ground plane specification

  read(file_unit,'(A)',IOSTAT=ierr)line
    if (ierr.NE.0) then 
      run_status='ERROR reading ground plane present/ absent information'
      CALL write_program_status()
      STOP 1
    end if 
  CALL convert_to_lower_case(line,line_length)
  
  bundle%ground_plane_present=(line.eq.'ground_plane')
  
  if (bundle%ground_plane_present) then
  
    read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_x,bundle%ground_plane_y,bundle%ground_plane_angle
    read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_nx,bundle%ground_plane_ny
    read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_cable_side
    
    if (ierr.NE.0) then 
      run_status='ERROR reading ground plane position and angle'
      CALL write_program_status()
      STOP 1
    end if 
! convert ground plane angle to radians
    bundle%ground_plane_angle=bundle%ground_plane_angle*pi/180d0
! calculate offset
    bundle%ground_plane_offset=-bundle%ground_plane_x*sin(bundle%ground_plane_angle)+  &
                                bundle%ground_plane_y*cos(bundle%ground_plane_angle)

    read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_sigma
    read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_w  
    read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_h      
    read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_Rdc
    read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_nsegx
    read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_nsegz
    read(file_unit,*,IOSTAT=ierr)bundle%ground_plane_nh
    
    if (ierr.NE.0) then 
      run_status='ERROR reading ground plane sigma, w h and Rdc'
      CALL write_program_status()
      STOP 1
    end if 
                                
! set the geometric data for the last cable in the list i.e. the ground plane

    cable=bundle%n_cables
    bundle%cable_x_offset(cable)=bundle%ground_plane_x
    bundle%cable_y_offset(cable)=bundle%ground_plane_y
    bundle%cable_angle(cable)=bundle%ground_plane_angle-pi/2d0  ! note there are two ground plane spec formats at the moment...
 
  else
! set x, y, angle and offset to 0
    bundle%ground_plane_x=0d0
    bundle%ground_plane_y=0d0
    bundle%ground_plane_nx=0d0
    bundle%ground_plane_ny=0d0
    bundle%ground_plane_cable_side=0
    bundle%ground_plane_angle=0d0
    bundle%ground_plane_offset=0d0
  end if
  
! read system dimension information
  read(file_unit,*)bundle%tot_n_conductors
  read(file_unit,*)bundle%tot_n_external_conductors
  read(file_unit,*)bundle%system_dimension
  
! read domain information
  read(file_unit,*,ERR=9000)bundle%tot_n_domains
     
  ALLOCATE( bundle%n_conductors(1:bundle%tot_n_domains) ) 
  
  ALLOCATE( bundle%L(1:bundle%tot_n_domains) ) 
  ALLOCATE( bundle%C(1:bundle%tot_n_domains) )
  ALLOCATE( bundle%Z(1:bundle%tot_n_domains))
  ALLOCATE( bundle%Y(1:bundle%tot_n_domains))
  
  ALLOCATE( bundle%terminal_conductor_list(1:bundle%tot_n_domains))

  n_domains=bundle%tot_n_domains

  do domain=1,bundle%tot_n_domains
  
    read(file_unit,*,ERR=9000)   ! comment line
    
    read(file_unit,*)bundle%n_conductors(domain)
    
    read(file_unit,*,ERR=9000)  ! comment line
    read(file_unit,*,ERR=9000)bundle%terminal_conductor_list(domain)%n_elements    
    ALLOCATE( bundle%terminal_conductor_list(domain)%element(1:bundle%terminal_conductor_list(domain)%n_elements) )
    do conductor=1,bundle%terminal_conductor_list(domain)%n_elements
      read(file_unit,*,ERR=9000)bundle%terminal_conductor_list(domain)%element(conductor)
    end do
        
    read(file_unit,*,ERR=9000)   ! comment line
    
    read(file_unit,*)matrix_dimension
    bundle%L(domain)%dim=matrix_dimension
    ALLOCATE( bundle%L(domain)%mat(1:matrix_dimension,1:matrix_dimension) ) 
    CALL dread_matrix(bundle%L(domain)%mat,matrix_dimension,matrix_dimension,matrix_dimension,file_unit)

    write(*,*)'L='
    CALL dwrite_matrix(bundle%L(domain)%mat,matrix_dimension,matrix_dimension,matrix_dimension,0) ! ******
    
    read(file_unit,*,ERR=9000)   ! comment line
    read(file_unit,*)matrix_dimension
    bundle%C(domain)%dim=matrix_dimension
    ALLOCATE( bundle%C(domain)%mat(1:matrix_dimension,1:matrix_dimension) ) 
    CALL dread_matrix(bundle%C(domain)%mat,matrix_dimension,matrix_dimension,matrix_dimension,file_unit)

    write(*,*)'C='
    CALL dwrite_matrix(bundle%C(domain)%mat,matrix_dimension,matrix_dimension,matrix_dimension,0) ! ******

    read(file_unit,*,ERR=9000) ! comment line 
    CALL read_Sfilter_matrix( bundle%Z(domain),file_unit )
      
    read(file_unit,*,ERR=9000) ! comment line 
    CALL read_Sfilter_matrix( bundle%Y(domain),file_unit )

  end do
  
  read(file_unit,*,ERR=9000)    ! comment line
  read(file_unit,*)matrix_dimension
  bundle%global_MI%dim=matrix_dimension
  ALLOCATE( bundle%global_MI%mat(1:matrix_dimension,1:matrix_dimension) ) 
  CALL dread_matrix(bundle%global_MI%mat,matrix_dimension,matrix_dimension,matrix_dimension,file_unit)
  
  read(file_unit,*,ERR=9000)    ! comment line
  read(file_unit,*)matrix_dimension
  bundle%global_MV%dim=matrix_dimension
  ALLOCATE( bundle%global_MV%mat(1:matrix_dimension,1:matrix_dimension) ) 
  CALL dread_matrix(bundle%global_MV%mat,matrix_dimension,matrix_dimension,matrix_dimension,file_unit)
  
  read(file_unit,*,ERR=9000)    ! comment line
  read(file_unit,*)matrix_dimension
  bundle%global_L%dim=matrix_dimension
  ALLOCATE( bundle%global_L%mat(1:matrix_dimension,1:matrix_dimension) ) 
  CALL dread_matrix(bundle%global_L%mat,matrix_dimension,matrix_dimension,matrix_dimension,file_unit)
  
  read(file_unit,*,ERR=9000)    ! comment line
  read(file_unit,*)matrix_dimension
  bundle%global_C%dim=matrix_dimension
  ALLOCATE( bundle%global_C%mat(1:matrix_dimension,1:matrix_dimension) ) 
  CALL dread_matrix(bundle%global_C%mat,matrix_dimension,matrix_dimension,matrix_dimension,file_unit)

  read(file_unit,*,ERR=9000) ! comment line 
  CALL read_Sfilter_matrix( bundle%global_Z,file_unit )
    
  read(file_unit,*,ERR=9000) ! comment line 
  CALL read_Sfilter_matrix( bundle%global_Y,file_unit )
  
! read the loss model information for each conductor

  ALLOCATE( bundle%conductor_impedance(1:bundle%tot_n_conductors) )
  
  read(file_unit,*)  ! comment line
  do conductor=1,bundle%tot_n_conductors
    CALL read_conductor_impedance_model(bundle%conductor_impedance(conductor),file_unit)
  end do
  
! read the x and y position for each conductor: this is needed for incident field excitation

  ALLOCATE( bundle%conductor_x_offset(1:bundle%tot_n_conductors) )
  ALLOCATE( bundle%conductor_y_offset(1:bundle%tot_n_conductors) )

  read(file_unit,*)  ! comment line
  do conductor=1,bundle%tot_n_conductors
    read(file_unit,*)bundle%conductor_x_offset(conductor),bundle%conductor_y_offset(conductor)
  end do
  
 ! numbering information
  read(file_unit,*)  ! comment line
  
  ALLOCATE( bundle%terminal_conductor_is_shield_flag(1:bundle%tot_n_conductors) )
  ALLOCATE( bundle%terminal_conductor_to_inner_domain(1:bundle%tot_n_conductors) )
  ALLOCATE( bundle%terminal_conductor_to_outer_domain(1:bundle%tot_n_conductors) )
  ALLOCATE( bundle%terminal_conductor_to_global_domain_conductor(1:bundle%tot_n_conductors) )
  ALLOCATE( bundle%terminal_conductor_to_local_domain_conductor(1:bundle%tot_n_conductors) )
  ALLOCATE( bundle%terminal_conductor_to_reference_terminal_conductor(1:bundle%tot_n_conductors) )
  
  do conductor=1,bundle%tot_n_conductors
  
    read(file_unit,*)bundle%terminal_conductor_is_shield_flag(conductor),      &
                     bundle%terminal_conductor_to_inner_domain(conductor),     &
                     bundle%terminal_conductor_to_outer_domain(conductor),     &
                     bundle%terminal_conductor_to_global_domain_conductor(conductor), &
                     bundle%terminal_conductor_to_local_domain_conductor(conductor), &
                     bundle%terminal_conductor_to_reference_terminal_conductor(conductor)
  end do
  
! read the conductor labels
  ALLOCATE( bundle%conductor_label(1:bundle%tot_n_conductors) )
  read(file_unit,*)  ! comment line
  do conductor=1,bundle%tot_n_conductors
    read(file_unit,'(A)')bundle%conductor_label(conductor)
  end do

! close .bundle file

  CLOSE(unit=file_unit)

  if (verbose) write(*,*)'Closed file:',trim(filename)

  RETURN
  
9000  run_status='ERROR reading the bundle file:'//trim(filename)
      CALL write_program_status()
      STOP 1
  
  END SUBROUTINE read_cable_bundle

! NAME
!     SUBROUTINE write_cable_bundle
!
!     write the cable bundle structure to a specified unit
!     
! COMMENTS
!    
!
! HISTORY
!    started 2/12/2015 CJS
!    stage 3 developments started 21/3/2015 CJS
!    generalise conductor impedance for transfer impedance model (stage 5) CJS 12/5/2016
!

  SUBROUTINE write_cable_bundle(bundle,file_unit)

  USE type_specifications
  USE constants
  USE general_module
  USE cable_module
  USE maths

  IMPLICIT NONE

! variables passed to subroutine

  type(bundle_specification_type),intent(INOUT) :: bundle
  integer,intent(IN)                            :: file_unit
  
! local variables

  character(len=filename_length)    :: filename
  character(len=line_length)        :: line
  integer                :: cable
  integer                :: domain
  integer                :: matrix_dimension
  integer                :: conductor
  integer                :: n_cables_to_write

! START

! Open the output (.bundle) file

  filename=trim(MOD_bundle_lib_dir)//trim(bundle%bundle_name)//bundle_file_extn
  open(unit=file_unit,file=trim(filename))

! write .bundle file
  write(file_unit,'(A)')trim(bundle%version)
  write(file_unit,'(A)')trim(bundle%bundle_name)
  
  write(file_unit,*)bundle%n_cables_without_ground_plane,'  ! number of cables not including ground plane'
  write(file_unit,*)bundle%n_cables,'  ! number of cables, cable name and x y coordinates follow...'

! write cable names and coordinates
  
  do cable=1,bundle%n_cables_without_ground_plane
    write(file_unit,'(A)')trim(bundle%cable(cable)%cable_name)
    write(file_unit,*)bundle%cable_x_offset(cable),bundle%cable_y_offset(cable),    &
                      bundle%cable_angle(cable),' x y coordinates and angle of cable '
  end do ! next cable
   
! write ground plane specification  
  
  if (bundle%ground_plane_present) then
    write(file_unit,'(A)')'ground_plane'
! note: convert angle back to degrees
    write(file_unit,*)bundle%ground_plane_x,bundle%ground_plane_y,bundle%ground_plane_angle*180d0/pi, &
                      ' x y coordinates and angle of ground plane '
    write(file_unit,*)bundle%ground_plane_nx,bundle%ground_plane_ny,' ground plane normal direction'
    write(file_unit,*)bundle%ground_plane_cable_side,' orientation of cables wrt ground plane'

    write(file_unit,*)bundle%ground_plane_sigma,' ground plane conductivity'
    write(file_unit,*)bundle%ground_plane_w,' ground plane width '  
    write(file_unit,*)bundle%ground_plane_h,' ground plane height (thickness) '      
    write(file_unit,*)bundle%ground_plane_Rdc,' ground plane Rdc '
    write(file_unit,*)bundle%ground_plane_nsegx,' ground plane nsegx '
    write(file_unit,*)bundle%ground_plane_nsegz,' ground plane nsegz '
    write(file_unit,*)bundle%ground_plane_nh,' ground plane nh '

  else
    write(file_unit,'(A)')'no_ground_plane'
  end if  
  
  write(file_unit,*)bundle%tot_n_conductors,'  # total number of conductors'
  write(file_unit,*)bundle%tot_n_external_conductors,'  # total number of external conductors'
  write(file_unit,*)bundle%system_dimension,'  # dimension of the matrix system characterising the MTL propagation'
  
  write(file_unit,*)bundle%tot_n_domains,'  # number of domains'

  do domain=1,bundle%tot_n_domains
  
    write(file_unit,*)'Domain number',domain
    write(file_unit,*)bundle%n_conductors(domain),' number of conductors in this domain'   
    
    write(file_unit,*)'terminal_conductor_list'
    write(file_unit,*)bundle%terminal_conductor_list(domain)%n_elements,' ! number of elements '
    do conductor=1,bundle%terminal_conductor_list(domain)%n_elements
      write(file_unit,*)bundle%terminal_conductor_list(domain)%element(conductor)
    end do
    
    matrix_dimension=bundle%L(domain)%dim
    write(file_unit,*)'Per-Unit-length Inductance Matrix, [L]'
    write(file_unit,*)matrix_dimension,' Dimension of [L]'
    CALL dwrite_matrix(bundle%L(domain)%mat,matrix_dimension,matrix_dimension,matrix_dimension,file_unit)
  
    write(file_unit,*)'Per-Unit-length Capacitance Matrix, [C]'
    write(file_unit,*)matrix_dimension,' Dimension of [C]'
    CALL dwrite_matrix(bundle%C(domain)%mat,matrix_dimension,matrix_dimension,matrix_dimension,file_unit)

    write(file_unit,*)'Per-Unit-length Impedance Matrix, Z'
    CALL write_Sfilter_matrix( bundle%Z(domain),file_unit )
      
    write(file_unit,*)'Per-Unit-length Admittance Matrix, Y'
    CALL write_Sfilter_matrix( bundle%Y(domain),file_unit )

  end do
  
  write(file_unit,*)'# Global current domain transformation matrix, [MI]'
  matrix_dimension=bundle%global_MI%dim-1 
  write(file_unit,*)matrix_dimension,' Dimension of [MI]'
  CALL dwrite_matrix(bundle%global_MI%mat,matrix_dimension,matrix_dimension,bundle%global_MI%dim,file_unit)
  
  write(file_unit,*)'# Global voltage domain transformation matrix, [MV]'
  matrix_dimension=bundle%global_MV%dim-1
  write(file_unit,*)matrix_dimension,' Dimension of [MV]'
  CALL dwrite_matrix(bundle%global_MV%mat,matrix_dimension,matrix_dimension,bundle%global_MI%dim,file_unit)
  
  write(file_unit,*)'# Global domain based inductance matrix, [L]'
  matrix_dimension=bundle%global_L%dim
  write(file_unit,*)matrix_dimension,' Dimension of [L]'
  CALL dwrite_matrix(bundle%global_L%mat,matrix_dimension,matrix_dimension,matrix_dimension,file_unit)
  
  write(file_unit,*)'# Global domain based capacitance matrix, [C]'
  matrix_dimension=bundle%global_C%dim
  write(file_unit,*)matrix_dimension,' Dimension of [C]'
  CALL dwrite_matrix(bundle%global_C%mat,matrix_dimension,matrix_dimension,matrix_dimension,file_unit)
 
  write(file_unit,*)'Global domain based Per-Unit-length Impedance Matrix, Z'
  CALL write_Sfilter_matrix( bundle%global_Z,file_unit )
      
  write(file_unit,*)'Global domain based Per-Unit-length Admittance Matrix, Y'
  CALL write_Sfilter_matrix( bundle%global_Y,file_unit )
  
! write the loss model information for each conductor
  write(file_unit,*)' # conductor impedance models' 
  do conductor=1,bundle%tot_n_conductors
  
    CALL write_conductor_impedance_model(bundle%conductor_impedance(conductor),file_unit)
                      
  end do
  
! write the x and y position for each conductor: this is needed for incident field excitation
  write(file_unit,*)' # conductor x y positions' 
  do conductor=1,bundle%tot_n_conductors
    write(file_unit,*)bundle%conductor_x_offset(conductor),bundle%conductor_y_offset(conductor)
  end do

 ! numbering information
  write(file_unit,'(A)')'is_shield  tc_to_in_domain  tc_to_out_domain  tc_to_gdc    tc_to_ldc    tc_to_ref_tc  ' 

  do conductor=1,bundle%tot_n_conductors
  
    write(file_unit,8000)bundle%terminal_conductor_is_shield_flag(conductor),      &
                         bundle%terminal_conductor_to_inner_domain(conductor),     &
                         bundle%terminal_conductor_to_outer_domain(conductor),     &
                         bundle%terminal_conductor_to_global_domain_conductor(conductor), &
                         bundle%terminal_conductor_to_local_domain_conductor(conductor), &
                         bundle%terminal_conductor_to_reference_terminal_conductor(conductor)
8000 format(L5,5I15)
  end do
 
! write the conductor labels
  write(file_unit,*)' # Conductor labels' 
  do conductor=1,bundle%tot_n_conductors
    write(file_unit,*)trim(bundle%conductor_label(conductor))
  end do

! close .bundle file

  CLOSE(unit=file_unit)
  
  RETURN

  END SUBROUTINE write_cable_bundle

! NAME
!     SUBROUTINE deallocate_cable_bundle
!
!     deallocate cable_bundle structure data
!     
! COMMENTS
!     
!
! HISTORY
!    started 2/12/2015 CJS
!    generalise conductor impedance for transfer impedance model (stage 5) CJS 12/5/2016
!

  SUBROUTINE deallocate_cable_bundle(bundle)

  USE type_specifications
  USE general_module
  USE cable_module

  IMPLICIT NONE

! variables passed to subroutine

  type(bundle_specification_type),intent(INOUT)    :: bundle
  
! local variables

  integer    :: cable
  integer    :: domain
  integer    :: conductor

! START

  do cable=1,bundle%n_cables
    CALL deallocate_cable(bundle%cable(cable))
  end do
  
  
  if (allocated(bundle%cable)) DEALLOCATE( bundle%cable )
  if (allocated(bundle%cable_x_offset )) DEALLOCATE( bundle%cable_x_offset )
  if (allocated(bundle%cable_y_offset )) DEALLOCATE( bundle%cable_y_offset )
  if (allocated(bundle%cable_angle ))    DEALLOCATE( bundle%cable_angle )
  
  if (allocated(bundle%n_conductors )) DEALLOCATE( bundle%n_conductors ) 
  
  if (allocated(bundle%terminal_conductor_list )) then  
    do domain=1,bundle%tot_n_domains    
      if (allocated(bundle%terminal_conductor_list(domain)%element )) &
        DEALLOCATE( bundle%terminal_conductor_list(domain)%element )
    end do 
    DEALLOCATE( bundle%terminal_conductor_list ) 
  end if
 
  if (allocated(bundle%L )) then  
    do domain=1,bundle%tot_n_domains    
      if (allocated(bundle%L(domain)%mat )) DEALLOCATE( bundle%L(domain)%mat )
    end do 
    DEALLOCATE( bundle%L ) 
  end if
 
  if (allocated(bundle%C )) then  
    do domain=1,bundle%tot_n_domains    
      if (allocated(bundle%C(domain)%mat )) DEALLOCATE( bundle%C(domain)%mat )
    end do 
    DEALLOCATE( bundle%C ) 
  end if
  
  if (ALLOCATED(bundle%Z)) then  
    do domain=1,bundle%tot_n_domains
      CALL deallocate_Sfilter_matrix( bundle%Z(domain) )
    end do  
    DEALLOCATE(bundle%Z)
  end if
  
  if (ALLOCATED(bundle%Y)) then  
    do domain=1,bundle%tot_n_domains
      CALL deallocate_Sfilter_matrix( bundle%Y(domain) )
    end do  
    DEALLOCATE(bundle%Y)
  end if
   
  if (allocated(bundle%global_MI%mat)) DEALLOCATE( bundle%global_MI%mat )
  if (allocated(bundle%global_MV%mat)) DEALLOCATE( bundle%global_MV%mat )
  
  if (allocated(bundle%global_L%mat)) DEALLOCATE( bundle%global_L%mat )
  if (allocated(bundle%global_C%mat)) DEALLOCATE( bundle%global_C%mat )
  
  CALL deallocate_Sfilter_matrix( bundle%global_Z )
  CALL deallocate_Sfilter_matrix( bundle%global_Y )

  if (ALLOCATED(bundle%conductor_impedance)) then  
    do conductor=1,bundle%tot_n_conductors
      CALL deallocate_conductor_impedance_model(bundle%conductor_impedance(conductor))
    end do
    DEALLOCATE(bundle%conductor_impedance)
  end if
  
  if (ALLOCATED( bundle%conductor_x_offset )) DEALLOCATE( bundle%conductor_x_offset )
  if (ALLOCATED( bundle%conductor_y_offset )) DEALLOCATE( bundle%conductor_y_offset )

! Numbering information required for transfer impedance calculation

  if(ALLOCATED( bundle%terminal_conductor_is_shield_flag ))  DEALLOCATE( bundle%terminal_conductor_is_shield_flag )
  if(ALLOCATED( bundle%terminal_conductor_to_inner_domain )) DEALLOCATE( bundle%terminal_conductor_to_inner_domain )
  if(ALLOCATED( bundle%terminal_conductor_to_outer_domain )) DEALLOCATE( bundle%terminal_conductor_to_outer_domain )
  if(ALLOCATED( bundle%terminal_conductor_to_global_domain_conductor ))  &
    DEALLOCATE( bundle%terminal_conductor_to_global_domain_conductor )
  if(ALLOCATED( bundle%terminal_conductor_to_local_domain_conductor ))  &
    DEALLOCATE( bundle%terminal_conductor_to_local_domain_conductor )
  if(ALLOCATED( bundle%terminal_conductor_to_reference_terminal_conductor )) &
    DEALLOCATE( bundle%terminal_conductor_to_reference_terminal_conductor )
  
  if (ALLOCATED( bundle%conductor_label )) then
    DEALLOCATE(bundle%conductor_label)
  end if

  RETURN

  END SUBROUTINE deallocate_cable_bundle
!
! NAME
!     SUBROUTINE check_cables_wrt_ground_plane
!
!     check that the cables are all on one side of the ground plane to ensure consistency
!     and indicate which side the cables are on by setting bundle%ground_plane_cable_side
!     
! COMMENTS
!     
!
! HISTORY
!    started 29/06/2016 CJS
!

  SUBROUTINE check_cables_wrt_ground_plane(bundle)

  USE type_specifications
  USE general_module
  USE cable_module

  IMPLICIT NONE

! variables passed to subroutine

  type(bundle_specification_type),intent(INOUT)    :: bundle
  
! local variables

  integer    :: cable
  
  real(dp)   :: norm_x,norm_y
  real(dp)   :: p
  real(dp)   :: gp_offset
  real(dp)   :: gp_offset_cable_1
  
! START

! calculate the ground plane normal direction

  norm_x=bundle%ground_plane_nx
  norm_y=bundle%ground_plane_ny

! The equation of the ground plane is r.norm-P=0, calculate P here  
  P=bundle%ground_plane_x*norm_x+bundle%ground_plane_y*norm_y

! loop over all of the cables, calculating the offset from the ground plane in the normal direction
! offset=r.norm-P, and check that all offsets have the same sign

 if (verbose) then
   write(*,*)'Checking the orientation of cables wrt the ground plane'
   write(*,*)'GP normal:',norm_x,norm_y
   write(*,*)'GP offset:',bundle%ground_plane_x+bundle%ground_plane_y
   write(*,*)'P=',P
 end if

  do cable=1,bundle%n_cables_without_ground_plane
  
    gp_offset=bundle%cable_x_offset(cable)*norm_x+bundle%cable_y_offset(cable)*norm_y-P
    
    if (verbose) then
      write(*,*)'Cable',cable,' offset',gp_offset
    end if
    
    if (cable.EQ.1) then
    
      gp_offset_cable_1=gp_offset
      
    else
! check consistency of ground plane offsets

      if (gp_offset_cable_1*gp_offset.LE.0d0) then
        run_status='ERROR there are cables both sides of the ground plane'
        CALL write_program_status()
        STOP 1
      end if
          
    end if  ! cable.NE.1
    
  end do  ! next cable to check
  
  if (gp_offset_cable_1.GT.0d0) then
    bundle%ground_plane_cable_side=1
  else
    bundle%ground_plane_cable_side=-1
  end if
  
  RETURN

  END SUBROUTINE check_cables_wrt_ground_plane
!
! NAME
!     SUBROUTINE check_cable_intersections
!
!     check that the cables do not intersect
!     the intersection test checks intersection of the outer shape of each cable
!     which is usually the outer dielectric layer apart from overshields or D connectors which
!     are assumed not to have dielectric
!     
! COMMENTS
!     Revised to work with all cable types
!
! HISTORY
!    started 8/10/2016 CJS
!
!
  SUBROUTINE check_cable_intersection(bundle)

  USE type_specifications
  USE general_module
  USE cable_module

  IMPLICIT NONE

! variables passed to subroutine

  type(bundle_specification_type),intent(IN)    :: bundle
  
! local variables

  integer    :: cable1,cable2
  integer    :: nec1,nec2,ec1,ec2
  
  integer    :: shape1,shape2
  integer    :: type1,type2
  
  real(dp)   :: ox,oy,theta
  real(dp)   :: r

  real(dp)   :: wd,hd
  integer    :: nc,ncrow(2)
  real(dp)   :: rw,p,s,o,W(2)
  
  integer :: npts1
  real(dp),allocatable :: shape1_x(:)
  real(dp),allocatable :: shape1_y(:)
   
  integer :: npts2
  real(dp),allocatable :: shape2_x(:)
  real(dp),allocatable :: shape2_y(:)

  logical :: intersect
  logical :: intersect2
  logical :: nested_1_in_2
  logical :: nested_2_in_1
  logical :: gp_intersect
  
  logical :: intersection_found
  
! START

  intersection_found=.FALSE.

! loop over cable 1
  do cable1=1,bundle%n_cables_without_ground_plane
  
    ox=bundle%cable_x_offset(cable1)
    oy=bundle%cable_y_offset(cable1)
    theta=bundle%cable_angle(cable1)
    type1=bundle%cable(cable1)%cable_type
    
    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
    end if
    
! loop over the external conductors of cable 1
    do ec1=1,nec1
    
      shape1=bundle%cable(cable1)%external_model(ec1)%conductor_type
         
! generate a list of points on cable 1 outer surface
      if (shape1.EQ.rectangle) then
    
        wd=bundle%cable(cable1)%external_model(1)%dielectric_width
        hd=bundle%cable(cable1)%external_model(1)%dielectric_height
        CALL generate_rectangle_points(npts1,shape1_x,shape1_y,ox,oy,theta,wd,hd)
      
      else if (shape1.EQ.circle) then
    
        r=bundle%cable(cable1)%external_model(ec1)%dielectric_radius
        CALL generate_circle_points(npts1,shape1_x,shape1_y,ox,oy,r)
     
      else if (shape1.EQ.Dshape) then
    
        nc=bundle%cable(cable1)%tot_n_conductors 
        rw=bundle%cable(cable1)%parameters(1)
        p=bundle%cable(cable1)%parameters(2)
        s=bundle%cable(cable1)%parameters(3)
        o=bundle%cable(cable1)%parameters(4)
        ncrow(1)=nc/2
        ncrow(2)=(nc-1)-ncrow(1)
        w(1)=(ncrow(1)-1)*p
        w(2)=(ncrow(2)-1)*p
        CALL generate_Dshape_points(npts1,shape1_x,shape1_y,ox,oy,w(1)/2d0,w(2)/2d0,s/2d0,rw+o,theta)
      
      end if
          
      if (bundle%ground_plane_present) then
! check for intersection of the cables with the ground plane
        gp_intersect=.FALSE.
        
        if (verbose) write(*,*)'TESTING: intersection of cable',cable1,' with the ground plane'

        CALL gptest(npts1,shape1_x,shape1_y,gp_intersect)

        if (gp_intersect) then
          if (verbose) write(*,*)'Cable ',cable1,' intersects the ground plane'
          intersection_found=.TRUE.
        end if
              
      end if
      
! loop over cable 2
      do cable2=cable1+1,bundle%n_cables_without_ground_plane
  
        ox=bundle%cable_x_offset(cable2)
        oy=bundle%cable_y_offset(cable2)
        theta=bundle%cable_angle(cable2)
        type2=bundle%cable(cable2)%cable_type
    
        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
        end if
    
! loop over the external conductors of cable 2
        do ec2=1,nec2

          shape2=bundle%cable(cable2)%external_model(ec2)%conductor_type
         
! generate a list of points on cable 1 outer surface
          if (shape2.EQ.rectangle) then
    
            wd=bundle%cable(cable2)%external_model(1)%dielectric_width
            hd=bundle%cable(cable2)%external_model(1)%dielectric_height
            CALL generate_rectangle_points(npts2,shape2_x,shape2_y,ox,oy,theta,wd,hd)
      
          else if (shape2.EQ.circle) then
    
            r=bundle%cable(cable2)%external_model(ec2)%dielectric_radius
            CALL generate_circle_points(npts2,shape2_x,shape2_y,ox,oy,r)
     
          else if (shape2.EQ.Dshape) then
    
            nc=bundle%cable(cable2)%tot_n_conductors 
            rw=bundle%cable(cable2)%parameters(1)
            p=bundle%cable(cable2)%parameters(2)
            s=bundle%cable(cable2)%parameters(3)
            o=bundle%cable(cable2)%parameters(4)
            ncrow(1)=nc/2
            ncrow(2)=(nc-1)-ncrow(1)
            w(1)=(ncrow(1)-1)*p
            w(2)=(ncrow(2)-1)*p
            CALL generate_Dshape_points(npts2,shape2_x,shape2_y,ox,oy,w(1)/2d0,w(2)/2d0,s/2d0,rw+o,theta)
      
          end if

! assume that there is no intersection and the shapes are not nested

          intersect=.FALSE.
          nested_1_in_2=.FALSE.
          nested_2_in_1=.FALSE.
          
! check for the intersection of cable 1 with cable 2 and flag if cable 1 is nested within cable 2.
        
          if (verbose) write(*,*)'TESTING: intersection of cable',cable1,' with cable ',cable2
          CALL shape_test(npts1,shape1_x,shape1_y,npts2,shape2_x,shape2_y,intersect,nested_1_in_2)

          if (intersect) then
            if (verbose) write(*,*)'Cable ',cable1,' intersects Cable',cable2
            intersection_found=.TRUE.
          end if

! Cables can only be nested if the outer cable is an overshield
          if ( nested_1_in_2 ) then
            if (type1.NE.cable_geometry_type_overshield) then
              if (verbose) write(*,*)'Cable ',cable2,' is inside Cable',cable1
              intersection_found=.TRUE.
            else
               if (verbose) write(*,*)'Cable ',cable2,' is OK inside Overshield Cable',cable1
            end if
          end if

! check for the intersection of cable 1 with cable 2 and flag if cable 1 is nested within cable 2.

          if (verbose) write(*,*)'TESTING: intersection of cable',cable2,' with cable ',cable1
          CALL shape_test(npts2,shape2_x,shape2_y,npts1,shape1_x,shape1_y,intersect2,nested_2_in_1)

          if (intersect) then
            if (verbose) write(*,*)'Cable ',cable2,' intersects Cable',cable1
            intersection_found=.TRUE.
          end if

! Cables can only be nested if the outer cable is an overshield          
          if ( nested_2_in_1 ) then
            if (type2.NE.cable_geometry_type_overshield) then
              if (verbose) write(*,*)'Cable ',cable1,' is inside Cable',cable2
              intersection_found=.TRUE.
            else
              if (verbose) write(*,*)'Cable ',cable1,' is OK inside Overshield Cable',cable2
            end if
          end if
          
          DEALLOCATE( shape2_x )
          DEALLOCATE( shape2_y )
                           
        end do ! next external conductor of cable 2
               
      end do ! next cable 2  
          
      DEALLOCATE( shape1_x )
      DEALLOCATE( shape1_y )
    
    end do  ! next external conductor of cable1
  
  end do  ! next cable1
  
  if (intersection_found) then
    run_status='ERROR there are cables in the bundle which intersect'
    CALL write_program_status()
    STOP 1
  end if
    
  RETURN

  END SUBROUTINE check_cable_intersection
!
! NAME
!     SUBROUTINE gptest
!
!     check that a cable does not intersect the ground plane
!     
! COMMENTS
!     
!
! HISTORY
!    started 20/4/2017 CJS
!
!
SUBROUTINE gptest(npts,shape_x,shape_y,gp_intersect)

USE type_specifications

IMPLICIT NONE

integer,intent(IN) :: npts
real(dp),intent(IN) :: shape_x(npts)
real(dp),intent(IN) :: shape_y(npts)

logical,intent(OUT) :: gp_intersect

! local variables

integer :: i

! START

do i=1,npts

  if (shape_y(i).LE.0D0) then
    gp_intersect=.TRUE.
    RETURN
  end if

end do ! next point

RETURN

END SUBROUTINE gptest
!
! NAME
!     SUBROUTINE shape_test
!
!     check that cables do not intersect, also flag whether the first cable is nested 
!     i.e. lies entirely within the second
!     
! COMMENTS
!     We assume that the shape of the cable is convex and that the shape is closed
!
! HISTORY
!    started 20/4/2017 CJS
!
!
SUBROUTINE shape_test(npts1,shape1_x,shape1_y,npts2,shape2_x,shape2_y,intersect,nested)

USE type_specifications

IMPLICIT NONE

integer,intent(IN)  :: npts1
real(dp),intent(IN) :: shape1_x(npts1)
real(dp),intent(IN) :: shape1_y(npts1)

integer,intent(IN)  :: npts2
real(dp),intent(IN) :: shape2_x(npts2)
real(dp),intent(IN) :: shape2_y(npts2)

logical,intent(OUT) :: intersect            ! flag to indicate that shapes 1 and 2 intersect
logical,intent(OUT) :: nested               ! flag to indicate shape 2 lies within shape 1

! local variables

integer :: i,j

real(dp) :: vx1,vy1
real(dp) :: vx2,vy2
real(dp) :: vp
real(dp) :: dirn

logical :: inside_point_found
logical :: outside_point_found

logical :: local_outside_point_found

! START

intersect=.FALSE.
nested=.FALSE.

inside_point_found=.FALSE.
outside_point_found=.FALSE.

! loop around test points in shape 2
do j=1,npts2

  local_outside_point_found=.FALSE.

  dirn=0.0  ! set to zero initially to indicate that it is not known
  
! loop around line segments in shape 1. Note there is one less line segment than the number of points
  do i=1,npts1-1

! vector along line segment in shape 1
    vx1=shape1_x(i+1)-shape1_x(i)
    vy1=shape1_y(i+1)-shape1_y(i)
    
! vector from first point in line segment to test point
    vx2=shape2_x(j)-shape1_x(i)
    vy2=shape2_y(j)-shape1_y(i)
    
! calculate the z component of the vector product of the two vectors
    vp=vx1*vy2-vx2*vy1
    
    if (dirn.EQ.0D0) then    ! maybe we need to use abs(dirn).LT.small here...
    
! this is the first value so this sets the direction to test for all other sets of points

      dirn=vp
      
    else
    
! test whether the sign of the direction has changed

      if (dirn*vp.LT.0.0) then
      
! the signs are different then this indicates a test point outside shape 1
  
        local_outside_point_found=.TRUE.
      
      end if  ! change in sign of direction
      
    end if  ! first test line segment.

  end do  ! next line segment in shape 1
  
  if (local_outside_point_found) then
    outside_point_found=.TRUE.
  else
    inside_point_found=.TRUE.  
  end if

end do ! next test point in shape 2

! work out whether we have an intersection

if (inside_point_found.AND.outside_point_found) intersect=.TRUE.
  
! work out whether we have a nested shape

if (inside_point_found.AND.(.NOT.outside_point_found)) nested=.TRUE.

RETURN

END SUBROUTINE 

END MODULE cable_bundle_module