!
! 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
!
!
! SUBROUTINE write_ladder_network(Hname,gain,CFterm,CFtype,CF_dim,max_order,R_add,nw,wmin,wmax, &
!                                node_1_in,node_2_in,node_1_out,node_2_out,vref_node,next_free_node,unit)
!
!
! NAME
!     write_ladder_network
!
! DESCRIPTION
!       write spice circuit file for the ladder network. 
!       See chapter 7 of the Theory manual. Figure 7.2 shows the ladder network for a simple impedance
!       Figure 7.10 shows the circuit for a non positive-real transfer function. 
!       The individual branch types are illustrated in figure 7.6 for the series impedance branches 
!       and figure 7.7 for the parallel admittance branches.
!
! SEE ALSO
!
!
! HISTORY
!
!     started 09/2017 CJS
!     16/11/2017 CJS Include network synthesis process to replace s-domain transfer functions
!     7/3/2018 CJS fix bug in which the next_free_node was not incremented if Radd.GT.0 i.e. for non PR functions
!
SUBROUTINE write_ladder_network(Hname,gain,CFterm,CFtype,CF_dim,max_order,R_add,nw,wmin,wmax, &
                                node_1_in,node_2_in,node_1_out,node_2_out,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

integer CF_dim,max_order
character(LEN=spice_name_length),intent(IN) :: Hname
real(dp),intent(IN)            :: gain

real(dp) :: CFterm(1:CF_dim,1:5)
integer  :: CFtype(1:CF_dim)
real(dp) :: R_add
integer  :: nw
real(dp) :: wmin,wmax
integer  :: node_1_in,node_2_in,node_1_out,node_2_out,unit
integer,intent(IN) :: vref_node                  ! spice subcircuit reference node number
integer,intent(INOUT) :: next_free_node          ! next free node number

! local variables

real(dp) :: R,L,C,L2,K
integer  :: type

integer  :: i,loop

integer  :: node1,node2,last_series_node,first_series_node,inode1,inode2,R_add_node

!START

! write the continued fraction to the screen
    
  write(*,*)'Write the spice circuit file for the ladder network'
  write(*,*)'      component     type          R           L           C '
  do loop=1,max_order 
    write(*,'(2I12,A5,4ES12.4)')loop,CFtype(loop),'     ',CFterm(loop,1),CFterm(loop,2),CFterm(loop,3)
  end do
  
  write(*,*)'R_add=',R_add
  
  open(unit=40,file='ngspice_circuit.cir')
  
  write(unit,'(A)')'*****Equivalent circuit model for s-domain transfer function *****'
  write(unit,'(A)')'* generated by network_syntheis.F90'
  
! Voltage Controlled current source
  write(unit,'(A)')'*'
  next_free_node=next_free_node+1
  write(spice_model_file_unit,'(A,2I6,2I6,ES16.6)')'GS1_'//trim(Hname),  &
                               vref_node,next_free_node,node_1_in,node_2_in,1.0 
  write(unit,'(A)')'*'
  
  first_series_node=next_free_node
  last_series_node=next_free_node
  next_free_node=next_free_node+1
  
  do loop=1,max_order 
  
    type=CFtype(loop)   

    if (type.GT.0) then
! a series component so connect to a new node or the last component connects to the reference
      node1=last_series_node
      if (loop.NE.max_order) then
        node2=next_free_node
        next_free_node=next_free_node+1
      else
        node2=vref_node  
      end if
      last_series_node=node2
    else if (type.LT.0) then
      node1=last_series_node
      node2=vref_node   
    end if
    
    R=CFterm(loop,1)
    L=CFterm(loop,2)
    C=CFterm(loop,3)
    L2=CFterm(loop,4)
    K=CFterm(loop,5)
    
    if (type.EQ.series_RLC) then   

      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'L'//trim(Hname),loop,' ',node1,node2,L
      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'C'//trim(Hname),loop,' ',node1,node2,C
      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'R'//trim(Hname),loop,' ',node1,node2,R
    
    else if (type.EQ.series_LC) then   

      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'L'//trim(Hname),loop,' ',node1,node2,L
      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'C'//trim(Hname),loop,' ',node1,node2,C
    
    else if (type.EQ.series_RC) then   

      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'C'//trim(Hname),loop,' ',node1,node2,C
      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'R'//trim(Hname),loop,' ',node1,node2,R
    
    else if (type.EQ.series_RL) then   

      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'L'//trim(Hname),loop,' ',node1,node2,L
      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'R'//trim(Hname),loop,' ',node1,node2,R
    
    else if (type.EQ.series_C) then   

      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'C'//trim(Hname),loop,' ',node1,node2,C
      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'R'//trim(Hname),loop,' ',node1,node2,1E8
    
    else if (type.EQ.series_L) then   

      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'L'//trim(Hname),loop,' ',node1,node2,L
    
    else if (type.EQ.series_R) then   

      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'R'//trim(Hname),loop,' ',node1,node2,R
    
    else if (type.EQ.series_BRUNE) then   

      next_free_node=next_free_node+1
      inode1=next_free_node
      next_free_node=next_free_node+1
      inode2=next_free_node
      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'R'//trim(Hname),loop   ,' ',node1,inode1,R
      next_free_node=next_free_node+1
      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'L1'//trim(Hname),loop  ,' ',inode1,inode2,L
      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'C'//trim(Hname),loop   ,' ',inode2,0 ,C
      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'Rbig'//trim(Hname),loop,' ',inode2,0 ,1E8
      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'L2'//trim(Hname),loop  ,' ',node2,inode2,L2
! note directions of inductor specifications      
      write(unit,'(A,I2.2,A,I2.2,A,I2.2,A,ES12.4)')'K'//trim(Hname),loop ,' L1'//trim(Hname),  &
                                             loop,' L2'//trim(Hname),loop,' ',K
      
    
    else if (type.EQ.shunt_RLC) then   

      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'L'//trim(Hname),loop,' ',node1      ,next_free_node ,L
      next_free_node=next_free_node+1
      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'C'//trim(Hname),loop,' ',next_free_node-1,next_free_node ,C
      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'R'//trim(Hname),loop,' ',next_free_node,node2     ,R
      next_free_node=next_free_node+1
    
    else if (type.EQ.shunt_LC) then   

      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'L'//trim(Hname),loop,' ',node1      ,next_free_node ,L
      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'C'//trim(Hname),loop,' ',next_free_node,node2     ,C
      next_free_node=next_free_node+1
    
    else if (type.EQ.shunt_RC) then   

      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'R'//trim(Hname),loop,' ',node1      ,next_free_node ,R
      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'C'//trim(Hname),loop,' ',next_free_node,node2     ,C
      next_free_node=next_free_node+1
      
    else if (type.EQ.shunt_RL) then   

      next_free_node=next_free_node+1
      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'R'//trim(Hname),loop,' ',node1      ,next_free_node ,R
      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'L'//trim(Hname),loop,' ',next_free_node,node2     ,L
      next_free_node=next_free_node+1
      
    else if (type.EQ.shunt_L) then   

      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'L'//trim(Hname),loop,' ',node1,node2     ,L
      
    else if (type.EQ.shunt_C) then   

      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'C'//trim(Hname),loop,' ',node1,node2     ,C
      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'R'//trim(Hname),loop,' ',node1,node2,1E8
      
    else if (type.EQ.shunt_R) then   

      write(unit,'(A,I2.2,A,I6,I6,ES12.4)')'R'//trim(Hname),loop,' ',node1,node2     ,R
   
    end if
           
    write(unit,'(A)')'*'
           
  end do  ! next element of the network
  write(unit,'(A)')'*'
  
  if (R_add.GT.0d0) then
! we need an additional branch for R_add
    write(*,*)'Adding additional resistance branch'
    R_add_node=next_free_node
    next_free_node=next_free_node+1
    write(unit,'(A)')'*'
    write(spice_model_file_unit,'(A,2I6,2I6,ES16.6)')'GS2_'//trim(Hname),  &
                               vref_node,R_add_node,node_1_in,node_2_in,1.0 
    write(unit,'(A)')'*'      
    write(unit,'(A,I6,I6,ES12.4)')'R_add'//trim(Hname),vref_node,R_add_node,R_add
  else
    R_add_node=Vref_node
  end if
  
! controlled source including the gain term
  if (R_add.GT.0d0) then
    write(spice_model_file_unit,'(A,2I6,2I6,ES16.6)')'ES2_'//trim(Hname),             &
                                                     node_1_out,node_2_out,   & 
                                                     first_series_node,R_add_node,gain
  else
    write(spice_model_file_unit,'(A,2I6,2I6,ES16.6)')'ES2_'//trim(Hname),             &
                                                     node_1_out,node_2_out,   & 
                                                     first_series_node,vref_node,gain
  
  end if

  write(unit,'(A)')'*'
  
  RETURN
  END SUBROUTINE write_ladder_network