!
! 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 PUL_LC_calc_wide_separation_approximation
!     SUBROUTINE PUL_LC_calc_overshield_wide_separation_approximation
!     SUBROUTINE calculate_height_over_ground_plane(x,y,theta,d,h)
!
! NAME
!     SUBROUTINE PUL_LC_calc_wide_separation_approximation
!
!     Calculation of PUL_parameters for cylindrical conductors in free space
!     using a wide separation approximation
!     The calculation allows the inclusion of a ground plane in the system of conductors
! 
! The calculation is based on the theory in the Theory Manual however more detailed analysis and derivations
! can be found in
! C.R.Paul, "Analysis of Multiconductor Transmission Lines" John Wiley and Sons 
!     
! COMMENTS
!     
!
! HISTORY
!    started 2/12/2015 CJS
!     8/5/2017         CJS: Include references to Theory_Manual
!

SUBROUTINE PUL_LC_calc_wide_separation_approximation(PUL)

  USE type_specifications
  USE constants
  USE general_module
  USE maths

  IMPLICIT NONE

! variables passed to subroutine

  type(PUL_type),intent(INOUT)    :: PUL    ! per-unit-length parameter calculation structure
  
! local variables

  integer    :: nc                 ! number of conductors
  integer    :: matrix_dimension
  
  integer    :: row,col
  
  real(dp)    :: x_0,y_0       ! x and y of reference conductor
  real(dp)    :: x_row,y_row   ! x and y of row conductor
  real(dp)    :: x_col,y_col   ! x and y of col conductor
  
  real(dp)    :: d_row_0       ! distance from row conductor to reference conductor
  real(dp)    :: d_row_col     ! distance from row conductor to col conductor
  real(dp)    :: d_col_0       ! distance from col conductor to reference conductor
  
  real(dp)    :: r_row         ! radius of row conductor 
  real(dp)    :: r_col         ! radius of col conductor 
  real(dp)    :: r_0           ! radius of reference conductor
  
  real(dp)    :: h_row         ! height of row conductor above the ground plane
  real(dp)    :: h_col         ! height of col conductor above the ground plane
  
  integer :: ierr
  
! START

! Work out the dimension of the per-unit-length matrices  

  nc=PUL%n_conductors
  matrix_dimension=nc-1
  
  if (verbose) then
    write(*,*)'CALLED PUL_LC_calc_wide_separation_approximation'
    write(*,*)'Number of conductors=',nc
    write(*,*)'matrix dimension=',matrix_dimension
  end if
  
! check we have a valid system
  if (matrix_dimension.LT.1) then
    run_status='ERROR in PUL_LC_calc_wide_separation_approximation, Matrix dimension is less than 1'
    CALL write_program_status()
    STOP 1
  end if
  
! Allocate the per-unit-length matrices  
  PUL%L%dim=matrix_dimension
  ALLOCATE( PUL%L%mat(1:PUL%L%dim,1:PUL%L%dim) )
  PUL%C%dim=matrix_dimension
  ALLOCATE( PUL%C%mat(1:PUL%C%dim,1:PUL%C%dim) )
  PUL%G%dim=matrix_dimension
  ALLOCATE( PUL%G%mat(1:PUL%G%dim,1:PUL%G%dim) )

! Allocate the impedance and admittance filter matrices  
  PUL%Zfilter%dim=matrix_dimension
  ALLOCATE(PUL%Zfilter%sfilter_mat(1:PUL%Zfilter%dim,1:PUL%Zfilter%dim))
  PUL%Yfilter%dim=matrix_dimension
  ALLOCATE(PUL%Yfilter%sfilter_mat(1:PUL%Yfilter%dim,1:PUL%Yfilter%dim))

! Calculate the inductance matrix first

  if (PUL%ground_plane_present) then
  
    do row=1,nc-1
    
! calculate the diagonal element here    
      x_row=PUL%x(row)
      y_row=PUL%y(row)
      r_row=PUL%r(row)
      
      CALL calculate_height_over_ground_plane(x_row,y_row,r_row,    &
                                              PUL%ground_plane_angle,PUL%ground_plane_offset,h_row)
                                
      PUL%L%mat(row,row)=(mu0/(2d0*pi))*log(2d0*h_row/r_row)  ! Paul equation 3.66a  Theory_Manual_Eqn 2.24

      do col=row+1,nc-1
! calculate the off diagonal elements here    

        x_col=PUL%x(col)
        y_col=PUL%y(col)
        r_col=PUL%r(col)
        d_row_col=sqrt((x_row-x_col)**2+(y_row-y_col)**2)
      
        CALL calculate_height_over_ground_plane(x_col,y_col,r_col,   &
                                            PUL%ground_plane_angle,PUL%ground_plane_offset,h_col)

        PUL%L%mat(row,col)=(mu0/(4d0*pi))*log(1d0+4d0*h_row*h_col/(d_row_col**2))   ! C.R.Paul equation 3.66b  Theory_Manual_Eqn 2.25
    
! The matrix is symmetric so we can set L[col,row]=L[row,col]
        PUL%L%mat(col,row)=PUL%L%mat(row,col)
      
      end do ! next col
      
    end do ! next row
  
  else
! There is no ground plane so choose the last conductor specified as the reference
    
    x_0=PUL%x(nc)
    y_0=PUL%y(nc)
    r_0=PUL%r(nc)
  
    do row=1,nc-1
    
! calculate the diagonal element here    
      x_row=PUL%x(row)
      y_row=PUL%y(row)
      r_row=PUL%r(row)
      d_row_0=sqrt((x_row-x_0)**2+(y_row-y_0)**2)
      PUL%L%mat(row,row)=(mu0/(2d0*pi))*log(d_row_0**2/(r_0*r_row))  ! Paul equation 3.63a Theory_Manual_Eqn 2.21
    
      do col=row+1,nc-1
! calculate the off diagonal elements here    

        x_col=PUL%x(col)
        y_col=PUL%y(col)
        d_col_0=sqrt((x_col-x_0)**2+(y_col-y_0)**2)
        d_row_col=sqrt((x_row-x_col)**2+(y_row-y_col)**2)

        PUL%L%mat(row,col)=(mu0/(2d0*pi))*log(d_row_0*d_col_0/(d_row_col*r_0))   ! C.R.Paul equation 3.63b Theory_Manual_Eqn 2.22
    
! The matrix is symmetric so we can set L[col,row]=L[row,col]
        PUL%L%mat(col,row)=PUL%L%mat(row,col)
      
      end do ! next col
      
    end do ! next row
  
  end if
  
! The capacitance matrix is calculated as mu0*eps0*[L_inverse] ! Theory_Manual_Eqn 2.23
  
  ierr=0   ! set ierr=0 on input to matrix inverse to cause the program to stop if we have a singular matrix
  CALL dinvert_Gauss_Jordan(PUL%L%mat,matrix_dimension,PUL%C%mat,matrix_dimension,ierr)
  
  PUL%C%mat(:,:)=mu0*eps0*PUL%C%mat(:,:)

! Conductance matrix is set to zero
  PUL%G%mat(:,:)=0d0
  
! calculate the impedance and admittance filters

  CALL Z_Y_from_L_C(PUL%L,PUL%C,PUL%Zfilter,PUL%Yfilter)
   
  if (verbose) then
    write(*,*)'Inductance matrix, L'
    CALL dwrite_matrix(PUL%L%mat,matrix_dimension,matrix_dimension,matrix_dimension,0)
    write(*,*)'Capacitance matrix, C'
    CALL dwrite_matrix(PUL%C%mat,matrix_dimension,matrix_dimension,matrix_dimension,0)
  end if
  
  END SUBROUTINE PUL_LC_calc_wide_separation_approximation

!
! NAME
!     SUBROUTINE PUL_LC_calc_overshield_wide_separation_approximation
!
!     Calculation of PUL_parameters for cylindrical conductors within a cylindrical return conductor (overshield)
!     using a wide separation approximation
! 
! The calculation is based on the theory in:
! C.R.Paul, "Analysis of Multiconductor Transmission Lines" John Wiley and Sons
!     
! COMMENTS
!     
!
! HISTORY
!    started 18/04/2016 CJS
!    24/04/2016 CJS Fix problem with off diagonal inductance matrix calculation due to problem with cos_theta_row_col calculation.
!    26/10/2016 CJS Fix 0/0 error when d_col=0 for off diagonal elements
!     8/5/2017         CJS: Include references to Theory_Manual
!
SUBROUTINE PUL_LC_calc_overshield_wide_separation_approximation(PUL)

  USE type_specifications
  USE constants
  USE general_module
  USE maths

  IMPLICIT NONE

! variables passed to subroutine

  type(PUL_type),intent(INOUT)    :: PUL    ! per-unit-length parameter calculation structure
  
! local variables

  integer    :: nc                 ! number of conductors
  integer    :: matrix_dimension
  
  integer    :: row,col
  
  real(dp)    :: xs,ys         ! x and y of reference conductor
  real(dp)    :: rs            ! radius of reference conductor
  
  real(dp)    :: x_row,y_row   ! x and y of row conductor
  real(dp)    :: d_row         ! distance from row conductor to centre of reference conductor
  real(dp)    :: r_row         ! radius of row conductor 

  real(dp)    :: x_col,y_col   ! x and y of col conductor
  real(dp)    :: d_col         ! distance from col conductor to centre of reference conductor
  real(dp)    :: r_col         ! radius of col conductor 
  
  real(dp)    :: d_row_col         ! distance from row conductor to col conductor
  real(dp)    :: theta_row         ! angular of row conductor
  real(dp)    :: theta_col         ! angular of row conductor
  real(dp)    :: cos_theta_row_col ! angular separation of conductors
   
  real(dp)    :: num
  real(dp)    :: den
  
  integer :: ierr
    
! START

  if (verbose) then
    write(*,*)'CALLED: PUL_LC_calc_overshield_wide_separation_approximation'
    write(*,*)'n_conductors (including overshield)=',PUL%n_conductors
  end if

! Work out the dimension of the per-unit-length matrices  

  nc=PUL%n_conductors
  matrix_dimension=nc-1  ! the conductor list includes the oversheild
  
! Allocate the per-unit-length matrices  
  PUL%L%dim=matrix_dimension
  ALLOCATE( PUL%L%mat(1:PUL%L%dim,1:PUL%L%dim) )
  PUL%C%dim=matrix_dimension
  ALLOCATE( PUL%C%mat(1:PUL%C%dim,1:PUL%C%dim) )
  PUL%G%dim=matrix_dimension
  ALLOCATE( PUL%G%mat(1:PUL%G%dim,1:PUL%G%dim) )

! Allocate the impedance and admittance filter matrices  
  PUL%Zfilter%dim=matrix_dimension
  ALLOCATE(PUL%Zfilter%sfilter_mat(1:PUL%Zfilter%dim,1:PUL%Zfilter%dim))
  PUL%Yfilter%dim=matrix_dimension
  ALLOCATE(PUL%Yfilter%sfilter_mat(1:PUL%Yfilter%dim,1:PUL%Yfilter%dim))

! Calculate the inductance matrix first
    
  xs=PUL%overshield_x
  ys=PUL%overshield_y
  rs=PUL%overshield_r
  
  do row=1,nc-1
    
! calculate the diagonal element here    
    x_row=PUL%x(row)
    y_row=PUL%y(row)
    r_row=PUL%r(row)
    d_row=sqrt((x_row-xs)**2+(y_row-ys)**2)
    theta_row=atan2((y_row-ys),(x_row-xs))

! Check we have a valid conductor system   
    if ((d_row+r_row).GE.rs) then
      if (verbose) then
        write(*,*)'ERROR in PUL_LC_calc_overshield_wide_separation_approximation'
        write(*,*)'Conductor lies outside or intersects the reference shield conductor'
        write(*,*)'Shield radius=',rs
        write(*,*)'di =',d_row
        write(*,*)'rwi=',r_row
      end if
      run_status='ERROR in PUL_LC_calc_overshield_wide_separation_approximation,'// &
                 ' Conductor lies outside or intersects the reference shield conductor'
      CALL write_program_status()
      STOP 1
    end if
    
    PUL%L%mat(row,row)=(mu0/(2d0*pi))*log( (rs**2-d_row**2)/(rs*r_row) )  ! Paul equation 3.67a Theory_Manual_Eqn 2.27
    
    if(verbose) then
      write(*,*)'row=',row
      write(*,*)'xs=',xs
      write(*,*)'ys=',ys
      write(*,*)'rs=',rs
      write(*,*)'x_row=',x_row
      write(*,*)'y_row=',y_row     
      write(*,*)'r_row=',r_row     
      write(*,*)'d_row=',d_row  
      write(*,*)'L=', PUL%L%mat(row,row)
    end if
    
    do col=row+1,nc-1
! calculate the off diagonal elements here    

      x_col=PUL%x(col)
      y_col=PUL%y(col)
      d_col=sqrt((x_col-xs)**2+(y_col-ys)**2)
      d_row_col=sqrt((x_row-x_col)**2+(y_row-y_col)**2)
      
      theta_col=atan2((y_col-ys),(x_col-xs))

      cos_theta_row_col=cos(theta_row-theta_col)
 
! Check we have a valid conductor system   
      if ((d_col+r_col).GE.rs) then
        if (verbose) then
          write(*,*)'ERROR in PUL_LC_calc_overshield_wide_separation_approximation'
          write(*,*)'Conductor lies outside or intersects the reference shield conductor'
          write(*,*)'Shield radius=',rs
          write(*,*)'di =',d_col
          write(*,*)'rwi=',r_col
        end if
        
        run_status='ERROR in PUL_LC_calc_overshield_wide_separation_approximation,'// &
                   ' Conductor lies outside or intersects the reference shield conductor'
        CALL write_program_status()
        STOP 1
 
      end if
    
      num=(d_row*d_col)**2+rs**4-2d0*d_row*d_col*(rs**2)*cos_theta_row_col

! divide den by d_col**2
      den=(d_row)**2+d_col**2-2d0*d_row*d_col*cos_theta_row_col

! Divide by d_col
      PUL%L%mat(row,col)=(mu0/(2d0*pi))*log( (1d0/rs)*sqrt(num/den) )   ! C.R.Paul equation 367b (revised) ! Theory_Manual_Eqn 2.28
    
! The matrix is symmetric so we can set L[col,row]=L[row,col]
      PUL%L%mat(col,row)=PUL%L%mat(row,col)
      
    end do ! next col
      
  end do ! next row
  
! The capacitance matrix is calculated as mu0*eps0*[L_inverse] Theory_Manual_Eqn 2.29

  if (verbose) then
    write(*,*)'L:'
    CALL dwrite_matrix(PUL%L%mat,matrix_dimension,matrix_dimension,matrix_dimension,0)
  end if
  
  ierr=0   ! set ierr=0 on input to matrix inverse to cause the program to stop if we have a singular matrix
  CALL dinvert_Gauss_Jordan(PUL%L%mat,matrix_dimension,PUL%C%mat,matrix_dimension,ierr)
  
  PUL%C%mat(:,:)=mu0*eps0*PUL%C%mat(:,:)

! Conductance matrix is set to zero
  PUL%G%mat(:,:)=0d0
  
! calculate the impedance and admittance filters

  CALL Z_Y_from_L_C(PUL%L,PUL%C,PUL%Zfilter,PUL%Yfilter)
  
  if (verbose) then
    write(*,*)'Inductance matrix, L'
    CALL dwrite_matrix(PUL%L%mat,matrix_dimension,matrix_dimension,matrix_dimension,0)
    write(*,*)'Capacitance matrix, C'
    CALL dwrite_matrix(PUL%C%mat,matrix_dimension,matrix_dimension,matrix_dimension,0)
  end if
  
  END SUBROUTINE PUL_LC_calc_overshield_wide_separation_approximation
!
! NAME
!     SUBROUTINE calculate_height_over_ground_plane
!
!     
! COMMENTS
!     
!
! HISTORY
!    started 2/12/2015 CJS
!
  
  SUBROUTINE calculate_height_over_ground_plane(x,y,r,theta,d,h)

  USE type_specifications
  USE general_module

  IMPLICIT NONE

! variables passed to subroutine

  real(dp),intent(IN)    :: x,y     ! coordinates of point in the x-y plane
  real(dp),intent(IN)    :: r       ! radius of the wire (to check that wire doesn't intersect the ground)
  real(dp),intent(IN)    :: theta   ! angle of the ground plane normal to the x axis (radians)
  real(dp),intent(IN)    :: d       ! distance from the origin to the ground plane along the normal direction

  real(dp),intent(OUT)   :: h       ! distance from the point to the ground plane 
  
! local variables

  real(dp)    :: nx,ny  ! unit normal to ground plane in x-y plane

! START

  nx=cos(theta)
  ny=sin(theta)
  
! projecting a vector from the origin to the wire position gives 
! the distance from the origin in the normal direction. Subtracting
! the distance from the origin to the ground plane in this direction 
! gives the height of the point above the ground
  
  h=x*nx+y*ny-d
  
! Do some checks to make sure that we have a physical situation (wires above ground, no wire-ground intersections)
  
  if (h.LT.0d0) then
    
    if (verbose) then
      write(*,*)'wire x,y,r',x,y,r
      write(*,*)'ground plane theta,d',theta*57.2957795d0,d
    end if
  
    run_status='ERROR in calculate_height_over_ground_plane, the wire is below the ground plane'
    CALL write_program_status()
    STOP 1
    
  else if (h.LE.r) then
    
    if (verbose) then
      write(*,*)'wire x,y,r',x,y,r
      write(*,*)'ground plane theta,d',theta*57.2957795d0,d
    end if
    
    run_status='ERROR in calculate_height_over_ground_plane, the wire intersects the ground plane'
    CALL write_program_status()
    STOP 1
    
  end if

! passed the checks so continue...

  RETURN

  END SUBROUTINE calculate_height_over_ground_plane