!
! 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 :: first_conductor,last_conductor,conductor
  
  character(LEN=filename_length) :: command                  ! string used for running external commands 
  integer :: exit_stat
  
! Process_Zc variables:

  integer  :: matrix_dimension
  complex(dp) :: Zc_temp
  complex(dp),allocatable :: function_to_fit(:)    ! complex function to be fitted using Sfilter_fit
  
  real(dp) :: ferr
  
  real(dp) :: gp_half_width

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,row,col

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

real(dp),allocatable    :: Rdc_domain(:,:)

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

logical :: subtract_dc_resistance

integer :: loop,floop

integer :: ierr

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

! work out the local frequency sampling: FastHenry uses a specific form based on log frequency
! and number of samples per decade
  
  write(*,*)'Filter fit order=',fit_order
  write(*,*)'fmin=',freq_spec%fmin
  write(*,*)'fmax=',freq_spec%fmax
  write(*,*)'n_frequencies=',freq_spec%n_frequencies  
  write(*,*)'ndec=',freq_spec%ndec 
  
  if (freq_spec%freq_range_type.EQ.'lin') then
    write(run_status,*)'ERROR in PUL_RL_FastHenry2: We must have a logarithmic frequency range specified'
    CALL write_program_status()
    STOP 1

  end if

! 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 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 1:  NO OVERSHIELD, WITH GROUND PLANE

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

! ********* STILL SOME HARD WIRED GROUND PLANE SPEC: TO BE SORTED OUT ************  

! ALSO NEED TO CHECK THAT EVERYTHING WE NEED IS SET...

    write(*,*)'Ground plane parameters:'
    write(*,*)'PUL%ground_plane_w=',PUL%ground_plane_w
    write(*,*)'PUL%ground_plane_h=',PUL%ground_plane_h
    write(*,*)'PUL%ground_plane_sigma=',PUL%ground_plane_sigma
    write(*,*)'PUL%ground_plane_nsegx=',PUL%ground_plane_nsegx
    write(*,*)'PUL%ground_plane_nsegz=',PUL%ground_plane_nsegz
    write(*,*)'PUL%ground_plane_nh=',PUL%ground_plane_nh

    if ( (PUL%ground_plane_w.EQ.0d0).OR.(PUL%ground_plane_h.EQ.0d0).OR.(PUL%ground_plane_sigma.EQ.0d0) &
     .OR.(PUL%ground_plane_nsegx.EQ.0).OR.(PUL%ground_plane_nsegz.EQ.0).OR.(PUL%ground_plane_nh.EQ.0) )then
      write(run_status,*)'ERROR in PUL_RL_FastHenry2: ground plane parameters not defined'
      CALL write_program_status()
      STOP 1
    end if

    gp_half_width=PUL%ground_plane_w/2.0
    
    write(fh2_input_file_unit,'(2ES12.4,A)')-gp_half_width,-PUL%ground_plane_h/2.0,' 0.0e-3	! x y z of ground plane point 1'
    write(fh2_input_file_unit,'(2ES12.4,A)') gp_half_width,-PUL%ground_plane_h/2.0,' 0.0e-3	! x y z of ground plane point 2'
    write(fh2_input_file_unit,'(2ES12.4,A)') gp_half_width,-PUL%ground_plane_h/2.0,' 1000.0e-3	! x y z of ground plane point 3'
    write(fh2_input_file_unit,'(ES12.4,A)')PUL%ground_plane_h,'	 ! thickness				      '
    write(fh2_input_file_unit,'(2I6,A)')PUL%ground_plane_nsegx,PUL%ground_plane_nsegz,  &
                                        '  ! discretisation in p1-p2 and p2=p3 directions'
    write(fh2_input_file_unit,'(ES12.4,A)')PUL%ground_plane_sigma,' ! conductivity, sigma			      '
    write(fh2_input_file_unit,'(I4,A)')PUL%ground_plane_nh,'	! 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	      '
! ********* STILL SOME HARD WIRED GROUND PLANE SPEC: TO BE SORTED OUT ************    
  
  else if ((PUL%overshield_present).AND.(.NOT.PUL%ground_plane_present)) then 

! SOLUTION TYPE 2:  OVERSHIELD, NO GROUND PLANE

    first_conductor=1
    last_conductor=PUL%n_conductors
  
  else if ((.NOT.PUL%overshield_present).AND.(.NOT.PUL%ground_plane_present)) then 

! SOLUTION TYPE 3:  NO OVERSHIELD, NO GROUND PLANE

    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 ( (conductor.EQ.last_conductor).AND.PUL%overshield_present) then

      write(fh2_input_file_unit,'(A)')'annulus    ! conductor 1 shape (cylinder, rectangle, annulus)'
      write(fh2_input_file_unit,'(2ES12.4,A)')PUL%overshield_x,PUL%overshield_y,' ! centre xc yc'
      write(fh2_input_file_unit,'(ES12.4,A)')PUL%overshield_r,' ! inner radius rc'
      write(fh2_input_file_unit,'(ES12.4,A)')PUL%overshield_r*1.01,' ! ***** outer radius rc: 1.01*rinner *****'
      write(fh2_input_file_unit,'(ES12.4,A)')PUL%overshield_r/2,' ! ***** conductor discretisation size, dl=r/2 *****'
      write(fh2_input_file_unit,'(ES12.4,A)')1e20,'  ! **** conductivity, sigma: high conductivity ****'
      
!      write(fh2_input_file_unit,'(2I4,A)')FH2_nw,FH2_nh,' ! segment discretisation in x and y, nwinc nhinc '
!      write(fh2_input_file_unit,'(2ES12.4,A)')FH2_rw,FH2_rh,'  ! segment discretisation ratio in x and y, rw rh '

      write(fh2_input_file_unit,'(A)')' 3 1 ! segment discretisation in x and y, nwinc nhinc '
      write(fh2_input_file_unit,'(A)')'1.0 1.0  ! segment discretisation ratio in x and y, rw rh '
     
    else  ! not an overshield
    
      if (PUL%shape(conductor).EQ.circle) then
      
        if (PUL%conductor_is_shield(conductor)) then
      
! ********* USED FOR SHIELDS IN THE EXTERNAL DOMAIN: STILL SOME HARD WIRED PARAMETERS ************    
          write(fh2_input_file_unit,'(A)')'annulus    ! conductor 1 shape (cylinder, rectangle, annulus)'    
          write(fh2_input_file_unit,'(2ES12.4,A)')PUL%x(conductor)+PUL%ox(conductor),&
                                                PUL%y(conductor)+PUL%oy(conductor),' ! centre xc yc'              
          write(fh2_input_file_unit,'(ES12.4,A)')PUL%r(conductor)/1.01d0,' ! inner radius rs/1.01'
          write(fh2_input_file_unit,'(ES12.4,A)')PUL%r(conductor),' ! outer radius rs'
          write(fh2_input_file_unit,'(ES12.4,A)')PUL%overshield_r/2,' ! ***** conductor discretisation size, dl=r/2 *****'
          write(fh2_input_file_unit,'(ES12.4,A)')1e20,'  ! **** conductivity, sigma: high conductivity ****'     
          write(fh2_input_file_unit,'(A)')' 3 1 ! segment discretisation in x and y, nwinc nhinc '
          write(fh2_input_file_unit,'(A)')'1.0 1.0  ! segment discretisation ratio in x and y, rw rh '
        
        else
        
          if(auto_cyl_grid) then
            write(fh2_input_file_unit,'(A)')'cylindrical  grid_block auto ! conductor 1 shape (cylinder, rectangle, annulus)'
          else
            write(fh2_input_file_unit,'(A)')'cylindrical  ! conductor 1 shape (cylinder, rectangle, annulus)'
          end if
          
          write(fh2_input_file_unit,'(2ES12.4,A)')PUL%x(conductor)+PUL%ox(conductor),&
                                                PUL%y(conductor)+PUL%oy(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)/FH2_nlayers_radius, &
                                               ' ! cylindrical conductor discretisation size, dl'
          write(fh2_input_file_unit,'(ES12.4,A)')PUL%sigma(conductor),'  ! conductivity, sigma '
          write(fh2_input_file_unit,'(2I4,A)')FH2_nw,FH2_nh,' ! segment discretisation in x and y, nwinc nhinc '
          write(fh2_input_file_unit,'(2ES12.4,A)')FH2_rw,FH2_rh,'  ! segment discretisation ratio in x and y, rw rh '
        
        end if    ! conductor_is_shield
        
      else if (PUL%shape(conductor).EQ.rectangle) then
    
        write(fh2_input_file_unit,'(A)')'rectangle	! conductor 1 shape (cylinder, rectangle, annulus)'
        write(fh2_input_file_unit,'(2ES12.4,A)')PUL%x(conductor)+PUL%ox(conductor),&
                                                PUL%y(conductor)+PUL%oy(conductor),' ! centre xc yc'
        write(fh2_input_file_unit,'(2ES12.4,A)')PUL%rw(conductor),PUL%rh(conductor),' ! width and height'
        write(fh2_input_file_unit,'(ES12.4,A)')PUL%sigma(conductor),'  ! conductivity, sigma '
        write(fh2_input_file_unit,'(2I4,A)')FH2_nw,FH2_nh,' ! segment discretisation in x and y, nwinc nhinc '
        write(fh2_input_file_unit,'(2ES12.4,A)')FH2_rw,FH2_rh,'  ! segment discretisation ratio in x and y, rw rh '
 
      else
      
        write(run_status,*)'ERROR in PUL_RL_FastHenry2 cannot handle shape:',PUL%shape(conductor)
        CALL write_program_status()
        STOP 1
  
      end if   ! shape
      
    end if ! overshield
    
  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,'(ES12.4,A)')freq_spec%ndec,'   ! Number of points per decade, ndec'

  CLOSE(unit=fh2_input_file_unit)
  
  write(*,*)'CALLING write_FH_input_file'
  command='write_FH_input_file  < fh2.txt'
  write(*,*)'COMMAND:'
  write(*,'(A)')command
  
  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'
  
  if(fh2_no_refinement) then
    command='fasthenry -a off fh2.inp'
  else
    command='fasthenry fh2.inp'
  end if
  
  write(*,*)'COMMAND:'
  write(*,'(A)')command
  
  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
  
  command='cp fh2.txt fh2.txt_'//trim(name)
  CALL execute_command_line(command,EXITSTAT=exit_stat)
  
  command='cp fh2.inp fh2.inp_'//trim(name)
  CALL execute_command_line(command,EXITSTAT=exit_stat)
  
  command='cp Zc.mat Zc.mat_'//trim(name)
  CALL execute_command_line(command,EXITSTAT=exit_stat)
  
  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( Rdc_domain(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)
  
  matrix_dimension=PUL%n_conductors-1
  
  if (matrix_dimension.NE.nrows) then
    write(run_status,*)'ERROR in PUL_RL_FastHenry2 matrix dimension',matrix_dimension,nrows
    CALL write_program_status()
    STOP 1
  end if

! Make the impedance matrix explicitly symmetric
  do floop=1,n_freq_in
    do r=1,nrows
      do c=r+1,ncols
    
        Zc_temp=(Zc_in(floop,r,c)+Zc_in(floop,c,r))/2d0
      
        Zc_in(floop,r,c)=Zc_temp
        Zc_in(floop,c,r)=Zc_temp
      
      end do
    end do
  end do

  if (L_from_fasthenry) then
! For the modal decomposition, use the inductance from the highest frequency calculated 

!  ******************* TO BE USED WITH CARE DUE TO OBSERVED FASTER THAN LIGHT MODES AND *************
!  ********* MODE BEATING IN SOME PLANE PROBLEMS E.G. TWO WIRES OVER GND IN FREE SPACE      *********
  
    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))      
      
      end do
    end do
  
  end if

! Save the resistance matrix from the lowest frequency calculated 
  do r=1,nrows
    do c=1,ncols
    
      Rdc_domain(r,c)=real(Zc_in(1,r,c))      
      
    end do
  end do
  
! We calculate the inductance and resistance matrix at a number of frequencies before fitting filter functions
! to each of the frequency dependent impedance matrix entries.
  
  write(*,*)'**************************'
  write(*,*)'Process FastHenry2 results'
  write(*,*)'**************************'
  write(*,*)'Number of frequencies (freq_spec) =',freq_spec%n_frequencies
  write(*,*)'Number of frequencies (FastHenry2)=',n_freq_in

! Frequency check  
  if (n_freq_in.NE.freq_spec%n_frequencies) then
    write(run_status,*)'ERROR in PUL_RL_FastHenry2 n_frequencies specification',n_freq_in,freq_spec%n_frequencies
    CALL write_program_status()
    STOP 1
  end if
  
  do floop=1,n_freq_in
    ferr=abs(freq_spec%freq_list(floop)-freq_in(floop))/(freq_spec%freq_list(floop)+freq_in(floop))
    if (ferr.GT.0.001) then
      write(run_status,*)'ERROR in PUL_RL_FastHenry2 frequency samples, floop=',floop, &
                         ' freq_spec=',freq_spec%freq_list(floop),' freq_FH2=',freq_in(floop)
      CALL write_program_status()
      STOP 1   
    end if
  end do
 
! loop over frequency    
  if (verbose) then
    do floop=1,n_freq_in
      write(*,*)freq_spec%freq_list(floop),freq_in(floop)
    end do
  end if
  
! we have frequency domain impedance matrix, Zc

! loop over the elements of Zc and the Zfilter matrix by filter fitting 
  ALLOCATE( function_to_fit(1:freq_spec%n_frequencies) )
    
  do row=1,matrix_dimension
    do col=row,matrix_dimension

! get the function values for this matrix element function_to_fit=R-Rdc+jwL
      do floop=1,freq_spec%n_frequencies
        function_to_fit(floop)=Zc_in(floop,row,col)-Rdc_domain(row,col)
      end do
        
! calculate the Zfilter matrix using the filter fitting process
! note aorder=border and no restrictions are put on the function

      CALL Calculate_Sfilter(function_to_fit,freq_spec%freq_list,freq_spec%n_frequencies,  &
                             PUL%Zfilter%sfilter_mat(row,col),fit_order+1,-1,2)    ! note type 2 filter fit=> a0=0.0

      if (col.NE.row) then       
! make the matrix symmetrical
        PUL%Zfilter%sfilter_mat(col,row)=PUL%Zfilter%sfilter_mat(row,col)
      end if
        
    end do ! next col
      
  end do ! next row
    
  DEALLOCATE( function_to_fit )  
 
  DEALLOCATE( freq_in )
  DEALLOCATE( Zc_in )
  DEALLOCATE( Rdc_domain )
  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