! 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 write_incident_field_circuit
!
! NAME
!     write_incident_field_circuit
!
! AUTHORS
!     Chris Smartt
!
! DESCRIPTION
!     This code writes the circuit components required to implement the incident field excitation model
!     The Spice model is described in Theory_Manual_Section 3.8 
!     and the Spice circuit is seen in Theory_Manual_Figures 3.17, and 3.18
!
!     INPUTS REQUIRED
!     1.    nodes for connection of the incident field sources
!     2.    time delays for the viictim domain modes
!     3.    bundle length
!     4.    modal decomposition matrices
!     5.    victim domain propagation correction filter functions
!     6.    pre-calculated coefficients related to the incident plane wave field excitation 
!     7.    nodes for the specification of the incident field function
!     8.    next free node and reference node for the spice sub-circuit
!
!
!     OUTPUTS
!     The components required to implement the model are written to
!     the subcircuit file
!     
! COMMENTS
!     Write the circuit elements required for the incident field excitation model
!
!     The excitation function comes from the voltage between two subcircuit termination nodes.
!
! HISTORY
!
!     STAGE 6 developments started 16/06/2016 CJS.
!     Include ground plane 28/6/2016 CJS
!     Replace Tz transmission lines with a very small delay by a very small resistance - causes long runtimes in Pspice. 6/12/2016.
!     9/5/2017 CJS Document software with reference to Theory_Manual 

!     15/10/2017 CJS Include source scaling to keep the voltages and currents in a sensible range
!     20/10/2017 CJS call subroutine to write delay lines so we can choose to use LTRA or T elements
!     16/11/2017 CJS Include network synthesis process to replace s-domain transfer functions

!
SUBROUTINE write_incident_field_circuit(n_victim_domain_modes, &
                             Vv_end1_node,Vv_end2_node, &
                             T_victim, &
                             length, &
                             TVI_victim,TI_victim, &
                             Hpv_filter, &
                             Tz,alpha0,alphaL, &
                             next_free_node,Einc_node1,Einc_node2,vref_node,victim_domain )

USE type_specifications
USE general_module
USE constants
USE spice_cable_bundle_module
USE filter_module
USE maths

IMPLICIT NONE

! variables passed to the subroutine

integer,intent(IN) :: n_victim_domain_modes   ! number of modes in the victim domain

! victim domain node lists for connection of the incident field source terms at each end.
integer,intent(IN) :: Vv_end1_node(1:n_victim_domain_modes)
integer,intent(IN) :: Vv_end2_node(1:n_victim_domain_modes)

! mode delays for the victim domain transmission line modes
real(dp),intent(IN) :: T_victim(1:n_victim_domain_modes)

! modal decomposition matrices for the victim domain
real(dp),intent(IN) :: TVI_victim(1:n_victim_domain_modes,1:n_victim_domain_modes)
real(dp),intent(IN) :: TI_victim(1:n_victim_domain_modes,1:n_victim_domain_modes)
 
real(dp),intent(IN) :: length                                   ! bundle length
TYPE(Sfilter),intent(IN) :: Hpv_filter(1:n_victim_domain_modes) ! victim mode proopagation correction terms

real(dp),intent(IN) :: Tz                                    ! incident field excitation delay time
real(dp),intent(IN) :: alpha0(1:n_victim_domain_modes)       ! pre-calculated model constants 
real(dp),intent(IN) :: alphaL(1:n_victim_domain_modes)       ! pre-calculated model constants 

integer,intent(INOUT) :: next_free_node      ! the next free node number for the spice sub-circuit 
integer,intent(IN) :: Einc_node1             ! excitation function sub-circuit terminal node 1
integer,intent(IN) :: Einc_node2             ! excitation function sub-circuit terminal node 2 
integer,intent(IN) :: vref_node              ! spice sub-circuit reference node number
integer,intent(IN) :: victim_domain          ! victim domain number

! local variables

integer :: v_mode

real(dp) :: Z0_delay=50d0   ! transmission line impedance for pure delay lines

! Names for delay lines

character(len=spice_name_length) :: delay_line_Tz_name
character(len=spice_name_length) :: delay_line_Tv_name(1:n_victim_domain_modes)
character(len=spice_name_length) :: delay_line_TzPTv_name(1:n_victim_domain_modes)

character(len=spice_name_length) :: delay_line_ZC_Tz_name
character(len=spice_name_length) :: delay_line_ZC_Tv_name(1:n_victim_domain_modes)
character(len=spice_name_length) :: delay_line_ZC_TzPTv_name(1:n_victim_domain_modes)

character(len=spice_name_length) :: delay_line_E1_Tz_name
character(len=spice_name_length) :: delay_line_E1_Tv_name(1:n_victim_domain_modes)
character(len=spice_name_length) :: delay_line_E1_TzPTv_name(1:n_victim_domain_modes)

character(len=spice_name_length) :: E_Einc_l_name(1:n_victim_domain_modes)
character(len=spice_name_length) :: E_Einc_s_name(1:n_victim_domain_modes)

character(len=spice_name_length) :: E_Einc_filter_name(1:n_victim_domain_modes)

! names for special case Tz=Tvictim components
character(len=spice_name_length) :: G_Einc_derivative_name(1:n_victim_domain_modes)
character(len=spice_name_length) :: L_Einc_derivative_name(1:n_victim_domain_modes)

! working strings

character(len=spice_name_length) :: name1,name2,name3
character(len=spice_name_length) :: Einc_string

! String used if a transmission line delay is very small and the T element is replaced by a resistance
character(len=spice_name_length) :: Rtempstring

! circuit to combine incident field terms

character(len=spice_name_length) :: combine_delays_s_E_name(1:n_victim_domain_modes,1:2)
character(len=spice_name_length) :: R_combine_delays_s_name(1:n_victim_domain_modes)

character(len=spice_name_length) :: combine_delays_l_E_name(1:n_victim_domain_modes,1:2)
character(len=spice_name_length) :: R_combine_delays_l_name(1:n_victim_domain_modes)

real(dp)       :: Evalue

! Nodes for delay lines 

! delay line Tz 
integer :: delay_line_Tz_s_nodes
integer :: delay_line_Tz_l_nodes

! delay line Tv 
integer :: delay_line_Tv_s_nodes(1:n_victim_domain_modes)
integer :: delay_line_Tv_l_nodes(1:n_victim_domain_modes)

! delay line Tz+Tv 
integer :: delay_line_TzPTv_s_nodes(1:n_victim_domain_modes)
integer :: delay_line_TzPTv_l_nodes(1:n_victim_domain_modes)

integer :: last_Hjw_source_s_node(1:n_victim_domain_modes)
integer :: last_Hjw_source_l_node(1:n_victim_domain_modes)

integer :: combine_delays_s_Enode
integer :: combine_delays_l_Enode

real(dp) :: gain            ! gain value for s-domain transfer functions

real(dp) :: Tz_minus_Tv     ! difference between source and victim mode velocities
logical  :: Tz_equal_Tv     ! if incident field delay and victim mode delay are the same then we have a different circuit topology

! loop variables
integer :: row,i

! nodes for special case Tz=Tvictim components
integer :: Einc_derivative_node(1:n_victim_domain_modes)

real(dp) :: source_scale     ! scaling factor for controlled sources to ensure that 
                             ! voltage/ currents stay in a sensible range

! START

if (verbose) write(*,*)'CALLED write_incient_field_excitation_circuit'

source_scale=c0  

! check for the case where Tz<0 and stop with an error for now. 
! Eventually we will cope with this situation by reversing the model.
if (Tz.LT.0d0) then
  run_status='ERROR in write_incient_field_excitation_circuit. Tz is less than zero'
  CALL write_program_status()
  STOP 1
end if

CALL write_spice_comment('START OF INCIDENT FIELD EXCITATION MODELS')

! Set Tz delay line node numbers

! delay line nodes for positive z propagation, source end
next_free_node=next_free_node+1
delay_line_Tz_s_nodes=next_free_node

! delay line nodes for positive z propagation, load end
next_free_node=next_free_node+1
delay_line_Tz_l_nodes=next_free_node
      
! Set Tz delay line component names

delay_line_Tz_name='T_Tz_Einc'
  
! mode impedance
  
delay_line_ZC_Tz_name='RZC_Tz_Einc'

! source terms
  
delay_line_E1_Tz_name='E1_Tz_Einc'
      
! Write Tz delay line components
! Theory_Manual_Fig 3.17, Theory_Manual_Eqn 3.168b

if (Tz.GT.Tz_min_delay) then
! significant delay so use a T element
  
  CALL write_spice_comment('Incident field delay lines, Tz delay') 
  

!  write(spice_model_file_unit,'(A30,4I6,A4,E16.6,A4,E16.6)')delay_line_Tz_name,  &
!                                           delay_line_Tz_s_nodes,vref_node, &
!                                           delay_line_Tz_l_nodes,vref_node, &
!                                           ' Z0=',Z0_delay,' TD=',Tz
  CALL write_delay_line(delay_line_Tz_name,  &
                        delay_line_Tz_s_nodes,vref_node, &
                        delay_line_Tz_l_nodes,vref_node, &
                        Z0_delay,Tz,length)
                                           

else
! the delay is very small so just use a very small series resistance instead of a delay line
  
  CALL write_spice_comment('Incident field delay lines, Tz delay=0. Delay line replaced by small resistance') 

  Rtempstring='R'//trim(delay_line_Tz_name)
  write(spice_model_file_unit,'(A30,2I6,E16.6)')Rtempstring,  &
                                           delay_line_Tz_s_nodes,delay_line_Tz_l_nodes,Rsmall
end if
 
! modal impedances on modal delay lines, Tz
 
CALL write_spice_comment('Matched impedance: Tz delay')
  
write(spice_model_file_unit,'(A30,2I6,E16.6)')delay_line_ZC_Tz_name, &
                             delay_line_Tz_l_nodes,vref_node,Z0_delay

! delay line controlled source for positive z propagation, Ts
  
CALL write_spice_comment('Delay line controlled sources Tz delay') 

write(spice_model_file_unit,'(A30,4I6,E16.6)')delay_line_E1_Tz_name,&
                            delay_line_Tz_s_nodes,vref_node, &
                            Einc_node1,Einc_node2,1.0 
 
! End of Tz delay line

do v_mode=1,n_victim_domain_modes ! loop over victim domain modes
  
! create Einc_string which labels the transfer impedance model number plus the source mode and victim mode numbers
  name1='EINC_vm_'
  CALL add_integer_to_string(name1,v_mode,Einc_string)

! we always need the Tz+T_victim(v_mode) delay lines so write these components now
! Theory_Manual_Fig 3.17, Theory_Manual_Eqn 3.168a

! Set TzPTv delay line nodes

! delay line nodes, source end
  next_free_node=next_free_node+1
  delay_line_TzPTv_s_nodes(v_mode)=next_free_node

! delay line nodes, load end
  next_free_node=next_free_node+1
  delay_line_TzPTv_l_nodes(v_mode)=next_free_node

! Set TzPTv delay line component names

  delay_line_TzPTv_name(v_mode)='T_TzPTv_'//trim(Einc_string)
  
! mode impedance
  
  delay_line_ZC_TzPTv_name(v_mode)='RZC_TzPTv_'//trim(Einc_string)

! source terms
  
  delay_line_E1_TzPTv_name(v_mode)='E1_TzPTv_'//trim(Einc_string)  
    
! Write TzPTv delay lines
! Theory_Manual_Fig 3.17, Theory_Manual_Eqn 3.168a

  CALL write_spice_comment('Delay line Tz+T_victim(v_mode)') 
  

!  write(spice_model_file_unit,'(A30,4I6,A4,E16.6,A4,E16.6)')delay_line_TzPTv_name(v_mode),&
!                                       delay_line_TzPTv_s_nodes(v_mode),vref_node, &
!                                       delay_line_TzPTv_l_nodes(v_mode),vref_node, &
!                                       ' Z0=',Z0_delay,' TD=',Tz+T_victim(v_mode)
  CALL write_delay_line(delay_line_TzPTv_name(v_mode),&
                        delay_line_TzPTv_s_nodes(v_mode),vref_node, &
                        delay_line_TzPTv_l_nodes(v_mode),vref_node, &
                        Z0_delay,Tz+T_victim(v_mode),length)

! modal impedances on modal delay lines, TzPTv
  
  CALL write_spice_comment('Matched impedance:  Tz+T_victim(v_mode) delay') 
  
  write(spice_model_file_unit,'(A30,2I6,E16.6)')delay_line_ZC_TzPTv_name(v_mode),        &
                               delay_line_TzPTv_l_nodes(v_mode),vref_node,Z0_delay

! delay line controlled source, TzPTv

  CALL write_spice_comment('Incident field delay line controlled source:  Tz+T_victim(v_mode) delay')
  
  write(spice_model_file_unit,'(A30,4I6,E16.6)')delay_line_E1_TzPTv_name(v_mode),&
                               delay_line_TzPTv_s_nodes(v_mode),vref_node, &
                               Einc_node1,Einc_node2,1.0

! End of TzPTv delay lines      

! Test for the special case when the source and victim mode delays are the same (or very close)
! In this case we need to use a different model to avoid a singularity in the normal model
! Theory_Manual_Section 3.8.1, Theory_Manual_Eqn 3.169

  Tz_minus_Tv=Tz-T_victim(v_mode)

  if (abs(Tz_minus_Tv).GT.Einc_min_delay) then
  
    Tz_equal_Tv=.FALSE.
    
! The whole incident field excitation circuit is implemented using delay lines
! Theory_Manual_Fig 3.17

! Set Tv delay line nodes

! delay line nodes source end
    next_free_node=next_free_node+1
    delay_line_Tv_s_nodes(v_mode)=next_free_node

! delay line nodes load end
    next_free_node=next_free_node+1
    delay_line_Tv_l_nodes(v_mode)=next_free_node

! Set Tv delay line component names

    delay_line_Tv_name(v_mode)='T_Tv_'//trim(Einc_string) 
  
! mode impedance
  
    delay_line_ZC_Tv_name(v_mode)='RZC_Tv_'//trim(Einc_string)

! source terms
  
    delay_line_E1_Tv_name(v_mode)='E1_Tv_'//trim(Einc_string)
  
! Write Tv delay lines

    CALL write_spice_comment('Incident field delay line T_victim(v_mode)')
    

!    write(spice_model_file_unit,'(A30,4I6,A4,E16.6,A4,E16.6)')delay_line_Tv_name(v_mode), &
!                                           delay_line_Tv_s_nodes(v_mode),vref_node, &
!                                           delay_line_Tv_l_nodes(v_mode),vref_node, &
!                                           ' Z0=',Z0_delay,' TD=',T_victim(v_mode)
    CALL write_delay_line(delay_line_Tv_name(v_mode), &
                          delay_line_Tv_s_nodes(v_mode),vref_node, &
                          delay_line_Tv_l_nodes(v_mode),vref_node, &
                          Z0_delay,T_victim(v_mode),length)

! modal impedances on modal delay lines, T_victim(v_mode)

  
    CALL write_spice_comment('Matched impedance T_victim(v_mode) delay')
  
    write(spice_model_file_unit,'(A30,2I6,E16.6)')delay_line_ZC_Tv_name(v_mode),   &
                             delay_line_Tv_l_nodes(v_mode),vref_node,Z0_delay
 
! delay line controlled source for positive z propagation, T_victim(v_mode)
  
    CALL write_spice_comment('Controlled source for T_victim(v_mode) delay') 
  
    write(spice_model_file_unit,'(A30,4I6,E16.6)')delay_line_E1_Tv_name(v_mode)      &
                              ,delay_line_Tv_s_nodes(v_mode),vref_node, &
                               Einc_node1,Einc_node2,1.0 
       
  else
! Special case Tz=T_victim(v_mode)as see in Theory_Manual_Eqn 3.169  Theory_Manual_Fig 3.18 

    Tz_equal_Tv=.TRUE.
      
! **** The special case required time derivative circuits operating on the delayed incident field excitation****
! new node for time derivative delayed incident field function
    next_free_node=next_free_node+1
    Einc_derivative_node(v_mode)=next_free_node
  
! ****** Names for the special case circuit for Tsource=Tvictim  
    G_Einc_derivative_name(v_mode)='G_Vp_ddt_'//trim(Einc_string)
  
    L_Einc_derivative_name(v_mode)='L_Vp_ddt_'//trim(Einc_string)

! ***** Inductive circuit to calculate the time derivative of delayed Einc *****
  
    CALL write_spice_comment('Controlled source for derivative of delayed Einc i.e. Einc(0,t-Tz)') 

    write(spice_model_file_unit,'(A30,4I6,E16.6)')G_Einc_derivative_name(v_mode), &
                             Einc_derivative_node(v_mode),vref_node, &
                             delay_line_Tz_l_nodes,vref_node,1.0 
                       

    CALL write_spice_comment('Inductor for derivative  of delayed Einc i.e. Einc(0,t-Tz)') 

  
    write(spice_model_file_unit,'(A30,2I6,E16.6)')L_Einc_derivative_name(v_mode), &

                                 Einc_derivative_node(v_mode),vref_node,1d0/source_scale

  end if ! Special case Tz-Tv =0
          
! The remaining part of the circuit combines all the contributions to the victim mode voltage source
! So as to evaluate the incident field terms in Theory_Manual_Eqn 3.168a,b
! We create the nodes for the summation circuit as we go
    
! create Einc_string which labels the transfer impedance model number plus victim mode number
  name1='EINC_vm_'
  CALL add_integer_to_string(name1,v_mode,Einc_string)

! first add the contributions which are common to both forms of circuit  

! Incident field excitation voltage source names      

  combine_delays_s_E_name(v_mode,1)='E_zt_dsum_s_'//trim(Einc_string)//'_E1'
  combine_delays_s_E_name(v_mode,2)='E_zt_dsum_s_'//trim(Einc_string)//'_E2'
  
 ! START OF CIRCUIT TO COMBINE INCIDENT FIELD EXCITATION TERMS

!  calculation of V_victim at z=0. Theory_Manual_Eqn 3.168a
 
  CALL write_spice_comment('Circuit to combine incident field excitation terms') 

! Einc, no delay:   Theory_Manual_Eqn 3.168a, line 2, term 1.
  next_free_node=next_free_node+1
  Evalue=+alpha0(v_mode)/(Tz+T_victim(v_mode))
  write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_s_E_name(v_mode,1) &
                                ,next_free_node,vref_node &
                                ,Einc_node1,Einc_node2 &
                                ,Evalue
  combine_delays_s_Enode=next_free_node
    
! Einc, delay=T_victim+Tz:  Theory_Manual_Eqn 3.168a, line 2, term 2
  next_free_node=next_free_node+1
  Evalue=-alpha0(v_mode)/(Tz+T_victim(v_mode))
  write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_s_E_name(v_mode,2) &
                                 ,next_free_node,combine_delays_s_Enode &
                                 ,delay_line_TzPTv_l_nodes(v_mode),vref_node &
                                 ,Evalue 
  combine_delays_s_Enode=next_free_node
    
!  calculation of V_victim at z=L. Theory_Manual_Eqn 3.168 a
                             			      			      
  if (.NOT.Tz_equal_Tv) then

! normal form based on delay lines,     Theory_Manual_Eqn 3.168b
    combine_delays_l_E_name(v_mode,1)='E_zt_dsum_l_'//trim(Einc_string)//'_E1'
    combine_delays_l_E_name(v_mode,2)='E_zt_dsum_l_'//trim(Einc_string)//'_E2'

! Einc, delay=T_victim: Theory_Manual_Eqn 3.168b, line 2, term 1.
    next_free_node=next_free_node+1
    Evalue=+alphaL(v_mode)/(T_victim(v_mode)-Tz)
    write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_l_E_name(v_mode,1) &
                                ,next_free_node,vref_node &
                                ,delay_line_Tv_l_nodes(v_mode),vref_node &
                                ,Evalue
    combine_delays_l_Enode=next_free_node
    
! Einc, delay=Tz: Theory_Manual_Eqn 3.168b, line 2, term 2.
    next_free_node=next_free_node+1
    Evalue=-alphaL(v_mode)/(T_victim(v_mode)-Tz)
    write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_l_E_name(v_mode,2) &
                                 ,next_free_node,combine_delays_l_Enode &
                                 ,delay_line_Tz_l_nodes,vref_node &
                                 ,Evalue 
    combine_delays_l_Enode=next_free_node

  else
! we need to calcuate the contribution to z=L using the time derivative of the delayed incident field  
! Theory_Manual_Eqn 3.169, Theory_Manual_Fig 3.18

    combine_delays_l_E_name(v_mode,1)='E_zt_dsum_l_'//trim(Einc_string)//'_E1'
     
    next_free_node=next_free_node+1
    Evalue=-alphaL(v_mode)

    Evalue=Evalue*source_scale

    write(spice_model_file_unit,'(A30,4I6,E16.6)')combine_delays_l_E_name(v_mode,1) &
                                 ,next_free_node,combine_delays_l_Enode  &
                                 ,Einc_derivative_node(v_mode),vref_node &
                                 ,Evalue 
    combine_delays_l_Enode=next_free_node
    
  end if  ! special case T_victim=Tz
    
  CALL write_spice_comment('Large resistance to complete the circuit for the series voltage sources') 
  

! Resistance to complete the circuit for the series voltage sources

  R_combine_delays_l_name(v_mode)='R_Einc_dsum_l_'//trim(Einc_string)

  write(spice_model_file_unit,'(A30,2I6,E16.6)')R_combine_delays_l_name(v_mode),combine_delays_l_Enode,vref_node,Rcombine_sources

! Resistance to complete the circuit for the series voltage sources

  R_combine_delays_s_name(v_mode)='R_Einc_dsum_s_'//trim(Einc_string)

  write(spice_model_file_unit,'(A30,2I6,E16.6)')R_combine_delays_s_name(v_mode),combine_delays_s_Enode,vref_node,Rcombine_sources 

  
! Write filter function for propagation correction filter 
! Theory_Manual_Equation 3.170

  CALL write_spice_comment('Incident field excitation sources, end 1')
  
  E_Einc_s_name(v_mode)='Einc_s_'//trim(Einc_string)
  CALL write_s_domain_controlled_voltage_source(E_Einc_s_name(v_mode),                          &
                                                combine_delays_s_Enode,vref_node,               &
                                                Vv_end1_node(v_mode),Vref_node,                 &
                                                Hpv_filter(v_mode),1d0,vref_node,next_free_node) ! note: gain set to 1.0


  CALL write_spice_comment('Incident field excitation sources, end 2')
  
  E_Einc_l_name(v_mode)='Einc_l_'//trim(Einc_string)
  CALL write_s_domain_controlled_voltage_source(E_Einc_l_name(v_mode),                          &
                                                combine_delays_l_Enode,vref_node,               &
                                                Vv_end2_node(v_mode),Vref_node,                 &
                                                Hpv_filter(v_mode),1d0,vref_node,next_free_node) ! note: gain set to 1.0

end do ! next victim mode

RETURN

END SUBROUTINE write_incident_field_circuit