!
! 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 .
!
! 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 .
!
! 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