!
! 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 calculate_propagation_correction_filters
!
! NAME
!     calculate_propagation_correction_filters
!
! DESCRIPTION
!     Work out the complex correction filter to be applied to the dispersionless transmission line in order
!     to correct for frequency dependent losses
!     The code works out the filter frequency response over the frequency range of interest then
!     a vector fitting process requesting the appropriate filter order
!     The filter produced by the vector fitting process is then read back
!     
! COMMENTS
!     
!
! HISTORY
!
!     started 22/04/2016 CJS based on phase 1 work in 2015
!     8/5/2017         CJS: Include references to Theory_Manual
!     23/10/2017       CJS: New approach to mode identification based on propagation constants rather than mode shapes. This
!                           may improve matters for degenerate modes
!     16/11/2017 CJS Include network synthesis process to replace s-domain transfer functions
!
!
SUBROUTINE calculate_propagation_correction_filters(domain,dim,model_order,Z_filter,Y_filter,L,C,TI_REF, &
               conductor_impedance,Hfilter,length,delay,freq_spec)

USE type_specifications
USE general_module
USE constants
USE cable_module
USE maths
USE filter_module
USE sfilter_fit_module
USE frequency_spec

IMPLICIT NONE

! variables passed to subroutine
integer,intent(IN)  :: domain                  ! domain number
integer,intent(IN)  :: dim                     ! matrix dimension
integer,intent(IN)  :: model_order             ! order of propagation correction model to be calculated

real(dp),intent(IN) :: length                  ! transmission line length
real(dp),intent(IN) :: delay(dim)              ! mode delays for the lossless idealised transmission line

type(Sfilter),intent(IN) :: Z_filter(dim,dim)  ! frequency dependent impedance matrix
type(Sfilter),intent(IN) :: Y_filter(dim,dim)  ! frequency dependent admittance matrix

real(dp),intent(IN) :: L(dim,dim)  ! 'high frequency' inductance matrix (that used to calculate modal decomposition for Spice model)
real(dp),intent(IN) :: C(dim,dim  )! 'high frequency' capacitance matrix (that used to calculate modal decomposition for Spice model)

real(dp),intent(IN) :: TI_REF(dim,dim)     ! modal decomposition matrix used in the spice model

! data for the loss model
type(conductor_impedance_model),intent(IN) :: conductor_impedance(dim+1)  ! conductor based impedance model (impedance due to finite conductivity)

type(Sfilter),intent(OUT) :: Hfilter(1:dim)  ! propagation correction filters for each mode to be returned from the subroutine

type(frequency_specification),intent(IN) :: freq_spec  ! frequency range specificatiobn for propagation filter fitting process

! local variables

real(dp)    :: Rdc(dim+1)   ! dc resistance on conductors
complex(dp) :: Zint(dim+1)  ! complex internal impedance on conductors

real(dp)    :: Rdc_temp       ! temporary variable for dc resistance on conductors
complex(dp) :: Zint_0(dim+1)  ! complex internal impedance on conductors at zero frequency

real(dp)    :: Rdc_t   ! dc transfer resistance 
complex(dp) :: Zint_t  ! complex transfer impedance

complex(dp) :: Z(dim,dim)  ! impedance matrix at a particular frequency
complex(dp) :: Y(dim,dim)  ! admittance matrix at a particular frequency

complex(dp) :: YZ(dim,dim) ! [Y][Z] matrix product

complex(dp) :: TI(dim,dim) ! Modal decomposition matrices of [YZ] product
complex(dp) :: TII(dim,dim)
complex(dp) :: GAMMA_SQR(dim) ! eigenvalues of modal decomposition 

complex(dp) :: TEST_MAT(dim,dim)

complex(dp) :: GAMMA_C(dim)      ! propagation constants for the frequency dependent model
complex(dp) :: last_GAMMA_C(dim) ! propagation constants for the frequency dependent model at the previous frequency
complex(dp) :: GAMMA_C_hf(dim)   ! propagation constants for the lossless high frequency model

! propagation factors for the length of transmission line for frequency dependent model
complex(dp) :: Exp_m_gamma_l(dim)  

! propagation factors for the length of transmission line for the high frequency lossless model
complex(dp) :: Exp_m_gamma_l_hf(dim)

complex(dp),allocatable :: Correction(:,:)  ! propagation correction values for each mode and frequency of interest
complex(dp),allocatable :: Hp(:)            ! propagation correction value for a particular mode for all frequencies of interest 
                                            ! this is an input to the filter_fit process

real(dp) :: f             ! frequency (Hz)
real(dp) :: w             ! angular frequency (rad/s)
integer :: frequency_loop ! frequency loop variable

! strings to construct filenames
character(LEN=line_length)::name1,name2,filename(dim)

! temporary loop variables
integer :: row,col,i,ii

real(dp) :: oparray(dim*2)

! variables used to reorder modes from the frequency dependent modal decomposition process
integer :: ref_mode_to_local_mode_list(dim)
integer :: mode_chosen_list(dim)
real(dp):: largest_value
integer :: largest_value_row
integer :: largest_value_col
real(dp):: smallest_value
integer :: smallest_value_row
logical :: value_picked(dim) ! flag to indicate that this array element has already been picked
integer :: local_mode

! error status for matrix inverse calculation.
integer  :: ierr

! switch for different mode identification algorithms
 integer :: mode_test_type
 integer,parameter :: mode_test_type_eigenvector=1
 integer,parameter :: mode_test_type_eigenvalue=2
 

! START 

  if (verbose) write(*,*)'CALLED: calculate_propagation_correction_filters'
  if (verbose) then
    write(*,*)'TYPE:',freq_spec%freq_range_type
    write(*,*)'fmin=',freq_spec%fmin,' fmax=',freq_spec%fmax
    write(*,*)'Number of frequencies=',freq_spec%n_frequencies
  end if
  
  
  mode_test_type=mode_test_type_eigenvector

 
! alllocate memory for the propagation correction values  
  ALLOCATE( Correction(1:dim,1:freq_spec%n_frequencies) )
  ALLOCATE( Hp(1:freq_spec%n_frequencies) )

!  if (plot_propagation_correction_filter_fit_data) then

! open files for plotting the propagation correction functions

! write a gnuplot file to plot this correction function    
    CALL add_integer_to_string('plot_correction_filters_domain_',domain,name1)
    open( unit=80,file=trim(name1)//'.plt' )
    if (freq_spec%freq_range_type.EQ.'log') write(80,*)'set logscale x'

    do i=1,dim
  
      CALL add_integer_to_string('domain_',domain,name1)
      CALL add_integer_to_string('_correction_for_mode_',i,name2)
      filename(i)=trim(name1)//trim(name2)//'.Z.fdmat'
      open(unit=80+i,file=filename(i))
    
      if (verbose) write(*,*)'Opening file:',trim(filename(i))

! write gnuplot information    
      write(80,'(A)')'set autoscale x'
      write(80,'(A)')'set autoscale y'
      write(80,'(A)')'set xlabel "frequency(Hz)"'
      write(80,'(A)')'set ylabel "Propagation correction"'
      write(80,'(A,A,A)')'plot "',trim(filename(i))//'_fit','" u 1:2 title "Re{Corr_fit}" w p,\'
      write(80,'(A,A,A)')'     "',trim(filename(i))//'_fit','" u 1:3 title "Im{Corr_fit}" w p,\'
      write(80,'(A,A,A)')'     "',trim(filename(i)),'" u 1:2 title "Re{Corr}" w l,\'
      write(80,'(A,A,A)')'     "',trim(filename(i)),'" u 1:3 title "Im{Corr}" w l'
      write(80,'(A)')'pause -1'
    
    end do
    
    close(unit=80)

!  end if ! plot_propagation_correction_filter_fit_data

! loop over frequency working out the correction required for each mode at each frequency
! Note we work downwards in frequency as we use high frequency values of L and C to define the
! reference modal decomposition and propagation constants
  
  do frequency_loop=freq_spec%n_frequencies,1,-1
  
    f=freq_spec%freq_list(frequency_loop)
  
! shift the frequency from d.c. as jwL, jwC are zero at d.c. and this messes up the process
    if (f.EQ.0d0) f=1d0  
  
    if (verbose) write(*,*)'frequency',frequency_loop,' of',freq_spec%n_frequencies,' f=',f

    w=2d0*pi*f
    
! calculate the conductor impedance contribution for all the conductors in the domain
    do row=1,dim+1
    
! calculate Zint_0 i.e. the impedance at zero frequency. This is the quantity to subtract from
! the conductor impedance rather than Rdc as it is different for common mode/ differential mode domains.
                                              
      CALL evaluate_conductor_impedance_model(conductor_impedance(row),0d0,Zint_0(row),Rdc_temp,Zint_t,Rdc_t)
                                                  
! calculate the conductor impedance contribution at this frequency
      CALL evaluate_conductor_impedance_model(conductor_impedance(row),f,Zint(row),Rdc(row),Zint_t,Rdc_t)
            
! Apply multiplication factor to the conductor impedance to correct for common mode/ differential modes in twisted pairs
! See note at the end of Theory_Manual_Section 3.5.4
      Zint_0(row)=Zint_0(row)*conductor_impedance(row)%Resistance_multiplication_factor
                                                                                           
! Apply multiplication factor to the conductor impedance to correct for common mode/ differential modes in twisted pairs
      Zint(row)=Zint(row)*conductor_impedance(row)%Resistance_multiplication_factor
            
    end do ! next row

! Use the pre_calulcated delay (input to the subroutine) to obtain GAMMA_C_hf 
! Theory_Manual_Eqns 3.81, 3.82, 3.83

    GAMMA_C_hf(:)=j*w*delay(:)/length
    
    if (frequency_loop.EQ.freq_spec%n_frequencies) then
      last_GAMMA_C(1:dim)=GAMMA_C_hf(dim)
    else
      last_GAMMA_C(1:dim)=GAMMA_C(dim)
    end if
    
! Do the same with the frequency dependent model
! This is the full frequency dependent model including frequency dependent dielectrics
! and the 'skin effect' loss model
! Note that we subtract the d.c. resistance from Zint 
  
! Fill the MTL impedance and admittance matrices  with the frequency dependent 
! values calculated from the rational function form of the Z and Y matrix elements
! evaluated at the current frequency

    do row=1,dim
      do col=row,dim

        Z(row,col)=evaluate_Sfilter_frequency_response(Z_filter(row,col),f)      
	
!  Add the internal impedance contribution from the reference conductor
        Z(row,col)=Z(row,col)+Zint(dim+1)-Zint_0(dim+1)
        	
!  Add the internal impedance contribution from this conductor
        if (row.eq.col) then
          Z(row,col)=Z(row,col)+Zint(row)-Zint_0(row)
        end if

!  The impedance matrix is symmetrical
        Z(col,row)=Z(row,col)

! There is a contribution to the admittance from the Y_filter due to frequency dependent dielectrics
        Y(row,col)=evaluate_Sfilter_frequency_response(Y_filter(row,col),f)
        
!  The admittance matrix is symmetrical
        Y(col,row)=Y(row,col)
     
      end do ! next col
    end do  ! next row

! calculate the modal decomposition of the YZ product (Theory_Manual_Eqn 2.33)
    YZ=matmul(Y,Z)
    CALL calc_eigenvectors(YZ,GAMMA_SQR,TI,dim) 

! Get the vector, gamma, with the frequency dependent propagation constant for each mode    
! we require gamma to have positive real and imaginary parts

    GAMMA_C(:)=sqrt(GAMMA_SQR(:))
    do row=1,dim
    
      if      ( (real(GAMMA_C(row)).LT.0d0).AND.(aimag(GAMMA_C(row)).LT.0d0) ) then
      
        GAMMA_C(row)=-GAMMA_C(row)
        
      else if ( (real(GAMMA_C(row)).GT.0d0).AND.(aimag(GAMMA_C(row)).LT.0d0) ) then
      
        if (real(GAMMA_C(row)).GT.1d-10) then
          write(*,*)'Problem with complex propagation constant, gamma=',GAMMA_C(row)        
          STOP
        else
          GAMMA_C(row)=-j*aimag(GAMMA_C(row))
        end if
        
      else if ( (real(GAMMA_C(row)).LT.0d0).AND.(aimag(GAMMA_C(row)).GT.0d0) ) then
      
        if (real(GAMMA_C(row)).LT.1d-10) then
          write(*,*)'Problem with complex propagation constant, gamma=',GAMMA_C(row)        
          STOP
        else
          GAMMA_C(row)=j*aimag(GAMMA_C(row))
        end if
        
      end if
      
    end do
    
! Work out how the modes calculated here relate to those calculated in the lossless case. 
! as there may be some re-ordering required see Theory_Manual_Section 3.5.4 below equation 3.97

    ierr=1   ! set ierr=1 on input to matrix inverse to cause the program to return here if we have a singular matrix

    TEST_MAT=MATMUL(TRANSPOSE(TI),TI_REF)

    if (mode_test_type.EQ.mode_test_type_eigenvector) then

      do i=1,dim      ! loop over modes
    
        largest_value=0d0
        largest_value_col=0
        largest_value_row=0

! find the largest value still in the matrix     
        do row=1,dim    
          do col=1,dim
     
            if ( abs(TEST_MAT(row,col)).GT.largest_value)then
! this is the largest un-chosen column value so far
              largest_value=abs(TEST_MAT(row,col))
              largest_value_row=row
              largest_value_col=col
            end if
       
          end do
        end do
      
        ref_mode_to_local_mode_list(largest_value_col)=largest_value_row
        
 ! set the values in this row and col to zero so we cannot pair these modes again
        do ii=1,dim  
          TEST_MAT(largest_value_row,ii)=(0d0,0d0)
          TEST_MAT(ii,largest_value_col)=(0d0,0d0)
        end do
     
      end do ! next mode
  
! check that we have identified a unique set of modes
      do row=1,dim
    
        if (ref_mode_to_local_mode_list(row).EQ.0) then
          run_status='ERROR in propagation_correction_filters: unable to identify modes uniquely'
          CALL write_program_status()
          STOP 1
        end if
        
      end do  ! next row
      
    else if (mode_test_type.EQ.mode_test_type_eigenvalue) then
    
      value_picked(:)=.FALSE.
      
      do row=1,dim
      
        smallest_value_row=0
        smallest_value=1d30
        
        do ii=1,dim
        
          if (.NOT.value_picked(ii)) then
            if ( abs(last_GAMMA_C(ii)-GAMMA_C(ii)).LT.smallest_value) then
! this is the smallest un-chosen column value so far
              smallest_value=abs(last_GAMMA_C(ii)-GAMMA_C(ii))
              smallest_value_row=ii
            end if          
          end if
          
        end do

        ref_mode_to_local_mode_list(row)=smallest_value_row
        value_picked(smallest_value_row)=.TRUE.
        
      end do
      
    end if
    
! work out the propagation correction for each of the modes  

! loop over reference modes i.e. those used in the spice model modal decomposition
    do row=1,dim
    
! the correction is applied at both ends of the transmission line so 
! each correction corrects for half the propagation distance hence use length/2 below
! Theory_manual_Eqn 3.97

       Exp_m_gamma_l_hf(row)=exp(-GAMMA_C_hf(row)*length/2d0)

! get the coresponding local mode number       
       local_mode=ref_mode_to_local_mode_list(row)
       
       Exp_m_gamma_l(row)   =exp(-GAMMA_C(local_mode)*length/2d0) 
       
       Correction(row,frequency_loop)= (Exp_m_gamma_l(row)/Exp_m_gamma_l_hf(row))
       
       if (plot_propagation_correction_filter_fit_data) then
! write propagation correction data for plotting
!         write(80+row,8000)f,real(Correction(row,frequency_loop)), &
!                           aimag(Correction(row,frequency_loop)),abs(Correction(row,frequency_loop))
         write(80+row,8000)f,real(GAMMA_C(local_mode)),aimag(GAMMA_C(local_mode)),real(GAMMA_C_hf(row)),aimag(GAMMA_C_hf(row))
       end if
8000 format(5E16.6)

    end do ! next row
  
  end do ! next frequency
  
  if (plot_propagation_correction_filter_fit_data) then
    do i=1,dim
      close(unit=80+i)
    end do
  end if
  
! for each of the modes calculate the best fit rationalfilter function for the propagation correction using 
! the Calculate_Sfilter process described in Theory_Manual_Section 5.2

  do i=1,dim

! Copy the complex frequency domain data to the local array  
    do frequency_loop=1,freq_spec%n_frequencies
    
      Hp(frequency_loop)=Correction(i,frequency_loop)
      
    end do

! Call the filter fitting process with this mode propagation correction data 
! with border=aorder and fit_type=1 i.e. Hfilter(f=0)=1 is imposed
    CALL Calculate_Sfilter(Hp,freq_spec%freq_list,freq_spec%n_frequencies,Hfilter(i),model_order,0,1) 
    
    if (plot_propagation_correction_filter_fit_data) then
    
      open(unit=80,file=trim(filename(i))//'_fit')
      
      do frequency_loop=1,freq_spec%n_frequencies
        
        write(80,8010)f,evaluate_Sfilter_frequency_response(Hfilter(i),f)
8010 format(3E16.6)
      
      end do ! next frequency
      
      close(unit=80)
      
    end if ! plot_propagation_correction_filter_fit_data
    
  end do ! next mode

! deallocate the local storage  
  DEALLOCATE( Correction )
  DEALLOCATE( Hp )
  
END SUBROUTINE calculate_propagation_correction_filters

