!
! 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:
! SUBROUTINE create_spice_subcircuit_model
!
! NAME
!     create_spice_subcircuit_model
!
! AUTHORS
!     Chris Smartt
!
! DESCRIPTION
!     This subroutine creates the Spice subcircuit model for the cable bundle
!     including transfer impedance coupling terms when requested
!     and incident field excitation model when requested
!     The input to the subroutine is the spice_bundle_model structure which contains 
!     all the information required to build the spice model.
!
!     The processs is split into the following stages:
! 
! STAGE 1: Initialisation: open the file for the spice subcircuit model
! STAGE 2: Allocate memory for node lists, propagation correction filter lists etc
! STAGE 3: Work out where the transfer impedance models and incident field sources need to go
! STAGE 4: Evaluate and write the d.c. resistance on each conductor
! STAGE_5: Implement the Domain decomposition in the Spice model
! For each domain:
!     STAGE_6, Within the current domain work out the modal transformation matrices,
!              mode impedances and mode delays to implement the modal decomposition in spice
!     STAGE_7, Within the current domain work out the modal propagation correction filter functions
!     STAGE_8, Write the modal decomposition spice model and modal propagation spice model
!
! STAGE_9, generate the information required for and then write the transfer impedance coupling model(s) as required
! STAGE_10, generate the information required for and then write the incident field excitation model as required
! STAGE 11: Deallocate the domain based information
!     
! COMMENTS
!     1. The notation needs to be organised a bit better.
!     2. We will need to split the node numbering, component naming and file writing processes into subroutines
!     otherwise it will become a complete mess... STARTED: writing processes have been separated out now.
!     3. There is locally allocated (domain based) memory for modal decomposition etc which is now saved globally
!        for the transfer impedance model. The local data arrays should be removed. 
!     4. Following the substitution of the revised modal decomposition process we need the check for unused variables and remove them
!
! HISTORY
!
!     started 9/12/2015 CJS: STAGE_1 developments
!     STAGE 2 developments started  2/02/2016 CJS 
!     STAGE_3 developments started 24/03/2016 CJS 
!     STAGE_4 developments started 22/04/2016 CJS 
!     STAGE_5 developments started 16/05/2016 CJS
!     STAGE_6 developments started 13/06/2016 CJS
!     20/06/2016 CJS: Implement a new, robust method for modal decomposition for lossless propagation in an inhomogeneous medium 
!                     i.e. what is required for the modal decomposition and method of characteristics model in Spice.
!     27/9/2016 CJS: Implement Xie's model for incident field excitaiton of shielded cables
!     15/12/2016 CJS: use the frequency_spec module for all frequency dependent processes
!     17/5/2017 CJS : fix problem with incident field excitation with impedance on the reference conductor for the revised termination model.
!     12/6/2017 CJS : fix problem which occurs if a domain is the victim domain for more than one transfer impedance model
!     16/11/2017 CJS Include network synthesis process to replace s-domain transfer functions
!     October 2023 CJS Proximity_effects branch: Embed FastHenry2 frequency dependent proximity effects
!
SUBROUTINE create_spice_subcircuit_model(spice_bundle_model)

USE type_specifications
USE general_module
USE constants
USE cable_module
USE cable_bundle_module
USE spice_cable_bundle_module
USE MTL_analytic_solution
USE maths
USE filter_module
USE frequency_spec

IMPLICIT NONE

! variables passed to the subroutine

TYPE(spice_model_specification_type),INTENT(IN)    :: spice_bundle_model

! local variables

! structure for modal decomposition information
TYPE::mode_decomp_list_type
  integer              :: n_modes
  real(dp),allocatable :: TVI(:,:)
  real(dp),allocatable :: TI(:,:)
  real(dp),allocatable :: delay(:)
  real(dp),allocatable :: impedance(:)
END TYPE mode_decomp_list_type

! structure to hold a list of filter functions
TYPE::Sfilter_list_type
  integer             :: n_modes
  Type(Sfilter),allocatable     :: Hp(:)
END TYPE Sfilter_list_type

! structure to hold a list of nodes
TYPE::node_list_type
  integer             :: n_nodes
  integer,allocatable :: node(:)
END TYPE node_list_type

character(len=filename_length)    :: filename

integer :: tot_n_conductors
integer :: n_domains
integer :: n_conductors
integer :: n_modes
integer :: domain
integer :: mode
integer :: end

integer :: global_domain_conductor         ! conductor number in the global domain nunbering
integer :: local_domain_conductor          ! conductor number in the local domain
integer :: terminal_conductor              ! conductor number at the sub-circuit terminals 
integer :: shield_conductor                ! terminal conductor number for a shield for transfer impedance model

integer :: conductor        ! conductor loop variable

! node numbering and component name stuff

integer :: next_free_node  ! this is the number of the next spice node to be used

integer :: vref_node       ! this is the reference node used within the transmission line sub-circuit

! node numbers for external connection of subcircuit at both ends of the transmission line
integer,allocatable :: external_end1_nodes(:)
integer,allocatable :: external_end2_nodes(:)

! node numbers for d.c. resistance at both ends of the transmission line
integer,allocatable :: rdc_end1_nodes(:)
integer,allocatable :: rdc_end2_nodes(:)

! node numbers for domain based transmission lines for all domains together
integer,allocatable :: all_domain_end1_nodes(:)
integer,allocatable :: all_domain_end2_nodes(:)

! node numbers for domain based transmission lines within a single domain
integer,allocatable :: domain_end1_nodes(:)
integer,allocatable :: domain_end2_nodes(:)

! node numbers for modal transmission lines
integer,allocatable :: modal_end1_nodes(:)
integer,allocatable :: modal_end2_nodes(:)

integer        :: dim     ! matrix dimension

! modal decomposition matrices
real(dp),allocatable     :: real_MI(:,:)
real(dp),allocatable     :: real_MV(:,:)
real(dp),allocatable     :: real_MVI(:,:)

! domain decomposition matrices
real(dp),allocatable     :: real_TII(:,:)
real(dp),allocatable     :: real_TV(:,:)
real(dp),allocatable     :: real_TI(:,:)
real(dp),allocatable     :: real_TVI(:,:)

! mode velocity and impedance within a domain. 
real(dp),allocatable     :: mode_velocity(:)
real(dp),allocatable     :: mode_impedance(:)

! Filter list used in the calculation of propagation correction filter functions
Type(Sfilter_list_type),allocatable     :: domain_Hp_list(:)

real(dp),allocatable     :: Rdc(:)           ! d.c. resistance of a conductor
complex(dp)              :: Zint_c           ! complex surface impedance of a conductor

real(dp)                 :: Rdc_t            ! d.c. resistance of a shield evaluated from the transfer impedance model
complex(dp)              :: Zint_t           ! complex transfer impedance of a shield conductor

! domain based list of conductor impedance models, used for propagation correction calcuulation
type(conductor_impedance_model),allocatable :: conductor_impedance_domain(:)
  
type(node_list_type),allocatable      :: domain_MOC_reference_end1_nodes(:) ! list of reference nodes for method of characteristics
type(node_list_type),allocatable      :: domain_MOC_reference_end2_nodes(:) ! list of reference nodes for method of characteristics

integer :: impedance_model_type       ! local type for a conductor impedance model

! Transfer impedance model stuff...

integer :: Zt_model                                     ! loop variable for the transfer impedance model number to be included

! transfer impedance model source domain information
integer :: source_domain
integer :: n_source_domain_modes

! transfer impedance model victim domain information
integer :: victim_domain              
integer :: n_victim_domain_modes
 
integer,allocatable :: ZT_source_domain(:)
integer,allocatable :: ZT_victim_domain(:)
  
type(node_list_type),allocatable      :: domain_MOC_ZT_reference_end1_nodes(:) ! list of reference nodes for ZT models: method of characteristics
type(node_list_type),allocatable      :: domain_MOC_ZT_reference_end2_nodes(:) ! list of reference nodes for ZT models: method of characteristics
  
type(node_list_type),allocatable      :: domain_ZT_internal_end1_nodes(:) ! list of internal victim domain nodes for ZT models
type(node_list_type),allocatable      :: domain_ZT_internal_end2_nodes(:) ! list of internal victim domain nodes for ZT models

type(node_list_type),allocatable      :: domain_MOC_V_plus_nodes(:) ! list of nodes for the +z travelling characteristic mode voltage in each domain
type(node_list_type),allocatable      :: domain_MOC_V_minus_nodes(:) ! list of nodes for the +z travelling characteristic mode voltage in each domain

type(mode_decomp_list_type),allocatable :: domain_mode_decomp(:)  ! domain base modal decompositiion data

logical,allocatable :: domain_is_ZT_victim(:)       ! list of flags to indicate transfer impedance victim domains
integer,allocatable :: n_ZT_victim_in_domain(:)     ! number of transfer impedance models with this victim domain
integer,allocatable :: count_ZT_victim_in_domain(:) ! counter for transfer impedance models within a domain

logical,allocatable :: domain_is_Einc_ZT_victim(:)  ! list of flags to indicate victim domains for incident field excitation and transfer impedance coupling (Xie model)

type(node_list_type) :: ZT_node1_end1,ZT_node1_end2              ! nodes used for the transfer impedance source terms in the victim domain
type(node_list_type) :: ZT_node2_end1,ZT_node2_end2

integer :: source_domain_shield_conductor           ! conductor number for shield in source domain
integer :: victim_domain_shield_conductor           ! conductor number for shield in victim domain

integer :: terminal_shield_conductor                ! shield conductor number according to the terminal conductor numbering system
integer :: domain_shield_conductor                  ! shield conductor number according to the domain conductor numbering system

! Information for incident field excitation coupling via transfer impedance

type(node_list_type),allocatable      :: domain_MOC_EINC_ZT_reference_end1_nodes(:) ! list of reference nodes for Einc_ZT models: method of characteristics
type(node_list_type),allocatable      :: domain_MOC_EINC_ZT_reference_end2_nodes(:) ! list of reference nodes for Einc_ZT models: method of characteristics

! Information for incident field excitation

logical,allocatable :: domain_is_Einc_victim(:)     
integer :: n_einc_domain_modes

integer :: Einc_node1,Einc_node2        ! subcircuit nodes for connecting the external incident field source

real(dp) :: Tz                          ! incident field excitation delay time

! coordinates and constants required for the incident field excitation model. These variables follow the notation of  #EQN_REFERENCE_REQUIRED
real(dp),allocatable :: xcoord(:)
real(dp),allocatable :: ycoord(:)
real(dp),allocatable :: alpha0(:)
real(dp),allocatable :: alphaL(:)
real(dp),allocatable :: c1(:)
real(dp),allocatable :: c2(:)
real(dp),allocatable :: c1b(:)

! Additional general variables
integer :: i
integer,parameter :: write_matrices=0   ! used for debugging only

! string used to generate comments
character(len=max_spice_line_length)    :: comment

! error integer used in matrix inverse calculation
integer :: ierr

! START

! STAGE 1: Initialisation: open the file for the spice subcircuit model

  filename=trim(MOD_spice_bundle_lib_dir)//  &
           trim(spice_bundle_model%spice_model_name)//trim(spice_model_file_extn)
  
  OPEN(unit=spice_model_file_unit,file=filename)

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

! Copy some of the bundle information into local variables for clarity  

  n_domains=spice_bundle_model%bundle%tot_n_domains
  tot_n_conductors=spice_bundle_model%bundle%tot_n_conductors

! STAGE 2: Allocate memory for node lists, propagation correction filter lists, 
! Allocate the node number lists for both ends of the transmission line

  ALLOCATE( external_end1_nodes(1:tot_n_conductors) )
  ALLOCATE( external_end2_nodes(1:tot_n_conductors) )

  ALLOCATE( rdc_end1_nodes(1:tot_n_conductors) )
  ALLOCATE( rdc_end2_nodes(1:tot_n_conductors) )

  ALLOCATE( all_domain_end1_nodes(1:tot_n_conductors) )
  ALLOCATE( all_domain_end2_nodes(1:tot_n_conductors) )
  
! set the reference node for the transmission line subcircuit 
  next_free_node=first_subcircuit_node_number
    
  CALL create_new_node(vref_node,next_free_node)
  
! allocate the domain based information
  ALLOCATE( domain_HP_list(1:n_domains) )  
  
! domain based information for transfer impedance models, including those with incident field excitation (Xie's model)
  ALLOCATE( domain_mode_decomp(1:n_domains) )
  
  ALLOCATE( domain_MOC_V_plus_nodes(1:n_domains) )
  ALLOCATE( domain_MOC_V_minus_nodes(1:n_domains) )
  
  ALLOCATE( domain_MOC_reference_end1_nodes(1:n_domains) )
  ALLOCATE( domain_MOC_reference_end2_nodes(1:n_domains) )
  
  ALLOCATE( domain_MOC_ZT_reference_end1_nodes(1:n_domains) )
  ALLOCATE( domain_MOC_ZT_reference_end2_nodes(1:n_domains) )
  
  ALLOCATE( domain_ZT_internal_end1_nodes(1:n_domains) )
  ALLOCATE( domain_ZT_internal_end2_nodes(1:n_domains) )
  
  ALLOCATE( domain_MOC_EINC_ZT_reference_end1_nodes(1:n_domains) )
  ALLOCATE( domain_MOC_EINC_ZT_reference_end2_nodes(1:n_domains) )
  
  ALLOCATE( domain_is_ZT_victim(1:n_domains) )
  ALLOCATE( domain_is_Einc_victim(1:n_domains) )
  ALLOCATE( domain_is_Einc_ZT_victim(1:n_domains) )
  
  ALLOCATE( n_ZT_victim_in_domain(1:n_domains) )
  ALLOCATE( count_ZT_victim_in_domain(1:n_domains) )
 
  domain_is_ZT_victim(1:n_domains)=.FALSE.  
  domain_is_Einc_victim(1:n_domains)=.FALSE.
  domain_is_Einc_ZT_victim(1:n_domains)=.FALSE.
  
  n_ZT_victim_in_domain(1:n_domains)=0
  
  Einc_node1=0
  Einc_node2=0
  
  if (spice_bundle_model%include_incident_field) then
  
    domain_is_Einc_victim(n_domains)=.TRUE.
    
! create nodes for the sub-circuit header which will be used as the incident field excitation source
    
    CALL create_new_node(Einc_node1,next_free_node)
    
! create new node for Einc_node2
    CALL create_new_node(Einc_node2,next_free_node)
! OLD - caused problems if there were impedanced on reference conductor
!    Einc_node2=vref_node
    
  end if
  
! STAGE 3: Work out where the transfer impedance models and incident field sources need to go

  if (verbose) then
    write(*,*)''
    write(*,*)'Number of transfer impedance models:',spice_bundle_model%n_transfer_impedances
  end if
  
  ALLOCATE( ZT_source_domain(1:spice_bundle_model%n_transfer_impedances) )
  ALLOCATE( ZT_victim_domain(1:spice_bundle_model%n_transfer_impedances) )
  
  do ZT_model=1,spice_bundle_model%n_transfer_impedances
  
    conductor=spice_bundle_model%Zt_conductor(ZT_model)

! check that we have identified a shield conductor    
    if (.NOT.spice_bundle_model%bundle%terminal_conductor_is_shield_flag(conductor)) then
      write(run_status,*)'ERROR in create_spice_subcircuit_model.Conductor number',conductor,' is not a shield '
      CALL write_program_status()
      STOP 1
    end if

! check that we have a viable transfer impedance model for this shield

    impedance_model_type=spice_bundle_model%bundle%conductor_impedance(conductor)%impedance_model_type
    
    if (     (impedance_model_type.NE.impedance_model_type_cylindrical_shield) &
        .AND.(impedance_model_type.NE.impedance_model_type_filter)              ) then
    
      if (verbose) then
        write(*,*)'Error building transfer impedance models'
        write(*,*)'Conductor number',conductor,' does not have a transfer impedance model specified'
        write(*,*)'Condutor impedance model type:',impedance_model_type
      end if
      write(run_status,*)'ERROR in create_spice_subcircuit_model.Conductor number',conductor,&
                         ' does not have a transfer impedance model specified '
      CALL write_program_status()
      STOP 1
      
    end if
    
! get the source and victim domain for the transfer impedance model
    
    if (spice_bundle_model%Zt_direction(ZT_model).EQ.1) then
      source_domain=spice_bundle_model%bundle%terminal_conductor_to_inner_domain(conductor)
      victim_domain=spice_bundle_model%bundle%terminal_conductor_to_outer_domain(conductor)
    elseif (spice_bundle_model%Zt_direction(ZT_model).EQ.-1) then
      source_domain=spice_bundle_model%bundle%terminal_conductor_to_outer_domain(conductor)
      victim_domain=spice_bundle_model%bundle%terminal_conductor_to_inner_domain(conductor)
    else
      run_status='Error building transfer impedance models. Transfer impedance coupling direction should be +1 or -1'
      CALL write_program_status()
      STOP 1
    end if
    
    ZT_source_domain(ZT_model)=source_domain
    ZT_victim_domain(ZT_model)=victim_domain
    domain_is_ZT_victim(victim_domain)=.TRUE.
    
    n_ZT_victim_in_domain(victim_domain)=n_ZT_victim_in_domain(victim_domain)+1
    
    if (verbose) then
      write(*,*)'Transfer impedance model:',ZT_model,' source domain:',source_domain,' victim domain:',victim_domain
    end if
    
! check whether the source domain has an incident field excitation. If it does then the victim domain needs
! will have additional terms due to incident field excitation via transfer impedance (Xie's model)

    if (domain_is_Einc_victim(source_domain)) then
    
      if (use_xie) then
        domain_is_Einc_ZT_victim(victim_domain)=.TRUE.
      
        if (verbose) then
          write(*,*)'Transfer impedance model with incident field: source domain:',source_domain,' victim domain:',victim_domain
        end if
        
      end if ! use the Xie model for incident field excitation of shielded cables
      
    end if ! source domain has an incident field excitation 
   
  end do ! next transfer impedance model

! Write the subcircuit header, this also creates the node numbers 
! for the external connections to the transmission line model

  CALL write_spice_subcircuit_header( spice_bundle_model,       &
                                      next_free_node,tot_n_conductors,external_end1_nodes,external_end2_nodes,         &
                                      spice_bundle_model%include_incident_field,Einc_node1,Einc_node2 )

! STAGE 4: Evaluate and write the d.c. resistance on each conductor. The extraction of the d.c. resistance
! to the conductor terminations makes the modal decomposition a weak function of frequency. See Theory_Manual_Section 3.6

  ALLOCATE( Rdc(1:tot_n_conductors) )
  
  do conductor=1,tot_n_conductors
  
! evaluate conductor impedance. Note that only the conductor d.c. resistance is required at the moment
    CALL evaluate_conductor_impedance_model(spice_bundle_model%bundle%conductor_impedance(conductor), &
                                             0d0,Zint_c,Rdc(conductor),Zint_t,Rdc_t)
                                             
    if (high_freq_Zt_model) then  ! if this is a shield then don't lump the d.c. resistance
      if (spice_bundle_model%bundle%terminal_conductor_is_shield_flag(conductor)) then
        Rdc(conductor)=Rsmall
      end if
    end if

! don't include zero impedances in the Spice model. Replace with something very small instead.
    if (Rdc(conductor).LT.Rsmall) Rdc(conductor)=Rsmall

  end do

  CALL write_spice_comment('D.C. RESISTANCE END 1')
  end=1
  CALL write_spice_dc_resistances( vref_node,next_free_node,Rdc,tot_n_conductors,end,  &
                                   external_end1_nodes,rdc_end1_nodes,spice_bundle_model%length )
  
  CALL write_spice_comment('D.C. RESISTANCE END 2')
  end=2
  CALL write_spice_dc_resistances( vref_node,next_free_node,Rdc,tot_n_conductors,end,  &
                                   external_end2_nodes,rdc_end2_nodes,spice_bundle_model%length )
                                                          
! End of d.c. resistances                                                               

! STAGE_5 Implement the Domain decomposition in the Spice model
! call with the external node numbering and return with the domain decomposition node numbering
! See Theory_Manual_Section 3.2
      
! Allocate variables for the domain decomposition
  dim=spice_bundle_model%bundle%system_dimension
  ALLOCATE( real_MI(1:dim,1:dim) )
  ALLOCATE( real_MV(1:dim,1:dim) )
  ALLOCATE( real_MVI(1:dim,1:dim) )

! Copy the domain transformation matrices and calculate the inverse of MV
  real_MI(:,:)=spice_bundle_model%bundle%global_MI%mat(:,:)
  real_MV(:,:)=spice_bundle_model%bundle%global_MV%mat(:,:)
  
  ierr=0   ! set ierr=0 on input to matrix inverse to cause the program to stop if we have a singular matrix
  CALL dinvert_Gauss_Jordan(real_MV,dim,real_MVI,dim,ierr)
  
  CALL write_spice_comment('DOMAIN TRANSFORMATION END 1')
  end=1
  CALL write_spice_domain_decomposition_equivalent_circuit(vref_node,next_free_node,end,tot_n_conductors,dim,   &
                                                           rdc_end1_nodes,real_MVI,real_MI,all_domain_end1_nodes)  
  
  CALL write_spice_comment('DOMAIN TRANSFORMATION END 2')
  end=2
  CALL write_spice_domain_decomposition_equivalent_circuit(vref_node,next_free_node,end,tot_n_conductors,dim,   &
                                                           rdc_end2_nodes,real_MVI,real_MI,all_domain_end2_nodes)  
                                                                
  DEALLOCATE( real_MI )
  DEALLOCATE( real_MV )
  DEALLOCATE( real_MVI )
                                                          
! End of domain decomposition                                                                

! reset the counter for domain conductor numbers                                   
  global_domain_conductor=0
        
! Begin the main loop over domains
  do domain=1,n_domains
  
    write(comment,'(A,I6)')'DOMAIN ',domain
    CALL write_spice_comment(comment)
    
! STAGE_6, Within the current domain work out the modal transformation matrices,
!          mode impedances and mode delays to implement the modal decomposition in spice
! See Theory_Manual_Section 3.4
    
! Allocate variables for the modal decomposition
    n_conductors=spice_bundle_model%bundle%n_conductors(domain)
    dim=spice_bundle_model%bundle%L(domain)%dim
    n_modes=dim
        
    ALLOCATE( real_TII(1:dim,1:dim) )
    ALLOCATE( real_TV(1:dim,1:dim) )
    ALLOCATE( real_TI(1:dim,1:dim) )
    ALLOCATE( real_TVI(1:dim,1:dim) )
    ALLOCATE( mode_velocity(1:dim) )
    ALLOCATE( mode_impedance(1:dim) )
  
  
    CALL write_spice_comment('Modal Decomposition')
    domain_mode_decomp(domain)%n_modes=dim
    ALLOCATE( domain_mode_decomp(domain)%TVI(1:dim,1:dim) )
    ALLOCATE( domain_mode_decomp(domain)%TI(1:dim,1:dim) )
    ALLOCATE( domain_mode_decomp(domain)%delay(1:dim) )
    ALLOCATE( domain_mode_decomp(domain)%impedance(1:dim) )

! OPTION NEVER IMPLEMENTED...
!    if (use_fasthenry.AND.use_lossy_hf_modes) then 
    
! Use the 'high frequency' L and C matrices to derive the modal decomposition for the spice model

      CALL modal_decomposition_LC(dim,spice_bundle_model%bundle%L(domain)%mat, &
                                      spice_bundle_model%bundle%C(domain)%mat, &
                                      real_TI,real_TII,real_TV,real_TVI,mode_velocity,mode_impedance)
  
! save the mode information on the idealised transmission line    
    
      domain_mode_decomp(domain)%TVI(1:dim,1:dim)=real_TVI(1:dim,1:dim)
      domain_mode_decomp(domain)%TI(1:dim,1:dim)=real_TI(1:dim,1:dim)
        
      do mode=1,n_modes
        domain_mode_decomp(domain)%delay(mode)=spice_bundle_model%length/mode_velocity(mode)
        domain_mode_decomp(domain)%impedance(mode)=mode_impedance(mode)
      end do
      
!    else
!      
!! Use the 'high frequency' L and C matrices to derive the modal decomposition for the spice model
!
!      CALL modal_decomposition_LC(dim,spice_bundle_model%bundle%L(domain)%mat, &
!                                      spice_bundle_model%bundle%C(domain)%mat, &
!                                      real_TI,real_TII,real_TV,real_TVI,mode_velocity,mode_impedance)
!  
!! save the mode information on the idealised transmission line    
!    
!      domain_mode_decomp(domain)%TVI(1:dim,1:dim)=real_TVI(1:dim,1:dim)
!      domain_mode_decomp(domain)%TI(1:dim,1:dim)=real_TI(1:dim,1:dim)
!        
!      do mode=1,n_modes
!        domain_mode_decomp(domain)%delay(mode)=spice_bundle_model%length/mode_velocity(mode)
!        domain_mode_decomp(domain)%impedance(mode)=mode_impedance(mode)
!      end do
!    
!    end if
     
! STAGE7:  work out the propagation correction filters for this domain
! See Theory_Manual_Section 3.6

! Allocate the conductor based loss model data
    
    dim=spice_bundle_model%bundle%n_conductors(domain)
    ALLOCATE( conductor_impedance_domain(1:dim+1) )

! Allocate the propagation correction filter
    dim=n_modes
    domain_HP_list(domain)%n_modes=n_modes
    ALLOCATE( domain_HP_list(domain)%Hp(1:n_modes) )  

! Retreive the conductor based impedance model data    
    do conductor=1,spice_bundle_model%bundle%n_conductors(domain)

      terminal_conductor=spice_bundle_model%bundle%terminal_conductor_list(domain)%element(conductor)
      
! copy the conductor impedance model for this cable conductor to the bundle structure        
        
      conductor_impedance_domain(conductor)%impedance_model_type=  &
                  spice_bundle_model%bundle%conductor_impedance(terminal_conductor)%impedance_model_type
      
      conductor_impedance_domain(conductor)%radius=  &
                  spice_bundle_model%bundle%conductor_impedance(terminal_conductor)%radius
      
      conductor_impedance_domain(conductor)%thickness=  &
                  spice_bundle_model%bundle%conductor_impedance(terminal_conductor)%thickness
                  
      conductor_impedance_domain(conductor)%width=  &
                  spice_bundle_model%bundle%conductor_impedance(terminal_conductor)%width
                  
      conductor_impedance_domain(conductor)%height=  &
                  spice_bundle_model%bundle%conductor_impedance(terminal_conductor)%height
                  
      conductor_impedance_domain(conductor)%conductivity=  &
                  spice_bundle_model%bundle%conductor_impedance(terminal_conductor)%conductivity
            
      conductor_impedance_domain(conductor)%Resistance_multiplication_factor=  &
                  spice_bundle_model%bundle%conductor_impedance(terminal_conductor)%Resistance_multiplication_factor
                  
      if ( (conductor_impedance_domain(conductor)%impedance_model_type.EQ.impedance_model_type_filter).OR. &
           (conductor_impedance_domain(conductor)%impedance_model_type.EQ.impedance_model_type_cylindrical_shield) ) then
                              
        conductor_impedance_domain(conductor)%ZT_filter=   &
               spice_bundle_model%bundle%conductor_impedance(terminal_conductor)%ZT_filter
               
      end if
         
    end do

    if (spice_bundle_model%Fit_model_order.NE.0) then
    
      CALL calculate_propagation_correction_filters(domain,dim,spice_bundle_model%FIT_model_order,          &
           spice_bundle_model%bundle%Z(domain)%Sfilter_mat,spice_bundle_model%bundle%Y(domain)%Sfilter_mat, &
           spice_bundle_model%bundle%L(domain)%mat,spice_bundle_model%bundle%C(domain)%mat,                 &
           domain_mode_decomp(domain)%TI(1:dim,1:dim),                                                      &
           conductor_impedance_domain,domain_HP_list(domain)%Hp,spice_bundle_model%length,                  &
           domain_mode_decomp(domain)%delay,                                                                &
           spice_bundle_model%prop_corr_fit_freq_spec)
           
     else
! the filter order is 0 and then therefore the propagation correction filter is equal to 1

       do mode=1,n_modes
         domain_HP_list(domain)%Hp(mode)=1d0
        end do
        
     end if
     
! STAGE 8: write the modal decomposition circuit elements at both ends of the transmission line
! The conductor based node numbers are known and the mode based node numbers are created

! Allocate the node number lists for both ends of the transmission line within this domain

    ALLOCATE( domain_end1_nodes(1:n_conductors) )
    ALLOCATE( domain_end2_nodes(1:n_conductors) )

! The domain based node list within a single domain comes from a subset of the
! domain based transmission lines for all domains together

    do local_domain_conductor=1,n_conductors-1

      global_domain_conductor=global_domain_conductor+1
      
      domain_end1_nodes(local_domain_conductor)=all_domain_end1_nodes(global_domain_conductor)
      domain_end2_nodes(local_domain_conductor)=all_domain_end2_nodes(global_domain_conductor)
      
    end do
    domain_end1_nodes(n_conductors)=vref_node
    domain_end2_nodes(n_conductors)=vref_node

! Allocate the modal node number lists for both ends of the transmission line for this domain
    ALLOCATE( modal_end1_nodes(1:n_conductors) )
    ALLOCATE( modal_end2_nodes(1:n_conductors) )

    end=1

    CALL write_spice_modal_decomposition_equivalent_circuit(vref_node,next_free_node,domain,end,n_conductors,dim,   &
                                                            domain_end1_nodes,real_TV,real_TII,modal_end1_nodes)

    end=2
    CALL write_spice_modal_decomposition_equivalent_circuit(vref_node,next_free_node,domain,end,n_conductors,dim,   &
                                                            domain_end2_nodes,real_TV,real_TII,modal_end2_nodes)

! allocate the method of characteristics reference node lists 
    
    ALLOCATE( domain_MOC_V_plus_nodes(domain)%node(1:n_modes) )
    ALLOCATE( domain_MOC_V_minus_nodes(domain)%node(1:n_modes) )
    
    ALLOCATE( domain_MOC_reference_end1_nodes(domain)%node(1:n_modes) )
    ALLOCATE( domain_MOC_reference_end2_nodes(domain)%node(1:n_modes) )
    
    ALLOCATE( domain_MOC_ZT_reference_end1_nodes(domain)%node(1:n_modes) )
    ALLOCATE( domain_MOC_ZT_reference_end2_nodes(domain)%node(1:n_modes) )
    
    ALLOCATE( domain_ZT_internal_end1_nodes(domain)%node(1:n_modes) )
    ALLOCATE( domain_ZT_internal_end2_nodes(domain)%node(1:n_modes) )
    
    ALLOCATE( domain_MOC_EINC_ZT_reference_end1_nodes(domain)%node(1:n_modes) )
    ALLOCATE( domain_MOC_EINC_ZT_reference_end2_nodes(domain)%node(1:n_modes) )

! Allow for the presence of additional voltage sources related to incident field coupling 
! and transfer impedance coupling sources as seen in Theory_Manual_Figure 3.2

! Set nodes to the overall reference node initially i.e. assume no transfer impedance or incident field models
    domain_MOC_reference_end1_nodes(domain)%node(1:n_modes)=vref_node
    domain_MOC_reference_end2_nodes(domain)%node(1:n_modes)=vref_node
    domain_MOC_ZT_reference_end1_nodes(domain)%node(1:n_modes)=vref_node
    domain_MOC_ZT_reference_end2_nodes(domain)%node(1:n_modes)=vref_node
    domain_MOC_EINC_ZT_reference_end1_nodes(domain)%node(1:n_modes)=vref_node
    domain_MOC_EINC_ZT_reference_end2_nodes(domain)%node(1:n_modes)=vref_node

! The order of sources is: MOC sources, ZT sources, Einc_ZT sources, Einc sources
    
! create nodes for the incident field excitation coupled via transfer impedance equivalent voltage sources 
! in the method of characteristics sub-circuit for the external domain if required
    
! create nodes for the incident field excitation equivalent voltage sources in the method of characteristics 
! sub-circuit for the external domain if required

    if (domain_is_Einc_victim(domain)) then
          
      do mode=1,n_modes
        domain_MOC_EINC_ZT_reference_end1_nodes(domain)%node(mode)=next_free_node     ! set to the next free node but don't increase next_free_node
        domain_MOC_ZT_reference_end1_nodes(domain)%node(mode)=next_free_node          ! set to the next free node but don't increase next_free_node
        CALL create_new_node(domain_MOC_reference_end1_nodes(domain)%node(mode),next_free_node)
        
        domain_MOC_EINC_ZT_reference_end2_nodes(domain)%node(mode)=next_free_node     ! set to the next free node but don't increase next_free_node
        domain_MOC_ZT_reference_end2_nodes(domain)%node(mode)=next_free_node          ! set to the next free node but don't increase next_free_node
        CALL create_new_node(domain_MOC_reference_end2_nodes(domain)%node(mode),next_free_node)
        
      end do
             
    end if ! this is the external domain and there is an incident field excitation

! create nodes for the incident field excitation coupling via ZT equivalent voltage sources 
! in the method of characteristics sub-circuit for the external domain if required

    if (domain_is_Einc_ZT_victim(domain)) then
          
      do mode=1,n_modes
        domain_MOC_ZT_reference_end1_nodes(domain)%node(mode)=next_free_node      ! set to the next free node but don't increase next_free_node
        CALL create_new_node(domain_MOC_reference_end1_nodes(domain)%node(mode),next_free_node)
        
        domain_MOC_ZT_reference_end2_nodes(domain)%node(mode)=next_free_node      ! set to the next free node but don't increase next_free_node
               
        CALL create_new_node(domain_MOC_reference_end2_nodes(domain)%node(mode),next_free_node)

      end do
             
    end if ! this domain is the victim of the incident field excitation via transfer impedance coupling (Xie's model)
     
! create nodes for the transfer impedance equivalent voltage sources in the method of characteristics 
! sub-circuit for the victim domain if required

    if ( domain_is_ZT_victim(domain) ) then
          
      do mode=1,n_modes
        
        CALL create_new_node(domain_MOC_reference_end1_nodes(domain)%node(mode),next_free_node)
        CALL create_new_node(domain_MOC_reference_end2_nodes(domain)%node(mode),next_free_node)
        
      end do
       
    end if ! this is a victim domain for a transfer impedance model
    
! loop over the modes in this domain and write the spice circuit elements for 
! the modal propagation. See Theory_Manual_Section 3.5

    do mode=1,n_modes
      
      CALL write_spice_method_of_characteristics_equivalent_circuit(modal_end1_nodes(mode),modal_end2_nodes(mode),vref_node, &
                                                                    domain_MOC_reference_end1_nodes(domain)%node(mode), &
                                                                    domain_MOC_reference_end2_nodes(domain)%node(mode), &
                                                                    domain_MOC_V_plus_nodes(domain)%node(mode), &
                                                                    domain_MOC_V_minus_nodes(domain)%node(mode), &
                                                                    next_free_node,domain,mode, &
                                                                    domain_mode_decomp(domain)%impedance(mode), &
                                                                    domain_mode_decomp(domain)%delay(mode),      &
                                                                    domain_HP_list(domain)%Hp(mode),spice_bundle_model%length)

    end do ! next mode in this domain

! deallocate modal decomposition stuff for this domain
   
    DEALLOCATE( real_TII )
    DEALLOCATE( real_TV )
    DEALLOCATE( real_TI )
    DEALLOCATE( real_TVI )
    DEALLOCATE( mode_velocity )
    DEALLOCATE( mode_impedance )
    
    do conductor=1,spice_bundle_model%bundle%n_conductors(domain)+1
      CALL deallocate_conductor_impedance_model(conductor_impedance_domain(conductor))
    end do
    DEALLOCATE( conductor_impedance_domain )
    
    DEALLOCATE( modal_end1_nodes )
    DEALLOCATE( modal_end2_nodes )

    DEALLOCATE( domain_end1_nodes )
    DEALLOCATE( domain_end2_nodes )

  end do ! next domain

! STAGE 9: Transfer impedance coupling model
! See Theory_Manual_Section 3.7

  count_ZT_victim_in_domain(1:n_domains)=0
 
  do ZT_model=1,spice_bundle_model%n_transfer_impedances

    source_domain=ZT_source_domain(ZT_model)
    victim_domain=ZT_victim_domain(ZT_model)
 
    source_domain=ZT_source_domain(ZT_model)
    n_source_domain_modes=spice_bundle_model%bundle%n_conductors(source_domain)-1
    
    victim_domain=ZT_victim_domain(ZT_model)
    n_victim_domain_modes=spice_bundle_model%bundle%n_conductors(victim_domain)-1
    
    if (verbose) then
    
      write(*,*)'Source domain:',source_domain
      write(*,*)'Number of source domain modes:',n_source_domain_modes
      
      write(*,*)'    Mode delay(s)           Mode Impedance (ohms)'
      do mode=1,n_source_domain_modes
        write(*,*)domain_mode_decomp(source_domain)%delay(mode),domain_mode_decomp(source_domain)%impedance(mode)
      end do
      
      write(*,*)'Source domain:',Victim_domain
      write(*,*)'Number of victim domain modes:',n_victim_domain_modes

      write(*,*)'    Mode delay(s)           Mode Impedance (ohms)'
      do mode=1,n_victim_domain_modes
        write(*,*)domain_mode_decomp(victim_domain)%delay(mode),domain_mode_decomp(victim_domain)%impedance(mode)
      end do
    
    end if
    
! get the information required to generate the transfer impedance model
     
    shield_conductor=spice_bundle_model%Zt_conductor(ZT_model)  ! this is the conductor number for the shield whose transfer impedance is being included

! get the local conductor number (i.e. within the domain) in the source and victim domains

    terminal_shield_conductor=spice_bundle_model%Zt_conductor(ZT_model)
    domain_shield_conductor=spice_bundle_model%bundle%terminal_conductor_to_local_domain_conductor(terminal_shield_conductor)
    
    if (spice_bundle_model%Zt_direction(ZT_model).EQ.1) then
! coupling from inside the shield to outside
    
      source_domain_shield_conductor=spice_bundle_model%bundle%n_conductors(source_domain)
      victim_domain_shield_conductor=domain_shield_conductor
      
    else
! coupling from inside the shield to outside
    
      source_domain_shield_conductor=domain_shield_conductor
      victim_domain_shield_conductor=spice_bundle_model%bundle%n_conductors(victim_domain)
    
    end if

! If this is not the final transfer impedance model for this victim domain then we
! must create a new node for the refrence voltage    
    count_ZT_victim_in_domain(victim_domain)=count_ZT_victim_in_domain(victim_domain)+1
    
    ALLOCATE( ZT_node1_end1%node(1:n_victim_domain_modes) )
    ALLOCATE( ZT_node1_end2%node(1:n_victim_domain_modes) )
    ALLOCATE( ZT_node2_end1%node(1:n_victim_domain_modes) )
    ALLOCATE( ZT_node2_end2%node(1:n_victim_domain_modes) )
    
    if (count_ZT_victim_in_domain(victim_domain).EQ.1) then
! we use the already created reference node for MOC sources
      ZT_node1_end1%node(1:n_victim_domain_modes)=domain_MOC_reference_end1_nodes(victim_domain)%node(1:n_victim_domain_modes)
      ZT_node1_end2%node(1:n_victim_domain_modes)=domain_MOC_reference_end2_nodes(victim_domain)%node(1:n_victim_domain_modes)
    else
! we use the previously created internal node
      ZT_node1_end1%node(1:n_victim_domain_modes)=domain_ZT_internal_end1_nodes(victim_domain)%node(1:n_victim_domain_modes)
      ZT_node1_end2%node(1:n_victim_domain_modes)=domain_ZT_internal_end2_nodes(victim_domain)%node(1:n_victim_domain_modes)
    end if
       
    if (count_ZT_victim_in_domain(victim_domain).NE.n_ZT_victim_in_domain(victim_domain)) then
! we create a new internal node 
      do mode=1,n_victim_domain_modes
        CALL create_new_node(ZT_node2_end1%node(mode),next_free_node)
        CALL create_new_node(ZT_node2_end2%node(mode),next_free_node)
      end do
! save the internal nodes 
      domain_ZT_internal_end1_nodes(victim_domain)%node(1:n_victim_domain_modes)=ZT_node2_end1%node(1:n_victim_domain_modes)
      domain_ZT_internal_end2_nodes(victim_domain)%node(1:n_victim_domain_modes)=ZT_node2_end2%node(1:n_victim_domain_modes)
    else
! we use the already created reference node for transfer impedance sources
      ZT_node2_end1%node(1:n_victim_domain_modes)=domain_MOC_ZT_reference_end1_nodes(victim_domain)%node(1:n_victim_domain_modes)
      ZT_node2_end2%node(1:n_victim_domain_modes)=domain_MOC_ZT_reference_end2_nodes(victim_domain)%node(1:n_victim_domain_modes)
    end if

    CALL write_transfer_impedance_circuit(n_source_domain_modes,n_victim_domain_modes, &
                                          ZT_node1_end1%node,ZT_node1_end2%node, &
                                          ZT_node2_end1%node,ZT_node2_end2%node, &
                                          domain_MOC_V_minus_nodes(source_domain)%node, &
                                          domain_MOC_V_plus_nodes(source_domain)%node, &
                                          domain_mode_decomp(source_domain)%delay, &
                                          domain_mode_decomp(source_domain)%impedance, &
                                          domain_mode_decomp(victim_domain)%delay, &
                                          source_domain_shield_conductor, &
                                          victim_domain_shield_conductor, &
                                          spice_bundle_model%length, &
                                          domain_mode_decomp(source_domain)%TVI,domain_mode_decomp(source_domain)%TI, &
                                          domain_mode_decomp(victim_domain)%TVI,domain_mode_decomp(victim_domain)%TI, &
                                          domain_HP_list(victim_domain)%Hp, &
                                          spice_bundle_model%bundle%conductor_impedance(shield_conductor)%ZT_filter, &
                                          next_free_node,vref_node,ZT_model,source_domain,victim_domain )

    
    DEALLOCATE( ZT_node1_end1%node )
    DEALLOCATE( ZT_node1_end2%node )
    DEALLOCATE( ZT_node2_end1%node )
    DEALLOCATE( ZT_node2_end2%node )
    
  end do ! next transfer impedance model

! STAGE 10: generate the information required for and then write the incident field excitation model as required
! See Theory_Manual_Sections 3.8, 3.9

  if (spice_bundle_model%include_incident_field) then

    domain=n_domains
    n_einc_domain_modes=spice_bundle_model%bundle%n_conductors(domain)-1
    
    write(*,*)'n_einc_domain_modes=',n_einc_domain_modes
    write(*,*)'Domain',domain
    
! calculate the coefficients required for the incident field excitation model (Paul, 1st edition equations 7.241a,b)

! assemble the data required for the incident field calculation.

! get the central coordinates of all the conductors in the external domain.
    ALLOCATE( xcoord(1:n_einc_domain_modes+1) )
    ALLOCATE( ycoord(1:n_einc_domain_modes+1) )

    write(*,*)'Incident field excitation, number of modes in external domain=',n_einc_domain_modes

! get the coordinates of the conductors in the external domain
    conductor=0
    do i=1,tot_n_conductors
    
      write(*,*)i,'outer_domain',spice_bundle_model%bundle%terminal_conductor_to_outer_domain(i),&
                  'inner_domain',spice_bundle_model%bundle%terminal_conductor_to_inner_domain(i)
                  
    
      if (spice_bundle_model%bundle%terminal_conductor_to_outer_domain(i).EQ.domain) then
! this conductor is in the outer domain so must be included
    
        conductor=conductor+1
        xcoord(conductor)=spice_bundle_model%bundle%conductor_x_offset(i)
        ycoord(conductor)=spice_bundle_model%bundle%conductor_y_offset(i)
 
       end if
      
    end do ! next conductor
    
    ALLOCATE( alpha0(1:n_einc_domain_modes) )
    ALLOCATE( alphaL(1:n_einc_domain_modes) )
    ALLOCATE( c1(1:n_einc_domain_modes) )
    ALLOCATE( c2(1:n_einc_domain_modes) )
    ALLOCATE( c1b(1:n_einc_domain_modes) )

! See Theory_Manual_Section 3.8
    CALL calculate_incident_field_sources(xcoord,ycoord,spice_bundle_model%Eamplitude,                       &
                                          spice_bundle_model%Ex,spice_bundle_model%Ey,spice_bundle_model%Ez, &
                                          spice_bundle_model%Hx,spice_bundle_model%Hy,spice_bundle_model%Hz, &
                                          spice_bundle_model%kx,spice_bundle_model%ky,spice_bundle_model%kz, &
                                          spice_bundle_model%bundle%ground_plane_present,                    &
                                          spice_bundle_model%bundle%ground_plane_x,                          &
                                          spice_bundle_model%bundle%ground_plane_y,                          &
                                          spice_bundle_model%bundle%ground_plane_theta,                      &
                                          spice_bundle_model%length,n_einc_domain_modes,                     &
                                          domain_mode_decomp(domain)%TI,domain_mode_decomp(domain)%delay,Tz,alpha0,alphaL)

    CALL write_incident_field_circuit(n_einc_domain_modes, &
                                      domain_MOC_EINC_ZT_reference_end1_nodes(domain)%node, &
                                      domain_MOC_EINC_ZT_reference_end2_nodes(domain)%node, &
                                      domain_mode_decomp(domain)%delay, &
                                      spice_bundle_model%length, &
                                      domain_mode_decomp(domain)%TVI,domain_mode_decomp(domain)%TI, &
                                      domain_HP_list(domain)%Hp, &
                                      Tz,alpha0,alphaL, &
                                      next_free_node,Einc_node1,Einc_node2,vref_node,domain )

! look for any victim domains for incident field coupling via transfer impedance                                 
! See Theory_Manual_Section 3.9
    
    do ZT_model=1,spice_bundle_model%n_transfer_impedances

      source_domain=ZT_source_domain(ZT_model)
      victim_domain=ZT_victim_domain(ZT_model)
         
! get the information required to generate the transfer impedance model

! Check whether the source domain for this transfer impedance model is the external domain
   
      if (     (source_domain.EQ.n_domains)              &
          .AND.(domain_is_Einc_ZT_victim(victim_domain))   ) then

! we must include additional terms within shielded cables for this victim domain
         
        shield_conductor=spice_bundle_model%Zt_conductor(ZT_model)  ! this is the conductor number for the shield whose transfer impedance is being included
        
        source_domain_shield_conductor=spice_bundle_model%bundle%terminal_conductor_to_local_domain_conductor(shield_conductor)
        
        n_source_domain_modes=spice_bundle_model%bundle%n_conductors(source_domain)-1
    
        n_victim_domain_modes=spice_bundle_model%bundle%n_conductors(victim_domain)-1
        
        victim_domain_shield_conductor=spice_bundle_model%bundle%n_conductors(victim_domain) ! always the reference conductor
       
        if (verbose) then
          write(*,*)'******Calculate ZT_incident field sources******'
          write(*,*)'Terminal conductor number of shield=',shield_conductor
          write(*,*)'Source domain=',source_domain
          write(*,*)'n_conductors in source domain',spice_bundle_model%bundle%n_conductors(domain)
          write(*,*)'n_source domain modes=',n_source_domain_modes
          write(*,*)'Source domain shield conductor=',domain_shield_conductor
          write(*,*)'n_einc_domain_modes (source domain modes)=',n_einc_domain_modes
          write(*,*)'Victim domain=',victim_domain
          write(*,*)'n_victim domain modes=',n_victim_domain_modes
          write(*,*)'Victim domain shield conductor=',victim_domain_shield_conductor
        end if
        
! See Theory_Manual_Section 3.9
        CALL calculate_ZT_incident_field_sources(xcoord,ycoord,spice_bundle_model%Eamplitude,                &
                                          spice_bundle_model%Ex,spice_bundle_model%Ey,spice_bundle_model%Ez, &
                                          spice_bundle_model%Hx,spice_bundle_model%Hy,spice_bundle_model%Hz, &
                                          spice_bundle_model%kx,spice_bundle_model%ky,spice_bundle_model%kz, &
                                          spice_bundle_model%bundle%ground_plane_present,                    &
                                          spice_bundle_model%bundle%ground_plane_x,                          &
                                          spice_bundle_model%bundle%ground_plane_y,                          &
                                          spice_bundle_model%bundle%ground_plane_theta,                      &
                                          spice_bundle_model%length,n_einc_domain_modes,                     &
                                          domain_shield_conductor,                                           &
                                          domain_mode_decomp(domain)%TI,domain_mode_decomp(domain)%delay,Tz,C1,C2,C1b)
       
        CALL write_ZT_incident_field_circuit(ZT_model,source_domain,victim_domain,n_source_domain_modes,n_victim_domain_modes, &
                                             domain_MOC_ZT_reference_end1_nodes(victim_domain)%node, &
                                             domain_MOC_ZT_reference_end2_nodes(victim_domain)%node, &
                                             domain_MOC_EINC_ZT_reference_end1_nodes(victim_domain)%node, &
                                             domain_MOC_EINC_ZT_reference_end2_nodes(victim_domain)%node, &
                                             Einc_node1,Einc_node2, &
                                             domain_MOC_V_minus_nodes(source_domain)%node, &
                                             domain_MOC_V_plus_nodes(source_domain)%node, &
                                             domain_mode_decomp(source_domain)%TVI,domain_mode_decomp(source_domain)%TI, &
                                             domain_mode_decomp(source_domain)%delay, &
                                             domain_mode_decomp(source_domain)%impedance, &
                                             domain_mode_decomp(victim_domain)%TVI,domain_mode_decomp(victim_domain)%TI, &
                                             source_domain_shield_conductor, &
                                             victim_domain_shield_conductor, &
                                             domain_mode_decomp(victim_domain)%delay, &
                                             domain_mode_decomp(victim_domain)%impedance, &
                                             spice_bundle_model%bundle%conductor_impedance(shield_conductor)%ZT_filter, &
                                             domain_HP_list(victim_domain)%Hp, &
                                             Tz, &
                                             C1,C2,C1b, &
                                             spice_bundle_model%length, &
                                             next_free_node,vref_node )      
    
      end if ! is this a victim domain for incident field coupling via transfer impedance   
    
    end do ! next transfer impedance model to check
                                      
    DEALLOCATE( xcoord )
    DEALLOCATE( ycoord )
    DEALLOCATE( alpha0 )
    DEALLOCATE( alphaL )
    DEALLOCATE( c1 )
    DEALLOCATE( c2 )
    DEALLOCATE( c1b )

  end if ! incident field excitation model is to be included
  
! End subcircuit definition
  write(spice_model_file_unit,'(A)')'*'
  write(spice_model_file_unit,'(A)')'.ends'
  write(spice_model_file_unit,'(A)')'*'
    
! STAGE 11: Dallocate the domain based information
  do domain=1,n_domains
    DEALLOCATE( domain_mode_decomp(domain)%TVI )
    DEALLOCATE( domain_mode_decomp(domain)%TI )
    DEALLOCATE( domain_mode_decomp(domain)%delay )
    DEALLOCATE( domain_mode_decomp(domain)%impedance )
  end do
  DEALLOCATE( domain_mode_decomp )
  
  DEALLOCATE( ZT_source_domain )
  DEALLOCATE( ZT_victim_domain )
  
  DEALLOCATE( domain_is_ZT_victim )
  DEALLOCATE( domain_is_Einc_victim )
  DEALLOCATE( domain_is_Einc_ZT_victim )
  
  DEALLOCATE( n_ZT_victim_in_domain )
  DEALLOCATE( count_ZT_victim_in_domain )
  
  do domain=1,n_domains
    
    DEALLOCATE( domain_MOC_V_plus_nodes(domain)%node )
    DEALLOCATE( domain_MOC_V_minus_nodes(domain)%node )
  
    DEALLOCATE( domain_MOC_reference_end1_nodes(domain)%node )
    DEALLOCATE( domain_MOC_reference_end2_nodes(domain)%node )
    
    DEALLOCATE( domain_MOC_ZT_reference_end1_nodes(domain)%node )
    DEALLOCATE( domain_MOC_ZT_reference_end2_nodes(domain)%node )
    
    if ( ALLOCATED(domain_ZT_internal_end1_nodes(domain)%node) ) DEALLOCATE( domain_ZT_internal_end1_nodes(domain)%node )
    if ( ALLOCATED(domain_ZT_internal_end2_nodes(domain)%node) ) DEALLOCATE( domain_ZT_internal_end2_nodes(domain)%node )

    DEALLOCATE( domain_MOC_EINC_ZT_reference_end1_nodes(domain)%node )
    DEALLOCATE( domain_MOC_EINC_ZT_reference_end2_nodes(domain)%node )
    
    do mode=1,domain_HP_list(domain)%n_modes
      CALL deallocate_Sfilter(domain_HP_list(domain)%Hp(mode))
    end do
    DEALLOCATE( domain_HP_list(domain)%Hp )
    
  end do
    
  DEALLOCATE( domain_MOC_V_plus_nodes )
  DEALLOCATE( domain_MOC_V_minus_nodes )
  
  DEALLOCATE( domain_MOC_reference_end1_nodes )
  DEALLOCATE( domain_MOC_reference_end2_nodes )
  
  DEALLOCATE( domain_MOC_ZT_reference_end1_nodes )
  DEALLOCATE( domain_MOC_ZT_reference_end2_nodes )
  
  DEALLOCATE( domain_ZT_internal_end1_nodes )
  DEALLOCATE( domain_ZT_internal_end2_nodes )
  
  DEALLOCATE( domain_MOC_EINC_ZT_reference_end1_nodes )
  DEALLOCATE( domain_MOC_EINC_ZT_reference_end2_nodes )
  DEALLOCATE( domain_HP_list )
  
! Deallocate the external node number lists for both ends of the transmission line

  DEALLOCATE( external_end1_nodes )
  DEALLOCATE( external_end2_nodes )

  DEALLOCATE( rdc_end1_nodes )
  DEALLOCATE( rdc_end2_nodes )

  DEALLOCATE( all_domain_end1_nodes )
  DEALLOCATE( all_domain_end2_nodes )
    
  DEALLOCATE( Rdc )

  CLOSE(unit=spice_model_file_unit)

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

  RETURN

END SUBROUTINE create_spice_subcircuit_model