! ! 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 . ! ! 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 . ! ! 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 ! impedance_model_type_FH2 ! 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 if (conductor_impedance%impedance_model_type.EQ.impedance_model_type_FH2) then read(unit,*,ERR=9000)conductor_impedance%Rdc 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' else if (conductor_impedance%impedance_model_type.EQ.impedance_model_type_FH2) then write(unit,*)conductor_impedance%Rdc,' # conductor dc resistance' 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) else if (conductor_impedance%impedance_model_type.EQ.impedance_model_type_FH2) then Z_c=conductor_impedance%Rdc Rdc_c=conductor_impedance%Rdc Z_t=(0d0,0d0) Rdc_t=0d0 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