!
! 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:
! PROGRAM rational_function_fit
!
! NAME
!     rational_function_fit
!
! AUTHORS
!     Chris Smartt
!
! DESCRIPTION
!     Read in complex frequency domain data from a file and calculate a 'best fit' rational function
!     to represent the data. The numerator and denominator orders of the rational function should be specified explicitly.
!    
!     The rational function is returned in a file.
!
!     This process can be used to generate models for frequency dependent dielectric relative permittivity functions
!     or for frequency dependent transfer impedance functions which can subsequently be used to develop cable bundle
!     models for Spice.
!
!
! COMMENTS
!     The rational function fitting process is done with no constraints applied at the moment.
!
! HISTORY
!
!     started 22/08/2016 CJS 
!
!
PROGRAM rational_function_fit

USE type_specifications
USE general_module
USE constants
USE frequency_spec
USE filter_module
USE Sfilter_fit_module

IMPLICIT NONE

! local variables

character(len=filename_length)    :: ipfilename     ! name of the input datafile
character(len=filename_length)    :: opfilename     ! filename for the output rational function
    
logical        :: file_exists

integer                  :: n_frequency_samples ! number of frequencies for which data is available
real(dp),allocatable     :: f(:)                ! Frequency data
complex(dp),allocatable  :: function(:)         ! Frequency dependent complex function data

integer        :: aorder                        ! numerator order of rational function model
integer        :: border                        ! denominator order of rational function model
integer        :: relative_border               ! denominator order relative to numerator order

type(Sfilter)  :: function_model                ! Frequency dependent rational function model of the input data

logical        :: stable                        ! flag to indicate whether this model has stable poles (LHS of s-plane) or not

character(len=line_length) :: line              ! line read from the input file

real(dp)       :: f_in                          ! frequency read from file
real(dp)       :: re                            ! real part of function
real(dp)       :: im                            ! imaginary part of function

complex(dp)    :: function_fit                  ! function value from the rational function model

integer        :: loop                          ! data reading loop variable
integer        :: frequency_sample              ! frequency loop variable
 
integer        :: ierr  ! integer to return error codes from file reads

! START

! Open the input file containing the complex frequency domain data

  verbose=.TRUE.

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

  write(*,*)'Enter the name of the complex frequency domain data'

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

  write(*,*)'Opened file:',trim(ipfilename)
  
  do loop=1,2     ! read the data file twice, the first time to work out the amount of data, the second to read the data into arrays

    frequency_sample=0       ! reset the counter for the frequency samples

100 CONTINUE  
      read(local_file_unit,'(A)',END=200)line   ! read the next line of the file into the string, line; exit at the end of the file
    
! try to read data in the expected format from the line just read, if the line does not have data in the correct format then skip it
      read(line,*,END=110,ERR=110)f_in,re,im
    
      frequency_sample=frequency_sample+1
    
      if (loop.EQ.2) then                            ! put the values in the frequency and function value arrays
        f(frequency_sample)=f_in
        function(frequency_sample)=cmplx(re,im)
      end if
  
110   CONTINUE    ! jump here if the line read doesn't have data in the correct format
      
      GOTO 100

200 CONTINUE    ! jump here when we reach the end of the file

    if (loop.EQ.1) then                             ! allocate space for all of the data points
    
      n_frequency_samples=frequency_sample
      ALLOCATE( f(1:n_frequency_samples) )
      ALLOCATE( function(1:n_frequency_samples) )
      rewind(local_file_unit)                       ! go back to the beginning of the file to read it again
      
    end if
  
  end do ! next read loop

! close the file with the input data
  CLOSE(unit=local_file_unit)
  
  write(*,*)'Number of frequency samples read:',n_frequency_samples
  write(*,*)
  
! Create a rational function model of the function data 
  
  write(*,*)'Enter the numerator order for the model, this should be no more than the denominator order+1'
  read(*,*)aorder
  
  write(*,*)'Enter the denominator order for the model'
  read(*,*)border
  
  if (aorder.GT.border+1) then
    run_status='ERROR: The  numerator order should be no more than the denominator order+1'
    CALL write_program_status()
    STOP 1
  end if
  
  relative_border=border-aorder
  
! call calculate_Sfilter for the specified model order

  write(*,*)'Calculating the rational function model...'
  CALL Calculate_Sfilter(function,f,n_frequency_samples,function_model,aorder,relative_border,0) ! call with fit_type=0 i.e. no constraints

! check the stability of the model i.e. the poles are on the LHS of the s-plane  
  CALL test_filter_pole_stability(function_model,stable)
  
  if (stable) then
    write(*,*)'This model is stable'
  else
    write(*,*)'********** WARNING ************'
    write(*,*)'*** This model is unstable ****'  
    write(*,*)'*******************************'
  end if
  
! Write vector fit models to file
  open(unit=local_file_unit,file='function_fit.fout')
  do frequency_sample=1,n_frequency_samples
  
    function_fit=evaluate_Sfilter_frequency_response(function_model,f(frequency_sample))
    
    write(local_file_unit,8000)f(frequency_sample),real(function(frequency_sample)),aimag(function(frequency_sample)),   &
                                                   real(function_fit),aimag(function_fit)
8000 format(5ES16.6)

  end do
  
  close(unit=local_file_unit)

! Open a file for the best fit model

!  write(*,*)'Enter the name of the best fit model'
!  read(*,'(A)')opfilename

  opfilename=trim(ipfilename)//'.function_fit'
  
  open(unit=local_file_unit,file=opfilename)

! Write the rational function model to the file

  write(local_file_unit,*)'# Rational function model fit to data from file:',trim(ipfilename)
  CALL Write_Sfilter(function_model,local_file_unit)

  close(unit=local_file_unit)

! deallocate memory and finish up
  CALL deallocate_Sfilter(function_model)
  
  DEALLOCATE(f)
  DEALLOCATE(function)

  run_status='Finished_Correctly'
  CALL write_program_status()
  
END PROGRAM rational_function_fit