!
! 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 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