!
! This file is part of SACAMOS, State of the Art CAble MOdels in Spice. 
! It was developed by the University of Nottingham and the Netherlands Aerospace 
! Centre (NLR) for ESA under contract number 4000112765/14/NL/HK.
! 
! Copyright (C) 2016-2017 University of Nottingham
! 
! SACAMOS is free software: you can redistribute it and/or modify it under the 
! terms of the GNU General Public License as published by the Free Software 
! Foundation, either version 3 of the License, or (at your option) any later 
! version.
! 
! SACAMOS is distributed in the hope that it will be useful, but 
! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 
! or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
! for more details.
! 
! A copy of the GNU General Public License version 3 can be found in the 
! file GNU_GPL_v3 in the root or at <http://www.gnu.org/licenses/>.
! 
! SACAMOS uses the EISPACK library (in /SRC/EISPACK). EISPACK is subject to 
! the GNU Lesser General Public License. A copy of the GNU Lesser General Public 
! License version can be found in the file GNU_LGPL in the root of EISPACK 
! (/SRC/EISPACK ) or at <http://www.gnu.org/licenses/>.
! 
! The University of Nottingham can be contacted at: ggiemr@nottingham.ac.uk
!
! File Contents:
! SUBROUTINE twisted_pair_set_parameters
! SUBROUTINE twisted_pair_set_internal_domain_information
!
! NAME
!     twisted_pair_set_parameters
!
! AUTHORS
!     Chris Smartt
!
! DESCRIPTION
!     Set the overall parameters for a twisted_pair cable
!
! COMMENTS
!      
!
! HISTORY
!
!     started 12/4/201 CJS 
!
!
SUBROUTINE twisted_pair_set_parameters(cable)

USE type_specifications

IMPLICIT NONE

! variables passed to subroutine

  type(cable_specification_type),intent(INOUT)    :: cable

! local variables

! START

  cable%cable_type=cable_geometry_type_twisted_pair
  cable%tot_n_conductors=2
  cable%tot_n_domains=2
  cable%n_external_conductors=1
  cable%n_internal_conductors=1
  cable%n_internal_domains=1
  cable%n_parameters=4
  cable%n_dielectric_filters=1
  cable%n_transfer_impedance_models=0

END SUBROUTINE twisted_pair_set_parameters
!
! NAME
!     twisted_pair_set_internal_domain_information
!
! AUTHORS
!     Chris Smartt
!
! DESCRIPTION
!     Set the overall parameters for a twisted_pair cable
!
! COMMENTS
!     Set the dimension of the domain transformation matrices to include an external reference conductor for the cable 
!     We need to set a radius for the 'common mode equivalent conductor'. Set to 1.5* wire radius for now...
!
! HISTORY
!
!     started 12/4/201 CJS 
!     8/9/2016 CJS common mode/ differential mode loss correction
!     19/9/2016 CJS frequency dependent dielectric in Laplace solver
!     8/5/2017         CJS: Include references to Theory_Manual
!
!
SUBROUTINE twisted_pair_set_internal_domain_information(cable)

USE type_specifications
USE constants
USE general_module
USE PUL_parameter_module

IMPLICIT NONE

! variables passed to subroutine

  type(cable_specification_type),intent(INOUT)    :: cable

! local variables

  integer :: n_conductors
  integer :: dim

  type(PUL_type)    :: PUL

  integer :: domain
  
! variables for cable parameter checks 
  logical :: cable_spec_error
  real(dp) :: rw
  real(dp) :: s
  real(dp) :: rd
  real(dp) :: sigma
 
  type(Sfilter) :: epsr
  
  character(LEN=error_message_length) :: message 

! START
    
  write(*,*)'CALLED   twisted_pair_set_internal_domain_information'

! Check the cable parameters

  rw=cable%parameters(1) 
  s=cable%parameters(2) 
  rd=cable%parameters(3) 
  sigma=cable%parameters(4) 
  
  epsr=cable%dielectric_filter(1)
  
  cable_spec_error=.FALSE.    ! assume no errors initially
  message=''
  CALL twisted_pair_check(rw,rd,s,cable_spec_error,cable%cable_name,message)
  CALL conductivity_check(sigma,cable_spec_error,cable%cable_name,message)
  CALL dielectric_check(epsr,cable_spec_error,cable%cable_name,message)
   
  if (cable_spec_error) then       
    run_status='ERROR in cable_model_builder, error on parameters for cable:'//trim(cable%cable_name)//'. '//trim(message)
    CALL write_program_status()
    STOP 1
  end if
    
! The differential mode is treated as an internal domain which doesn't couple to other conductors
  domain=1
    
  cable%n_internal_conductors_in_domain(domain)=2
    
! The number of modes in the internal domain is 1
  dim=1
  cable%L_domain(domain)%dim=dim
  ALLOCATE(cable%L_domain(domain)%mat(dim,dim))
  cable%C_domain(domain)%dim=dim
  ALLOCATE(cable%C_domain(domain)%mat(dim,dim))

! the L and C matrix elements are for the differential mode
! this is calculated as if the conductors were in free space with no other conductors

  if (use_laplace) then
  
! allocate memory for the PUL parameter solver interface

    if(verbose) write(*,*)'Domain:',domain
    if(verbose) write(*,*)'Allocating PUL data structure for shielded twisted pairs'
    n_conductors=2
    
    CALL allocate_and_reset_PUL_data(PUL,n_conductors)  
    
    PUL%shape(1:n_conductors)=circle
    
    PUL%x(1)=-s/2d0
    PUL%y(1)=0.0
    PUL%r(1)=rw
    PUL%rd(1)=rd
    PUL%epsr(1)=epsr 
    
    PUL%x(2)=s/2d0
    PUL%y(2)=0.0
    PUL%r(2)=rw
    PUL%rd(2)=rd
    PUL%epsr(2)=epsr

    PUL%epsr_background = 1d0  ! permittivity of homogeneous dielectric medium surrounding conductors (air)
  
! no ground plane
    PUL%ground_plane_present=.FALSE.
      
! no overshield
    PUL%overshield_present=.FALSE.
    
    CALL PUL_LC_Laplace(PUL,cable%cable_name,cable%Y_fit_model_order,cable%Y_fit_freq_spec,domain) 
    
    cable%L_domain(domain)%mat(1,1)=PUL%L%mat(1,1)
    cable%C_domain(domain)%mat(1,1)=PUL%C%mat(1,1)
    
    cable%Z_domain(domain)%dim=dim
    ALLOCATE(cable%Z_domain(domain)%sfilter_mat(dim,dim))
    cable%Y_domain(domain)%dim=dim
    ALLOCATE(cable%Y_domain(domain)%sfilter_mat(dim,dim))
    
    cable%Z_domain(domain)%sfilter_mat(:,:)=PUL%Zfilter%sfilter_mat(:,:)
    cable%Y_domain(domain)%sfilter_mat(:,:)=PUL%Yfilter%sfilter_mat(:,:)
    
  else

! See C.R. Paul, 1st edition, equation 3.46 ! Theory_Manual_Eqn 2.27, 2.28

    cable%C_domain(domain)%mat(1,1)=pi*eps0/log( s/(2d0*rw)+sqrt( (s/(2d0*rw))**2 -1) )
    cable%L_domain(domain)%mat(1,1)=(mu0*eps0)/cable%C_domain(domain)%mat(1,1)
    
    CALL Z_Y_from_L_C(cable%L_domain(domain),cable%C_domain(domain),cable%Z_domain(domain),cable%Y_domain(domain))
    
  end if
  
  if (use_laplace) CALL deallocate_PUL_data(PUL)  ! deallocate the PUL data structure
    
! Set the domain decomposition matrices ! Theory_Manual_Eqn 6.9, 6.10

! The dimension of the domain transformation matrices is 3
  dim=3
  cable%MI%dim=dim
  ALLOCATE(cable%MI%mat(dim,dim))
  cable%MV%dim=dim
  ALLOCATE(cable%MV%mat(dim,dim))

  cable%MI%mat(1,1)=0.5D0
  cable%MI%mat(1,2)=-0.5d0
  cable%MI%mat(1,3)=0d0
  
  cable%MI%mat(2,1)=1d0
  cable%MI%mat(2,2)=1d0
  cable%MI%mat(2,3)=0d0
  
  cable%MI%mat(3,1)=1d0
  cable%MI%mat(3,2)=1d0
  cable%MI%mat(3,3)=1d0


  cable%MV%mat(1,1)=1D0
  cable%MV%mat(1,2)=-1d0
  cable%MV%mat(1,3)=0d0
  
  cable%MV%mat(2,1)=0.5d0
  cable%MV%mat(2,2)=0.5d0
  cable%MV%mat(2,3)=-1d0
  
  cable%MV%mat(3,1)=0d0
  cable%MV%mat(3,2)=0d0
  cable%MV%mat(3,3)=1d0

! Set the local reference conductor numbering  
  ALLOCATE( cable%local_reference_conductor(2) )
  cable%local_reference_conductor(1)=2              ! differential mode, reference is the second conductor
  cable%local_reference_conductor(2)=0              ! common mode, reference is in the external domain and not known

! Set the local domain information: include a reference conductor in the count
  ALLOCATE( cable%local_domain_n_conductors(1:cable%tot_n_domains) )
  cable%local_domain_n_conductors(1)=2              ! differential mode domain 
  cable%local_domain_n_conductors(2)=2              ! common mode: reference in external domain 
  
! Set the external domain conductor and dielectric information
   
  ALLOCATE( cable%external_model(cable%n_external_conductors) )
  CALL reset_external_conductor_model(cable%external_model(1))
  cable%external_model(1)%conductor_type=circle
  cable%external_model(1)%conductor_radius =rw*Twisted_pair_equivalent_radius ! equivalent radius model for common mode
  cable%external_model(1)%dielectric_radius=rw*Twisted_pair_equivalent_radius+(rd-rw) ! add dielectric layer of the same thickness as specified
  cable%external_model(1)%dielectric_epsr=epsr
    
! set the conductor impedance model for the conductors
  cable%conductor_impedance(1)%impedance_model_type=impedance_model_type_cylindrical_with_conductivity
  cable%conductor_impedance(1)%radius=rw
  cable%conductor_impedance(1)%conductivity=sigma
  cable%conductor_impedance(1)%Resistance_multiplication_factor=1.5d0
  
  cable%conductor_impedance(2)%impedance_model_type=impedance_model_type_cylindrical_with_conductivity
  cable%conductor_impedance(2)%radius=rw
  cable%conductor_impedance(2)%conductivity=sigma
  cable%conductor_impedance(2)%Resistance_multiplication_factor=1.5d0
  
  CALL deallocate_Sfilter(epsr)

END SUBROUTINE twisted_pair_set_internal_domain_information
!
! NAME
!     twisted_pair_plot
!
! AUTHORS
!     Chris Smartt
!
! DESCRIPTION
!     plot twisted pair cable 
!
! COMMENTS
!     the angle has NO impact here due to the twisting
!
! HISTORY
!
!     started 14/4/2016 CJS 
!
!
SUBROUTINE twisted_pair_plot(cable,x_offset,y_offset,theta,xmin,xmax,ymin,ymax)

USE type_specifications
USE general_module

IMPLICIT NONE

! variables passed to subroutine

  type(cable_specification_type),intent(IN)    :: cable
  
  real(dp),intent(IN) :: x_offset,y_offset,theta
  real(dp),intent(INOUT) ::  xmin,xmax,ymin,ymax

! local variables

  real(dp) :: x,y,r
  real(dp) :: s

! START

! plot inner conductor, 1
  r=cable%parameters(1)   ! conductor radius
  s=cable%parameters(2)   ! conductor separation
  x=x_offset+(s/2d0)
  y=y_offset

  CALL write_circle(x,y,r,conductor_geometry_file_unit,xmin,xmax,ymin,ymax)

! plot inner conductor, 2
  x=x_offset-(s/2d0)
  y=y_offset

  CALL write_circle(x,y,r,conductor_geometry_file_unit,xmin,xmax,ymin,ymax)

! plot conductor1 dielectric
  r=cable%parameters(3)   ! dielectric radius, conductor 1
  x=x_offset+(s/2d0)
  y=y_offset

  CALL write_circle(x,y,r,dielectric_geometry_file_unit,xmin,xmax,ymin,ymax)

! plot conductor2 dielectric
  x=x_offset-(s/2d0)
  y=y_offset

  CALL write_circle(x,y,r,dielectric_geometry_file_unit,xmin,xmax,ymin,ymax)

  
  RETURN
  
END SUBROUTINE twisted_pair_plot