!
! 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:
! 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