!

! 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:
! PROGRAM cable_model_builder
!
! NAME
!     cable_model_builder
!
! AUTHORS
!     Chris Smartt
!
! DESCRIPTION
!     cable_model_builder collects all the information required to characterise
!     an individual cable from a .cable_spec file. The code performs any calculations necessary to 
!     derive the per-unit-length parameters of any shielded (internal) domains as required. 
!
!     The software supports the following cable types:
!
!     Frequency Dependent Cylindrical conductor with dielectric
!     Frequency Dependent Coaxial cable with transfer impedance and shield surface impedance loss
!     Frequency Dependent Twinax cable with transfer impedance and shield surface impedance loss
!     Frequency Dependent Twisted pair
!     Frequency Dependent Shielded twisted pair with transfer impedance and shield surface impedance loss
!     Frequency Dependent Spacewire with transfer impedance and shield surface impedance loss
!     Frequency Dependent Overshield with transfer impedance and shield surface impedance loss
!     Frequency Dependent flex cable
!     D connector
!
!     The input to the program is the name of a cable specification file. A file name.cable_spec must exist, containing
!     all the data required to specify a cable.
!
!     The output of the cable_model_builder code is a file name.cable which can be used as an input to the 
!     cable_bundle_model_builder software
!
!     The .cable file is placed in a directory specified in the .cable_spec file. This may be the local directory (./) 
!     or another specified path. In this way the software can interact with a library of cable models (MOD).
!
!     The program may be run with the cable name specified in the command line i.e. 'cable_model_builder name'
!     or it the name is absent from the command line, the user is prompted for the name.
!
!
! COMMENTS
!     Updated to V2
!
! HISTORY
!
!     started 25/11/2015 CJS Started STAGE_1 developments 
!     started 21/03/2016 CJS Started STAGE_3 developments 
!     7/4/2016  CJS add twinax, twisted_pair and shielded_twisted_pair cable types
!     27/4/2016 CJS add the conductor impedance (loss) model
!     12/5/2016 CJS Generalise the conductor impedance model to include transfer impedance
!     21/7/2016 CJS add models for FD_twisted_pair ZT_FD_shielded_twisted_pair ZT_FD_spacewire ZT_FD_twinax
!     5/9/2016  CJS Include the shielded cable models with surface impedance loss:
!                   FD_coax2, ZT_FD_coax2, FD_twisted_pair2, ZT_FD_shielded_twisted_pair2, ZT_FD_spacewire2 ZT_FD_twinax2
!     23/9/2016 CJS Include rectangular conductor type
!     16/11/2016 CJS Include Dconnector type
!     December 2016 CJS Version 2: Rationalise cable types so that there is only a single version of each type of cable
!     24/2/2017 CJS Allow the input name to include a path i.e. the _spec file does not need to be local.

!     23/10/2017 CJS Check the pole stability of the dielectric and transfer impedance filters.
!     16/11/2017 CJS Include network synthesis process to replace s-domain transfer functions

!     13/3/2018 CJS Add flag for direct/ iterative matrix solver in Laplace solution

!     October 2023 CJS Proximity_effects branch: Embed FastHenry2 frequency dependent impedance matrix calculation

!
!
PROGRAM cable_model_builder

USE type_specifications
USE general_module
USE constants
USE cable_module
USE filter_module

IMPLICIT NONE

! local variables

! command line argument value and length
character(len=filename_length)    :: argument1
integer                           :: argument1_length

character(len=filename_length)    :: cable_name_with_path  ! name of the cable including the path
character(len=filename_length)    :: cable_path            ! path to the cable_spec file
character(len=filename_length)    :: cable_name            ! name of the cable
character(len=filename_length)    :: filename              ! filename of the .cable_spec file
    
logical        :: file_exists

! Structure to hold the cable specification (see CABLE_MODULES/cable_module.F90)
type(cable_specification_type)    :: cable_spec

! string for reading flags
character(len=line_length)    :: line

character(len=line_length)    :: stripped_line

integer        :: ierr   ! integer to return error codes from file reads

! temporary variables for reading number of conductors, parameters, frequency dependent parameters, transfer impedance functions.
integer        :: nc_in,np_in,nFDp_in,nFDZT_in

logical        :: check_passed                   ! flag for stability testing of S domain transfer function filters

! loop variable
integer        :: i

! START

! Open the input file describing the cable parameters
! This file could be created by the associated GUI or otherwise generated

  program_name="cable_model_builder"
  run_status='Started'
  CALL write_program_status()
  
  CALL read_version()
    
  CALL write_license()

! get the first command line argument. If set then this is the cable name, if it is not set then
! it must be read

  CALL get_command_argument(1 , argument1, argument1_length)

  if (argument1_length.NE.0) then
  
    cable_name_with_path =trim(argument1)
  
  else

    write(*,*)'Enter the name of the cable specification data (without .cable_spec extension)'

    read(*,'(A)')cable_name_with_path 

  end if
  
  CALL strip_path(cable_name_with_path,cable_path,cable_name)

  filename=trim(cable_name_with_path)//cable_spec_file_extn

  inquire(file=trim(filename),exist=file_exists)
  if (.NOT.file_exists) then
    run_status='ERROR in cable_model_builder, Cannot find the file:'//trim(filename)
    CALL write_program_status()
    STOP 1
  end if 
  
! open and read the file
  
  OPEN(unit=cable_spec_file_unit,file=filename)

  write(*,*)'Opened file:',trim(filename)
  
! set the cable name to be the same as the name of the cable specification data
  cable_spec%cable_name=cable_name

! set the version tag in the cable structure
  cable_spec%version=SPICE_CABLE_MODEL_BUILDER_version

! read the directory for the cable models
  read(cable_spec_file_unit,*)  ! comment line
  read(cable_spec_file_unit,'(A)')MOD_cable_lib_dir
  CALL path_format(MOD_cable_lib_dir)
  
! ensure that the paths exist
  CALL check_and_make_path(MOD_cable_lib_dir)

! read the cable type
  read(cable_spec_file_unit,'(A)')cable_spec%cable_type_string
  
  CALL convert_to_lower_case(cable_spec%cable_type_string,line_length)
  
  if (cable_spec%cable_type_string.eq.'cylindrical') then

    CALL cylindrical_set_parameters(cable_spec)
    
  else if (cable_spec%cable_type_string.eq.'twisted_pair') then

    CALL twisted_pair_set_parameters(cable_spec)
    
  else if (cable_spec%cable_type_string.eq.'coax') then

    CALL coax_set_parameters(cable_spec)
        
  else if (cable_spec%cable_type_string.eq.'twinax') then

    CALL twinax_set_parameters(cable_spec)
    
  else if (cable_spec%cable_type_string.eq.'shielded_twisted_pair') then

    CALL shielded_twisted_pair_set_parameters(cable_spec)

  else if (cable_spec%cable_type_string.eq.'spacewire') then

    CALL spacewire_set_parameters(cable_spec)

  else if (cable_spec%cable_type_string.eq.'original_flex_cable') then

    CALL flex_cable_set_parameters(cable_spec)

  else if (cable_spec%cable_type_string.eq.'flex_cable') then

    CALL ML_flex_cable_set_parameters(cable_spec)

  else if (cable_spec%cable_type_string.eq.'dconnector') then

    CALL Dconnector_set_parameters(cable_spec)

  else if (cable_spec%cable_type_string.eq.'overshield') then

    CALL overshield_set_parameters(cable_spec)

  else
   
    run_status='ERROR in cable_model_builder, Unknown cable type:'//trim(cable_spec%cable_type_string)
    CALL write_program_status()
    STOP 1
      
  end if

! read the number of conductors from the .cable_spec file 

  read(cable_spec_file_unit,*,IOSTAT=ierr)nc_in
  if (ierr.NE.0) then 
    run_status='ERROR reading number of conductors'
    CALL write_program_status()
    STOP 1
  end if
  
  if (cable_spec%tot_n_conductors.NE.0) then ! This cable type has a specific number of conductors so
! check that the number of conductors is consistent with this cable type

    if (nc_in.NE.cable_spec%tot_n_conductors) then
      run_status='ERROR, Wrong number of conductors'
      CALL write_program_status()
      STOP 1
    end if
    
  else
! the number of conductors is set in the cable specification so set the 
! paramter values accordingly

    cable_spec%tot_n_conductors=nc_in
    
! if the number of external conductors is not set then we assume that all condcutors are external (e.g. flex cable)
    if (cable_spec%n_external_conductors.NE.0) then
      cable_spec%n_external_conductors=nc_in
    end if
    
  end if
  
! read the number of parameters from the .cable_spec file 

  read(cable_spec_file_unit,*,IOSTAT=ierr)np_in
  if (ierr.NE.0) then 
    run_status='ERROR reading number of parameters'
    CALL write_program_status()
    STOP 1
  end if
  

  if (cable_spec%n_parameters.NE.0) then
  

! Check that the number of parameters is consistent with this cable type

    if (np_in.NE.cable_spec%n_parameters) then
      run_status='ERROR, Wrong number of parameters'
      CALL write_program_status()
      STOP 1
    end if
  
  else

! for a ML_flex cable the number of parameters depends on the cable specification  
    cable_spec%n_parameters=np_in
    

  end if

! allocate and read parameters

  ALLOCATE( cable_spec%parameters(1:cable_spec%n_parameters) )
  
  do i=1,cable_spec%n_parameters

    read(cable_spec_file_unit,*,IOSTAT=ierr)cable_spec%parameters(i)
    if (ierr.NE.0) then 
      run_status='ERROR reading cable_spec%parameters'
      CALL write_program_status()
      STOP 1
    end if
  
  end do

! read the number of dielectric filters from the .cable_spec file or assume 0 if the information is not there

  read(cable_spec_file_unit,*,IOSTAT=ierr)nFDp_in
  if (ierr.NE.0) then 
! Assume that there are no frequency dependent dielectric models
    backspace(cable_spec_file_unit)
    nFDp_in=0
  end if
  
! Check that the number of dielectric filters is consistent with this cable type
  if (nFDp_in.NE.cable_spec%n_dielectric_filters) then
    run_status='ERROR, Wrong number of dielectric models for this cable type'
    CALL write_program_status()
    STOP 1
  end if
  
  if (cable_spec%n_dielectric_filters.GT.0) then
    ALLOCATE( cable_spec%dielectric_filter(1:cable_spec%n_dielectric_filters) )
    do i=1,cable_spec%n_dielectric_filters
      read(cable_spec_file_unit,*,IOSTAT=ierr)    ! comment line
      CALL read_Sfilter(cable_spec%dielectric_filter(i),cable_spec_file_unit)

      
! test pole stability
      CALL test_filter_pole_stability(cable_spec%dielectric_filter(i),check_passed)
  
      if (.NOT.check_passed) then
        write(*,*)'ERROR in cable_model_builder'
        write(*,*)'The dielectric filter function has unstable poles'
        write(*,*)'Dielectric filter function number:',i
        CALL write_Sfilter(cable_spec%dielectric_filter(i),0)
        write(run_status,'(A,I2)')'ERROR, Unstable poles in dielectric filter function: ',i
        CALL write_program_status()
        STOP 1
      else
        if (verbose) write(*,*)'The dielectric filter function has stable poles: function number ',i    
      end if
          

    end do
  end if
   
! read the number of transfer impedance filters from the .cable_spec file or assume 0 if the information is not there

  read(cable_spec_file_unit,*,IOSTAT=ierr)nFDZT_in
  if (ierr.NE.0) then 
! Assume that there are no transfer impedance filters
    backspace(cable_spec_file_unit)
    nFDZT_in=0
  end if
  
! Check that the number of transfer impedance filters is consistent with this cable type
  if (nFDZT_in.NE.cable_spec%n_transfer_impedance_models) then
    run_status='ERROR, Wrong number of transfer impedance models for this cable type'
    CALL write_program_status()
    STOP 1
  end if
  
  if (cable_spec%n_transfer_impedance_models.GT.0) then
    ALLOCATE( cable_spec%transfer_impedance(1:cable_spec%n_transfer_impedance_models) )
    do i=1,cable_spec%n_transfer_impedance_models
      read(cable_spec_file_unit,*,IOSTAT=ierr)    ! comment line
      CALL read_Sfilter(cable_spec%transfer_impedance(i),cable_spec_file_unit)

      
! test pole stability
      CALL test_filter_pole_stability(cable_spec%transfer_impedance(i),check_passed)
  
      if (.NOT.check_passed) then
        write(*,*)'ERROR in cable_rmpedance filter function has unstable poles'
        write(*,*)'Transfer impedance filter function number:',i
        CALL write_Sfilter(cable_spec%transfer_impedance(i),0)
        write(run_status,'(A,I2)')'ERROR, Unstable poles in transfer impedance filter function: ',i
        CALL write_program_status()
        STOP 1
      else
        if (verbose) write(*,*)'The transfer impedance filter function has stable poles: function number ',i    
      end if

    end do
  end if
   
! Set deafult propagation correction transfer function fit information
  cable_spec%Y_fit_model_order=0  
  CALL reset_frequency_specification(cable_spec%Y_fit_freq_spec)
  CALL set_up_frequency_specification(cable_spec%Y_fit_freq_spec)
  
! Read the optional propagation correction transfer function fit information
  read(cable_spec_file_unit,*,IOSTAT=ierr)cable_spec%Y_fit_model_order
  if (ierr.NE.0) then
! Assume there is no filter fit information specified so move on to the next stage
    backspace(cable_spec_file_unit)
    goto 100
  end if
  
  write(*,*)'Reading the filter fit frequency range'
  
  CALL read_and_set_up_frequency_specification(cable_spec%Y_fit_freq_spec,cable_spec_file_unit)

100 continue
  
! the end of the file can contain flags to control the running of the software and the output

  do

    read(cable_spec_file_unit,'(A)',END=110,ERR=110)line

    CALL convert_to_lower_case(line,line_length)

! Set flags according to the information at the end of the .spice_model_spec file
  
    if (INDEX(line,'verbose').NE.0) verbose=.TRUE.
    if (INDEX(line,'use_s_xfer').NE.0) use_s_xfer=.TRUE.
    if (INDEX(line,'no_s_xfer').NE.0) use_s_xfer=.FALSE.
    if (INDEX(line,'use_laplace').NE.0) use_Laplace=.TRUE.
    if (INDEX(line,'no_laplace').NE.0) use_Laplace=.FALSE.

    if (INDEX(line,'plot_potential').NE.0) plot_potential=.TRUE.
    if (INDEX(line,'no_plot_potential').NE.0) plot_potential=.FALSE.
    if (INDEX(line,'plot_mesh').NE.0) plot_mesh=.TRUE.
    if (INDEX(line,'no_plot_mesh').NE.0) plot_mesh=.FALSE.

    
    if (INDEX(line,'direct_solver').NE.0) direct_solver=.TRUE.
    if (INDEX(line,'iterative_solver').NE.0) direct_solver=.FALSE.

! redefine mesh generation parameters if required
    if (INDEX(line,'laplace_boundary_constant').NE.0) then

      read(cable_spec_file_unit,*,END=9000,ERR=9000)Laplace_boundary_constant

    end if
    
    if (INDEX(line,'laplace_surface_mesh_constant').NE.0) then

      read(cable_spec_file_unit,*,END=9000,ERR=9000)Laplace_surface_mesh_constant

    end if
    
    if (INDEX(line,'twisted_pair_equivalent_radius').NE.0) then

      read(cable_spec_file_unit,*,END=9000,ERR=9000)Twisted_pair_equivalent_radius

    end if

    
    if (INDEX(line,'max_mesh_edge_length').NE.0) then

      read(cable_spec_file_unit,*,END=9000,ERR=9000)max_mesh_edge_length
    end if
    
    if (INDEX(line,'cg_tol').NE.0) then
      read(cable_spec_file_unit,*,END=9000,ERR=9000)cg_tol

    end if

    
    if (INDEX(line,'no_fasthenry').NE.0) use_FastHenry=.FALSE.
    
    if (INDEX(line,'use_fasthenry').EQ.1) then
      use_FastHenry=.TRUE.
      
      stripped_line=line(14:line_length)
      read(stripped_line,*,ERR=305,END=305)FH2_nlayers_radius,FH2_nw,FH2_nh,FH2_rw,FH2_rh
      
      if(INDEX(line,'auto_cyl_grid').NE.0) then
        auto_cyl_grid=.TRUE.
      end if
      if(INDEX(line,'no_refinement').NE.0) then
        fh2_no_refinement=.TRUE.
      end if
      
305   CONTINUE
    
      write(*,*)'use_FastHenry=.TRUE.'
      write(*,*)'Parameters:',FH2_nlayers_radius,FH2_nw,FH2_nh,FH2_rw,FH2_rh
      write(*,*)'auto_cyl_grid    :',auto_cyl_grid
      write(*,*)'fh2_no_refinement:',fh2_no_refinement
      
    end if
    
    if (INDEX(line,'l_from_fasthenry').NE.0) L_from_fasthenry=.TRUE.

  end do  ! continue until all flags are read - indicated by an end of file.

110 CONTINUE
  
! close the file with the cable data
  CLOSE(unit=cable_spec_file_unit)
  

  if (use_FastHenry) use_Laplace=.TRUE.  ! We must use the Laplace solver too before we run FastHenry2 so set this flag
  

  if (cable_spec%n_internal_domains.GT.0) then
! construct the information for the internal domains

! allocate internal domain information
    ALLOCATE(cable_spec%n_internal_conductors_in_domain(1:cable_spec%n_internal_domains))
    ALLOCATE(cable_spec%L_domain(1:cable_spec%n_internal_domains))
    ALLOCATE(cable_spec%C_domain(1:cable_spec%n_internal_domains))
    ALLOCATE(cable_spec%Z_domain(1:cable_spec%n_internal_domains))
    ALLOCATE(cable_spec%Y_domain(1:cable_spec%n_internal_domains))

  end if ! There are internal domains requiring parameter specification

! allocate the conductor based impedance model and reset all parameters assuming PEC for now
  ALLOCATE( cable_spec%conductor_impedance(1:cable_spec%tot_n_conductors) )
  do i=1,cable_spec%tot_n_conductors
    cable_spec%conductor_impedance(i)%impedance_model_type=impedance_model_type_PEC
    cable_spec%conductor_impedance(i)%radius=0d0
    cable_spec%conductor_impedance(i)%conductivity=0d0
    cable_spec%conductor_impedance(i)%Resistance_multiplication_factor=1d0
  end do
      
! Set the cable structure information according to each cable type
  
  if (cable_spec%cable_type.EQ.cable_geometry_type_cylindrical) then

    CALL cylindrical_set_internal_domain_information(cable_spec)
  
  else if (cable_spec%cable_type.EQ.cable_geometry_type_twisted_pair) then

    CALL twisted_pair_set_internal_domain_information(cable_spec)

  else if (cable_spec%cable_type.EQ.cable_geometry_type_coax) then

    CALL coax_set_internal_domain_information(cable_spec)

  else if (cable_spec%cable_type.EQ.cable_geometry_type_twinax) then

    CALL twinax_set_internal_domain_information(cable_spec)

  else if (cable_spec%cable_type.EQ.cable_geometry_type_shielded_twisted_pair) then

    CALL shielded_twisted_pair_set_internal_domain_information(cable_spec)
    
  else if (cable_spec%cable_type.EQ.cable_geometry_type_spacewire) then

    CALL spacewire_set_internal_domain_information(cable_spec)

  else if (cable_spec%cable_type.EQ.cable_geometry_type_flex_cable) then

    CALL flex_cable_set_internal_domain_information(cable_spec)

  else if (cable_spec%cable_type.EQ.cable_geometry_type_ML_flex_cable) then

    CALL ML_flex_cable_set_internal_domain_information(cable_spec)

  else if (cable_spec%cable_type.EQ.cable_geometry_type_dconnector) then

    CALL dconnector_set_internal_domain_information(cable_spec)
  
  else if (cable_spec%cable_type.EQ.cable_geometry_type_overshield) then

    CALL overshield_set_internal_domain_information(cable_spec)
  
  end if

! Write the cable information to a .cable file

  CAll write_cable(cable_spec,cable_file_unit)

! deallocate memory used
  
  CALL deallocate_frequency_specification(cable_spec%Y_fit_freq_spec)

  CAll deallocate_cable(cable_spec)

! finish up

  run_status='Finished_Correctly'
  CALL write_program_status()
  

  STOP
 
9000    run_status='ERROR reading control parameter from the cable_spec file'
        CALL write_program_status()
        STOP 1
  

END PROGRAM cable_model_builder