!
! 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_bundle_model_builder
!
! NAME
!     cable_bundle_model_builder
!
! AUTHORS
!     Chris Smartt
!
! DESCRIPTION
!     The cable_bundle_model_builder brings together a number of cable models, ground plane and 
!     overshield models with their position in the bundle cross section as required to build a model of a cable bundle. 
!     The output of the bundle_model_builder code is a file name.bundle which can be used to generate Spice cable bundle models
!     for a particular modellling scenario using the program spice_cable_bundle_model_builder
!
!     The input to the program is the name of a cablebundle specification file. A file name.bundle_spec must exist, containing
!     all the data required to specify a cable bundle.
!
!     The .cable files referred to in the .bundle_spec file are looked for in a directory specified in the .bundle_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 output of the cable_bundle_model_builder code is a file name.bundle which can be used as an input to the 
!     spice_cable_bundle_model_builder software
!
!     The .bundle file is placed in a directory specified in the .bundle_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 bundle name specified in the command line i.e. 'cable_bundle_model_builder name'
!     or it the name is absent from the command line, the user is prompted for the name.

!
! COMMENTS
!      Update to V2
!
!
! HISTORY
!
!     started 17/11/2015 CJS: STAGE_1 developments
!     started 3/2/2016 CJS:   STAGE_2 developments - multi-conductor bundles
!     started 22/3/2016 CJS:  STAGE_3 developments - multi-domain bundles (shielded cables)
!     started 13/4/2016 CJS:  Simplify the code and improve the notation. At this stage the format of the .bundle_spec file changed
!                             Shift the complex bundle building processes into BUNDLE_DOMAIN_CREATION/create_global_domain_structure.F90
!     20/4/2016 CJS: Allow the code to read x,y  OR x, y, theta formats for cable position and orientation. Similar for ground plane
!     27/4/2016 CJS: Include a conductor based impedance (loss) model
!     29/6/2016 CJS: Check that the cables are all on one side of the ground plane if it is present
!     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.
!     13/3/2018 CJS Add flag for direct/ iterative matrix solver in Laplace solution and inf/finite ground plane
!     19/6/2018 CJS Add flag for Neumann/ Asymptotic boundary condition in Laplace solver. Default is Neumann
!     21/12/2023: Trap the special case of a tristed pair in free space which needs a special process...
!
PROGRAM cable_bundle_model_builder

USE type_specifications
USE general_module
USE constants
USE cable_module
USE cable_bundle_module
USE PUL_parameter_module
USE filter_module
USE maths
USE bundle_domain_creation

IMPLICIT NONE

! local variables

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

character(len=filename_length)    :: bundle_name_with_path  ! name of the bundle including the path
character(len=filename_length)    :: bundle_path            ! path to the bundle_spec file
character(len=filename_length)    :: bundle_name            ! name of the bundle
character(len=filename_length)    :: filename               ! filename of the .bundle_spec file

! structure for the bundle specification (see CABLE_BUNDLE_MODULES/cable_bundle_module.F90)
type(bundle_specification_type)    ::bundle_spec
    
logical    :: file_exists

character(len=line_length)    :: line

integer        :: ierr,ierr2   ! integers to return error codes from file reads

! loop variables
integer        :: cable
integer        :: domain
integer        :: dim
    
logical    :: must_use_laplace

! START

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

! Open the input file describing the cable bundle parameters
! This file could be created by the associated GUI or otherwise generated
  CALL get_command_argument(1 , argument1, argument1_length)

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

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

    read(*,'(A)')bundle_name_with_path

  end if
    
  CALL strip_path(bundle_name_with_path,bundle_path,bundle_name)

  filename=trim(bundle_name_with_path)//bundle_spec_file_extn

  inquire(file=trim(filename),exist=file_exists)
  if (.NOT.file_exists) then
    run_status='ERROR Cannot find the file:'//trim(filename)
    CALL write_program_status()
    STOP 1
  end if 

! set the version tag in the cable_bundle structure
  bundle_spec%version=SPICE_CABLE_MODEL_BUILDER_version
  
! open and read the .bundle_spec file
  
  OPEN(unit=bundle_spec_file_unit,file=filename)

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

! read .bundle_spec file. 

! read the MOD directory information for the cable and bundle models
  read(bundle_spec_file_unit,*)  ! comment line
  read(bundle_spec_file_unit,'(A)')MOD_cable_lib_dir
  CALL path_format(MOD_cable_lib_dir)
  read(bundle_spec_file_unit,*)  ! comment line
  read(bundle_spec_file_unit,'(A)')MOD_bundle_lib_dir
  CALL path_format(MOD_bundle_lib_dir)

! ensure that the paths exist
  if (.NOT.path_exists(MOD_cable_lib_dir)) then
    run_status='ERROR MOD_cable_lib_dir does not exist '//trim(MOD_cable_lib_dir)
    CALL write_program_status()
    STOP 1
  end if
  CALL check_and_make_path(MOD_bundle_lib_dir)

  read(bundle_spec_file_unit,*,IOSTAT=ierr)bundle_spec%n_cables
  if (ierr.NE.0) then 
    run_status='ERROR reading bundle_spec%n_cables '
    CALL write_program_status()
    STOP 1
  end if

! Allocate the cable based data structures including an additional cable for the ground plane if required
  ALLOCATE( bundle_spec%cable(1:bundle_spec%n_cables+1) )
  ALLOCATE( bundle_spec%cable_x_offset(1:bundle_spec%n_cables+1) )
  ALLOCATE( bundle_spec%cable_y_offset(1:bundle_spec%n_cables+1) )
  ALLOCATE( bundle_spec%cable_angle(1:bundle_spec%n_cables+1) )
  
! read the name (which indicates the name.cable file specifying the cable) 
! and position in the x-y bundle cross section of each cable

  do cable=1,bundle_spec%n_cables

    read(bundle_spec_file_unit,'(A)',IOSTAT=ierr)bundle_spec%cable(cable)%cable_name
    if (ierr.NE.0) then 
      run_status='ERROR reading bundle_spec%cable(cable)%cable_name '
      CALL write_program_status()
      STOP 1
    end if 
    
    read(bundle_spec_file_unit,'(A)',IOSTAT=ierr)line   
    read(line,*,IOSTAT=ierr)bundle_spec%cable_x_offset(cable),    &
                            bundle_spec%cable_y_offset(cable),    &
                            bundle_spec%cable_angle(cable)
    if (ierr.NE.0) then 
    
! Try reading the old format, set angle to zero and read x and y only
      bundle_spec%cable_angle(cable)=0d0
      read(line,*,IOSTAT=ierr2)bundle_spec%cable_x_offset(cable),    &
                               bundle_spec%cable_y_offset(cable)      
      if (ierr2.NE.0) then 
        run_status='ERROR reading bundle_spec%cable_x_offset(cable),bundle_spec%cable_y_offset(cable)'
        CALL write_program_status()
        STOP 1
      end if
    end if 
    
    bundle_spec%cable_angle(cable)=bundle_spec%cable_angle(cable)*pi/180d0   ! convert angle to radians
    
  end do  ! read the informaton for the next cable 
  
! read ground plane specification (if it exists, if it is absent then assume we have no ground plane)

  read(bundle_spec_file_unit,'(A)',IOSTAT=ierr)line
  if (ierr.NE.0) then 
    run_status='ERROR reading ground plane present/ absent information'
    CALL write_program_status()
    STOP 1
  end if 
  CALL convert_to_lower_case(line,line_length)
  
  bundle_spec%n_cables_without_ground_plane=bundle_spec%n_cables
  
  if (line.eq.'ground_plane') then
  
    bundle_spec%ground_plane_present=.TRUE.
    
! include an additional cable in the cable list for the ground plane (space for this is already allocated)

    bundle_spec%n_cables=bundle_spec%n_cables+1
    cable=bundle_spec%n_cables
    
! We now assume that the ground plane is along the x axis in the bundle cross section

! set x, y, theta, angle and offset to 0
    bundle_spec%ground_plane_x=0d0
    bundle_spec%ground_plane_y=0d0
    bundle_spec%ground_plane_theta=0d0
    
    bundle_spec%ground_plane_angle=pi/2d0    
    bundle_spec%ground_plane_nx=cos(bundle_spec%ground_plane_angle)
    bundle_spec%ground_plane_ny=sin(bundle_spec%ground_plane_angle)          
    bundle_spec%ground_plane_offset=bundle_spec%ground_plane_nx*bundle_spec%ground_plane_x+ &
                                    bundle_spec%ground_plane_ny*bundle_spec%ground_plane_y 
    
! copy to cable structure
    bundle_spec%cable_x_offset(cable)=bundle_spec%ground_plane_x
    bundle_spec%cable_y_offset(cable)=bundle_spec%ground_plane_y
    bundle_spec%cable_angle(cable)=bundle_spec%ground_plane_theta

! create conductor impedance model for the new cable    
    ALLOCATE(bundle_spec%cable(cable)%conductor_impedance(1:1))
    bundle_spec%cable(cable)%conductor_impedance(1)%impedance_model_type=impedance_model_type_PEC
    bundle_spec%cable(cable)%conductor_impedance(1)%Resistance_multiplication_factor=1d0
    bundle_spec%cable(cable)%conductor_impedance(1)%radius=0d0
    bundle_spec%cable(cable)%conductor_impedance(1)%conductivity=0d0
        
  else if (line.eq.'no_ground_plane') then
  
    bundle_spec%ground_plane_present=.FALSE. ! no ground plane present
  
! set x, y, theta, angle and offset to 0
    bundle_spec%ground_plane_x=0d0
    bundle_spec%ground_plane_y=0d0
    bundle_spec%ground_plane_theta=0d0
    bundle_spec%ground_plane_angle=0d0
    bundle_spec%ground_plane_offset=0d0
    
  else   ! no ground plane information was given
  
    bundle_spec%ground_plane_present=.FALSE. ! no ground plane present
    
! go back one line of the input file before continuing to read the .bundle_spec file

    backspace(unit=bundle_spec_file_unit)
  
  end if ! ground plane
     
! Set deafult propagation correction transfer function fit information
  bundle_spec%Y_fit_model_order=0  
  CALL reset_frequency_specification(bundle_spec%Y_fit_freq_spec)
  CALL set_up_frequency_specification(bundle_spec%Y_fit_freq_spec)
  
! Read the optional propagation correction transfer function fit information
  read(bundle_spec_file_unit,*,IOSTAT=ierr)bundle_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(bundle_spec_file_unit)
    goto 100
  end if
  
  write(*,*)'Reading the filter fit frequency range'
  
  CALL read_and_set_up_frequency_specification(bundle_spec%Y_fit_freq_spec,bundle_spec_file_unit)

100 continue
  
! the file can contain flags to control the running of the software and the output
  rewind(bundle_spec_file_unit)

  write(*,*)'Processing flags'

  do
    read(bundle_spec_file_unit,*,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.
    
    if (INDEX(line,'inf_gnd').NE.0) inf_gnd=.TRUE.
    if (INDEX(line,'finite_gnd').NE.0) inf_gnd=.FALSE.
    
    if (INDEX(line,'abc').NE.0) use_ABC=.TRUE.
    if (INDEX(line,'neumann').NE.0) use_ABC=.FALSE.

! redefine mesh generation parameters if required
    if (INDEX(line,'laplace_boundary_constant').NE.0) then
      read(bundle_spec_file_unit,*,END=9000,ERR=9000)Laplace_boundary_constant
    end if
    
    if (INDEX(line,'laplace_surface_mesh_constant').NE.0) then
      read(bundle_spec_file_unit,*,END=9000,ERR=9000)Laplace_surface_mesh_constant
    end if
    
    if (INDEX(line,'twisted_pair_equivalent_radius').NE.0) then
      read(bundle_spec_file_unit,*,END=9000,ERR=9000)Twisted_pair_equivalent_radius
    end if
    
    if (INDEX(line,'max_mesh_edge_length').NE.0) then
      read(bundle_spec_file_unit,*,END=9000,ERR=9000)max_mesh_edge_length
    end if
    
    if (INDEX(line,'gp_edge_length').NE.0) then
      read(bundle_spec_file_unit,*,END=9000,ERR=9000)gp_edge_length
    end if
    
    if (INDEX(line,'cg_tol').NE.0) then
      read(bundle_spec_file_unit,*,END=9000,ERR=9000)cg_tol
    end if

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

110 CONTINUE

! close the bundle_spec file

  CLOSE(unit=bundle_spec_file_unit)

  write(*,*)'Closed file:',trim(filename)
  
  write(*,*)'Total number of cables (including ground plane):',bundle_spec%n_cables
  write(*,*)'Number of cables (not including ground plane)  :',bundle_spec%n_cables_without_ground_plane

! read the .cable files required to build the bundle
! Also check which cable types need to use the Laplace solver in the external domain
  must_use_laplace=.FALSE.
  do cable=1,bundle_spec%n_cables_without_ground_plane
  
    CALL read_cable( bundle_spec%cable(cable),cable_file_unit)
    if (    (bundle_spec%cable(cable)%cable_type.EQ.cable_geometry_type_flex_cable)    &
        .OR.(bundle_spec%cable(cable)%cable_type.EQ.cable_geometry_type_ML_flex_cable) ) then
      must_use_laplace=.TRUE.
    end if
    
    if ( (bundle_spec%cable(cable)%cable_type.EQ.cable_geometry_type_dconnector)   &
         .AND.(bundle_spec%n_cables.GT.1) ) then
      must_use_laplace=.TRUE.
    end if
    
  end do ! read the next cable file in the bundle

! 21/12/2023: Check for the special case of a tristed pair in free space which needs
! a special process...
  
  if (bundle_spec%n_cables.Eq.1) then
    cable=1
    if (bundle_spec%cable(cable)%cable_type.EQ.cable_geometry_type_twisted_pair) then
    run_status='ERROR: twisted pair cable in free space with no other conductors. Model this using two cylindrical wires'
    CALL write_program_status()
    STOP 1
    end if 
  end if


! Check whether Laplace solver must be used

  if (must_use_laplace.AND.(.NOT.use_Laplace)) then
    run_status='ERROR: The laplace solver must be used for this cable bundle.'
    CALL write_program_status()
    STOP 1
  end if
  
! set the ground plane cable data (if it is present)

  if (bundle_spec%ground_plane_present) then
  
    cable=bundle_spec%n_cables
    bundle_spec%cable(cable)%version=SPICE_CABLE_MODEL_BUILDER_version
    bundle_spec%cable(cable)%cable_name='ground plane'
    bundle_spec%cable(cable)%cable_type_string='ground_plane' 
    bundle_spec%cable(cable)%cable_type=cable_geometry_type_ground_plane
    
    CALL ground_plane_set_parameters(bundle_spec%cable(cable))
    CALL ground_plane_set_internal_domain_information(bundle_spec%cable(cable))
    
! check that all the conductors are on one side of the ground plane
    CALL check_cables_wrt_ground_plane(bundle_spec)

  end if
  
! Write the bundle cross section (condcutors and dielectrics) to files suitable for gnuplot
  CALL plot_bundle_cross_section(bundle_spec)

! check whether there is any intersection between cables  
  CALL check_cable_intersection(bundle_spec)

! now we have the complete bundle specification information we can go and build the bundle structure

  CALL create_global_domain_structure(bundle_spec)

! Create the arrays of x and y coordinates of each conductor which are used in the incident field excitation
! moved from before create_global_domain_structure

  CALL set_conductor_positions_for_Einc(bundle_spec)
  
  if (verbose) then
! Write the .bundle structure to screen

    write(*,*)'____________________________________________'
    write(*,*)''
    write(*,*)'Domain inductance matrices [L]'
    write(*,*)''
  
    do domain=1,bundle_spec%tot_n_domains
      write(*,*)'Domain=',domain,' matrix dimension=',bundle_spec%L(domain)%dim
      dim=bundle_spec%L(domain)%dim
      CALL dwrite_matrix(bundle_spec%L(domain)%mat,dim,dim,dim,0)
     end do

    write(*,*)'____________________________________________'
    write(*,*)''
    write(*,*)'Domain capacitance matrices [C]'
    write(*,*)''
  
    do domain=1,bundle_spec%tot_n_domains
      write(*,*)'Domain=',domain,' matrix dimension=',bundle_spec%C(domain)%dim
      dim=bundle_spec%C(domain)%dim
      CALL dwrite_matrix(bundle_spec%C(domain)%mat,dim,dim,dim,0)
    end do

    write(*,*)'____________________________________________'
    write(*,*)''
    write(*,*)'Global external to domain conductor current transformation matrix, [MI]'
    write(*,*)''
  
    dim=bundle_spec%global_MI%dim
    CALL dwrite_matrix(bundle_spec%global_MI%mat,dim,dim,dim,0)

    write(*,*)'____________________________________________'
    write(*,*)''
    write(*,*)'Global external to domain conductor voltage transformation matrix, [MV]'
    write(*,*)''
  
    dim=bundle_spec%global_MV%dim
    CALL dwrite_matrix(bundle_spec%global_MV%mat,dim,dim,dim,0)

    write(*,*)'____________________________________________'
    write(*,*)''
    write(*,*)'Global Inductance matrix, [L]'
    write(*,*)''
  
    dim=bundle_spec%global_L%dim
    CALL dwrite_matrix(bundle_spec%global_L%mat,dim,dim,dim,0)

    write(*,*)'____________________________________________'
    write(*,*)''
    write(*,*)'Global Capacitance matrix, [C]'
    write(*,*)''
  
    dim=bundle_spec%global_C%dim
    CALL dwrite_matrix(bundle_spec%global_C%mat,dim,dim,dim,0)
  
  end if ! verbose  
  
  CALL write_cable_bundle(bundle_spec,bundle_file_unit)
  
  CALL deallocate_frequency_specification(bundle_spec%Y_fit_freq_spec)

  CALL deallocate_cable_bundle(bundle_spec)

! finish up

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

END PROGRAM cable_bundle_model_builder