!
! 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:
! SUBROUTINE calc_incident_field_components
! SUBROUTINE calculate_incident_field_sources
! SUBROUTINE calculate_ZT_incident_field_sources
! SUBROUTINE calc_incident_field_FD_values
! SUBROUTINE calc_incident_field_FD_values_GP
! SUBROUTINE calculate_lumped_incident_field_sources
!
! NAME
!     calc_incident_field_components
!
! AUTHORS
!     Chris Smartt
!
! DESCRIPTION
!     This subroutine calculates the cartesian components of the incident field
!     excitation and the propagation (k) vector from the inputs which are in
!     a spherical coordinate system ! #FIGURE_REFERENCE_REQUIRED
!     
! COMMENTS
!     
!
! HISTORY
!
!     started 13/6/2016 CJS: STAGE_6 developments
!     8/5/2017         CJS: Include references to Theory_Manual
!
SUBROUTINE calc_incident_field_components(Ktheta,Kphi,Etheta,Ephi, &
                                          Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz)                                  

USE type_specifications
USE general_module
USE constants
USE eispack
USE maths

IMPLICIT NONE

! variables passed to the subroutine

real(dp),intent(IN)    :: ktheta        ! theta component of k vector
real(dp),intent(IN)    :: kphi          ! phi component of k vector
real(dp),intent(IN)    :: etheta        ! theta component of E field polarisation vector
real(dp),intent(IN)    :: ephi          ! phi component of E field polarisation vector

! Cartesian field components
real(dp),intent(OUT)   :: Ex
real(dp),intent(OUT)   :: Ey
real(dp),intent(OUT)   :: Ez
real(dp),intent(OUT)   :: Hx
real(dp),intent(OUT)   :: Hy
real(dp),intent(OUT)   :: Hz

! Cartesian propagation vector components
real(dp),intent(OUT)   :: kx
real(dp),intent(OUT)   :: ky
real(dp),intent(OUT)   :: kz  

! local variables

real(dp)    :: r,theta,phi

real(dp)    :: x,y,z

real(dp)    :: dxdt,dxdp
real(dp)    :: dydt,dydp
real(dp)    :: dzdt,dzdp
real(dp)    :: vx,vy,vz
real(dp)    :: norm

! START

! calculate unit length incident wave vector Theory_Manual_Eqn 3.135

  r=1d0
  theta=Ktheta
  phi=Kphi
    
  kx=r*sin(theta)*cos(phi)
  ky=r*sin(theta)*sin(phi)
  kz=r*cos(theta)
    
! Calculate cartesian unit vectors in the theta and phi directions.   Theory_Manual_Eqn 3.137, 3.138
    
  dxdt= r*cos(theta)*cos(phi)
  dxdp=-r*sin(phi)
  dydt= r*cos(theta)*sin(phi)
  dydp= r*cos(phi)
  dzdt=-r*sin(theta)
  dzdp=0d0

! Calculate the cartesian E field vector from the sum of Etheta and Ephi contributions.  Theory_Manual_Eqn 3.139
  vx=Etheta*dxdt+Ephi*dxdp
  vy=Etheta*dydt+Ephi*dydp
  vz=Etheta*dzdt+Ephi*dzdp

! Normalise the E field vector to a length of 1    
  norm=sqrt(vx*vx+vy*vy+vz*vz)
    
  if (norm.eq.0d0) then
    if (verbose) then
      write(*,*)'Error calculating incident field data'
      write(*,8000)'E vector(theta)=',dxdt,dydt,dzdt
      write(*,8000)'E vector(phi)=',dxdp,dydp,dzdp
      write(*,8000)'E vector=',Etheta,Ephi
    end if
    run_status='ERROR in calc_incident_field_components: Zero length incident E field vector'
    CALL write_program_status()
    STOP 1
  end if

  Ex=vx/norm
  Ey=vy/norm
  Ez=vz/norm   
     
! H = (K x E)/Z0.  Theory_Manual_Eqn 3.142

  Hx=(ky*Ez-kz*Ey)/Z0
  Hy=(kz*Ex-kx*Ez)/Z0
  Hz=(kx*Ey-ky*Ex)/Z0
  
  if (verbose) then
  
    write(*,*)'Incident_field_components:'
    write(*,8000)'K=',kx,ky,kz
    write(*,8000)'E=',Ex,Ey,Ez
    write(*,8000)'H=',Hx,Hy,Hz

8000 format(A,3E16.6)
  end if
    
  RETURN

END SUBROUTINE calc_incident_field_components
! 
! NAME
!     calculate_incident_field_sources
!
! DESCRIPTION
!     calculate constants required to determine the controlled sources 
!     for Spice models of transient incident field illumination using the
!     method of characterisitcs based solution
!     
!
! COMMENTS
!     Refer to Theory_Manual_Section 3.7
!     Also discussed in more depth in C.R. Paul, "Analysis of Multiconductor Transmission Lines" 
!     1st edition, section 7.3.2 
! 
!
! HISTORY
!
!     started 16/06/2016 CJS
!     28/6/2016 CJS :Include the presence of a ground plane if required.
!     8/5/2017         CJS: Include references to Theory_Manual
!
!
SUBROUTINE calculate_incident_field_sources(xcoord,ycoord,Eamplitude,Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz,                   &
                                            ground_plane_present,ground_plane_x,ground_plane_y,ground_plane_theta, &
                                            length,n_modes,TI,Tdi,Tdz,alpha_0,alpha_L)
USE type_specifications
USE constants
USE general_module

IMPLICIT NONE

! variables passed to subroutine

integer,intent(IN)   :: n_modes                               ! input: number of modes in the illuminated domain
real(dp),intent(IN)  :: xcoord(n_modes+1),ycoord(n_modes+1)   ! input: centre coordinates of conductors
real(dp),intent(IN)  :: Eamplitude,Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz ! input: incident field specification
logical ,intent(IN)  :: ground_plane_present                  ! input: flag indicating the presence of a ground plane
real(dp),intent(IN)  :: ground_plane_x,ground_plane_y         ! input: x and y coordinates of a point on the ground plane
real(dp),intent(IN)  :: ground_plane_theta                    ! input: angle of the ground plane from the x axis
real(dp),intent(IN)  :: length                                ! input: length of transmission line
real(dp),intent(IN)  :: TI(n_modes,n_modes)                   ! input: mode decomposition matrix
real(dp),intent(IN)  :: Tdi(n_modes)                          ! input: mode delays
real(dp),intent(OUT)  :: alpha_0(n_modes),alpha_L(n_modes)    ! output: constants for controlled sources at z=0 and z=L
real(dp),intent(OUT)  :: Tdz                                  ! output: incident field delay time


! local variables

integer  :: n_conductors

real(dp) :: xk,yk
real(dp) :: rvx,rvy,rvz
real(dp) :: Tdxyk(n_modes)         

real(dp) :: vi

integer :: i,k

! START
  
  n_conductors=n_modes+1

  alpha_0(:)=0d0
  alpha_L(:)=0d0

  if (Eamplitude.eq.0d0) then  ! there is no incident field so return zero constants
    return
  end if

! calculate 1/velocity in each direction (avoids divide by zeros using reciprocal values)
! obtained from Theory_Manual_Eqn 3.136
  rvx=kx/c0
  rvy=ky/c0
  rvz=kz/c0

! Time delay in Z direction
  Tdz=length*rvz
  
! loop over all conductors except the reference conductor

  do k=1,n_conductors-1
    
! Assume that the final (reference) conductor lies at the origin of the xy plane
! and calculate the vector xk,yk from the reference conductor to conductor k.
! If there is a ground plane then calculate the vector to the closest point on the ground plane
! i.e. in the y direction since the ground plane is along the x axis.

    if (ground_plane_present) then
      xk=0d0
    else
      xk=xcoord(k)-xcoord(n_modes+1)
    end if
    yk=ycoord(k)-ycoord(n_modes+1)
    

! work out the delay propagating from this conductor to the reference      
    Tdxyk(k)=xk*rvx+yk*rvy
        
  end do
  
! loop over modes, i

  do i=1,n_modes
      
! loop over all conductors except the reference conductor

    do k=1,n_conductors-1
    
! Again, assume that the final (reference) conductor lies at the origin of the xy plane
! and calculate the vector xk,yk from the reference conductor to conductor k.
! If there is a ground plane then calculate the vector to the closest point on the ground plane
! i.e. in the y direction since the ground plane is along the x axis.
      
      if (ground_plane_present) then
! only the height over the ground plane matters so set x coordinate to zero and y is the height over ground
        xk=0d0
      else
        xk=xcoord(k)-xcoord(n_modes+1)
      end if
      
      yk=ycoord(k)-ycoord(n_modes+1)

! note alpha_0,alpha_L are used to calculate the response of mode i due to excitations on all conductors, k
! hence the summation with over k with coefficients TI(k,i)
! Theory_Manual_Eqn 3.167
      
      alpha_0(i)=alpha_0(i)+( Ez*Tdxyk(k)*length - (Ex*xk+Ey*yk)*(Tdi(i)+Tdz) )*TI(k,i)  ! note, TI(k,i) = TI(conductor,mode)
      alpha_L(i)=alpha_L(i)+( Ez*Tdxyk(k)*length + (Ex*xk+Ey*yk)*(Tdi(i)-Tdz) )*TI(k,i) 
               
    end do ! next conductor
  
  end do ! next mode
  
  if (ground_plane_present) then ! See Theory_Manual_Section 3.7.2
! the alpha terms are doubled
    alpha_0(:)=2d0*alpha_0(:)
    alpha_L(:)=2d0*alpha_L(:)
  end if
  
  
  if (.NOT.verbose) RETURN
  
  write(*,*)'calculate_incident_field_sources'
  write(*,*)'ex=',ex
  write(*,*)'ey=',ey
  write(*,*)'ez=',ez
  write(*,*)'kx=',kx
  write(*,*)'ky=',ky
  write(*,*)'kz=',kz
  write(*,*)'rvx=',rvx
  write(*,*)'rvy=',rvy
  write(*,*)'rvz=',rvz
  
  write(*,*)'Tdz=',Tdz
   
  RETURN
   
END SUBROUTINE calculate_incident_field_sources
! 
! NAME
!     calculate_ZT_incident_field_sources
!
! DESCRIPTION
!     calculate the controlled sources for Spice modelling of transient incident field illumination 
!     of shielded cables (Xie's model), method of characterisitcs based solution
!     
!
! COMMENTS
!     Refer to Theory_Manual_Section 3.8
!
!
! HISTORY
!
!     started 3/10/2016 CJS 
!     8/5/2017         CJS: Include references to Theory_Manual
!
!
SUBROUTINE calculate_ZT_incident_field_sources(xcoord,ycoord,Eamplitude,Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz,                   &
                                              ground_plane_present,ground_plane_x,ground_plane_y,ground_plane_theta, &
                                              length,n_modes,shield_conductor,TI,Tdi,Tdz,c1,c2,c1b)
USE type_specifications
USE constants
USE general_module

IMPLICIT NONE

! variables passed to subroutine

integer,intent(IN)   :: n_modes                               ! input: number of modes in the illuminated domain
real(dp),intent(IN)  :: xcoord(n_modes+1),ycoord(n_modes+1)   ! input: centre coordinates of conductors
real(dp),intent(IN)  :: Eamplitude,Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz ! input: incident field specification
logical,intent(IN)   :: ground_plane_present                  ! input: flag indicating the presence of a ground plane
real(dp),intent(IN)  :: ground_plane_x,ground_plane_y         ! input: x and y coordinates of a point on the ground plane
real(dp),intent(IN)  :: ground_plane_theta                    ! input: angle of the ground plane from the x axis
real(dp),intent(IN)  :: length                                ! input: length of transmission line
integer,intent(IN)   :: shield_conductor                      ! input: the number of the shield conductor
real(dp),intent(IN)  :: TI(n_modes,n_modes)                   ! input: mode decomposition matrix
real(dp),intent(IN)  :: Tdi(n_modes)                          ! input: mode delays
real(dp),intent(OUT)  :: c1(n_modes),c2(n_modes),c1b(n_modes) ! output: constants for controlled sources in shielded domains
real(dp),intent(OUT)  :: Tdz                                  ! output: incident field delay time

! local variables

integer  :: n_conductors

real(dp) :: xk,yk
real(dp) :: rvx,rvy,rvz

real(dp) :: Tx(n_modes),dk,kw_rx,kw_ry,kw_r,vw_r  

real(dp) :: vi

real(dp) :: Tsum

! loop variables
integer :: i,k

! START

  if(verbose) write(*,*)'CALLED: calculate_ZT_incident_field_sources'
  
  n_conductors=n_modes+1
 
  if(verbose) then
    write(*,*)'n_modes=',n_modes
    write(*,*)'shield_conductor=',shield_conductor
  end if

! Zero the coefficients initially 
  c1(:)=0d0
  c2(:)=0d0
  c1b(:)=0d0

  if (Eamplitude.eq.0d0) then  ! there is no incident field so return zero constants
    return
  end if

! calculate 1/velocity in each direction (avoids divide by zeros using reciprocal values)
! obtained from Theory_Manual_Eqn 3.136
  rvx=kx/c0
  rvy=ky/c0
  rvz=kz/c0

  Tdz=length*rvz
  
! loop over all conductors except the reference conductor

  do k=1,n_conductors-1
    
! Assume that the final (reference) conductor lies at the origin of the xy plane
! and calculate the vector xk,yk from the reference conductor to conductor k.
! If there is a ground plane then calculate the vector to the closest point on the ground plane
! i.e. in the y direction since the ground plane is along the x axis.

    if (ground_plane_present) then
! only the height over the ground plane matters so set x coordinate to zero and y is the height over ground
      xk=0d0
    else
      xk=xcoord(k)-xcoord(n_modes+1)
    end if
    yk=ycoord(k)-ycoord(n_modes+1)
        
! distance from this conductor to the reference conductor
    dk=sqrt(xk*xk+yk*yk)
    
! unit vector from this conductor to the reference conductor
    kw_rx=xk/dk
    kw_ry=yk/dk
    
! project the k vector along this direction
    kw_r=kx*kw_rx+ky*kw_ry
    
! reciprocal velocity in this direction
    vw_r=kw_r/c0
    
! Tx definition in Xie paper - see also Theory_Manual_Section 3.8 comment in the text below equation 3.198
! where here we have generalised to allow for general conductor offsets in the x-y plane of the bundle
    Tx(k)=length*vw_r  
    
  end do
  
! Calculate C1(i), C2(i) , C1b(i) where i is the source domain mode.
! We only have contributions relating to the shield conductor whose transfer impedance we are 
! including in the model. 

! loop over modes, i in the source domain

  do i=1,n_modes
    
! loop over all conductors except the reference conductor

    do k=1,n_conductors-1
    
    
! Again, assume that the final (reference) conductor lies at the origin of the xy plane
! and calculate the vector xk,yk from the reference conductor to conductor k.
! If there is a ground plane then calculate the vector to the closest point on the ground plane
! i.e. in the y direction since the ground plane is along the x axis.
      if (ground_plane_present) then
        xk=0d0
      else
        xk=xcoord(k)-xcoord(n_modes+1)
      end if
    
      yk=ycoord(k)-ycoord(n_modes+1)
    
      dk=sqrt(xk*xk+yk*yk)

! note c1,c2 are used to calculate the response of mode i due to excitations on all conductors, k
! hence the summation with over k with coefficients TI(k,i)
! See the two coefficients for the delayed incident field in Theory_Manual_Eqn 3.195 to 3.198
    
      c1(i)=c1(i)+dk*length*(-Ez*Tx(k) +(Ex*xk/dk+Ey*yk/dk)*(Tdz-Tdi(i)) )*TI(k,i) 
      c2(i)=c2(i)+dk*length*(-Ez*Tx(k) +(Ex*xk/dk+Ey*yk/dk)*(Tdz+Tdi(i)) )*TI(k,i) 

! C1b is for the special case when Tdz=Tdi, see the last paragraph in Theory_Manual_Section 3.8.3
      c1b(i)=c1b(i)+dk*length*( Ex*xk/dk+Ey*yk/dk )*TI(k,i) 
       
    end do ! next conductor
          
  end do ! next mode
  
  if (ground_plane_present) then  ! See Theory_Manual_Section 3.7.2
! the alpha terms are doubled
    c1(:)=2d0*c1(:)         
    c2(:)=2d0*c2(:)
    c1b(:)=2d0*c1b(:)         
  end if
  
  
  if (.NOT.verbose) RETURN
  
  write(*,*)'calculate_incident_field_sources'
  write(*,*)'ex=',ex
  write(*,*)'ey=',ey
  write(*,*)'ez=',ez
  write(*,*)'kx=',kx
  write(*,*)'ky=',ky
  write(*,*)'kz=',kz
  write(*,*)'rvx=',rvx
  write(*,*)'rvy=',rvy
  write(*,*)'rvz=',rvz
  
  write(*,*)'Tdz=',Tdz
  
! loop over conductors

  do k=1,n_modes
  
    xk=xcoord(k)-xcoord(n_modes+1)
    yk=ycoord(k)-ycoord(n_modes+1)
    
! distance from this conductor to the reference conductor
    dk=sqrt(xk*xk+yk*yk)
    
! unit vector from this conductor to the reference conductor
    kw_rx=xk/dk
    kw_ry=yk/dk
    
! project the k vector along this direction
    kw_r=kx*kw_rx+ky*kw_ry
          
    write(*,*)'conductor ',k
    write(*,*)'xk',xk
    write(*,*)'yk',yk
    write(*,*)'dk=',dk
    write(*,*)'kw_rx=',kw_rx
    write(*,*)'kw_ry=',kw_ry
    write(*,*)'kw_r=',kw_r
    write(*,*)'Tx=',Tx(k)
    
  end do
  
! loop over modes, i

  do i=1,n_modes
   
    write(*,*)'mode ',i
    write(*,*)'Tdi=',Tdi(i)
    write(*,*)'Tdi(i)+Tdz=',Tdi(i)+Tdz
    write(*,*)'c1=',c1(i)
    write(*,*)'c2=',c2(i)
    write(*,*)'c1b=',c1b(i)
   
  end do 
  
  RETURN
   
END SUBROUTINE calculate_ZT_incident_field_sources
!
! NAME
!     calc_incident_field_values
!
! DESCRIPTION
!     given coordinates in space and frequency calculate the incident E and H field component values
!     
! COMMENTS
!     
!
! HISTORY
!
!     started 29/04/2015 CJS
!     8/5/2017           CJS: Include references to Theory_Manual
!
!
SUBROUTINE calc_incident_field_FD_values(x,y,z,f,Eamplitude,Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz, &
                                         Ex_value,Ey_value,Ez_value,Hx_value,Hy_value,Hz_value)

USE type_specifications
USE constants

IMPLICIT NONE

complex(dp),intent(IN) :: Eamplitude                    ! E field amplitude
real(dp),intent(IN)    :: Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz    ! Cartesian components of field and propgation vector
real(dp),intent(IN)    :: x,y,z                         ! coordinates in 3D space of field evaluation point
real(dp),intent(IN)    :: f                             ! frequency

complex(dp),intent(OUT) :: Ex_value,Ey_value,Ez_value,Hx_value,Hy_value,Hz_value  ! output cartesian field values

! local variables

real(dp) :: w            ! angular frequency

complex(dp) :: phase     ! phase at field evaluation point

! START

! Theory_Manual_Eqn 3.140

w=2d0*pi*f

phase=exp( -j*(w/c0)*( kx*x+ky*y+kz*z ) )

Ex_value=Eamplitude*Ex*phase
Ey_value=Eamplitude*Ey*phase
Ez_value=Eamplitude*Ez*phase

! Theory_Manual_Eqn 3.140 and 3.142

Hx_value=Eamplitude*Hx*phase
Hy_value=Eamplitude*Hy*phase
Hz_value=Eamplitude*Hz*phase

RETURN

END SUBROUTINE calc_incident_field_FD_values
!
! NAME
!     calc_incident_field_FD_values_GP
!
! DESCRIPTION
!     given coordinates in space and frequency calculate the incident E and H field component values
!     with a ground plane present
!     The ground plane is assumed to lie along the x axis
!
! COMMENTS
!     
!
! HISTORY
!
!     started 28/06/2016 CJS
!
!
SUBROUTINE calc_incident_field_FD_values_GP(x,y,z,f,Eamplitude,Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz, &
                                         Ex_value,Ey_value,Ez_value,Hx_value,Hy_value,Hz_value)

USE type_specifications
USE constants

IMPLICIT NONE

complex(dp),intent(IN) :: Eamplitude                    ! E field amplitude
real(dp),intent(IN)    :: Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz    ! Cartesian components of field and propgation vector
real(dp),intent(IN)    :: x,y,z                         ! coordinates in 3D space of field evaluation point
real(dp),intent(IN)    :: f                             ! frequency

complex(dp),intent(OUT) :: Ex_value,Ey_value,Ez_value,Hx_value,Hy_value,Hz_value  ! output cartesian field values

! local variables

real(dp) :: w            ! angular frequency

complex(dp) :: phase     ! phase at field evaluation point

! START

w=2d0*pi*f

! Theory_Manual_Eqn 3.140 with the addition of a reflected plane wave from the ground plane (y=0 plane)
! in which ky_reflected=-ky_incident and
! Ex_reflected=-Ex_incident,  Ez_reflected=-Ez_incident,  Hy_reflected=-Hy_incident

! Incident plane wave
phase=exp( -j*(w/c0)*( kx*x+ky*y+kz*z ) )

Ex_value=Eamplitude*Ex*phase
Ey_value=Eamplitude*Ey*phase
Ez_value=Eamplitude*Ez*phase

Hx_value=Eamplitude*Hx*phase
Hy_value=Eamplitude*Hy*phase
Hz_value=Eamplitude*Hz*phase

! Add reflected plane wave
phase=exp( -j*(w/c0)*(kx*x-ky*y+kz*z ) )

Ex_value=Ex_value-Eamplitude*Ex*phase
Ey_value=Ey_value+Eamplitude*Ey*phase
Ez_value=Ez_value-Eamplitude*Ez*phase

Hx_value=Hx_value+Eamplitude*Hx*phase
Hy_value=Hy_value-Eamplitude*Hy*phase
Hz_value=Hz_value+Eamplitude*Hz*phase

RETURN

END SUBROUTINE calc_incident_field_FD_values_GP
! 
! NAME
!     calculate_lumped_incident_field_sources
!
! DESCRIPTION
!     calculate the lumped sources at z=l due to incident field excitation
!     
!
! COMMENTS
!     The code is based on Theory_Manual_Section 2.4
!     Detailed discussions and derivations are found in 
!     C.R. Paul, "Analysis of Multiconductor Transmission Lines" 1st edition, section 7.2.1. 
!     
!     The sources are calculated by numerical integration
!    
! HISTORY
!
!     started 29/04/2015 CJS ESA Spice cable modelling, Phase 1 
!     8/5/2017           CJS: Include references to Theory_Manual
!
!
SUBROUTINE calculate_lumped_incident_field_sources(xcoord,ycoord,is_shielded,Eamplitude,Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz, &
                                                   ground_plane_present,length,f,VFT,IFT,dim,  &
                                                   T,TI,Y,Z,ZC,YC,GAMMA_C)

USE type_specifications
USE constants
USE maths

IMPLICIT NONE
integer,intent(IN)     :: dim                              ! matrix system dimension
real(dp),intent(IN)    :: xcoord(dim+1),ycoord(dim+1)      ! note that we need the reference conductor position in space
logical,intent(IN)     :: is_shielded(dim)                 ! flag to indicated shielded conductors
complex(dp),intent(IN) :: Eamplitude                       ! incident field amplitude
real(dp),intent(IN)    :: Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz       ! incident field cartesian field components and propagation vector
real(dp),intent(IN)    :: length                           ! transmission line length 
real(dp),intent(IN)    :: f                                ! frequency
logical,intent(IN)     :: ground_plane_present             ! flag to indicate the presence of a ground plane
complex(dp),intent(IN) :: T(dim,dim)                       ! modal decomposition matrix
complex(dp),intent(IN) :: TI(dim,dim)                      ! modal decomposition matrix (inverse of [T])
complex(dp),intent(IN) :: Y(dim,dim)                       ! global admittance matrix
complex(dp),intent(IN) :: Z(dim,dim)                       ! global impedance matrix 
complex(dp),intent(IN) :: ZC(dim,dim)                      ! characteristic impedance matrix
complex(dp),intent(IN) :: YC(dim,dim)                      ! characteristic admittance matrix
complex(dp),intent(IN) :: GAMMA_C(dim)                     ! modal propagation constants

complex(dp),intent(OUT) :: VFT(dim),IFT(dim)               ! incident field voltage and current sources

! local variables

real(dp) :: xp,yp,zp
real(dp) :: vx,vy,norm

real(dp) :: w,lambda

real(dp) :: dz          ! space step for numerical integration
integer  :: nz,zloop    ! loop variables for integration along the transmission line

! Cartesian field values on a conductor
complex(dp) :: Ex_c,Ey_c,Ez_c
complex(dp) :: Hx_c,Hy_c,Hz_c

! Cartesian field values on the reference conductor
complex(dp) :: Ex_r,Ey_r,Ez_r
complex(dp) :: Hx_r,Hy_r,Hz_r

! Temporary matrices for calculating the chain parameter matrices
complex(dp) :: M1(dim,dim),YI(dim,dim)
complex(dp) :: sinh_gamma_z(dim,dim)
complex(dp) :: cosh_gamma_z(dim,dim)
complex(dp) :: TSTI(dim,dim),TCTI(dim,dim)

! chain parameter matrix elements
complex(dp) :: PHI_11(dim,dim)
complex(dp) :: PHI_12(dim,dim)
complex(dp) :: PHI_21(dim,dim)
complex(dp) :: PHI_22(dim,dim)

!temporary vectors used in the calculations
complex(dp) :: VF1(dim),VF2(dim),VF3(dim),VF(dim),IF1(dim),IF(dim)
complex(dp) :: TV1(dim),TV2(dim),TV3(dim),TV4(dim)

! error code for matrix inverse 
integer :: ierr   

! loop variables
integer :: conductor
integer :: i

! START

  IFT(:)=(0d0,0d0)
  VFT(:)=(0d0,0d0)

  if (Eamplitude.eq.0d0) then   ! there is no incident field so return with zero value source terms
    return
  end if

! calculate angular frequency(rad/s) and wavelength(m)
  w=2d0*pi*f
  lambda=c0/f
  
  ierr=0      ! exit on error with the matrix inverse
  CALL cinvert_Gauss_Jordan(Y,dim,YI,dim,ierr)
  
  dz=lambda/20d0   ! space step for numerical integration
  nz=NINT(length/dz)
  
! ensure a minimum number of integration points
  if (nz.LT.10) nz=10
  dz=length/nz
  
! calculate source integrals   
   
! loop over the transmission line length
  do zloop=1,nz
  
    zp=(zloop-1)*dz+dz/2d0    ! z coordinate at the current point. We assume the integrand is uniform over intervals, dz.
    
! Calculate the chain parameter matrices for the current z value Theory_Manual_Eqns 2.58

    sinh_gamma_z(:,:)=(0d0,0d0)
    cosh_gamma_z(:,:)=(0d0,0d0)

    do i=1,dim
      cosh_gamma_z(i,i)=0.5d0*( exp(  GAMMA_C(i)*(length-zp) )+exp( -GAMMA_C(i)*(length-zp) ) )
      sinh_gamma_z(i,i)=0.5d0*( exp(  GAMMA_C(i)*(length-zp) )-exp( -GAMMA_C(i)*(length-zp) ) )
    end do
       
    M1=matmul(T,cosh_gamma_z)
    TCTI=matmul(M1,TI)

    M1=matmul(T,sinh_gamma_z)
    TSTI=matmul(M1,TI)   

! Phi11
    M1=matmul(YI,TCTI)
    PHI_11=matmul(M1,Y)
    
! Phi12
    M1=matmul(ZC,TSTI)
    PHI_12=-M1
    
! Phi21
    M1=matmul(TSTI,YC)
    PHI_21=-M1
      
! Phi22
    PHI_22=TCTI  
            
! loop over conductors (other than the reference)
    do conductor=1,dim
          
! calculate the cartesian field values at the current conductor and the reference conductor
          
      if (.NOT.ground_plane_present) then
                                           
! get the field values at the reference conductor i.e. the (dim+1)th conductor 
        xp=xcoord(dim+1)
        yp=ycoord(dim+1)
    
        CALL calc_incident_field_FD_values(xp,yp,zp,f,Eamplitude,Ex,Ey,Ez,Hx,Hy,Hz, &
                                           kx,ky,kz,Ex_r,Ey_r,Ez_r,Hx_r,Hy_r,Hz_r)
      
! get the field values at the conductor 
        xp=xcoord(conductor)
        yp=ycoord(conductor)
      
        CALL calc_incident_field_FD_values(xp,yp,zp,f,Eamplitude,Ex,Ey,Ez,Hx,Hy,Hz, &
                                           kx,ky,kz,Ex_c,Ey_c,Ez_c,Hx_c,Hy_c,Hz_c)
                                           
! calculate a vector from the reference to the conductor
        vx=xp-xcoord(dim+1)
        vy=yp-ycoord(dim+1)
                                           
      else
      
! get the field values at the reference conductor i.e. the (dim+1)th conductor 
! under the conductor we are interested in
        xp=xcoord(dim+1)     ! The 'reference' is at the origin when a ground plane is present 
        yp=ycoord(dim+1)
        CALL calc_incident_field_FD_values_GP(xp,yp,zp,f,Eamplitude,Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz, &
                                              Ex_r,Ey_r,Ez_r,Hx_r,Hy_r,Hz_r)
      
! get the field values at the conductor 
        xp=xcoord(conductor)
        yp=ycoord(conductor)
        CALL calc_incident_field_FD_values_GP(xp,yp,zp,f,Eamplitude,Ex,Ey,Ez,Hx,Hy,Hz,kx,ky,kz, &
                                              Ex_c,Ey_c,Ez_c,Hx_c,Hy_c,Hz_c)
                                                         
! calculate a vector from the reference to the conductor
        vx=xp-xcoord(dim+1)  ! Made the vector consistent with the field evaluation points 16/8/2016
        vy=yp-ycoord(dim+1)  ! The 'reference' is at the origin when a ground plane is present 
    
      end if
            
! the first part is due to the z component of the E field

      VF1(conductor)=(Ez_c-Ez_r)    ! C. Paul, equation 7.24(e), Ez contribution !  Theory_Manual_Eqn 2.53(second term) using 2.54
            
! the second part is due to the z derivative of the incident field integrated from the reference to conductor
! note the field is assumed to be linearly varying between its endpoint values in the integration

      VF2(conductor)=j*(w/c0)*kz*(0.5d0*(Ex_c+Ex_r)*vx+0.5d0*(Ey_c+Ey_r)*vy) ! C. Paul, equation 7.24(e), Et contribution, Theory_Manual_Eqn 2.53, first term
      
! An alternative is the time derivative of the B field integrated from reference to conductor
!      VF3(conductor)=j*w*mu0*(-0.5d0*(Hx_c+Hx_r)*vy+0.5d0*(Hy_c+Hy_r)*vx) ! C. Paul, equation 7.24(c), B.dA contribution (H field form) 
      
! the IF source term is due to the incident field integrated from the reference to conductor
! note the field is assumed to be linearly varying between its endpoint values in the integration
! this needs to be pre-multiplied by G+jwC but this happens later...
 ! C. Paul, equation 7.24(d), Et contribution, Theory_Manual_Eqn 2.55 using a linear interpolation of the field
 ! between the reference conductor and the conductor k.
       
      IF1(conductor)=-(0.5d0*(Ex_c+Ex_r)*vx+0.5d0*(Ey_c+Ey_r)*vy)
                                
    end do ! next conductor
    
    VF(:)=VF1(:)+VF2(:)  ! E field form
!     VF(:)=VF3(:)       ! H field form


! multiply IF1 by Y, the MTL admittance matrix (G+jwC) ! Theory_Manual_Eqn 2.56 transformed to the time domain
     IF=matmul(Y,IF1)

! Special process for shielded conductors - set IF to zero
    do conductor=1,dim
    
      if (is_shielded(conductor)) then
      
! set the current sources to zero, note that voltage sources must remain even for shielded conductors. 
! Note also that provided shielded cables are centred on their respective shield then this step should
! not be necessary due to the structure of the [Y] matrix and (IF)=[Y](IF1)
        
        IF(conductor)=(0d0,0d0)
        
       end if
                                
     end do ! next conductor
     
! add the contribution to the integration in Theory_Manual_Eqns 2.60
! C. Paul, equation 7.29b,7.29c      
     TV1=matmul(PHI_11,VF)
     TV2=matmul(PHI_12,IF)
    
     TV3=matmul(PHI_21,VF)
     TV4=matmul(PHI_22,IF)

! Add contribution to the integral over z assuming constant value of integrand over this section, dz
    do conductor=1,dim
      VFT(conductor)=VFT(conductor)+dz*(TV1(conductor)+TV2(conductor))  
      IFT(conductor)=IFT(conductor)+dz*(TV3(conductor)+TV4(conductor)) 
    end do
        
  end do ! next z value
    
END SUBROUTINE calculate_lumped_incident_field_sources