!
! 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 network_synthesis
!
! NAME
!     network_synthesis
!
! AUTHORS
!     Chris Smartt
!
! DESCRIPTION
!     testing some ideas for the implementation of rational transfer functions as
!     passive equivalent circuits
!
!     PROCESS:
!
!  1. Test that the rational function does constitute a physical impedance function
!  2. Generate a ladder network which synthesises this impedance function where the
!     elements of the ladder network consist of RLC elements
!  3. Evaluate the frequency response of the original function
!  4. Evaluate the frequency response of the continued fraction form
!  5. Write a Spice file for the network
!
!
! COMMENTS
!     The network synthesis algorithm is described in the Theory Manual, section 7.2.1
!
! HISTORY
!
!     started 22/08/2017 CJS 
!     14/9/2017 CJS Generalise the circuit elements to RLC combinations 
!                   using partial fraction expansions
!     10/2017 CJS include Brune synthesis method to complete the process
!
!
SUBROUTINE network_synthesis(Hin,gain_in,Hname,in1,in2,on1,on2,vref_node,next_free_node,unit)

USE type_specifications
USE general_module
USE constants
USE frequency_spec
USE filter_module
USE Sfilter_fit_module

IMPLICIT NONE

type(Sfilter),intent(IN)       :: Hin                     ! rational transfer function model 
real(dp),intent(IN)            :: gain_in
character(LEN=spice_name_length),intent(IN) :: Hname

integer,intent(IN)   :: in1,in2    ! controlling nodes
integer,intent(IN)   :: on1,on2    ! controlled (output) nodes
integer,intent(IN)   :: vref_node  ! local reference node
integer,intent(INOUT):: next_free_node  ! next free node node
integer,intent(IN)   :: unit       ! file unit for Spice sub-circuit

! local variables

type(Sfilter)    :: H                     ! rational transfer function model 
type(Sfilter)    :: HR                    ! rational transfer function model 
type(Sfilter_PR) :: H_PR                  ! pole-residue transfer function model 
type(Sfilter)    :: T1                    ! temporary filter function
type(Sfilter)    :: Rdc_filter            ! temporary filter function

integer        :: aorder,border,max_order,CF_dim
integer        :: n_branches

integer        :: test

integer        :: i,loop,type_loop                          ! loop variables

type(Polynomial)  :: num
type(Polynomial)  :: den

real(dp) :: value

real(dp),allocatable :: CFterm(:,:)
integer,allocatable  :: CFtype(:)

real(dp)    :: wmin,wmax,wstep
integer     :: nw
complex(dp) :: s,sn,num_fs,den_fs,CF_term,last_CF_term
complex(dp) :: H_rational,H_CF
real(dp)    :: R,L,C

real(dp)    :: R_min     ! minimum value of the resistance of the input function
real(dp)    :: w_R_min     ! angular frequency for minimum value of the resistance of the input function
real(dp)    :: R_add     ! additional resistance required to make the function positive real

logical :: stable
logical :: found
logical :: remainder_OK
logical :: remainder_zero
logical :: multiple_poles

integer           :: type,last_type
real(dp)          :: term

integer :: on,od
logical :: pole_at_zero
logical :: zero_at_zero

real(dp) :: ascale,bscale,scale
real(dp) :: gain

! START

  if (verbose) then
    write(*,*)'******************************'
    write(*,*)'CALLED Equivalent_Circuit_Test'
    write(*,*)'******************************'
  end if
!  verbose=.TRUE.
 
  if (verbose) then
    write(*,*)'CALLED with function:'
    CALL write_Sfilter(Hin,0)
  end if
 
! copy the input filter and renormalise it 
! Set an appropriate normalisation for the filter function 
! so that the coefficients don't get out of hand.

  H=renormalise_Sfilter(Hin)
  
  if (verbose) then
    write(*,*)'Renormalised filter:'
    CALL write_Sfilter(H,0)
  end if

!  CALL get_min_order_poly(H%a)
!  CALL get_min_order_poly(H%b)

! get a 'scale' for the filter function and include this in the gain term
! to prevent component values getting too big or small

  ascale=abs(H%a%coeff(H%a%order))
  bscale=abs(H%b%coeff(H%b%order))
    
  H%a%coeff(:)=H%a%coeff(:)/ascale
  H%b%coeff(:)=H%b%coeff(:)/bscale
  
  scale=ascale/bscale
  
  gain=gain_in*scale
  
  if (verbose) then
    write(*,*)'Scaled function:'
    CALL write_Sfilter(H,0)
    write(*,*)'Calculate minimum resistance of function:'
  end if
  
  CALL calculate_min_resistance_value(H,R_min,w_R_min)
  
  if (verbose) then
    write(*,*)'Minimum Resistance value is',R_min
    write(*,*)'at w=',w_R_min,' f=',w_R_min/(2d0*pi)
  end if
  
  if (R_min.LT.0d0) then
  
    if (verbose) write(*,*)'Adding 1.5*abs(Minimum Resistance value) to H'
  
    R_add=1.5d0*abs(R_min)
    Rdc_filter=allocate_Sfilter(0,0)
    Rdc_filter%wnorm=Hin%wnorm
    Rdc_filter%a%coeff(0)=R_add
    Rdc_filter%b%coeff(0)=1d0
  
    CALL deallocate_Sfilter(T1)
    T1=H+Rdc_filter
    CALL deallocate_Sfilter(H)
    H=T1
    CALL deallocate_Sfilter(T1)
    CALL deallocate_Sfilter(Rdc_filter)

  else
  
    R_add=0d0

  end if
  
  if (verbose) then
    write(*,*)'Revised H:'
    CALL write_Sfilter(H,0)
  end if
  
! now we have ensured that Re(H)>=0 at all frequencies, do the checks again.
! Check the transfer funcion for stability and for whether it is positive real

  CALL check_transfer_function(H,stable) 
  
  if (stable) then  
    if (verbose) write(*,*)'INPUT FUNCTION IS A STABLE, PHYSICAL IMPEDANCE'
  else  
    if (verbose) write(*,*)'INPUT FUNCTION IS NOT STABLE'
    run_status='INPUT FUNCTION IS NOT STABLE, even after adding a stabilising d.c. resistance'
    CALL write_program_status()
    STOP   
  end if
    
! Max_order is the maximum number of components, including  
! resistive and reactive components

  max_order=2*max(H%a%order,H%b%order) +1
  
  if (verbose) then
    write(*,*)'Maximum order is estimated to be ',max_order
  end if
  
! allocate the continued fraction data  
  CF_dim=max_order
  allocate( CFterm(1:max_order,1:5) ) ! note 5 terms
  CFterm(:,:)=0d0
  allocate( CFtype(1:max_order) )
  CFtype(:)=0
  
  n_branches=0
  type=type_impedance  ! H(s) is an impedance function to start with
  
  do loop=1,max_order
  
    remainder_OK=.FALSE.
    remainder_zero=.FALSE.

! Loop for trying both impedance and admittance functions   
    do type_loop=1,2
    
      if (verbose) then
        write(*,*)'Stage ',loop,' of ',max_order       
        if (type.EQ.type_impedance) then
          write(*,*)'TRYING TO CALCULATE AN IMPEDANCE'        
        else
          write(*,*)'TRYING TO CALCULATE AN ADMITTANCE'
         end if
      end if

! check the function for multiple poles. If there are no multiple poles then we can
! go on and look for viable branches

      CALL check_for_multiple_roots(H%b,multiple_poles,.TRUE.)
      
      if (.NOT.multiple_poles) then

! Calculate the partial fraction expansion of the function H(s)
        H_PR=Convert_filter_S_to_S_PR(H)
        if (verbose) CALL write_S_PR_filter(H_PR)
      
        do test=1,8
           
! Test number 1: looking for RLC branch
          select case (test)
        
          case(1)
! Test number 1: looking for RLC branch
            CALL RLC_test(H_PR,type,CFtype(loop),CFterm(loop,1),CFterm(loop,2),CFterm(loop,3),found, &
                          HR,remainder_OK,remainder_zero)
                              
          case(2)
! Test number 2: looking for LC branch
            CALL LC_test(H_PR,type,CFtype(loop),CFterm(loop,1),CFterm(loop,2),CFterm(loop,3),found, &
                          HR,remainder_OK,remainder_zero)
       
          case(3)
! Test number 3: looking for RC branch
            CALL RC_test(H_PR,type,CFtype(loop),CFterm(loop,1),CFterm(loop,2),CFterm(loop,3),found, &
                          HR,remainder_OK,remainder_zero)
      
          case(4)
! Test number 4: looking for RL branch
            CALL RL_test(H_PR,type,CFtype(loop),CFterm(loop,1),CFterm(loop,2),CFterm(loop,3),found, &
                         HR,remainder_OK,remainder_zero)
     
          case(5)
! Test number 5: looking for C branch
            CALL C_test(H_PR,type,CFtype(loop),CFterm(loop,1),CFterm(loop,2),CFterm(loop,3),found, &
                         HR,remainder_OK,remainder_zero)
     
          case(6)
! Test number 6: looking for L branch
            CALL L_test(H_PR,type,CFtype(loop),CFterm(loop,1),CFterm(loop,2),CFterm(loop,3),found, &
                         HR,remainder_OK,remainder_zero)
     
          case(7)
! Test number 7: looking for R branch
            CALL R_test(H_PR,type,CFtype(loop),CFterm(loop,1),CFterm(loop,2),CFterm(loop,3),found, &
                         HR,remainder_OK,remainder_zero)
     
          case(8)
! Test number 8: looking for R2 branch
            CALL R2_test(H_PR,type,CFtype(loop),CFterm(loop,1),CFterm(loop,2),CFterm(loop,3),found, &
                         HR,remainder_OK,remainder_zero)

          end select
                    
          if (found.AND.remainder_OK) then ! Adopt this as a viable branch
            n_branches=n_branches+1
            if (remainder_zero) then
              GOTO 2000
            else
              GOTO 1000
            end if
          end if

        end do ! next test
      
      end if

! If we have been unsuccessful finding a circuit element
! change from impedance->admittance of vice vers    

      if (type.EQ.type_impedance) then
        type=type_admittance
      else
        type=type_impedance
      end if
      
! Calculate the reciprocal filter function

      num=H%a
      den=H%b
      CALL deallocate_poly(H%a)
      CALL deallocate_poly(H%b)
      H%a=den
      H%b=num
      CALL deallocate_poly(num)
      CALL deallocate_poly(den)
          
    end do ! next type_loop
    
! if we get here then we have not found a viable way to proceed with
! building the circuit
    if (verbose) write(*,*)'CANNOT FIND VIABLE COMPONENT'
    if (verbose) write(*,*)'TRYING BRUNE CYCLE'

! convert to an impedance function if required
    if (type.EQ.type_admittance) then
      type=type_impedance     
! Calculate the reciprocal filter function
      num=H%a
      den=H%b
      CALL deallocate_poly(H%a)
      CALL deallocate_poly(H%b)
      H%a=den
      H%b=num
      CALL deallocate_poly(num)
      CALL deallocate_poly(den)
    end if
    
    CALL BRUNE_test(H,type,CFtype(loop),CFterm(loop,1),CFterm(loop,2),CFterm(loop,3), &
                                           CFterm(loop,4),CFterm(loop,5),found,          &
                                           HR,remainder_OK,remainder_zero)
                                           
    if (found.AND.remainder_OK) then ! Adopt this as a viable branch
      n_branches=n_branches+1
      if (remainder_zero) then
        GOTO 2000
      else
        GOTO 1000
      end if
    end if
                       
    run_status='ERROR: CANNOT FIND VIABLE COMPONENT'
    CALL write_program_status()
    STOP
    
! jump here when we have found the next circuit element
1000 CONTINUE
        
    if (verbose) write(*,*)'Prepare for the next stage'
    
    CALL deallocate_Sfilter(H)
    H=HR
               
  end do  ! next circuit element
  
! If we end up here then there is a problem because the remainder is not zero
! and we are supposed to have worked out all the circuit elements by now.


! jump here when the continued fraction truncates with a zero remainder
2000 continue

!  CALL write_CF_local(CFterm,CFtype,CF_dim,n_branches)

  if (verbose) write(*,*)''
  
  if (verbose) then

INCLUDE 'include_write_frequency_response.F90'

  end if

! Write a spice circuit model for the ladder network derived from the
! continued fraction expansion

  CALL write_ladder_network(Hname,gain,CFterm,CFtype,CF_dim,n_branches,R_add,nw,wmin,wmax, &
                            in1,in2,on1,on2,vref_node,next_free_node,unit)

! deallocate memory and finish up
  
  CALL deallocate_poly(num)
  CALL deallocate_poly(den)

  CALL deallocate_Sfilter(H)
  CALL deallocate_Sfilter(HR)
  CALL deallocate_Sfilter(T1)
  CALL deallocate_Sfilter(Rdc_filter)
  deallocate( CFterm )
   
  RETURN
  
9000 write(*,*)'Error reading transfer function file'
  STOP
  
END SUBROUTINE network_synthesis