!
! 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
!
! OVERALL DESCRIPTION
!     These subroutines perform processes related to the conductor impedance models available
!
!     The subroutines operate on the conductor_impedance_model structure defined in cable_module.F90
! 
! Models are:
! 
! impedance_model_type_PEC                               
! impedance_model_type_cylindrical_with_conductivity     
! impedance_model_type_filter                            
! impedance_model_type_cylidrical_shell_with_conductivity
! impedance_model_type_cylindrical_shield                
! impedance_model_type_rectangular_with_conductivity    

! File Contents:
! SUBROUTINE read_conductor_impedance_model
! SUBROUTINE write_conductor_impedance_model
! SUBROUTINE evaluate_conductor_impedance_model
! SUBROUTINE calculate_Rdc
! SUBROUTINE calculate_internal_impedance
! SUBROUTINE ber_bei
! SUBROUTINE deallocate_conductor_impedance_model
!
! NAME
!     read_conductor_impedance_model
!
! AUTHORS
!     Chris Smartt
!
! DESCRIPTION
!  Read the conductor_impedance_model structure from file
!  The impedance model type is read first then the parameters relating to that particular model type
!
! COMMENTS
!      
!
! HISTORY
!
!     started 12/5/2016 CJS 
!     24/8/2016 CJS          Include cylindrical shell model
!
!
SUBROUTINE read_conductor_impedance_model(conductor_impedance,unit)

USE type_specifications

IMPLICIT NONE

! variables passed to subroutine

  type(conductor_impedance_model),intent(OUT) :: conductor_impedance   ! impedance model structure to read
  integer,intent(IN)                          :: unit                  ! file unit to read from

! local variables

! START

! reset to default values initially
  conductor_impedance%radius=0d0
  conductor_impedance%conductivity=0d0
  conductor_impedance%thickness=0d0
  conductor_impedance%width=0d0
  conductor_impedance%height=0d0
  conductor_impedance%Rdc=0d0
  conductor_impedance%Resistance_multiplication_factor=1d0

  read(unit,*,ERR=9000)conductor_impedance%impedance_model_type
  
  write(*,*)'MODEL TYPE ',conductor_impedance%impedance_model_type
  
  if (conductor_impedance%impedance_model_type.EQ.impedance_model_type_PEC) then

    RETURN   ! no additional information required for PEC type

  else if (conductor_impedance%impedance_model_type.EQ.impedance_model_type_cylindrical_with_conductivity) then
  
    read(unit,*,ERR=9000)conductor_impedance%radius
    read(unit,*,ERR=9000)conductor_impedance%conductivity
    read(unit,*,ERR=9000)conductor_impedance%Resistance_multiplication_factor

  else if (conductor_impedance%impedance_model_type.EQ.impedance_model_type_cylidrical_shell_with_conductivity) then
  
    read(unit,*,ERR=9000)conductor_impedance%radius
    read(unit,*,ERR=9000)conductor_impedance%conductivity
    read(unit,*,ERR=9000)conductor_impedance%Resistance_multiplication_factor
    read(unit,*,ERR=9000)conductor_impedance%thickness
    
    write(*,*)'THICKNESS ',conductor_impedance%thickness

  else if (conductor_impedance%impedance_model_type.EQ.impedance_model_type_filter) then

    CALL read_Sfilter(conductor_impedance%ZT_filter,unit)
    
  else if (conductor_impedance%impedance_model_type.EQ.impedance_model_type_cylindrical_shield) then
     
    read(unit,*,ERR=9000)conductor_impedance%radius
    read(unit,*,ERR=9000)conductor_impedance%conductivity
    read(unit,*,ERR=9000)conductor_impedance%Resistance_multiplication_factor
    read(unit,*,ERR=9000)conductor_impedance%thickness
    CALL read_Sfilter(conductor_impedance%ZT_filter,unit)
  
  else if (conductor_impedance%impedance_model_type.EQ.impedance_model_type_rectangular_with_conductivity) then
  
    read(unit,*,ERR=9000)conductor_impedance%width
    read(unit,*,ERR=9000)conductor_impedance%height
    read(unit,*,ERR=9000)conductor_impedance%conductivity
    read(unit,*,ERR=9000)conductor_impedance%Resistance_multiplication_factor
    
  else
  
    write(run_status,*)'ERROR in read_conductor_impedance_model: Unknown impedance model type:' &
                       ,conductor_impedance%impedance_model_type
    CALL write_program_status()
    STOP 1
    end if
  
  RETURN
  
9000 run_status='ERROR in read_conductor_impedance_model'
  CALL write_program_status()
  STOP 1

END SUBROUTINE read_conductor_impedance_model

!
! NAME
!     write_conductor_impedance_model
!
! AUTHORS
!     Chris Smartt
!
!
! DESCRIPTION
!  Write the conductor_impedance_model structure from file
!  The impedance model type is written first then the parameters relating to that particular model type
!
! COMMENTS
!      
!
! HISTORY
!
!     started 12/5/2016 CJS 
!     24/8/2016 CJS          Include cylindrical shell model
!
!

SUBROUTINE write_conductor_impedance_model(conductor_impedance,unit)

USE type_specifications

IMPLICIT NONE

! variables passed to subroutine

  type(conductor_impedance_model),intent(IN)    :: conductor_impedance  ! impedance model to write
  integer,intent(IN)                            :: unit                 ! unit to write to

! local variables

! START

  write(unit,*)conductor_impedance%impedance_model_type,' # Conductor impedance model type'
  
  if (conductor_impedance%impedance_model_type.EQ.impedance_model_type_PEC) then

    RETURN   ! no additional information required for PEC type

  else if (conductor_impedance%impedance_model_type.EQ.impedance_model_type_cylindrical_with_conductivity) then
  
    write(unit,*)conductor_impedance%radius,' # conductor radius'
    write(unit,*)conductor_impedance%conductivity,' # conductivity'
    write(unit,*)conductor_impedance%Resistance_multiplication_factor,' # Resistance_multiplication_factor'
    
  else if (conductor_impedance%impedance_model_type.EQ.impedance_model_type_cylidrical_shell_with_conductivity) then
  
    write(unit,*)conductor_impedance%radius,' # conductor radius'
    write(unit,*)conductor_impedance%conductivity,' # conductivity'
    write(unit,*)conductor_impedance%Resistance_multiplication_factor,' # Resistance_multiplication_factor'
    write(unit,*)conductor_impedance%thickness,' # shield thickness'

  else if (conductor_impedance%impedance_model_type.EQ.impedance_model_type_filter) then

    CALL write_Sfilter(conductor_impedance%ZT_filter,unit)
    
  else if (conductor_impedance%impedance_model_type.EQ.impedance_model_type_cylindrical_shield) then
     
    write(unit,*)conductor_impedance%radius,' # conductor radius'
    write(unit,*)conductor_impedance%conductivity,' # conductivity'
    write(unit,*)conductor_impedance%Resistance_multiplication_factor,' # Resistance_multiplication_factor'
    write(unit,*)conductor_impedance%thickness,' # shield thickness'
    CALL write_Sfilter(conductor_impedance%ZT_filter,unit)

  else if (conductor_impedance%impedance_model_type.EQ.impedance_model_type_rectangular_with_conductivity) then
  
    write(unit,*)conductor_impedance%width,' # conductor width'
    write(unit,*)conductor_impedance%height,' # conductor height'
    write(unit,*)conductor_impedance%conductivity,' # conductivity'
    write(unit,*)conductor_impedance%Resistance_multiplication_factor,' # Resistance_multiplication_factor'
    
  end if
  
  RETURN

END SUBROUTINE write_conductor_impedance_model
!
! NAME
!     evaluate_conductor_impedance_model
!
! AUTHORS
!     Chris Smartt
!
! DESCRIPTION
!  Given the parameters of a conductor impedance model and a frequency
!  evaluate: complex surface impedance, d.c. surface impedance
!          : complex transfer imedance, d.c. transfer impedance
!          as required.
!
! COMMENTS
!      
!
! HISTORY
!
!     started 12/5/2016 CJS 
!     24/8/2016 CJS          Include cylindrical shield model
!     24/8/2016 CJS          Include self (condcutor) and transfer impedances in the calculation
!
!

SUBROUTINE evaluate_conductor_impedance_model(conductor_impedance,f,Z_c,Rdc_c,Z_t,Rdc_t)

USE type_specifications
USE constants

IMPLICIT NONE

! variables passed to subroutine

  type(conductor_impedance_model),intent(IN) :: conductor_impedance  ! impedance model
  real(dp),intent(IN)      :: f                 ! frequency
  
  complex(dp),intent(OUT)   :: Z_c              ! Complex conductor surface impedance value to be returned
  real(dp),intent(OUT)      :: Rdc_c            ! d.c. conductor resistance to be returned
  complex(dp),intent(OUT)   :: Z_t              ! Complex transfer impedance value to be returned
  real(dp),intent(OUT)      :: Rdc_t            ! d.c. transfer resistance to be returned

! local variables

  real(dp)    :: sigma,rw,rs,t,w
  real(dp)    :: den
  
  real(dp)    :: f0

! START

  if (conductor_impedance%impedance_model_type.EQ.impedance_model_type_PEC) then

    Z_c=(0d0,0d0)
    Rdc_c=0d0
    Z_t=(0d0,0d0)
    Rdc_t=0d0

  else if (conductor_impedance%impedance_model_type.EQ.impedance_model_type_cylindrical_with_conductivity) then
  
    sigma=conductor_impedance%conductivity
    rw=conductor_impedance%radius        
    CALL calculate_internal_impedance(sigma,rw,f,Z_c,Rdc_c)
            
    Z_t=(0d0,0d0)
    Rdc_t=0d0

  else if (conductor_impedance%impedance_model_type.EQ.impedance_model_type_cylidrical_shell_with_conductivity) then
  
    sigma=conductor_impedance%conductivity
    rs=conductor_impedance%radius        
    t=conductor_impedance%thickness
    CALL calculate_internal_impedance_shell(sigma,rs,t,f,Z_c,Rdc_c)

! set the transfer impedance to be zero
    Z_t=(0d0,0d0)
    Rdc_t=0d0

  else if (conductor_impedance%impedance_model_type.EQ.impedance_model_type_filter) then

! transfer impedance calculation  
    Z_t=evaluate_Sfilter_frequency_response(conductor_impedance%ZT_filter,f)
    Rdc_t=conductor_impedance%ZT_filter%a%coeff(0)/conductor_impedance%ZT_filter%b%coeff(0)

! Set the conductor impedance to be equal to the transfer impedance. This model is no longer used and may be removed
    Z_c=Z_t
    Rdc_c=Rdc_t
   
  else if (conductor_impedance%impedance_model_type.EQ.impedance_model_type_cylindrical_shield) then

! transfer impedance calculation  
    Z_t=evaluate_Sfilter_frequency_response(conductor_impedance%ZT_filter,f)
    Rdc_t=conductor_impedance%ZT_filter%a%coeff(0)/conductor_impedance%ZT_filter%b%coeff(0)

! conductor impedance calculation
    sigma=conductor_impedance%conductivity
    rs=conductor_impedance%radius        

    t=conductor_impedance%thickness
    
    CALL calculate_internal_impedance_shell(sigma,rs,t,f,Z_c,Rdc_c)

  else if (conductor_impedance%impedance_model_type.EQ.impedance_model_type_rectangular_with_conductivity) then

! The rectangular conductor model assumes Z(jw)=(1+sqrt(f/f0)(1+j)) See C.R. Paul,
! "Analysis of Multiconductor Transmission Lines" 1st Edition, p 175-177 and p 320-323.

    Z_t=(0d0,0d0)
    Rdc_t=0d0
    
    sigma=conductor_impedance%conductivity
    w=conductor_impedance%width
    t=conductor_impedance%height
    
    CALL calculate_internal_impedance_rectangular(sigma,w,t,f,Z_c,Rdc_c)
    
  end if
  
  RETURN

END SUBROUTINE evaluate_conductor_impedance_model
!
! NAME
!     calculate_Rdc
!
! DESCRIPTION
!     Calculate the d.c.resistance of a cylindrical conductor given the
!     conductivity and radius 
!     If the conductivity is set to zero then return zero internal impedance i.e. for infinite conductivity
!     
! COMMENTS
!     
!
! HISTORY
!
!     started 27/04/2016 CJS based on phase 1 work in 2015
!
!
SUBROUTINE calculate_Rdc(sigma,r,Rdc)

USE type_specifications
USE constants

IMPLICIT NONE

! variables passed to subroutine

real(dp),intent(IN) :: sigma       ! conductivity
real(dp),intent(IN) :: r           ! shell radius

real(dp),intent(OUT) :: Rdc        ! d.c. resistance

! START

  if ((r.EQ.0d0).OR.(sigma.eq.0d0)) then
! zero parameters - indicates a perfect conductor so return zero d.c. resistance
    Rdc=0d0
    RETURN
    
  end if

! Evaluate the d.c. resistance
  Rdc=1d0/(sigma*pi*(r**2))       ! #EQN_REFERENCE_REQUIRED

  RETURN
  
END SUBROUTINE calculate_Rdc
!
! NAME
!     calculate_internal_impedance
!
! DESCRIPTION
!     Calculate the internal impedance of a cylindrical conductor given the
!     conductivity, radius and frequency
!     This software uses an exact formula using Kelvin functions if the 
!     radius if the wire is less than a specified number of skin depths (10), otherwise we use
!     an approximate formula. This is due to convergence problems when rw >> delta
!     If the conductivity is set to zero then return zero internal impedance i.e. for infinite conductivity
!     
! COMMENTS
!     The equation references are from the theory manual
!     but more detail is found in C. R. Paul,"Analysis of Multiconductor Transmission Lines" 1st edition
!     and S.A.Schelkunoff, 'The Electromagnetic Theory of Coaxial Transmission Lines and Cylindrical Shields', 
!     Bell System Technical Journal,Vol 13, No 4, pp 532-579, 1934.
!
! HISTORY
!
!     started 27/04/2016 CJS based on phase 1 work in 2015
!     8/5/2017         CJS: Include references to Theory_Manual
!
!
SUBROUTINE calculate_internal_impedance(sigma,r,f,Zint,Rdc)

USE type_specifications
USE constants

IMPLICIT NONE

! variables passed to subroutine

real(dp),intent(IN) :: sigma       ! conductivity
real(dp),intent(IN) :: r           ! shell radius
real(dp),intent(IN) :: f           ! frequency

complex(dp),intent(OUT) :: Zint    ! Condcutor impedance
real(dp),intent(OUT) :: Rdc        ! d.c. resistance

! local variables

real(dp) :: q,delta

! Kelvin function variables
real(dp) ber,bei,dber,dbei

! START 

  if (sigma.eq.0d0) then
! zero conductivity - indicates a perfect conductor so return zero internal impedance  
    zint=(0d0,0d0)
    Rdc=0d0
    RETURN
    
  end if
  
! special case for f=0. In this case return the d.c. resistance in both Rdc and zint
  if (f.eq.0d0) then
    CALL calculate_Rdc(sigma,r,Rdc)
    zint=Rdc
    RETURN
  end if

! skin depth
  delta=1d0/(sqrt(pi*f*mu0*sigma))    ! Theory_Manual_Eqn 3.63

  if (r.LT.10d0*delta) then
! use the Kelvin function form
! Q parameter equation          3.196b !  Theory_Manual_Eqn 3.64
    q=sqrt(2d0)*r/delta
  
    CALL ber_bei(q,ber,bei,dber,dbei)

! Evaluate the impedance, equation 3.196a !  Theory_Manual_Eqn 3.62
    zint=(1d0/(sqrt(2d0)*pi*r*sigma*delta))*(ber+j*bei)/(dbei-j*dber) 
  
  else
    
! Evaluate the impedance, equation 3.202b  !  C. R. Paul, Equation 3.202b
    zint=(1d0/(2d0*r))*sqrt(mu0/(pi*sigma))*sqrt(f)+j*2d0*pi*f*(1d0/(4d0*pi*r))*sqrt(mu0/(pi*sigma))/sqrt(f)
    
  end if  ! r greater than 10 skin depths
  
! Evaluate the d.c. resistance
  CALL calculate_Rdc(sigma,r,Rdc)

  RETURN
  
END SUBROUTINE calculate_internal_impedance
!
! NAME
!     calculate_internal_impedance_shell
!
! DESCRIPTION
!     Calculate the internal impedance of a cylindrical shell conductor given the
!     conductivity, radius, thickness and frequency
!     
! COMMENTS
!     The calculation is based on the Zaa calculation in equation 82 of:
! S. A. Schelkunoff, "The Electroomagnetic Theory of Coaxial Transmission 
! Lines and Cylindrical Shields" Bell Sys Tech J, vol 13, pp 532-579, Oct 1934.
!
! HISTORY
!
!     started 24/08/2016 CJS 
!     8/5/2017         CJS: Include references to Theory_Manual
!
!
SUBROUTINE calculate_internal_impedance_shell(sigma,r,t,f,Zint,Rdc)

USE type_specifications
USE constants

IMPLICIT NONE

! variables passed to subroutine

real(dp),intent(IN) :: sigma       ! conductivity
real(dp),intent(IN) :: r           ! shell radius
real(dp),intent(IN) :: t           ! shell thickness
real(dp),intent(IN) :: f           ! frequency

complex(dp),intent(OUT) :: Zint    ! Condcutor impedance
real(dp),intent(OUT) :: Rdc        ! d.c. resistance

! local variables

real(dp)    :: delta    ! skin depth
complex(dp) :: gamma    ! complex propagation constant in conductor
complex(dp) :: sinh_gt
complex(dp) :: cosh_gt

! START 

  if (sigma.eq.0d0) then
! zero conductivity - indicates a perfect conductor so return zero internal impedance  
    zint=(0d0,0d0)
    Rdc=0d0
    RETURN
    
  end if
  
  Rdc=1d0/(2d0*pi*r*t*sigma) !  Theory_Manual_Eqn 3.65
 
  if (f.NE.0d0) then
! skin depth calculation

    delta=1d0/(sqrt(pi*f*mu0*sigma))          ! skin depth in conductor Theory_Manual_Eqn 3.63
    gamma=cmplx(1d0,1d0)/cmplx(delta)         ! complex propagation constant in shield Theory_Manual_Eqn 3.66
        
    sinh_gt=(exp(gamma*t)-exp(-gamma*t))/(2d0,0d0)
    cosh_gt=(exp(gamma*t)+exp(-gamma*t))/(2d0,0d0)

! Theory_Manual_Eqn 3.67 
    zint=Rdc*gamma*t*cosh_gt/sinh_gt 
    
  else
! at zero frequency Zint=Rdc
  
    Zint=Rdc
 
  end if
  
  RETURN
  
END SUBROUTINE calculate_internal_impedance_shell
!
! NAME
!     calculate_internal_impedance_rectangular
!
! DESCRIPTION
!     Calculate the internal impedance of a rectangular conductor given the
!     conductivity, width, thickness and frequency
!     The equation references are from the theory manual
!     but more detail is found in C. R. Paul,"Analysis of Multiconductor Transmission Lines" 1st edition
!     
! COMMENTS
!     
!
! HISTORY
!
!     started 10/10/2016 CJS 
!     8/5/2017         CJS: Include references to Theory_Manual
!
!
SUBROUTINE calculate_internal_impedance_rectangular(sigma,w,t,f,Zint,Rdc)

USE type_specifications
USE constants

IMPLICIT NONE

! variables passed to subroutine

real(dp),intent(IN) :: sigma       ! conductivity
real(dp),intent(IN) :: w           ! conductor width
real(dp),intent(IN) :: t           ! conductor thickness
real(dp),intent(IN) :: f           ! frequency

complex(dp),intent(OUT) :: Zint    ! Condcutor impedance
real(dp),intent(OUT) :: Rdc        ! d.c. resistance

! local variables

! START 

  if (sigma.eq.0d0) then
! zero conductivity - indicates a perfect conductor so return zero internal impedance  
    Zint=(0d0,0d0)
    Rdc=0d0
    RETURN
    
  end if
    
  Rdc=1d0/(w*t*sigma) ! Theory_Manual_Eqn 3.69

  if (f.NE.0d0) then

    Zint=Rdc+(0.5d0/(w+t))*sqrt(mu0/sigma)*sqrt(j*2d0*pi*f) !  Theory_Manual_Eqn 3.68 with 3.70
    
  else
! at zero frequency Zint=Rdc
  
    Zint=Rdc
 
  end if
  
  RETURN
  
END SUBROUTINE calculate_internal_impedance_rectangular
!
! NAME
!     ber_bei
!
! DESCRIPTION
!     Calculate Kelvin functions and their derivatives
!     using a series expansion method 
! 
! COMMENTS
!     may not be the most efficient way to do this...
!
! HISTORY
!
!     started 27/04/2016 CJS based on phase 1 work in 2015
!
!
SUBROUTINE ber_bei(x,ber,bei,dber,dbei)

! calculate Kelvin functions and derivatives
! See Abramowitz and Stegun, "Handbook of Mathematical Functions" equation 9.9.10

USE type_specifications
USE constants

IMPLICIT NONE

! variables passed to subroutine

real(dp),intent(IN)  :: x ! input argument

real(dp),intent(OUT) :: ber,bei,dber,dbei  ! output kelvin fuction and derivative values

! local variables

real(dp)        :: arg,sign,fac,term,last_term,p
integer         :: k

real(dp),parameter :: cvg_test=1d-12
integer,parameter  :: kmax=1000

! START

arg=x*x/4d0

ber=1d0
dber=0d0
last_term=ber
sign=-1d0

do k=2,kmax,2

  fac=dble(k*(k-1))
  term=last_term*sign*(arg/fac)*(arg/fac)
  ber=ber+term
  
  p=dble(k)*2d0
  dber=dber+term*p/x
    
  if (abs(term).LT.cvg_test) GOTO 100  
  
  last_term=term

end do

100 CONTINUE

bei=arg
last_term=bei
dbei=x/2d0
sign=-1d0

do k=3,kmax,2

  fac=dble(k*(k-1))
  term=last_term*sign*(arg/fac)*(arg/fac)
  bei=bei+term
  
  p=dble(k)*2d0
  dbei=dbei+term*p/x
  
  if (abs(term).LT.cvg_test) GOTO 200  
  
  last_term=term

end do

200 CONTINUE

RETURN

END SUBROUTINE ber_bei

!
! NAME
!     deallocate_conductor_impedance_model
!
! AUTHORS
!     Chris Smartt
!
! DESCRIPTION
!  deallocate filters required for comnductor impedance models
!
! COMMENTS
!      
!
! HISTORY
!
!     started 12/5/2016 CJS 
!
!

SUBROUTINE deallocate_conductor_impedance_model(conductor_impedance)

USE type_specifications

IMPLICIT NONE

! variables passed to subroutine

  type(conductor_impedance_model),intent(INOUT)    :: conductor_impedance

! local variables

! START

  CALL deallocate_Sfilter(conductor_impedance%ZT_filter)

END SUBROUTINE deallocate_conductor_impedance_model