!
! 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
!
! HISTORY
!
!     started 22/08/2016 CJS 
!
!
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

! 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

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

! START

! Open the input file describing the braid parameters
! This file could be created by the associated GUI or otherwise generated

  verbose=.TRUE.

  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))
  
  R0=4d0/(pi*dw*dw*N*C*sigma*cos(alpha))    ! d.c. resistance. equation 5.63 of D1
  
  S=pi*wh*lh/4d0   ! hole area
  req=sqrt(S/pi)
  kappa=1.84d0/req  

! Calculate an equivalent shield thickness from R0, D and sigma

  T=1d0/(2d0*pi*sigma*(D/2d0)*R0)
  
  Ck=0.875d0*exp(-j*kappa*T)

  Dm=D+2d0*dw
  b=(2d0*pi*Dm/C)*cos(alpha)-N*dw
  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))*((1-K)**(3d0/3d0))*(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

! 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'
    write(*,*)'req=',req,'  hole equivalent radius'
    write(*,*)'k=',kappa,'  hole cutoff k value'
    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
  
  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)

  end do ! next frequency
  
  close(unit=83)

! 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

! call 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'
  CALL Calculate_Sfilter(Zd,frequency_data%freq_list,frequency_data%n_frequencies,Zd_filter,order,1,0) ! call with fit_type=0
  
! 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'
  CALL Calculate_Sfilter(Zc,frequency_data%freq_list,frequency_data%n_frequencies,Zc_filter,order,0,0) ! call with fit_type=0

! 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,*)'# Transfer impedance model'
  CALL Write_Sfilter(Zt_filter,shield_model_file_unit)

! Write the shield conductor model to the shield model file

  write(shield_model_file_unit,*)'Conductor surface impedance model'
  CALL Write_Sfilter(Zc_filter,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'
  write(shield_model_file_unit,*)'req=',req,'  hole equivalent radius'
  write(shield_model_file_unit,*)'k=',kappa,'  hole cutoff k value'
  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