!
! 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:
! PROGRAM shield_conductor_and_transfer_impedance_model_builder
!
! NAME
! shield_conductor_and_transfer_impedance_model_builder
!
! AUTHORS
! Chris Smartt
!
! DESCRIPTION
! Calculate the transfer impedance and conductor impedance of a shield
! based on the geometrical and electrical properties of the braid forming
! the shield. The calculations are based on the NLR theory presented in the Theory Manual
!
! Example input file:
!0.002 ! braid diameter, D (m)
!8 ! Number of carriers, C
!10 ! Number of wires in a carrier, N
!0.000025 ! diameter of a single wire, d (m)
!5E7 ! conductivity of wires (S/m)
!50.0 ! pitch angle of the braid (degrees)
!4 ! order of transfer impedance model
!log ! frequency range type
!1E5 1E9 1000 ! fmin, fmax, number of frequencies
!
! COMMENTS
! The braid equivalent thickness is calculated from the d.c. resistance of the braid and the braid conductivity
! See appendix 1 of the user guide for the theory
!
! HISTORY
!
! started 22/08/2016 CJS
! started 29/09/2017 CJS Include more rigorous hole inductance calculation
! using elliptic integrals
! 11/10/2017 CJS special case for 0 order models to ensure that the d.c. resistance is correct
!
PROGRAM shield_conductor_and_transfer_impedance_model_builder
USE type_specifications
USE general_module
USE constants
USE frequency_spec
USE filter_module
USE Sfilter_fit_module
IMPLICIT NONE
! local variables
character(len=filename_length) :: braid_name ! name of the braid model
character(len=filename_length) :: filename ! filename for the braid model specification
logical :: file_exists
real(dp) :: D ! braid diameter, D (m)
integer :: C ! Number of carriers, C
integer :: N ! Number of wires in a carrier, N
real(dp) :: dw ! diameter of a single wire, d (m)
real(dp) :: sigma ! conductivity of wires (S/m)
real(dp) :: alpha ! pitch angle of the braid (degrees)
integer :: order ! order of transfer impedance and conductor impedance model
type(frequency_specification) :: frequency_data ! frequency range data for transfer impedance and conductor impedance model
complex(dp),allocatable :: Zd(:) ! Frequency dependent diffusion impedance data
complex(dp),allocatable :: Zt(:) ! Frequency dependent transfer impedance data
complex(dp),allocatable :: Zc(:) ! Frequency dependent shield conductor impedance data
complex(dp) :: Zt_fit ! Vector fit model Frequency dependent transfer impedance data
complex(dp) :: Zc_fit ! Vector fit model Frequency dependent shield conductor impedance data
type(Sfilter) :: Zd_filter ! Frequency dependent diffusion impedance rational function model
type(Sfilter) :: Zt_filter ! Frequency dependent transfer impedance rational function model
type(Sfilter) :: Zc_filter ! Frequency dependent shield conductor impedance rational function model
complex(dp) :: M ! Total contribution originating from braid magnetic field leakage
complex(dp) :: M12 ! per-unit-length hole inductance
complex(dp) :: Mb ! braid inductance
complex(dp) :: Ms ! skin inductance
real(dp) :: w ! angular frequency
real(dp) :: R0 ! d.c. resistance of shield
complex(dp) :: gamma ! complex propagation constant in shield
real(dp) :: delta ! skin depth in shield
real(dp) :: T ! shield conductor thickness
real(dp) :: lh ! hole length
real(dp) :: wh ! hole width
real(dp) :: S ! hole area
real(dp) :: req ! equivalent hole radius
real(dp) :: b ! width between holes
real(dp) :: hh ! average height for braid inductance calculation
real(dp) :: Dm ! Mean diameter of braid for braid inductance calculation
real(dp) :: v ! Number of holes per unit length
real(dp) :: gc ! constant used in hole inductance calculation
real(dp) :: F ! Fill factor of braid
real(dp) :: K ! Optical coverage of braid
real(dp) :: Ck ! constant used in hole inductance calculation
real(dp) :: RT ! transfer resistance in ZT=R+jwL model
real(dp) :: LT ! transfer inductance in ZT=R+jwL model
! variables for intermediate quantities used in the calculations
real(dp) :: P
real(dp) :: kappa
complex(dp) :: sinh_gT
complex(dp) :: cosh_gT
complex(dp) :: u
complex(dp) :: nu
real(dp) :: k1 ! factor used in Kley's model for braid inductance
real(dp) :: mum ! hole polarizability factor
real(dp) :: a ! braid radius
real(dp) :: e
real(dp) :: Th,CkT ! attenuation factor
integer :: floop ! frequency loop variable
integer :: Zt_aorder,Zt_border ! order of transfer impedance model
integer :: i
integer :: ierr ! integer to return error codes from file reads
integer :: model_type
integer,parameter :: NLR=1
integer,parameter :: Kley=2
! function types
real(dp) :: Em
real(dp) :: Km
! START
! Open the input file describing the braid parameters
! This file could be created by the associated GUI or otherwise generated
verbose=.TRUE.
! model_type=NLR
model_type=Kley
program_name="shield_conductor_and_transfer_impedance_model_builder"
run_status='Started'
CALL write_program_status()
CALL read_version()
CALL write_license()
write(*,*)'Enter the name of the shield braid specification data (without .braid_spec extension)'
read(*,'(A)')braid_name
filename=trim(braid_name)//braid_spec_file_extn
inquire(file=trim(filename),exist=file_exists)
if (.NOT.file_exists) then
run_status='ERROR in shield_conductor_and_transfer_impedance_model_builder, Cannot find the file:'//trim(filename)
CALL write_program_status()
STOP 1
end if
! open and read the file
OPEN(unit=braid_spec_file_unit,file=filename)
if(verbose) write(*,*)'Opened file:',trim(filename)
read(braid_spec_file_unit,*,IOSTAT=ierr)D
if (ierr.NE.0) then
run_status='ERROR reading shield diameter'
CALL write_program_status()
STOP 1
end if
read(braid_spec_file_unit,*,IOSTAT=ierr)C
if (ierr.NE.0) then
run_status='ERROR reading number of carriers'
CALL write_program_status()
STOP 1
end if
read(braid_spec_file_unit,*,IOSTAT=ierr)N
if (ierr.NE.0) then
run_status='ERROR reading number of number of wires in a carrier'
CALL write_program_status()
STOP 1
end if
read(braid_spec_file_unit,*,IOSTAT=ierr)dw
if (ierr.NE.0) then
run_status='ERROR reading wire diameter'
CALL write_program_status()
STOP 1
end if
read(braid_spec_file_unit,*,IOSTAT=ierr)sigma
if (ierr.NE.0) then
run_status='ERROR reading wire conductivity'
CALL write_program_status()
STOP 1
end if
read(braid_spec_file_unit,*,IOSTAT=ierr)alpha
if (ierr.NE.0) then
run_status='ERROR reading pitch angle of the braid'
CALL write_program_status()
STOP 1
end if
! convert alpha to radians
alpha=alpha*pi/180d0
read(braid_spec_file_unit,*,IOSTAT=ierr)order
if (ierr.NE.0) then
run_status='ERROR reading the model order for the transfer impedance and conductor impedance models'
CALL write_program_status()
STOP 1
end if
CALL read_and_set_up_frequency_specification(frequency_data,braid_spec_file_unit)
! close the file with the cable data
CLOSE(unit=braid_spec_file_unit)
! Evaluate the shield transfer impedance and conductor impedance over the specified frequency range
ALLOCATE(Zd(1:frequency_data%n_frequencies))
ALLOCATE(Zt(1:frequency_data%n_frequencies))
ALLOCATE(Zc(1:frequency_data%n_frequencies))
! Calculate the solution parameters which are frequency independent
gc=(2d0/pi)**(3d0/2d0)
P=C*tan(alpha)/(2d0*pi*D)
v=P*C
F=N*C*dw/(2d0*pi*D*cos(alpha)) ! fill factor
K=2d0*F-F*F ! optical coverage
lh=(1d0-F)*N*dw/(F*sin(alpha))
wh=(1d0-F)*N*dw/(F*cos(alpha))
S=pi*wh*lh/4d0 ! hole area
if (lh.GT.wh) then
e=sqrt(1d0-(wh/lh)**2)
else
e=sqrt(1d0-(lh/wh)**2)
end if
write(*,*)lh,wh,e
Dm=D+2d0*dw
! b=(2d0*pi*Dm/C)*cos(alpha)-N*dw
b=N*dw*(1d0-F)/F
! Calculate an equivalent shield thickness from R0, D and sigma
R0=4d0/(pi*dw*dw*N*C*sigma*cos(alpha)) ! d.c. resistance. equation 5.63 of D1
T=1d0/(2d0*pi*sigma*(D/2d0)*R0)
if (model_type.EQ.NLR) then
! Chimney effect stuff
req=sqrt(S/pi)
kappa=1.84d0/req
Ck=0.875d0*exp(-kappa*T)
if (b.GT.dw) then
hh=2d0*dw/(1d0+b/dw)
else
hh=dw
end if
! Hole inductance term ! equation 5.82 of D1
M12=1.08D0*gc*(pi*mu0/(6d0*C))*((1d0-K)**(3d0/2d0))*(2d0-cos(alpha))*Ck
! Braid inductance term
Mb=-mu0*(hh/(4d0*pi*Dm))*(1d0-(tan(alpha))**2)
! Skin inductance term. Assumed to be zero here.
Ms=0d0
else if (model_type.EQ.Kley) then
! Chimney effect stuff
! Tesche eqn 9.51
a=D/2d0+dw
Th=9.6d0*F* ((K*K*dw/(2d0*a))**0.33333)
Ck=0.875*exp(-Th)
! Hole inductance term !
M12=(pi*4.0*pi*1E-7/(6.0*C))*((1.0-K)**(1.5))*e*e/(Em(e)-(1.0-e*e)*Km(e))*Ck
! Braid inductance term
hh=dw
k1=(pi/4d0)/(0.6667D0*F*cos(alpha)+pi/10d0)
Mb=-mu0*(hh/(4d0*pi*Dm))*(0.22d0/(F*cos(alpha)))*cos(2d0*k1*alpha) ! OK NLR and Tesche(9.54)
! Skin inductance term. Assumed to be zero here.
Ms=0d0
end if
! Total field leakage contributions
M=M12+Mb+Ms
if (verbose) then
write(*,*)'braid circummference, cb=',pi*D
write(*,*)'C=',C,' Number of carriers'
write(*,*)'N=',N,' Number of conductors in each carrier'
write(*,*)'W=',N*dw,' Width of each carrier'
write(*,*)'W=',N*dw/cos(alpha),' Width of each carrier in circumferential direction'
write(*,*)'cb/(C/2)=',2d0*pi*D/C,' circumferential dimension for each carrier (note overlap)'
write(*,*)'P=',P
write(*,*)'v=',v,' Number of holes per unit length in braid'
write(*,*)'F=',F,' Fill factor'
write(*,*)'K=',K,' Optical coverage'
write(*,*)'l=',lh,' hole length'
write(*,*)'w=',wh,' hole width'
write(*,*)'S=',S, ' hole area'
if (model_type.EQ.NLR) then
write(*,*)'req=',req,' hole equivalent radius'
write(*,*)'k=',kappa,' hole cutoff k value'
else if (model_type.EQ.Kley) then
write(*,*)'Th=',Th,' hole attenuation factor'
end if
write(*,*)'Ck=',Ck,' hole inductance factor'
write(*,*)'Ro',R0,' braid d.c. resistance'
write(*,*)'T ',T,' braid equivalent thickness'
write(*,*)'mean braid diameter, Dm=',Dm
write(*,*)'Width between holes, b=',b
write(*,*)'Average height for braid inductance, hh=',hh
write(*,*)'M12=',M12,' Hole inductance'
write(*,*)'Mb =',Mb,' Braid inductance'
write(*,*)'Ms=',Ms,' Skin inductance'
write(*,*)'M=',M,' Total transfer inductance'
end if ! verbose
CALL plot_braid_figure(C,N,dw,alpha,D,braid_name)
open(unit=83,file='Zt_Zc.dat')
do floop=1,frequency_data%n_frequencies
w=2d0*pi*frequency_data%freq_list(floop)
! Diffusion impedance term
delta=sqrt(2d0/(w*mu0*sigma)) ! skin depth in conductor
gamma=cmplx(1d0,1d0)/cmplx(delta) ! complex propagation constant in shield
sinh_gT=(exp(gamma*T)-exp(-gamma*T))/(2d0,0d0)
cosh_gT=(exp(gamma*T)+exp(-gamma*T))/(2d0,0d0)
Zd(floop)=R0*gamma*T/sinh_gT ! equation 5.62 of D1
! Terms for calculation in Schelkunoff's notation
nu=j*w*mu0/gamma
u=t*sqrt(2d0*sigma*w*mu0)
! Transfer impedance is the sum of the diffusion impedance and the transfer inductance
Zt(floop)=Zd(floop)+j*w*M
! Conductor impedance term
Zc(floop)=R0*gamma*T*cosh_gT/sinh_gT
write(83,8000)frequency_data%freq_list(floop),real(Zt(floop)),aimag(Zt(floop)),real(Zc(floop)),aimag(Zc(floop))
8000 format(5ES16.6)
if (floop.EQ.frequency_data%n_frequencies) then
! work out the R+jwL model
RT=R0
LT=aimag(Zt(floop))/w
end if
end do ! next frequency
close(unit=83)
write(*,*)'________________________________________________________________'
write(*,*)''
write(*,*)'ZT=RT+jwLT model:'
write(*,*)'RT=',RT,' ohms/m'
write(*,*)'LT=',LT,' H/m'
write(*,*)'________________________________________________________________'
write(*,*)''
! Create a rational function model of the transfer impedance data in two stages
! First, create a rational function model of the frequency dependent diffusion impedance
! Then add the transfer inductance term jwM
! calculate_Sfilter with border=aorder+1 and with fit_type=0 i.e. Zd->0 as f-> infinity
if (verbose) write(*,*)'Calculate_Sfilter for diffusion impedance, Zd'
if (order.NE.0) then
CALL Calculate_Sfilter(Zd,frequency_data%freq_list,frequency_data%n_frequencies, &
Zd_filter,order,1,0) ! call with fit_type=0
else
! for a zero order model, return the d.c. resistance of the sheild
Zd_filter=allocate_Sfilter(0,0)
Zd_filter%wnorm=1d0
Zd_filter%a%coeff(0)=R0
Zd_filter%b%coeff(0)=1d0
end if
! Add the transfer inductance term to the diffusion impedance filter i.e. Zt=Zd+j*w*M
Zt_aorder=max(Zd_filter%a%order,Zd_filter%b%order+1)
Zt_border=Zd_filter%b%order
Zt_filter=allocate_Sfilter(Zt_aorder,Zt_border)
Zt_filter%wnorm=Zd_filter%wnorm
if (verbose) write(*,*)'Zt_filter%wnorm=',Zt_filter%wnorm
! copy a coefficients from Zd_filter to Zt_filter
do i=0,Zd_filter%a%order
Zt_filter%a%coeff(i)=Zd_filter%a%coeff(i)
end do
! copy b coefficients from Zd_filter to Zt_filter
do i=0,Zd_filter%b%order
Zt_filter%b%coeff(i)=Zd_filter%b%coeff(i)
end do
! add the jwM term
do i=0,Zd_filter%b%order
Zt_filter%a%coeff(i+1)=Zt_filter%a%coeff(i+1)+(M*Zt_filter%wnorm)*Zd_filter%b%coeff(i)
end do
! Create a rational function model of the shield conductor impedance data
! call calculate_Sfilter with border=aorder and with fit_type=0
if (verbose) write(*,*)'Calculate_Sfilter for surface impedance, Zc'
if (order.NE.0) then
CALL Calculate_Sfilter(Zc,frequency_data%freq_list,frequency_data%n_frequencies, &
Zc_filter,order,0,0) ! call with fit_type=0
else
! for a zero order model, return the d.c. resistance of the sheild
Zc_filter=allocate_Sfilter(0,0)
Zc_filter%wnorm=1d0
Zc_filter%a%coeff(0)=R0
Zc_filter%b%coeff(0)=1d0
end if
! Write vector fit models to file
open(unit=84,file='Zt_fit.fout')
open(unit=85,file='Zc_fit.fout')
do floop=1,frequency_data%n_frequencies
Zt_fit=evaluate_Sfilter_frequency_response(Zt_filter,frequency_data%freq_list(floop))
Zc_fit=evaluate_Sfilter_frequency_response(Zc_filter,frequency_data%freq_list(floop))
write(84,8000)frequency_data%freq_list(floop),real(Zt_fit),aimag(Zt_fit)
write(85,8000)frequency_data%freq_list(floop),real(Zc_fit),aimag(Zc_fit)
end do
close(unit=84)
close(unit=85)
! Open a file for the shield model
filename=trim(braid_name)//shield_model_file_extn
open(unit=shield_model_file_unit,file=filename)
! Write the shield equivalent thickness and conductivity
write(shield_model_file_unit,*)D/2d0,' # Parameter Shield radius (m)'
write(shield_model_file_unit,*)T, ' # Parameter Equivalent shield thickness (m)'
write(shield_model_file_unit,*)sigma,' # Parameter Shield conductivity (S/m)'
! Write the transfer impedance model to the shield model file
write(shield_model_file_unit,*)' '
write(shield_model_file_unit,*)'# Transfer impedance model'
write(shield_model_file_unit,*)' '
CALL Write_Sfilter(Zt_filter,shield_model_file_unit)
! Write the shield conductor model to the shield model file
write(shield_model_file_unit,*)' '
write(shield_model_file_unit,*)'# Conductor surface impedance model'
write(shield_model_file_unit,*)' '
CALL Write_Sfilter(Zc_filter,shield_model_file_unit)
write(shield_model_file_unit,*)' '
! Write the solution parameters to the shield model file
write(shield_model_file_unit,*)'# Shield parameters used in the shield model calculation'
write(shield_model_file_unit,*)'braid circummference, cb=',pi*D
write(shield_model_file_unit,*)'C=',C,' Number of carriers'
write(shield_model_file_unit,*)'N=',N,' Number of conductors in each carrier'
write(shield_model_file_unit,*)'W=',N*dw,' Width of each carrier'
write(shield_model_file_unit,*)'W=',N*dw/cos(alpha),' Width of each carrier in circumferential direction'
write(shield_model_file_unit,*)'cb/(C/2)=',2d0*pi*D/C,' circumferential dimension for each carrier (note overlap)'
write(shield_model_file_unit,*)'P=',P
write(shield_model_file_unit,*)'v=',v,' Number of holes per unit length in braid'
write(shield_model_file_unit,*)'F=',F,' Fill factor'
write(shield_model_file_unit,*)'K=',K,' Optical coverage'
write(shield_model_file_unit,*)'l=',lh,' hole length'
write(shield_model_file_unit,*)'w=',wh,' hole width'
write(shield_model_file_unit,*)'S=',S, ' hole area'
if (model_type.EQ.NLR) then
write(shield_model_file_unit,*)'req=',req,' hole equivalent radius'
write(shield_model_file_unit,*)'k=',kappa,' hole cutoff k value'
else if (model_type.EQ.Kley) then
write(shield_model_file_unit,*)'Th=',Th,' hole attenuation factor'
end if
write(shield_model_file_unit,*)'Ck=',Ck,' hole inductance factor'
write(shield_model_file_unit,*)'Ro',R0,' braid d.c. resistance'
write(shield_model_file_unit,*)'T ',T,' braid equivalent thickness'
write(shield_model_file_unit,*)'mean braid diameter, Dm=',Dm
write(shield_model_file_unit,*)'Width between holes, b=',b
write(shield_model_file_unit,*)'Average height for braid inductance, hh=',hh
write(shield_model_file_unit,*)'M12=',M12,' Hole inductance'
write(shield_model_file_unit,*)'Mb =',Mb,' Braid inductance'
write(shield_model_file_unit,*)'Ms=',Ms,' Skin inductance'
write(shield_model_file_unit,*)'M=',M,' Total transfer inductance'
! Close the shield model file
close(unit=shield_model_file_unit)
! deallocate memory and finish up
DEALLOCATE(Zd)
DEALLOCATE(Zt)
DEALLOCATE(Zc)
run_status='Finished_Correctly'
CALL write_program_status()
END PROGRAM shield_conductor_and_transfer_impedance_model_builder
real(dp) FUNCTION Em(m)
USE type_specifications
IMPLICIT NONE
! Elliptic integral approximation, see Abramowitz and Stegun 17.3.35
real(dp) m
real(dp),parameter :: a1=0.4630151
real(dp),parameter :: a2=0.1077812
real(dp),parameter :: b1=0.2452727
real(dp),parameter :: b2=0.0412496
real(dp) m1
! START
if ((m.LT.0.0).OR.(m.GT.1.0)) then
write(*,*)'Error: m out of range 0<m<1 in Em(m)'
write(*,*)'m=',m
STOP
end if
m1=1.0-m
Em=1.0+a1*m1+a2*m1*m1+(b1*m1+b2*m1*m1)*log(1.0/m1)
END
real(dp) FUNCTION Km(m)
USE type_specifications
IMPLICIT NONE
! Elliptic integral approximation, see Abramowitz and Stegun 17.3.33
real(dp) m
real(dp),parameter :: a0=1.3862944
real(dp),parameter :: a1=0.1119723
real(dp),parameter :: a2=0.0725296
real(dp),parameter :: b0=0.5
real(dp),parameter :: b1=0.1213478
real(dp),parameter :: b2=0.0288729
real(dp) m1
! START
if ((m.LT.0.0).OR.(m.GT.1.0)) then
write(*,*)'Error: m out of range 0<m<1 in Km(m)'
write(*,*)'m=',m
STOP
end if
m1=1.0-m
Km=a0+a1*m1+a2*m1*m1+(b0+b1*m1+b2*m1*m1)*log(1.0/m1)
END
!
! ____________________________________________________________
!
!
SUBROUTINE plot_braid_figure(C,N,dw,alpha,D,braid_name)
USE type_specifications
USE general_module
USE constants
IMPLICIT NONE
! variables passed to subroutine
integer :: N ! Number of wires in a carrier, N
integer :: C ! Number of carriers, C
real(dp) :: dw ! diameter of a single wire, d (m)
real(dp) :: alpha ! pitch angle of the braid (degrees)
real(dp) :: D ! braid diameter, D (m)
character(len=filename_length) :: braid_name ! name of the braid model
! local variables
real(dp) :: lx,ly
real(dp) :: dwx,dwy,dydx
real(dp) :: dyc
real(dp) :: x,y
real(dp) :: x1,y1
real(dp) :: x2,y2
real(dp) :: dy,dx
real(dp) :: wC,tC,tH,yH,tW
real(dp) :: l1,l2,l3,l4
real(dp) :: ax,ay,ox,l,ox_C
integer :: i,ii
! START
ly=pi*D ! extent of the circumferential direction
if (alpha.NE.0d0) then
lx=ly/tan(alpha)
lx=max(lx,1.62d0*ly)
else
lx=1.62d0*ly
end if
dwy=dw/cos(alpha)
dwx=dw/sin(alpha)
dydx=tan(alpha)
dyc=ly/(C/2d0)
wC=N*dw ! width of a carrier
tC=wC/sin(2d0*alpha)
tW=dw/sin(2d0*alpha)! width of a conductor
yH=dyc-wC/cos(alpha)
tH=yH/(2d0*sin(alpha))
l1=3d0*tH+2d0*tC
l2=tC
l3=tH
l4=tC
ax=-cos(alpha) ! gradient of start point line
ay=sin(alpha)
ox=dyc/tan(alpha) ! offset along x
write(*,*)'ly=',ly
write(*,*)'space for each carrier in circumferential direction (dyc)',dyc
write(*,*)'carrier overlap length ',tC
write(*,*)'conductor overlap length ',tW
write(*,*)'hole width ',yH
write(*,*)'hole edge length ',tH
write(*,*)'offset along x ',ox
open(unit=local_file_unit,file='braid_figure.dat')
!GOTO 1234
! WIRES GOING IN FIRST DIRECTION (UP IN FIGURE)
! loop over carriers
do ii=0,C/2-1
! plot a single carrier
do i=0,N ! the number of lines to plot is N+1
! calculate the first end point of the line
l=i*tW
x1=0d0+ax*l-ii*ox
y1=0d0+ay*l
CALL write_braid_line(x1,y1,dydx,lx,ly,l1,l2,l3,l4,local_file_unit)
end do ! next line (wire in carrier)
end do ! next carrier
! WIRES GOING IN SECOND DIRECTION (DOWN IN FIGURE)
1234 CONTINUE
dydx=-tan(alpha)
ax=-cos(alpha) ! gradient of start point line
ay=-sin(alpha)
! loop over carriers
do ii=0,C/2-1
ox_C=-3d0*ox
if (mod(ii,4).EQ.1) then
ox_C=-2d0*ox
else if (mod(ii,4).EQ.2) then
ox_C=-1d0*ox
else if (mod(ii,4).EQ.3) then
ox_C=-0d0*ox
end if
! plot a single carrier
do i=0,N ! the number of lines to plot is N+1
! calculate the first end point of the line
l=i*tW
x1=0d0+ax*l-ii*ox-ox/2d0
y1=ly+ay*l+dyc/2d0
CALL write_braid_line(x1,y1,dydx,lx,ly,l1,l2,l3,l4,local_file_unit)
end do ! next line (wire in carrier)
end do ! next carrier
close(unit=local_file_unit)
! write a file to plot the braid figure
open(unit=local_file_unit,file='plot_braid.plt')
write(local_file_unit,'(A,A,A)')'set title "',trim(braid_name),'"'
write(local_file_unit,'(A)')'set xlabel "Longitudinal direction"'
write(local_file_unit,'(A)')'set ylabel "Circumferential"'
write(local_file_unit,'(A)')'plot "braid_figure.dat" u 1:2 w l'
write(local_file_unit,'(A)')'pause -1'
close(unit=local_file_unit)
END SUBROUTINE plot_braid_figure
!
! ____________________________________________________________
!
!
SUBROUTINE write_braid_line(x,y,m,lx,ly,l1,l2,l3,l4,unit)
USE type_specifications
IMPLICIT NONE
real(dp),intent(IN) :: x,y ! point on the line
real(dp),intent(IN) :: m ! gradient
real(dp),intent(IN) :: lx ! xlimit
real(dp),intent(IN) :: ly ! ylimit
real(dp),intent(IN) :: l1 ! length1
real(dp),intent(IN) :: l2 ! length2
real(dp),intent(IN) :: l3 ! length3
real(dp),intent(IN) :: l4 ! length4
integer,intent(IN) :: unit
! local variables
real(dp) :: d1
real(dp) :: d2
logical :: line_end
! START
d1=0d0 ! start of line
10 CONTINUE
d2=d1+l1
! write(*,*)'Write braid line'
! write(*,*)x,y
! write(*,*)d1,d2
CALL write_braid_line_l(x,y,m,lx,ly,d1,d2,line_end,unit)
if (line_end) RETURN
d1=d2+l2 ! allow blank space l1 to l2
d2=d1+l3
! write(*,*)'Write braid line'
! write(*,*)x,y
! write(*,*)d1,d2
CALL write_braid_line_l(x,y,m,lx,ly,d1,d2,line_end,unit)
if (line_end) RETURN
d1=d2+l4 ! allow blank space l3 to l4
GOTO 10
RETURN
END SUBROUTINE write_braid_line
!
! ____________________________________________________________
!
!
SUBROUTINE write_braid_line_l(x0,y0,m,lx,ly,l1,l2,line_end,unit)
USE type_specifications
IMPLICIT NONE
real(dp),intent(IN) :: x0,y0 ! intial point on the line
real(dp),intent(IN) :: m ! gradient
real(dp),intent(IN) :: lx ! xlimit
real(dp),intent(IN) :: ly ! ylimit
real(dp),intent(IN) :: l1 ! length1
real(dp),intent(IN) :: l2 ! length2
integer,intent(IN) :: unit
logical,intent(OUT) :: line_end ! flag the end of the line
! local variables
real(dp) :: c ! constant for line equation1
real(dp) :: ax,ay ! constants for parametric line equation
real(dp) :: x1,y1 ! first point on the line segment
real(dp) :: x2,y2 ! second point on the line segment
real(dp) :: l_line ! length of line segment
real(dp) :: l_segment ! length of split line segment
integer :: i
! START
line_end=.FALSE.
! get the equation of the line in the form y=mx+c
c=y0-m*x0
! get the parametric equation of the line in the form x=x+ax*l, y=y+ay*l
ax=sqrt(1d0/(1d0+m*m))
ay=m*ax
! calculate the first end point of the line
x1=x0+ax*l1
y1=y0+ay*l1
! return if both ends of the line are outsude the range[0:lx]
if (x1.GT.lx) then
line_end=.TRUE.
RETURN
end if
x2=x0+ax*l2
if (x2.LT.0d0) then
RETURN
end if
! the first point x value may be less than zero so move along the line to x1=0
if (x1.LT.0d0) then
c=y1-m*x1
x1=0d0
y1=m*x1+c
end if
! the first point may be out of the y plotting range so move it back in
! by adding (subtracting) an integer times the period ly
if (y1.GT.0d0) then
i=INT(y1/ly)
else
i=INT(y1/ly)-1
end if
y1=y1-i*ly
! we now have an initial point inside the plotting region
! calculate the second end point relative to the first end point
x2=x0+ax*l2
y2=y0+ay*l2
y2=y2-i*ly ! subtract the same number of periods in y as the first point
! limit the line to lx
if (x2.GT.lx) then
c=y2-m*x2
x2=lx
y2=m*x2+c
line_end=.TRUE.
end if
l_line=sqrt( (x2-x1)**2+(y2-y1)**2 ) ! length of the line to be plotted
! write the first point
write(unit,'(2ES16.4)')x1,y1
if ( (y2.LE.ly).AND.(y2.GE.0d0) ) then
! write second point of line
write(unit,'(2ES16.4)')x2,y2
write(unit,*)
else if (y2.GT.ly) then
y2=ly
c=y1-m*x1
x2=(y2-c)/m
l_segment=sqrt( (x2-x1)**2+(y2-y1)**2 )
! write second point of split line
write(unit,*)'# split line, point 1, l_segment=',l_segment
write(unit,'(2ES16.4)')x2,y2
write(unit,*)
x1=x2
y1=y2-ly
! calculate the second end point relative to the new starting point
x2=x1+ax*(l_line-l_segment)
y2=y1+ay*(l_line-l_segment)
write(unit,*)'# split line, point 2,3'
write(unit,'(2ES16.4)')x1,y1
write(unit,'(2ES16.4)')x2,y2
write(unit,*)
else if (y2.LT.0d0) then
y2=0d0
c=y1-m*x1
x2=(y2-c)/m
l_segment=sqrt( (x2-x1)**2+(y2-y1)**2 )
! write second point of split line
write(unit,*)'# -split line, point 1, l_segment=',l_segment
write(unit,'(2ES16.4)')x2,y2
write(unit,*)
x1=x2
y1=y2+ly
! calculate the second end point relative to the first end point
x2=x1+ax*(l_line-l_segment)
y2=y1+ay*(l_line-l_segment)
write(unit,*)'# -split line, point 2,3'
write(unit,'(2ES16.4)')x1,y1
write(unit,'(2ES16.4)')x2,y2
write(unit,*)
end if
RETURN
END SUBROUTINE write_braid_line_l