!
! 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
!
!
!
! Description
!     This is a fortran fragment which is included in network_synthesis.F90
!     It is used for testing that the equivalent circuit model reproduces the 
!     required transfer function.
!
! Comments:
!     Note that this doesn't include the Brune sub-circuit. 
!
! History
!
!     16/11/2017 CJS Include network synthesis process to replace s-domain transfer functions
!
!
! evaluate the frequency response of the continued fraction function

  open(unit=30,file='Continued_fraction.fout')

  wstep=(wmax-wmin)/(nw-1)
  
  do i=1,nw
  
    s=(0d0,1d0)*(wmin+wstep*(i-1))
      
    last_CF_term=(0d0,0d0)
    last_type=CFtype(n_branches)
    
    do loop=n_branches,1,-1    ! note evaluate the continued fraction from the bottom up
       
      R=CFterm(loop,1)
      L=CFterm(loop,2)
      C=CFterm(loop,3)
      
      type=CFtype(loop)
      
      if (type.EQ.series_RLC) then  
       
        CF_term=s*L/(s*s*L*C+s*L/R+1D0)
    
      else if (type.EQ.series_LC) then   
       
        CF_term=s*L/(s*s*L*C+1d0)
        
      else if (type.EQ.series_RC) then   
       
        CF_term=R/(s*C*R+1d0)
    
      else if (type.EQ.series_RL) then   
       
        CF_term=s*L/(s*L/R+1D0)
    
      else if (type.EQ.series_C) then   
       
        CF_term=1d0/(s*C)
    
      else if (type.EQ.series_L) then   
       
        CF_term=s*L
    
      else if (type.EQ.series_R) then   
       
        CF_term=R    
   
! admittance terms
    
      else if (type.EQ.shunt_RLC) then   
        
        CF_term=s*C/(s*s*L*C+s*R*C+1D0)
      
      else if (type.EQ.shunt_LC) then   
       
        CF_term=s*C/(s*s*L*C+1D0)
    
      else if (type.EQ.shunt_RC) then   
       
        CF_term=s*C/(s*C*R+1D0)
    
      else if (type.EQ.shunt_RL) then   
       
        CF_term=1d0/(s*L+R)
    
      else if (type.EQ.shunt_C) then   
       
        CF_term=(s*C)
    
      else if (type.EQ.shunt_L) then   
       
        CF_term=1d0/(s*L)
    
      else if (type.EQ.shunt_R) then   
    
        CF_term=1d0/R

      end if 
      
      if (loop.NE.max_order) then 
        
        if (last_type*type.LT.0) then
! we switch from impedance to admittance or vice versa in the network
          CF_term=CF_term+1d0/last_CF_term
        else
          CF_term=CF_term+last_CF_term         
        end if
        
      end if
      
      last_CF_term=CF_term
      last_type=type
   
    end do
    
    if (last_type.GT.0d0) then
      H_CF=last_CF_term
    else
      H_CF=1d0/last_CF_term
    end if
    write(30,8030)(wmin+wstep*(i-1))/6.28318530718,real(H_CF),aimag(H_CF)
8030  format(3ES16.6)
  end do
  
  close(unit=30)