!
! This file is part of SACAMOS, State of the Art CAble MOdels for Spice. 
! It was developed by the University of Nottingham and the Netherlands Aerospace 
! Centre (NLR) for ESA under contract number 4000112765/14/NL/HK.
! 
! Copyright (C) 2016-2018 University of Nottingham
! 
! SACAMOS is free software: you can redistribute it and/or modify it under the 
! terms of the GNU General Public License as published by the Free Software 
! Foundation, either version 3 of the License, or (at your option) any later 
! version.
! 
! SACAMOS is distributed in the hope that it will be useful, but 
! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 
! or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
! for more details.
! 
! A copy of the GNU General Public License version 3 can be found in the 
! file GNU_GPL_v3 in the root or at <http://www.gnu.org/licenses/>.
! 
! SACAMOS uses the EISPACK library (in /SRC/EISPACK). EISPACK is subject to 
! the GNU Lesser General Public License. A copy of the GNU Lesser General Public 
! License version can be found in the file GNU_LGPL in the root of EISPACK 
! (/SRC/EISPACK ) or at <http://www.gnu.org/licenses/>.
! 
! The University of Nottingham can be contacted at: ggiemr@nottingham.ac.uk
!
!
! File Contents:
! 
!     SUBROUTINE PUL_RL_FastHenry2
!
! NAME
!     SUBROUTINE PUL_RL_FastHenry2
!
! DESCRIPTION
!     Wrapping subroutine to control the calculation of frequency dependent impedance
!     using FastHenry2
!
!     The process is divided into the following stages:
! STAGE 1: work out the configuration for the calculation i.e. is there a ground plane, is the outer boundary free space or a conductor (overshield)
! STAGE 2: Create the input file for FastHenry2
! STAGE 4: Call FastHenry2
! STAGE 5: Read the frequency dependent impedance matrices and fit rational functions to these
!     
! COMMENTS
!   
!
! HISTORY
!    started 23/10/2023
!
SUBROUTINE PUL_RL_FastHenry2(PUL,name,fit_order,freq_spec,domain)
!
USE type_specifications
USE constants
USE general_module
USE maths
USE filter_module
USE filter_module
USE Sfilter_fit_module
USE frequency_spec
!
IMPLICIT NONE

! variables passed to subroutine

  type(PUL_type), intent(INOUT)            :: PUL        ! per-unit-length parameter calculation structure
  character(LEN=line_length),intent(IN)    :: name       ! string used as the base for filenames 
  integer, intent(IN)                      :: fit_order  ! filter fit_order for frequency dependent dielectrics
  type(frequency_specification),intent(IN) :: freq_spec  ! filter frequency range specification for frequency dependent dielectrics
  integer, intent(IN)                      :: domain     ! domain number used to label the mesh files
   
! parameters 
  integer,parameter :: inside =1 
  integer,parameter :: outside=2
 
! local variables

  integer :: ndec
  
  integer :: first_conductor,last_conductor,conductor
  
  character(LEN=filename_length) :: command                  ! string used for running external commands 
  integer :: exit_stat
  
! Process_Zc variables:


character(LEN=256) :: line
character(LEN=256) :: freq_and_dim_string

integer :: local_line_length

integer :: i

integer :: n_freq_in

real(dp) :: freq,w
integer :: nrows,ncols,r,c

complex(dp),allocatable :: Zc_in(:,:,:)
real(dp),allocatable    :: freq_in(:)

real(dp),allocatable    :: values(:),re,im

logical :: subtract_dc_resistance

integer :: loop,floop

integer :: ierr

  
! START
  if (verbose) write(*,*)'CALLED: PUL_LC_Laplace'

! STAGE 1: work out the configuration for the calculation i.e. is there a ground plane, is the outer boundary free space or a conductor (overshield)
  
  write(*,*)'Filter fit order=',fit_order
  write(*,*)'fmin=',freq_spec%fmin
  write(*,*)'fmax=',freq_spec%fmax
  write(*,*)'n_frequencies=',freq_spec%n_frequencies
  
  ndec=NINT(real(freq_spec%n_frequencies)/log10(freq_spec%fmax/freq_spec%fmin))
  
! write the file which is used to generate to FastHenry2 input file
  OPEN(unit=fh2_input_file_unit,file='fh2.txt')
  
  write(fh2_input_file_unit,'(A)')'fh2.inp'
  
  if ((.NOT.PUL%overshield_present).AND.(PUL%ground_plane_present)) then 

! SOLUTION TYPE 2:  NO OVERSHIELD, WITH GROUND PLANE

    first_conductor=1
    last_conductor=PUL%n_conductors-1
    
    write(fh2_input_file_unit,'(A)')'ground_plane'

! ********* HARD WIRED GROUND PLANE SPEC: TO BE SORTED OUT ************    
    write(fh2_input_file_unit,'(A)')'-25e-3 0.0e-3  0.0e-3	! x y z of ground plane point 1	              '
    write(fh2_input_file_unit,'(A)')' 25e-3 0.0e-3  0.0e-3	! x y z of ground plane point 2 	      '
    write(fh2_input_file_unit,'(A)')' 25e-3 0.0e-3 1000.0e-3	! x y z of ground plane point 3 	      '
    write(fh2_input_file_unit,'(A)')'0.01e-3			! thickness				      '
    write(fh2_input_file_unit,'(A)')'10 10			! discretisation in p1-p2 and p2=p3 directions'
    write(fh2_input_file_unit,'(A)')'4.461E+07  		! conductivity, sigma			      '
    write(fh2_input_file_unit,'(A)')'1  			! gp discretisation in thickness,  nhinc      '
    write(fh2_input_file_unit,'(A)')'2.0			! gp discretisation in thickness ratio, rh    '
    write(fh2_input_file_unit,'(A)')'0e-3 0.0e-3 0.0e-3 	! x y z of ground plane node at end 1	      '
    write(fh2_input_file_unit,'(A)')'0e-3 0.0e-3 1000.0e-3	! x y z of ground plane node at end 2	      '
! ********* HARD WIRED GROUND PLANE SPEC: TO BE SORTED OUT ************    
  
  else if ((.NOT.PUL%overshield_present).AND.(.NOT.PUL%ground_plane_present)) then 
    
    first_conductor=1
    last_conductor=PUL%n_conductors
  
  else
  
    write(run_status,*)'ERROR in PUL_RL_FastHenry2 cannot handle this configuration'
    CALL write_program_status()
    STOP 1
    
  end if
  
     
  write(fh2_input_file_unit,'(I2,A)')last_conductor-first_conductor+1,'       ! number of conductors'
  
  do conductor=first_conductor,last_conductor
  
    write(fh2_input_file_unit,'(I2,A)')conductor,'	 ! conductor number'
    
    if (PUL%shape(conductor).EQ.circle) then
    
      write(fh2_input_file_unit,'(A)')'cylindrical	! conductor 1 shape (cylinder, rectangle, annulus)'
      write(fh2_input_file_unit,'(2ES12.4,A)')PUL%x(conductor),PUL%y(conductor),' ! centre xc yc'
      write(fh2_input_file_unit,'(ES12.4,A)')PUL%r(conductor),' ! radius rc'
      write(fh2_input_file_unit,'(ES12.4,A)')PUL%r(conductor)/8,' ! cylindrical conductor discretisation size, dl'
      write(fh2_input_file_unit,'(A)')'4.461E+07	  ! conductivity, sigma '
      write(fh2_input_file_unit,'(A)')'3  3		  ! segment discretisation in x and y, nwinc nhinc '
      write(fh2_input_file_unit,'(A)')'2.0 2.0  	  ! segment discretisation ratio in x and y, rw rh '

    else if (PUL%shape(conductor).EQ.rectangle) then
    
      write(fh2_input_file_unit,'(A)')'rectangle	! conductor 1 shape (cylinder, rectangle, annulus)'
  
    else
      
      write(run_status,*)'ERROR in PUL_RL_FastHenry2 cannot handle shape:',PUL%shape(conductor)
      CALL write_program_status()
      STOP 1
  
    end if
    
  end do ! next conductor
  
  write(fh2_input_file_unit,'(ES12.4,A)')freq_spec%fmin,'   ! Minimum frequency, fmin'
  write(fh2_input_file_unit,'(ES12.4,A)')freq_spec%fmax,'   ! Maximum frequency, fmax'
  write(fh2_input_file_unit,'(I8,A)')ndec,'   ! Number of points per decade, ndec'

  CLOSE(unit=fh2_input_file_unit)
  
  write(*,*)'CALLING FastHenry2'
  command='/home/chris/SACAMOS_FASTHENRY2/FASTHENRY2_WRITE_INPUT_FILE/bin/write_FH_input_file < fh2.txt'
  CALL execute_command_line(command,EXITSTAT=exit_stat)
    
! Check that the FastHenry2 input file generated correctly
  if (exit_stat.NE.0) then  
!  write_FH_input_file returned with a non zdero exit code indicating an error
    write(run_status,*)'ERROR in PUL_RL_FastHenry2 running write_FH_input_file: exit_stat=',exit_stat
    CALL write_program_status()
    STOP 1
  end if
  
  write(*,*)'CALLING FastHenry2'
  command='fasthenry fh2.inp'
  CALL execute_command_line(command,EXITSTAT=exit_stat)
    
! Check that the fasthenry ran correctly
  if (exit_stat.NE.0) then  
!  fasthenry returned with a non zdero exit code indicating an error
    write(run_status,*)'ERROR in PUL_RL_FastHenry2 running fasthenry: exit_stat=',exit_stat
    CALL write_program_status()
    STOP 1
  end if
  
  if (verbose) then
    write(*,*)
    write(*,*)' Processing frequency dependent impedance file from FastHenry2'
    write(*,*)
  end if
  
  OPEN(unit=fh2_output_file_unit,file='Zc.mat')  
  
  do loop=1,2

    n_freq_in=0

10  CONTINUE
    read(fh2_output_file_unit,'(A)',END=1000,ERR=1000)line
    if (line(1:32).NE.'Impedance matrix for frequency =') GOTO 10

!     if we get here then we have some impedance matrix data to read
!      write(*,*)'Found impedance matrix:',trim(line)

! increase the number of frequencies for which we have an impedance matrix    
      n_freq_in=n_freq_in+1

      local_line_length=len(trim(line))
      freq_and_dim_string=line(33:local_line_length)
  
!      write(*,*)'reading matrix for:',trim(freq_and_dim_string)

! replace the 'x' by a space then read the frequency, nrows and ncols
      local_line_length=len(trim(freq_and_dim_string)) 
      do i=1,local_line_length
        if (freq_and_dim_string(i:i).EQ.'x') freq_and_dim_string(i:i)=' '
      end do
  
!      write(*,*)'data string:',trim(freq_and_dim_string)
  
      read(freq_and_dim_string,*)freq,nrows,ncols
      if (nrows.NE.ncols) then
        write(*,*)'****** ERROR: nrows.NE.ncols ******'
      end if
!      write(*,*)'frequency=',freq,' n=',nrows
    
      if (loop.EQ.2) then
        freq_in(n_freq_in)=freq
      end if
        
      do r=1,nrows
  
        read(fh2_output_file_unit,'(A)',END=1000,ERR=1000)line

        if (loop.EQ.2) then
      
!   replace the 'j's by a space then read the complex data
        local_line_length=len(trim(line)) 
          
        do i=1,local_line_length
          if (line(i:i).EQ.'j') line(i:i)=' '
        end do
  
        read(line,*)(values(i),i=1,2*ncols)
    
        do c=1,ncols
    
          re=values(c*2-1)
          im=values(c*2)
                
          Zc_in(n_freq_in,r,c)=cmplx(re,im)
	
        end do ! next column
	 
      end if
  
    end do ! next row
     
    GOTO 10 ! continue to read the file

! Jump here when the file has been read

1000 CONTINUE
 
    if (loop.Eq.1) then
      if (verbose) write(*,*)'Number of frequencies=',n_freq_in
      ALLOCATE( freq_in(n_freq_in) )
      ALLOCATE( Zc_in(n_freq_in,nrows,ncols) )
      ALLOCATE( values(2*ncols) )
      rewind(unit=fh2_output_file_unit)
    end if
  
  end do ! next file read loop
  
  CLOSE(unit=fh2_output_file_unit)

! For now use the inductance from the highest frequency calculated 
  do r=1,nrows
    do c=1,ncols
    
      freq=freq_in(n_freq_in)
      w=2.0*pi*freq
      
      PUL%L%mat(r,c)=real(Zc_in(n_freq_in,r,c)/(j*w))
      PUL%Zfilter%sfilter_mat(r,c)=jwA_filter( PUL%L%mat(r,c) )
      
      write(*,*)r,c,PUL%L%mat(r,c)
      
    end do
  end do
 
  DEALLOCATE( freq_in )
  DEALLOCATE( Zc_in )
  DEALLOCATE( values )
    
  if (verbose) then
    write(*,*)'FINISHED: PUL_RL_FastHenry2'
    write(*,*)'Inductance matrix, L'
    CALL dwrite_matrix(PUL%L%mat,nrows,ncols,nrows,0)
  end if
  
  if (verbose) write(*,*)'FINISHED PUL_RL_FastHenry2'

!  STOP 1
  
END SUBROUTINE PUL_RL_FastHenry2