diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..0ea3d3f92b974e57fbe8dc1a2f45355c690c9152
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,18 @@
+.pydevproject
+.project
+.cproject
+.settings
+obj
+bin
+core*
+*~
+*.pyc
+*.so
+*.so*
+.pylintrc
+.metadata
+.idea
+.cvsignore
+.nse_depinfo
+software
+oldsrc
diff --git a/Doxyfile b/Doxyfile
new file mode 100644
index 0000000000000000000000000000000000000000..e3b078396cf875f937a3807401b4cc4ca222a6d0
--- /dev/null
+++ b/Doxyfile
@@ -0,0 +1,1344 @@
+# Doxyfile 1.5.3
+
+# This file describes the settings to be used by the documentation system
+# doxygen (www.doxygen.org) for a project
+#
+# All text after a hash (#) is considered a comment and will be ignored
+# The format is:
+# TAG = value [value, ...]
+# For lists items can also be appended using:
+# TAG += value [value, ...]
+# Values that contain spaces should be placed between quotes (" ")
+
+#---------------------------------------------------------------------------
+# Project related configuration options
+#---------------------------------------------------------------------------
+
+# This tag specifies the encoding used for all characters in the config file that
+# follow. The default is UTF-8 which is also the encoding used for all text before
+# the first occurrence of this tag. Doxygen uses libiconv (or the iconv built into
+# libc) for the transcoding. See http://www.gnu.org/software/libiconv for the list of
+# possible encodings.
+
+DOXYFILE_ENCODING = UTF-8
+
+# The PROJECT_NAME tag is a single word (or a sequence of words surrounded
+# by quotes) that should identify the project.
+
+PROJECT_NAME = Interpolator
+
+# The PROJECT_NUMBER tag can be used to enter a project or revision number.
+# This could be handy for archiving the generated documentation or
+# if some version control system is used.
+
+PROJECT_NUMBER = 1
+
+# The OUTPUT_DIRECTORY tag is used to specify the (relative or absolute)
+# base path where the generated documentation will be put.
+# If a relative path is entered, it will be relative to the location
+# where doxygen was started. If left blank the current directory will be used.
+
+OUTPUT_DIRECTORY = doc
+
+# If the CREATE_SUBDIRS tag is set to YES, then doxygen will create
+# 4096 sub-directories (in 2 levels) under the output directory of each output
+# format and will distribute the generated files over these directories.
+# Enabling this option can be useful when feeding doxygen a huge amount of
+# source files, where putting all generated files in the same directory would
+# otherwise cause performance problems for the file system.
+
+CREATE_SUBDIRS = NO
+
+# The OUTPUT_LANGUAGE tag is used to specify the language in which all
+# documentation generated by doxygen is written. Doxygen will use this
+# information to generate all constant output in the proper language.
+# The default language is English, other supported languages are:
+# Afrikaans, Arabic, Brazilian, Catalan, Chinese, Chinese-Traditional,
+# Croatian, Czech, Danish, Dutch, Finnish, French, German, Greek, Hungarian,
+# Italian, Japanese, Japanese-en (Japanese with English messages), Korean,
+# Korean-en, Lithuanian, Norwegian, Polish, Portuguese, Romanian, Russian,
+# Serbian, Slovak, Slovene, Spanish, Swedish, and Ukrainian.
+
+OUTPUT_LANGUAGE = English
+
+# If the BRIEF_MEMBER_DESC tag is set to YES (the default) Doxygen will
+# include brief member descriptions after the members that are listed in
+# the file and class documentation (similar to JavaDoc).
+# Set to NO to disable this.
+
+BRIEF_MEMBER_DESC = YES
+
+# If the REPEAT_BRIEF tag is set to YES (the default) Doxygen will prepend
+# the brief description of a member or function before the detailed description.
+# Note: if both HIDE_UNDOC_MEMBERS and BRIEF_MEMBER_DESC are set to NO, the
+# brief descriptions will be completely suppressed.
+
+REPEAT_BRIEF = YES
+
+# This tag implements a quasi-intelligent brief description abbreviator
+# that is used to form the text in various listings. Each string
+# in this list, if found as the leading text of the brief description, will be
+# stripped from the text and the result after processing the whole list, is
+# used as the annotated text. Otherwise, the brief description is used as-is.
+# If left blank, the following values are used ("$name" is automatically
+# replaced with the name of the entity): "The $name class" "The $name widget"
+# "The $name file" "is" "provides" "specifies" "contains"
+# "represents" "a" "an" "the"
+
+ABBREVIATE_BRIEF = "The $name class " \
+ "The $name widget " \
+ "The $name file " \
+ is \
+ provides \
+ specifies \
+ contains \
+ represents \
+ a \
+ an \
+ the
+
+# If the ALWAYS_DETAILED_SEC and REPEAT_BRIEF tags are both set to YES then
+# Doxygen will generate a detailed section even if there is only a brief
+# description.
+
+ALWAYS_DETAILED_SEC = NO
+
+# If the INLINE_INHERITED_MEMB tag is set to YES, doxygen will show all
+# inherited members of a class in the documentation of that class as if those
+# members were ordinary class members. Constructors, destructors and assignment
+# operators of the base classes will not be shown.
+
+INLINE_INHERITED_MEMB = YES
+
+# If the FULL_PATH_NAMES tag is set to YES then Doxygen will prepend the full
+# path before files name in the file list and in the header files. If set
+# to NO the shortest path that makes the file name unique will be used.
+
+FULL_PATH_NAMES = NO
+
+# If the FULL_PATH_NAMES tag is set to YES then the STRIP_FROM_PATH tag
+# can be used to strip a user-defined part of the path. Stripping is
+# only done if one of the specified strings matches the left-hand part of
+# the path. The tag can be used to show relative paths in the file list.
+# If left blank the directory from which doxygen is run is used as the
+# path to strip.
+
+STRIP_FROM_PATH = /home/claudio/devel/utils/interpolator
+
+# The STRIP_FROM_INC_PATH tag can be used to strip a user-defined part of
+# the path mentioned in the documentation of a class, which tells
+# the reader which header file to include in order to use a class.
+# If left blank only the name of the header file containing the class
+# definition is used. Otherwise one should specify the include paths that
+# are normally passed to the compiler using the -I flag.
+
+STRIP_FROM_INC_PATH =
+
+# If the SHORT_NAMES tag is set to YES, doxygen will generate much shorter
+# (but less readable) file names. This can be useful is your file systems
+# doesn't support long names like on DOS, Mac, or CD-ROM.
+
+SHORT_NAMES = NO
+
+# If the JAVADOC_AUTOBRIEF tag is set to YES then Doxygen
+# will interpret the first line (until the first dot) of a JavaDoc-style
+# comment as the brief description. If set to NO, the JavaDoc
+# comments will behave just like regular Qt-style comments
+# (thus requiring an explicit @brief command for a brief description.)
+
+JAVADOC_AUTOBRIEF = NO
+
+# If the QT_AUTOBRIEF tag is set to YES then Doxygen will
+# interpret the first line (until the first dot) of a Qt-style
+# comment as the brief description. If set to NO, the comments
+# will behave just like regular Qt-style comments (thus requiring
+# an explicit \brief command for a brief description.)
+
+QT_AUTOBRIEF = NO
+
+# The MULTILINE_CPP_IS_BRIEF tag can be set to YES to make Doxygen
+# treat a multi-line C++ special comment block (i.e. a block of //! or ///
+# comments) as a brief description. This used to be the default behaviour.
+# The new default is to treat a multi-line C++ comment block as a detailed
+# description. Set this tag to YES if you prefer the old behaviour instead.
+
+MULTILINE_CPP_IS_BRIEF = NO
+
+# If the DETAILS_AT_TOP tag is set to YES then Doxygen
+# will output the detailed description near the top, like JavaDoc.
+# If set to NO, the detailed description appears after the member
+# documentation.
+
+DETAILS_AT_TOP = NO
+
+# If the INHERIT_DOCS tag is set to YES (the default) then an undocumented
+# member inherits the documentation from any documented member that it
+# re-implements.
+
+INHERIT_DOCS = YES
+
+# If the SEPARATE_MEMBER_PAGES tag is set to YES, then doxygen will produce
+# a new page for each member. If set to NO, the documentation of a member will
+# be part of the file/class/namespace that contains it.
+
+SEPARATE_MEMBER_PAGES = NO
+
+# The TAB_SIZE tag can be used to set the number of spaces in a tab.
+# Doxygen uses this value to replace tabs by spaces in code fragments.
+
+TAB_SIZE = 8
+
+# This tag can be used to specify a number of aliases that acts
+# as commands in the documentation. An alias has the form "name=value".
+# For example adding "sideeffect=\par Side Effects:\n" will allow you to
+# put the command \sideeffect (or @sideeffect) in the documentation, which
+# will result in a user-defined paragraph with heading "Side Effects:".
+# You can put \n's in the value part of an alias to insert newlines.
+
+ALIASES =
+
+# Set the OPTIMIZE_OUTPUT_FOR_C tag to YES if your project consists of C
+# sources only. Doxygen will then generate output that is more tailored for C.
+# For instance, some of the names that are used will be different. The list
+# of all members will be omitted, etc.
+
+OPTIMIZE_OUTPUT_FOR_C = NO
+
+# Set the OPTIMIZE_OUTPUT_JAVA tag to YES if your project consists of Java
+# sources only. Doxygen will then generate output that is more tailored for Java.
+# For instance, namespaces will be presented as packages, qualified scopes
+# will look different, etc.
+
+OPTIMIZE_OUTPUT_JAVA = NO
+
+# If you use STL classes (i.e. std::string, std::vector, etc.) but do not want to
+# include (a tag file for) the STL sources as input, then you should
+# set this tag to YES in order to let doxygen match functions declarations and
+# definitions whose arguments contain STL classes (e.g. func(std::string); v.s.
+# func(std::string) {}). This also make the inheritance and collaboration
+# diagrams that involve STL classes more complete and accurate.
+
+BUILTIN_STL_SUPPORT = YES
+
+# If you use Microsoft's C++/CLI language, you should set this option to YES to
+# enable parsing support.
+
+CPP_CLI_SUPPORT = NO
+
+# If member grouping is used in the documentation and the DISTRIBUTE_GROUP_DOC
+# tag is set to YES, then doxygen will reuse the documentation of the first
+# member in the group (if any) for the other members of the group. By default
+# all members of a group must be documented explicitly.
+
+DISTRIBUTE_GROUP_DOC = NO
+
+# Set the SUBGROUPING tag to YES (the default) to allow class member groups of
+# the same type (for instance a group of public functions) to be put as a
+# subgroup of that type (e.g. under the Public Functions section). Set it to
+# NO to prevent subgrouping. Alternatively, this can be done per class using
+# the \nosubgrouping command.
+
+SUBGROUPING = YES
+
+#---------------------------------------------------------------------------
+# Build related configuration options
+#---------------------------------------------------------------------------
+
+# If the EXTRACT_ALL tag is set to YES doxygen will assume all entities in
+# documentation are documented, even if no documentation was available.
+# Private class members and static file members will be hidden unless
+# the EXTRACT_PRIVATE and EXTRACT_STATIC tags are set to YES
+
+EXTRACT_ALL = NO
+
+# If the EXTRACT_PRIVATE tag is set to YES all private members of a class
+# will be included in the documentation.
+
+EXTRACT_PRIVATE = NO
+
+# If the EXTRACT_STATIC tag is set to YES all static members of a file
+# will be included in the documentation.
+
+EXTRACT_STATIC = NO
+
+# If the EXTRACT_LOCAL_CLASSES tag is set to YES classes (and structs)
+# defined locally in source files will be included in the documentation.
+# If set to NO only classes defined in header files are included.
+
+EXTRACT_LOCAL_CLASSES = YES
+
+# This flag is only useful for Objective-C code. When set to YES local
+# methods, which are defined in the implementation section but not in
+# the interface are included in the documentation.
+# If set to NO (the default) only methods in the interface are included.
+
+EXTRACT_LOCAL_METHODS = NO
+
+# If this flag is set to YES, the members of anonymous namespaces will be extracted
+# and appear in the documentation as a namespace called 'anonymous_namespace{file}',
+# where file will be replaced with the base name of the file that contains the anonymous
+# namespace. By default anonymous namespace are hidden.
+
+EXTRACT_ANON_NSPACES = NO
+
+# If the HIDE_UNDOC_MEMBERS tag is set to YES, Doxygen will hide all
+# undocumented members of documented classes, files or namespaces.
+# If set to NO (the default) these members will be included in the
+# various overviews, but no documentation section is generated.
+# This option has no effect if EXTRACT_ALL is enabled.
+
+HIDE_UNDOC_MEMBERS = YES
+
+# If the HIDE_UNDOC_CLASSES tag is set to YES, Doxygen will hide all
+# undocumented classes that are normally visible in the class hierarchy.
+# If set to NO (the default) these classes will be included in the various
+# overviews. This option has no effect if EXTRACT_ALL is enabled.
+
+HIDE_UNDOC_CLASSES = YES
+
+# If the HIDE_FRIEND_COMPOUNDS tag is set to YES, Doxygen will hide all
+# friend (class|struct|union) declarations.
+# If set to NO (the default) these declarations will be included in the
+# documentation.
+
+HIDE_FRIEND_COMPOUNDS = NO
+
+# If the HIDE_IN_BODY_DOCS tag is set to YES, Doxygen will hide any
+# documentation blocks found inside the body of a function.
+# If set to NO (the default) these blocks will be appended to the
+# function's detailed documentation block.
+
+HIDE_IN_BODY_DOCS = NO
+
+# The INTERNAL_DOCS tag determines if documentation
+# that is typed after a \internal command is included. If the tag is set
+# to NO (the default) then the documentation will be excluded.
+# Set it to YES to include the internal documentation.
+
+INTERNAL_DOCS = NO
+
+# If the CASE_SENSE_NAMES tag is set to NO then Doxygen will only generate
+# file names in lower-case letters. If set to YES upper-case letters are also
+# allowed. This is useful if you have classes or files whose names only differ
+# in case and if your file system supports case sensitive file names. Windows
+# and Mac users are advised to set this option to NO.
+
+CASE_SENSE_NAMES = YES
+
+# If the HIDE_SCOPE_NAMES tag is set to NO (the default) then Doxygen
+# will show members with their full class and namespace scopes in the
+# documentation. If set to YES the scope will be hidden.
+
+HIDE_SCOPE_NAMES = NO
+
+# If the SHOW_INCLUDE_FILES tag is set to YES (the default) then Doxygen
+# will put a list of the files that are included by a file in the documentation
+# of that file.
+
+SHOW_INCLUDE_FILES = YES
+
+# If the INLINE_INFO tag is set to YES (the default) then a tag [inline]
+# is inserted in the documentation for inline members.
+
+INLINE_INFO = YES
+
+# If the SORT_MEMBER_DOCS tag is set to YES (the default) then doxygen
+# will sort the (detailed) documentation of file and class members
+# alphabetically by member name. If set to NO the members will appear in
+# declaration order.
+
+SORT_MEMBER_DOCS = YES
+
+# If the SORT_BRIEF_DOCS tag is set to YES then doxygen will sort the
+# brief documentation of file, namespace and class members alphabetically
+# by member name. If set to NO (the default) the members will appear in
+# declaration order.
+
+SORT_BRIEF_DOCS = NO
+
+# If the SORT_BY_SCOPE_NAME tag is set to YES, the class list will be
+# sorted by fully-qualified names, including namespaces. If set to
+# NO (the default), the class list will be sorted only by class name,
+# not including the namespace part.
+# Note: This option is not very useful if HIDE_SCOPE_NAMES is set to YES.
+# Note: This option applies only to the class list, not to the
+# alphabetical list.
+
+SORT_BY_SCOPE_NAME = NO
+
+# The GENERATE_TODOLIST tag can be used to enable (YES) or
+# disable (NO) the todo list. This list is created by putting \todo
+# commands in the documentation.
+
+GENERATE_TODOLIST = YES
+
+# The GENERATE_TESTLIST tag can be used to enable (YES) or
+# disable (NO) the test list. This list is created by putting \test
+# commands in the documentation.
+
+GENERATE_TESTLIST = YES
+
+# The GENERATE_BUGLIST tag can be used to enable (YES) or
+# disable (NO) the bug list. This list is created by putting \bug
+# commands in the documentation.
+
+GENERATE_BUGLIST = YES
+
+# The GENERATE_DEPRECATEDLIST tag can be used to enable (YES) or
+# disable (NO) the deprecated list. This list is created by putting
+# \deprecated commands in the documentation.
+
+GENERATE_DEPRECATEDLIST= YES
+
+# The ENABLED_SECTIONS tag can be used to enable conditional
+# documentation sections, marked by \if sectionname ... \endif.
+
+ENABLED_SECTIONS =
+
+# The MAX_INITIALIZER_LINES tag determines the maximum number of lines
+# the initial value of a variable or define consists of for it to appear in
+# the documentation. If the initializer consists of more lines than specified
+# here it will be hidden. Use a value of 0 to hide initializers completely.
+# The appearance of the initializer of individual variables and defines in the
+# documentation can be controlled using \showinitializer or \hideinitializer
+# command in the documentation regardless of this setting.
+
+MAX_INITIALIZER_LINES = 30
+
+# Set the SHOW_USED_FILES tag to NO to disable the list of files generated
+# at the bottom of the documentation of classes and structs. If set to YES the
+# list will mention the files that were used to generate the documentation.
+
+SHOW_USED_FILES = YES
+
+# If the sources in your project are distributed over multiple directories
+# then setting the SHOW_DIRECTORIES tag to YES will show the directory hierarchy
+# in the documentation. The default is NO.
+
+SHOW_DIRECTORIES = NO
+
+# The FILE_VERSION_FILTER tag can be used to specify a program or script that
+# doxygen should invoke to get the current version for each file (typically from the
+# version control system). Doxygen will invoke the program by executing (via
+# popen()) the command <command> <input-file>, where <command> is the value of
+# the FILE_VERSION_FILTER tag, and <input-file> is the name of an input file
+# provided by doxygen. Whatever the program writes to standard output
+# is used as the file version. See the manual for examples.
+
+FILE_VERSION_FILTER =
+
+#---------------------------------------------------------------------------
+# configuration options related to warning and progress messages
+#---------------------------------------------------------------------------
+
+# The QUIET tag can be used to turn on/off the messages that are generated
+# by doxygen. Possible values are YES and NO. If left blank NO is used.
+
+QUIET = NO
+
+# The WARNINGS tag can be used to turn on/off the warning messages that are
+# generated by doxygen. Possible values are YES and NO. If left blank
+# NO is used.
+
+WARNINGS = YES
+
+# If WARN_IF_UNDOCUMENTED is set to YES, then doxygen will generate warnings
+# for undocumented members. If EXTRACT_ALL is set to YES then this flag will
+# automatically be disabled.
+
+WARN_IF_UNDOCUMENTED = YES
+
+# If WARN_IF_DOC_ERROR is set to YES, doxygen will generate warnings for
+# potential errors in the documentation, such as not documenting some
+# parameters in a documented function, or documenting parameters that
+# don't exist or using markup commands wrongly.
+
+WARN_IF_DOC_ERROR = YES
+
+# This WARN_NO_PARAMDOC option can be abled to get warnings for
+# functions that are documented, but have no documentation for their parameters
+# or return value. If set to NO (the default) doxygen will only warn about
+# wrong or incomplete parameter documentation, but not about the absence of
+# documentation.
+
+WARN_NO_PARAMDOC = NO
+
+# The WARN_FORMAT tag determines the format of the warning messages that
+# doxygen can produce. The string should contain the $file, $line, and $text
+# tags, which will be replaced by the file and line number from which the
+# warning originated and the warning text. Optionally the format may contain
+# $version, which will be replaced by the version of the file (if it could
+# be obtained via FILE_VERSION_FILTER)
+
+WARN_FORMAT = "$file:$line: $text "
+
+# The WARN_LOGFILE tag can be used to specify a file to which warning
+# and error messages should be written. If left blank the output is written
+# to stderr.
+
+WARN_LOGFILE =
+
+#---------------------------------------------------------------------------
+# configuration options related to the input files
+#---------------------------------------------------------------------------
+
+# The INPUT tag can be used to specify the files and/or directories that contain
+# documented source files. You may enter file names like "myfile.cpp" or
+# directories like "/usr/src/myproject". Separate the files or directories
+# with spaces.
+
+INPUT = .
+
+# This tag can be used to specify the character encoding of the source files that
+# doxygen parses. Internally doxygen uses the UTF-8 encoding, which is also the default
+# input encoding. Doxygen uses libiconv (or the iconv built into libc) for the transcoding.
+# See http://www.gnu.org/software/libiconv for the list of possible encodings.
+
+INPUT_ENCODING = UTF-8
+
+# If the value of the INPUT tag contains directories, you can use the
+# FILE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp
+# and *.h) to filter out the source-files in the directories. If left
+# blank the following patterns are tested:
+# *.c *.cc *.cxx *.cpp *.c++ *.java *.ii *.ixx *.ipp *.i++ *.inl *.h *.hh *.hxx
+# *.hpp *.h++ *.idl *.odl *.cs *.php *.php3 *.inc *.m *.mm *.py
+
+FILE_PATTERNS = *.c \
+ *.cc \
+ *.cxx \
+ *.cpp \
+ *.c++ \
+ *.d \
+ *.java \
+ *.ii \
+ *.ixx \
+ *.ipp \
+ *.i++ \
+ *.inl \
+ *.h \
+ *.hh \
+ *.hxx \
+ *.hpp \
+ *.h++ \
+ *.idl \
+ *.odl \
+ *.cs \
+ *.php \
+ *.php3 \
+ *.inc \
+ *.m \
+ *.mm \
+ *.dox \
+ *.py \
+ *.C \
+ *.CC \
+ *.C++ \
+ *.II \
+ *.I++ \
+ *.H \
+ *.HH \
+ *.H++ \
+ *.CS \
+ *.PHP \
+ *.PHP3 \
+ *.M \
+ *.MM \
+ *.PY
+
+# The RECURSIVE tag can be used to turn specify whether or not subdirectories
+# should be searched for input files as well. Possible values are YES and NO.
+# If left blank NO is used.
+
+RECURSIVE = NO
+
+# The EXCLUDE tag can be used to specify files and/or directories that should
+# excluded from the INPUT source files. This way you can easily exclude a
+# subdirectory from a directory tree whose root is specified with the INPUT tag.
+
+EXCLUDE =
+
+# The EXCLUDE_SYMLINKS tag can be used select whether or not files or
+# directories that are symbolic links (a Unix filesystem feature) are excluded
+# from the input.
+
+EXCLUDE_SYMLINKS = NO
+
+# If the value of the INPUT tag contains directories, you can use the
+# EXCLUDE_PATTERNS tag to specify one or more wildcard patterns to exclude
+# certain files from those directories. Note that the wildcards are matched
+# against the file with absolute path, so to exclude all test directories
+# for example use the pattern */test/*
+
+EXCLUDE_PATTERNS =
+
+# The EXCLUDE_SYMBOLS tag can be used to specify one or more symbol names
+# (namespaces, classes, functions, etc.) that should be excluded from the output.
+# The symbol name can be a fully qualified name, a word, or if the wildcard * is used,
+# a substring. Examples: ANamespace, AClass, AClass::ANamespace, ANamespace::*Test
+
+EXCLUDE_SYMBOLS =
+
+# The EXAMPLE_PATH tag can be used to specify one or more files or
+# directories that contain example code fragments that are included (see
+# the \include command).
+
+EXAMPLE_PATH =
+
+# If the value of the EXAMPLE_PATH tag contains directories, you can use the
+# EXAMPLE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp
+# and *.h) to filter out the source-files in the directories. If left
+# blank all files are included.
+
+EXAMPLE_PATTERNS = *
+
+# If the EXAMPLE_RECURSIVE tag is set to YES then subdirectories will be
+# searched for input files to be used with the \include or \dontinclude
+# commands irrespective of the value of the RECURSIVE tag.
+# Possible values are YES and NO. If left blank NO is used.
+
+EXAMPLE_RECURSIVE = NO
+
+# The IMAGE_PATH tag can be used to specify one or more files or
+# directories that contain image that are included in the documentation (see
+# the \image command).
+
+IMAGE_PATH =
+
+# The INPUT_FILTER tag can be used to specify a program that doxygen should
+# invoke to filter for each input file. Doxygen will invoke the filter program
+# by executing (via popen()) the command <filter> <input-file>, where <filter>
+# is the value of the INPUT_FILTER tag, and <input-file> is the name of an
+# input file. Doxygen will then use the output that the filter program writes
+# to standard output. If FILTER_PATTERNS is specified, this tag will be
+# ignored.
+
+INPUT_FILTER =
+
+# The FILTER_PATTERNS tag can be used to specify filters on a per file pattern
+# basis. Doxygen will compare the file name with each pattern and apply the
+# filter if there is a match. The filters are a list of the form:
+# pattern=filter (like *.cpp=my_cpp_filter). See INPUT_FILTER for further
+# info on how filters are used. If FILTER_PATTERNS is empty, INPUT_FILTER
+# is applied to all files.
+
+FILTER_PATTERNS =
+
+# If the FILTER_SOURCE_FILES tag is set to YES, the input filter (if set using
+# INPUT_FILTER) will be used to filter the input files when producing source
+# files to browse (i.e. when SOURCE_BROWSER is set to YES).
+
+FILTER_SOURCE_FILES = NO
+
+#---------------------------------------------------------------------------
+# configuration options related to source browsing
+#---------------------------------------------------------------------------
+
+# If the SOURCE_BROWSER tag is set to YES then a list of source files will
+# be generated. Documented entities will be cross-referenced with these sources.
+# Note: To get rid of all source code in the generated output, make sure also
+# VERBATIM_HEADERS is set to NO. If you have enabled CALL_GRAPH or CALLER_GRAPH
+# then you must also enable this option. If you don't then doxygen will produce
+# a warning and turn it on anyway
+
+SOURCE_BROWSER = NO
+
+# Setting the INLINE_SOURCES tag to YES will include the body
+# of functions and classes directly in the documentation.
+
+INLINE_SOURCES = NO
+
+# Setting the STRIP_CODE_COMMENTS tag to YES (the default) will instruct
+# doxygen to hide any special comment blocks from generated source code
+# fragments. Normal C and C++ comments will always remain visible.
+
+STRIP_CODE_COMMENTS = YES
+
+# If the REFERENCED_BY_RELATION tag is set to YES (the default)
+# then for each documented function all documented
+# functions referencing it will be listed.
+
+REFERENCED_BY_RELATION = NO
+
+# If the REFERENCES_RELATION tag is set to YES (the default)
+# then for each documented function all documented entities
+# called/used by that function will be listed.
+
+REFERENCES_RELATION = NO
+
+# If the REFERENCES_LINK_SOURCE tag is set to YES (the default)
+# and SOURCE_BROWSER tag is set to YES, then the hyperlinks from
+# functions in REFERENCES_RELATION and REFERENCED_BY_RELATION lists will
+# link to the source code. Otherwise they will link to the documentstion.
+
+REFERENCES_LINK_SOURCE = YES
+
+# If the USE_HTAGS tag is set to YES then the references to source code
+# will point to the HTML generated by the htags(1) tool instead of doxygen
+# built-in source browser. The htags tool is part of GNU's global source
+# tagging system (see http://www.gnu.org/software/global/global.html). You
+# will need version 4.8.6 or higher.
+
+USE_HTAGS = NO
+
+# If the VERBATIM_HEADERS tag is set to YES (the default) then Doxygen
+# will generate a verbatim copy of the header file for each class for
+# which an include is specified. Set to NO to disable this.
+
+VERBATIM_HEADERS = NO
+
+#---------------------------------------------------------------------------
+# configuration options related to the alphabetical class index
+#---------------------------------------------------------------------------
+
+# If the ALPHABETICAL_INDEX tag is set to YES, an alphabetical index
+# of all compounds will be generated. Enable this if the project
+# contains a lot of classes, structs, unions or interfaces.
+
+ALPHABETICAL_INDEX = NO
+
+# If the alphabetical index is enabled (see ALPHABETICAL_INDEX) then
+# the COLS_IN_ALPHA_INDEX tag can be used to specify the number of columns
+# in which this list will be split (can be a number in the range [1..20])
+
+COLS_IN_ALPHA_INDEX = 5
+
+# In case all classes in a project start with a common prefix, all
+# classes will be put under the same header in the alphabetical index.
+# The IGNORE_PREFIX tag can be used to specify one or more prefixes that
+# should be ignored while generating the index headers.
+
+IGNORE_PREFIX =
+
+#---------------------------------------------------------------------------
+# configuration options related to the HTML output
+#---------------------------------------------------------------------------
+
+# If the GENERATE_HTML tag is set to YES (the default) Doxygen will
+# generate HTML output.
+
+GENERATE_HTML = YES
+
+# The HTML_OUTPUT tag is used to specify where the HTML docs will be put.
+# If a relative path is entered the value of OUTPUT_DIRECTORY will be
+# put in front of it. If left blank `html' will be used as the default path.
+
+HTML_OUTPUT = html
+
+# The HTML_FILE_EXTENSION tag can be used to specify the file extension for
+# each generated HTML page (for example: .htm,.php,.asp). If it is left blank
+# doxygen will generate files with .html extension.
+
+HTML_FILE_EXTENSION = .html
+
+# The HTML_HEADER tag can be used to specify a personal HTML header for
+# each generated HTML page. If it is left blank doxygen will generate a
+# standard header.
+
+HTML_HEADER =
+
+# The HTML_FOOTER tag can be used to specify a personal HTML footer for
+# each generated HTML page. If it is left blank doxygen will generate a
+# standard footer.
+
+HTML_FOOTER =
+
+# The HTML_STYLESHEET tag can be used to specify a user-defined cascading
+# style sheet that is used by each HTML page. It can be used to
+# fine-tune the look of the HTML output. If the tag is left blank doxygen
+# will generate a default style sheet. Note that doxygen will try to copy
+# the style sheet file to the HTML output directory, so don't put your own
+# stylesheet in the HTML output directory as well, or it will be erased!
+
+HTML_STYLESHEET =
+
+# If the HTML_ALIGN_MEMBERS tag is set to YES, the members of classes,
+# files or namespaces will be aligned in HTML using tables. If set to
+# NO a bullet list will be used.
+
+HTML_ALIGN_MEMBERS = YES
+
+# If the GENERATE_HTMLHELP tag is set to YES, additional index files
+# will be generated that can be used as input for tools like the
+# Microsoft HTML help workshop to generate a compressed HTML help file (.chm)
+# of the generated HTML documentation.
+
+GENERATE_HTMLHELP = NO
+
+# If the HTML_DYNAMIC_SECTIONS tag is set to YES then the generated HTML
+# documentation will contain sections that can be hidden and shown after the
+# page has loaded. For this to work a browser that supports
+# JavaScript and DHTML is required (for instance Mozilla 1.0+, Firefox
+# Netscape 6.0+, Internet explorer 5.0+, Konqueror, or Safari).
+
+HTML_DYNAMIC_SECTIONS = NO
+
+# If the GENERATE_HTMLHELP tag is set to YES, the CHM_FILE tag can
+# be used to specify the file name of the resulting .chm file. You
+# can add a path in front of the file if the result should not be
+# written to the html output directory.
+
+CHM_FILE =
+
+# If the GENERATE_HTMLHELP tag is set to YES, the HHC_LOCATION tag can
+# be used to specify the location (absolute path including file name) of
+# the HTML help compiler (hhc.exe). If non-empty doxygen will try to run
+# the HTML help compiler on the generated index.hhp.
+
+HHC_LOCATION =
+
+# If the GENERATE_HTMLHELP tag is set to YES, the GENERATE_CHI flag
+# controls if a separate .chi index file is generated (YES) or that
+# it should be included in the master .chm file (NO).
+
+GENERATE_CHI = NO
+
+# If the GENERATE_HTMLHELP tag is set to YES, the BINARY_TOC flag
+# controls whether a binary table of contents is generated (YES) or a
+# normal table of contents (NO) in the .chm file.
+
+BINARY_TOC = NO
+
+# The TOC_EXPAND flag can be set to YES to add extra items for group members
+# to the contents of the HTML help documentation and to the tree view.
+
+TOC_EXPAND = NO
+
+# The DISABLE_INDEX tag can be used to turn on/off the condensed index at
+# top of each HTML page. The value NO (the default) enables the index and
+# the value YES disables it.
+
+DISABLE_INDEX = NO
+
+# This tag can be used to set the number of enum values (range [1..20])
+# that doxygen will group on one line in the generated HTML documentation.
+
+ENUM_VALUES_PER_LINE = 4
+
+# If the GENERATE_TREEVIEW tag is set to YES, a side panel will be
+# generated containing a tree-like index structure (just like the one that
+# is generated for HTML Help). For this to work a browser that supports
+# JavaScript, DHTML, CSS and frames is required (for instance Mozilla 1.0+,
+# Netscape 6.0+, Internet explorer 5.0+, or Konqueror). Windows users are
+# probably better off using the HTML help feature.
+
+GENERATE_TREEVIEW = NO
+
+# If the treeview is enabled (see GENERATE_TREEVIEW) then this tag can be
+# used to set the initial width (in pixels) of the frame in which the tree
+# is shown.
+
+TREEVIEW_WIDTH = 250
+
+#---------------------------------------------------------------------------
+# configuration options related to the LaTeX output
+#---------------------------------------------------------------------------
+
+# If the GENERATE_LATEX tag is set to YES (the default) Doxygen will
+# generate Latex output.
+
+GENERATE_LATEX = YES
+
+# The LATEX_OUTPUT tag is used to specify where the LaTeX docs will be put.
+# If a relative path is entered the value of OUTPUT_DIRECTORY will be
+# put in front of it. If left blank `latex' will be used as the default path.
+
+LATEX_OUTPUT = latex
+
+# The LATEX_CMD_NAME tag can be used to specify the LaTeX command name to be
+# invoked. If left blank `latex' will be used as the default command name.
+
+LATEX_CMD_NAME = latex
+
+# The MAKEINDEX_CMD_NAME tag can be used to specify the command name to
+# generate index for LaTeX. If left blank `makeindex' will be used as the
+# default command name.
+
+MAKEINDEX_CMD_NAME = makeindex
+
+# If the COMPACT_LATEX tag is set to YES Doxygen generates more compact
+# LaTeX documents. This may be useful for small projects and may help to
+# save some trees in general.
+
+COMPACT_LATEX = NO
+
+# The PAPER_TYPE tag can be used to set the paper type that is used
+# by the printer. Possible values are: a4, a4wide, letter, legal and
+# executive. If left blank a4wide will be used.
+
+PAPER_TYPE = a4wide
+
+# The EXTRA_PACKAGES tag can be to specify one or more names of LaTeX
+# packages that should be included in the LaTeX output.
+
+EXTRA_PACKAGES =
+
+# The LATEX_HEADER tag can be used to specify a personal LaTeX header for
+# the generated latex document. The header should contain everything until
+# the first chapter. If it is left blank doxygen will generate a
+# standard header. Notice: only use this tag if you know what you are doing!
+
+LATEX_HEADER =
+
+# If the PDF_HYPERLINKS tag is set to YES, the LaTeX that is generated
+# is prepared for conversion to pdf (using ps2pdf). The pdf file will
+# contain links (just like the HTML output) instead of page references
+# This makes the output suitable for online browsing using a pdf viewer.
+
+PDF_HYPERLINKS = NO
+
+# If the USE_PDFLATEX tag is set to YES, pdflatex will be used instead of
+# plain latex in the generated Makefile. Set this option to YES to get a
+# higher quality PDF documentation.
+
+USE_PDFLATEX = NO
+
+# If the LATEX_BATCHMODE tag is set to YES, doxygen will add the \\batchmode.
+# command to the generated LaTeX files. This will instruct LaTeX to keep
+# running if errors occur, instead of asking the user for help.
+# This option is also used when generating formulas in HTML.
+
+LATEX_BATCHMODE = NO
+
+# If LATEX_HIDE_INDICES is set to YES then doxygen will not
+# include the index chapters (such as File Index, Compound Index, etc.)
+# in the output.
+
+LATEX_HIDE_INDICES = NO
+
+#---------------------------------------------------------------------------
+# configuration options related to the RTF output
+#---------------------------------------------------------------------------
+
+# If the GENERATE_RTF tag is set to YES Doxygen will generate RTF output
+# The RTF output is optimized for Word 97 and may not look very pretty with
+# other RTF readers or editors.
+
+GENERATE_RTF = NO
+
+# The RTF_OUTPUT tag is used to specify where the RTF docs will be put.
+# If a relative path is entered the value of OUTPUT_DIRECTORY will be
+# put in front of it. If left blank `rtf' will be used as the default path.
+
+RTF_OUTPUT = rtf
+
+# If the COMPACT_RTF tag is set to YES Doxygen generates more compact
+# RTF documents. This may be useful for small projects and may help to
+# save some trees in general.
+
+COMPACT_RTF = NO
+
+# If the RTF_HYPERLINKS tag is set to YES, the RTF that is generated
+# will contain hyperlink fields. The RTF file will
+# contain links (just like the HTML output) instead of page references.
+# This makes the output suitable for online browsing using WORD or other
+# programs which support those fields.
+# Note: wordpad (write) and others do not support links.
+
+RTF_HYPERLINKS = NO
+
+# Load stylesheet definitions from file. Syntax is similar to doxygen's
+# config file, i.e. a series of assignments. You only have to provide
+# replacements, missing definitions are set to their default value.
+
+RTF_STYLESHEET_FILE =
+
+# Set optional variables used in the generation of an rtf document.
+# Syntax is similar to doxygen's config file.
+
+RTF_EXTENSIONS_FILE =
+
+#---------------------------------------------------------------------------
+# configuration options related to the man page output
+#---------------------------------------------------------------------------
+
+# If the GENERATE_MAN tag is set to YES (the default) Doxygen will
+# generate man pages
+
+GENERATE_MAN = NO
+
+# The MAN_OUTPUT tag is used to specify where the man pages will be put.
+# If a relative path is entered the value of OUTPUT_DIRECTORY will be
+# put in front of it. If left blank `man' will be used as the default path.
+
+MAN_OUTPUT = man
+
+# The MAN_EXTENSION tag determines the extension that is added to
+# the generated man pages (default is the subroutine's section .3)
+
+MAN_EXTENSION = .3
+
+# If the MAN_LINKS tag is set to YES and Doxygen generates man output,
+# then it will generate one additional man file for each entity
+# documented in the real man page(s). These additional files
+# only source the real man page, but without them the man command
+# would be unable to find the correct page. The default is NO.
+
+MAN_LINKS = NO
+
+#---------------------------------------------------------------------------
+# configuration options related to the XML output
+#---------------------------------------------------------------------------
+
+# If the GENERATE_XML tag is set to YES Doxygen will
+# generate an XML file that captures the structure of
+# the code including all documentation.
+
+GENERATE_XML = NO
+
+# The XML_OUTPUT tag is used to specify where the XML pages will be put.
+# If a relative path is entered the value of OUTPUT_DIRECTORY will be
+# put in front of it. If left blank `xml' will be used as the default path.
+
+XML_OUTPUT = xml
+
+# The XML_SCHEMA tag can be used to specify an XML schema,
+# which can be used by a validating XML parser to check the
+# syntax of the XML files.
+
+XML_SCHEMA =
+
+# The XML_DTD tag can be used to specify an XML DTD,
+# which can be used by a validating XML parser to check the
+# syntax of the XML files.
+
+XML_DTD =
+
+# If the XML_PROGRAMLISTING tag is set to YES Doxygen will
+# dump the program listings (including syntax highlighting
+# and cross-referencing information) to the XML output. Note that
+# enabling this will significantly increase the size of the XML output.
+
+XML_PROGRAMLISTING = YES
+
+#---------------------------------------------------------------------------
+# configuration options for the AutoGen Definitions output
+#---------------------------------------------------------------------------
+
+# If the GENERATE_AUTOGEN_DEF tag is set to YES Doxygen will
+# generate an AutoGen Definitions (see autogen.sf.net) file
+# that captures the structure of the code including all
+# documentation. Note that this feature is still experimental
+# and incomplete at the moment.
+
+GENERATE_AUTOGEN_DEF = NO
+
+#---------------------------------------------------------------------------
+# configuration options related to the Perl module output
+#---------------------------------------------------------------------------
+
+# If the GENERATE_PERLMOD tag is set to YES Doxygen will
+# generate a Perl module file that captures the structure of
+# the code including all documentation. Note that this
+# feature is still experimental and incomplete at the
+# moment.
+
+GENERATE_PERLMOD = NO
+
+# If the PERLMOD_LATEX tag is set to YES Doxygen will generate
+# the necessary Makefile rules, Perl scripts and LaTeX code to be able
+# to generate PDF and DVI output from the Perl module output.
+
+PERLMOD_LATEX = NO
+
+# If the PERLMOD_PRETTY tag is set to YES the Perl module output will be
+# nicely formatted so it can be parsed by a human reader. This is useful
+# if you want to understand what is going on. On the other hand, if this
+# tag is set to NO the size of the Perl module output will be much smaller
+# and Perl will parse it just the same.
+
+PERLMOD_PRETTY = YES
+
+# The names of the make variables in the generated doxyrules.make file
+# are prefixed with the string contained in PERLMOD_MAKEVAR_PREFIX.
+# This is useful so different doxyrules.make files included by the same
+# Makefile don't overwrite each other's variables.
+
+PERLMOD_MAKEVAR_PREFIX =
+
+#---------------------------------------------------------------------------
+# Configuration options related to the preprocessor
+#---------------------------------------------------------------------------
+
+# If the ENABLE_PREPROCESSING tag is set to YES (the default) Doxygen will
+# evaluate all C-preprocessor directives found in the sources and include
+# files.
+
+ENABLE_PREPROCESSING = YES
+
+# If the MACRO_EXPANSION tag is set to YES Doxygen will expand all macro
+# names in the source code. If set to NO (the default) only conditional
+# compilation will be performed. Macro expansion can be done in a controlled
+# way by setting EXPAND_ONLY_PREDEF to YES.
+
+MACRO_EXPANSION = NO
+
+# If the EXPAND_ONLY_PREDEF and MACRO_EXPANSION tags are both set to YES
+# then the macro expansion is limited to the macros specified with the
+# PREDEFINED and EXPAND_AS_DEFINED tags.
+
+EXPAND_ONLY_PREDEF = NO
+
+# If the SEARCH_INCLUDES tag is set to YES (the default) the includes files
+# in the INCLUDE_PATH (see below) will be search if a #include is found.
+
+SEARCH_INCLUDES = YES
+
+# The INCLUDE_PATH tag can be used to specify one or more directories that
+# contain include files that are not input files but should be processed by
+# the preprocessor.
+
+INCLUDE_PATH =
+
+# You can use the INCLUDE_FILE_PATTERNS tag to specify one or more wildcard
+# patterns (like *.h and *.hpp) to filter out the header-files in the
+# directories. If left blank, the patterns specified with FILE_PATTERNS will
+# be used.
+
+INCLUDE_FILE_PATTERNS =
+
+# The PREDEFINED tag can be used to specify one or more macro names that
+# are defined before the preprocessor is started (similar to the -D option of
+# gcc). The argument of the tag is a list of macros of the form: name
+# or name=definition (no spaces). If the definition and the = are
+# omitted =1 is assumed. To prevent a macro definition from being
+# undefined via #undef or recursively expanded use the := operator
+# instead of the = operator.
+
+PREDEFINED =
+
+# If the MACRO_EXPANSION and EXPAND_ONLY_PREDEF tags are set to YES then
+# this tag can be used to specify a list of macro names that should be expanded.
+# The macro definition that is found in the sources will be used.
+# Use the PREDEFINED tag if you want to use a different macro definition.
+
+EXPAND_AS_DEFINED =
+
+# If the SKIP_FUNCTION_MACROS tag is set to YES (the default) then
+# doxygen's preprocessor will remove all function-like macros that are alone
+# on a line, have an all uppercase name, and do not end with a semicolon. Such
+# function macros are typically used for boiler-plate code, and will confuse
+# the parser if not removed.
+
+SKIP_FUNCTION_MACROS = YES
+
+#---------------------------------------------------------------------------
+# Configuration::additions related to external references
+#---------------------------------------------------------------------------
+
+# The TAGFILES option can be used to specify one or more tagfiles.
+# Optionally an initial location of the external documentation
+# can be added for each tagfile. The format of a tag file without
+# this location is as follows:
+# TAGFILES = file1 file2 ...
+# Adding location for the tag files is done as follows:
+# TAGFILES = file1=loc1 "file2 = loc2" ...
+# where "loc1" and "loc2" can be relative or absolute paths or
+# URLs. If a location is present for each tag, the installdox tool
+# does not have to be run to correct the links.
+# Note that each tag file must have a unique name
+# (where the name does NOT include the path)
+# If a tag file is not located in the directory in which doxygen
+# is run, you must also specify the path to the tagfile here.
+
+TAGFILES =
+
+# When a file name is specified after GENERATE_TAGFILE, doxygen will create
+# a tag file that is based on the input files it reads.
+
+GENERATE_TAGFILE =
+
+# If the ALLEXTERNALS tag is set to YES all external classes will be listed
+# in the class index. If set to NO only the inherited external classes
+# will be listed.
+
+ALLEXTERNALS = NO
+
+# If the EXTERNAL_GROUPS tag is set to YES all external groups will be listed
+# in the modules index. If set to NO, only the current project's groups will
+# be listed.
+
+EXTERNAL_GROUPS = YES
+
+# The PERL_PATH should be the absolute path and name of the perl script
+# interpreter (i.e. the result of `which perl').
+
+PERL_PATH = /usr/bin/perl
+
+#---------------------------------------------------------------------------
+# Configuration options related to the dot tool
+#---------------------------------------------------------------------------
+
+# If the CLASS_DIAGRAMS tag is set to YES (the default) Doxygen will
+# generate a inheritance diagram (in HTML, RTF and LaTeX) for classes with base
+# or super classes. Setting the tag to NO turns the diagrams off. Note that
+# this option is superseded by the HAVE_DOT option below. This is only a
+# fallback. It is recommended to install and use dot, since it yields more
+# powerful graphs.
+
+CLASS_DIAGRAMS = YES
+
+# You can define message sequence charts within doxygen comments using the \msc
+# command. Doxygen will then run the mscgen tool (see http://www.mcternan.me.uk/mscgen/) to
+# produce the chart and insert it in the documentation. The MSCGEN_PATH tag allows you to
+# specify the directory where the mscgen tool resides. If left empty the tool is assumed to
+# be found in the default search path.
+
+MSCGEN_PATH =
+
+# If set to YES, the inheritance and collaboration graphs will hide
+# inheritance and usage relations if the target is undocumented
+# or is not a class.
+
+HIDE_UNDOC_RELATIONS = YES
+
+# If you set the HAVE_DOT tag to YES then doxygen will assume the dot tool is
+# available from the path. This tool is part of Graphviz, a graph visualization
+# toolkit from AT&T and Lucent Bell Labs. The other options in this section
+# have no effect if this option is set to NO (the default)
+
+HAVE_DOT = NO
+
+# If the CLASS_GRAPH and HAVE_DOT tags are set to YES then doxygen
+# will generate a graph for each documented class showing the direct and
+# indirect inheritance relations. Setting this tag to YES will force the
+# the CLASS_DIAGRAMS tag to NO.
+
+CLASS_GRAPH = YES
+
+# If the COLLABORATION_GRAPH and HAVE_DOT tags are set to YES then doxygen
+# will generate a graph for each documented class showing the direct and
+# indirect implementation dependencies (inheritance, containment, and
+# class references variables) of the class with other documented classes.
+
+COLLABORATION_GRAPH = YES
+
+# If the GROUP_GRAPHS and HAVE_DOT tags are set to YES then doxygen
+# will generate a graph for groups, showing the direct groups dependencies
+
+GROUP_GRAPHS = YES
+
+# If the UML_LOOK tag is set to YES doxygen will generate inheritance and
+# collaboration diagrams in a style similar to the OMG's Unified Modeling
+# Language.
+
+UML_LOOK = NO
+
+# If set to YES, the inheritance and collaboration graphs will show the
+# relations between templates and their instances.
+
+TEMPLATE_RELATIONS = NO
+
+# If the ENABLE_PREPROCESSING, SEARCH_INCLUDES, INCLUDE_GRAPH, and HAVE_DOT
+# tags are set to YES then doxygen will generate a graph for each documented
+# file showing the direct and indirect include dependencies of the file with
+# other documented files.
+
+INCLUDE_GRAPH = YES
+
+# If the ENABLE_PREPROCESSING, SEARCH_INCLUDES, INCLUDED_BY_GRAPH, and
+# HAVE_DOT tags are set to YES then doxygen will generate a graph for each
+# documented header file showing the documented files that directly or
+# indirectly include this file.
+
+INCLUDED_BY_GRAPH = YES
+
+# If the CALL_GRAPH, SOURCE_BROWSER and HAVE_DOT tags are set to YES then doxygen will
+# generate a call dependency graph for every global function or class method.
+# Note that enabling this option will significantly increase the time of a run.
+# So in most cases it will be better to enable call graphs for selected
+# functions only using the \callgraph command.
+
+CALL_GRAPH = NO
+
+# If the CALLER_GRAPH, SOURCE_BROWSER and HAVE_DOT tags are set to YES then doxygen will
+# generate a caller dependency graph for every global function or class method.
+# Note that enabling this option will significantly increase the time of a run.
+# So in most cases it will be better to enable caller graphs for selected
+# functions only using the \callergraph command.
+
+CALLER_GRAPH = NO
+
+# If the GRAPHICAL_HIERARCHY and HAVE_DOT tags are set to YES then doxygen
+# will graphical hierarchy of all classes instead of a textual one.
+
+GRAPHICAL_HIERARCHY = YES
+
+# If the DIRECTORY_GRAPH, SHOW_DIRECTORIES and HAVE_DOT tags are set to YES
+# then doxygen will show the dependencies a directory has on other directories
+# in a graphical way. The dependency relations are determined by the #include
+# relations between the files in the directories.
+
+DIRECTORY_GRAPH = YES
+
+# The DOT_IMAGE_FORMAT tag can be used to set the image format of the images
+# generated by dot. Possible values are png, jpg, or gif
+# If left blank png will be used.
+
+DOT_IMAGE_FORMAT = png
+
+# The tag DOT_PATH can be used to specify the path where the dot tool can be
+# found. If left blank, it is assumed the dot tool can be found in the path.
+
+DOT_PATH =
+
+# The DOTFILE_DIRS tag can be used to specify one or more directories that
+# contain dot files that are included in the documentation (see the
+# \dotfile command).
+
+DOTFILE_DIRS =
+
+# The MAX_DOT_GRAPH_MAX_NODES tag can be used to set the maximum number of
+# nodes that will be shown in the graph. If the number of nodes in a graph
+# becomes larger than this value, doxygen will truncate the graph, which is
+# visualized by representing a node as a red box. Note that doxygen if the number
+# of direct children of the root node in a graph is already larger than
+# MAX_DOT_GRAPH_NOTES then the graph will not be shown at all. Also note
+# that the size of a graph can be further restricted by MAX_DOT_GRAPH_DEPTH.
+
+DOT_GRAPH_MAX_NODES = 50
+
+# The MAX_DOT_GRAPH_DEPTH tag can be used to set the maximum depth of the
+# graphs generated by dot. A depth value of 3 means that only nodes reachable
+# from the root by following a path via at most 3 edges will be shown. Nodes
+# that lay further from the root node will be omitted. Note that setting this
+# option to 1 or 2 may greatly reduce the computation time needed for large
+# code bases. Also note that the size of a graph can be further restricted by
+# DOT_GRAPH_MAX_NODES. Using a depth of 0 means no depth restriction.
+
+MAX_DOT_GRAPH_DEPTH = 1000
+
+# Set the DOT_TRANSPARENT tag to YES to generate images with a transparent
+# background. This is disabled by default, which results in a white background.
+# Warning: Depending on the platform used, enabling this option may lead to
+# badly anti-aliased labels on the edges of a graph (i.e. they become hard to
+# read).
+
+DOT_TRANSPARENT = NO
+
+# Set the DOT_MULTI_TARGETS tag to YES allow dot to generate multiple output
+# files in one run (i.e. multiple -o and -T options on the command line). This
+# makes dot run faster, but since only newer versions of dot (>1.8.10)
+# support this, this feature is disabled by default.
+
+DOT_MULTI_TARGETS = NO
+
+# If the GENERATE_LEGEND tag is set to YES (the default) Doxygen will
+# generate a legend page explaining the meaning of the various boxes and
+# arrows in the dot generated graphs.
+
+GENERATE_LEGEND = YES
+
+# If the DOT_CLEANUP tag is set to YES (the default) Doxygen will
+# remove the intermediate dot files that are used to generate
+# the various graphs.
+
+DOT_CLEANUP = YES
+
+#---------------------------------------------------------------------------
+# Configuration::additions related to the search engine
+#---------------------------------------------------------------------------
+
+# The SEARCHENGINE tag specifies whether or not a search engine should be
+# used. If set to NO the values of all tags below this one will be ignored.
+
+SEARCHENGINE = NO
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..ba8c334bc025a5393e2423d08f7bff176bd6d335
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,177 @@
+##
+# Makefile for inerpolator classes and library
+# branched off libhlf release_02
+#
+# $Author: claudio $
+#
+# $Revision: 1.7 $
+#
+# $Log: Makefile,v $
+# Revision 1.7 2017-03-06 13:19:03 claudio
+# -O2 produces bugs...
+#
+# Revision 1.6 2017-03-06 13:13:55 claudio
+# optimize, cleaner check
+#
+# Revision 1.5 2017-03-06 13:02:02 claudio
+# test_inverse
+#
+# Revision 1.4 2017-03-03 15:39:58 claudio
+# advance minor release
+#
+# Revision 1.3 2017-03-03 15:32:26 claudio
+# correctetd spline::inverse_calculate if returned value is at range limit
+#
+# Revision 1.2 2012-09-19 14:21:49 claudio
+# compile as position independent code
+#
+# Revision 1.1.1.1 2012-09-05 08:17:22 claudio
+# Interpolator utility classes
+#
+#
+#
+
+
+
+CXX= g++
+CC = $(CXX)
+LD = $(CXX)
+DEBUG = -g
+RANLIB = ranlib
+#OPTIM = -O2
+
+INCLUDE = -I.
+
+CXXFLAGS += -fPIC -D_REENTRANT $(DEBUG) $(OPTIM) $(WARN) $(INCLUDE)
+CFLAGS = $(CXXFLAGS)
+PROJECTHOME = .
+##############################################
+# support for shared libray versioning
+#
+LFLAGS_SONAME = -Wl,-soname,
+LFLAGS_RPATH = -Wl,-rpath,
+SHLDFLAGS = -shared
+BASELIBNAME = libinterpolator
+SHLIB_SUFFIX = so
+
+# release numbers for libraries
+#
+ LIBVERSION = 1
+ LIBRELEASE = 1
+ LIBSUBRELEASE = 1
+#
+
+LIBRARY = $(BASELIBNAME).a
+DT_SONAME = $(BASELIBNAME).$(SHLIB_SUFFIX).$(LIBVERSION)
+DT_SHLIB = $(BASELIBNAME).$(SHLIB_SUFFIX).$(LIBVERSION).$(LIBRELEASE).$(LIBSUBRELEASE)
+SHLIB = $(BASELIBNAME).$(SHLIB_SUFFIX)
+
+################################################
+#installation root
+PREFIX = /runtime
+##########################
+# installation directories
+DESTDIRLIB = $(PREFIX)/lib
+DESTDIRHEADERS = $(PREFIX)/include
+# header files to be installed for distribution
+INSTHEADERS = periodicspline.h spline.h multipolynomial.h interpolatingpolynomial.h interpolator.h
+################################################
+
+OBJS = periodicspline.o spline.o multipolynomial.o interpolator.o
+
+
+########################################################################################
+# compiler warnings
+WARN = -Wall -Wextra -Wno-deprecated -Woverloaded-virtual -Wsign-promo -Wshadow
+#WARN += -Weffc++ -Wold-style-cast -Wuninitialized -Wmissing-braces -Wparentheses -Wsequence-point -Wreturn-type -Wswitch -Wswitch-default -Wswitch-enum -Winit-self -Wundef -Wmissing-field-initializers -Winline
+WARN+= -pedantic -Wno-long-long # test with mysql due to long long variables!
+########################################################################################
+
+#######################################
+#options for coverage and performance analysis
+#DEBUG += -fprofile-arcs -ftest-coverage -pg
+#LDFLAGS += -fprofile-arcs -ftest-coverage -pg
+#####################################
+
+
+all: lib libstatic test
+
+lib: $(PROJECTHOME)/lib/$(DT_SHLIB)
+
+libstatic: $(PROJECTHOME)/lib/$(LIBRARY)
+
+check: check_spline check_multipoly
+
+
+test: test_spline test_inverse test_freespline test_multipoly test_periodicspline spltest
+
+test_spline: test_spline.o spline.o interpolator.o
+
+test_inverse: test_inverse.o spline.o interpolator.o
+
+test_freespline: test_freespline.o spline.o interpolator.o
+
+test_periodicspline: test_periodicspline.o periodicspline.o interpolator.o
+
+spltest: spltest.o spline.o interpolator.o
+
+test_multipoly: test_multipoly.o multipolynomial.o interpolator.o
+
+
+check_spline: test_spline
+ ./test_spline
+ gcov test_spline.cpp
+ gcov spline.cpp
+ gcov interpolator.cpp
+ valgrind --num-callers=20 --leak-check=yes --leak-resolution=high --leak-check=full --show-reachable=yes --log-file=test_spline.grind ./test_spline
+
+check_multipoly: test_multipoly
+ gcov test_multipoly.cpp
+ gcov multipoly.cpp
+ gcov interpolator.cpp
+ valgrind --num-callers=20 --leak-check=yes --leak-resolution=high --leak-check=full --show-reachable=yes --log-file=test_multipoly.grind ./test_multipoly
+
+$(PROJECTHOME)/lib/$(LIBRARY):$(OBJS) $(MAKEFILE)
+# @echo "Archiving $(LIBRARY) ... $(OBJS)"
+ mkdir -p -m 755 $(PROJECTHOME)/lib
+ ar $(ARFLAGS) $(PROJECTHOME)/lib/$(LIBRARY) $(OBJS)
+ $(RANLIB) $(PROJECTHOME)/lib/$(LIBRARY)
+# @echo "done"
+
+$(PROJECTHOME)/lib/$(DT_SHLIB): $(OBJS) $(MAKEFILE)
+# @echo "Creating $(DT_SONAME) ...$(OBJS)"
+ mkdir -p -m 755 $(PROJECTHOME)/lib
+ $(LD) $(SHLDFLAGS) $(LFLAGS_SONAME)$(DT_SONAME) -o $(PROJECTHOME)/lib/$(DT_SHLIB) $(OBJS)
+ ln -sf $(PROJECTHOME)/lib/$(DT_SHLIB) $(PROJECTHOME)/lib/$(SHLIB)
+ ln -sf $(PROJECTHOME)/lib/$(DT_SHLIB) $(PROJECTHOME)/lib/$(DT_SONAME)
+# @echo "done"
+
+doc:
+ doxygen 2>doxy.log
+
+docclean:
+ rm -rf doc doxy.log
+clean:
+ rm -f *.o core* *.gcov *.gcno *.gcda *.grind *grind.core.* *.valgrind dump* gmon.out doxygen_log.txt dump* *.csv
+ rm -f $(PROJECTHOME)/lib/*.a $(PROJECTHOME)/lib/*.so*
+ rm -f test_spline test_inverse test_periodicspline test_freespline test_multipoly test_magnets test_dipole test_quadrupole test_sextupole test_kicker test_solenoid spltest
+ rm -f *.csv
+
+distclean: clean docclean
+ rm -f *~ .kdbgrc.*
+
+install: lib
+ install -d $(DESTDIRLIB)
+ install -d $(DESTDIRHEADERS)
+ install $(PROJECTHOME)/lib/$(DT_SHLIB) $(DESTDIRLIB)
+ install $(INSTHEADERS) $(DESTDIRHEADERS)
+ rm -f $(DESTDIRLIB)/$(SHLIB)
+ rm -f $(DESTDIRLIB)/$(DT_SONAME)
+ cd $(DESTDIRLIB); ln -s $(DT_SHLIB) $(SHLIB);ln -s $(DT_SHLIB) $(DT_SONAME)
+
+uninstall:
+ rm -f $(DESTDIRLIB)/$(LIBRARY)
+ rm -f $(DESTDIRLIB)/$(SHLIB)
+ rm -f $(DESTDIRLIB)/$(DT_SHLIB)
+ rm -f $(DESTDIRLIB)/$(DT_SONAME)
+ rm -rf $(DESTDIRHEADERS)
diff --git a/README.md b/README.md
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..eda029306ba503a58f20c083bd7538622ee20e1a 100644
--- a/README.md
+++ b/README.md
@@ -0,0 +1,5 @@
+# interpolator
+
+C++ library for interpolation of funziontions of real variable y=f(x).
+
+Various splines and mutipolynomial methods are implemented.
diff --git a/interpolatingpolynomial.h b/interpolatingpolynomial.h
new file mode 100644
index 0000000000000000000000000000000000000000..6e91dadf7d5bddac7c2b5f58361cee5eaedf9fce
--- /dev/null
+++ b/interpolatingpolynomial.h
@@ -0,0 +1,42 @@
+/* $Author: claudio $
+*
+* $Revision: 1.1.1.1 $
+*
+* $Log: interpolatingpolynomial.h,v $
+* Revision 1.1.1.1 2012-09-05 08:17:22 claudio
+* Interpolator utility classes
+*
+*/
+
+#ifndef INTERPOLATOR_INTERPOLATINGPOLYNOMIAL_H
+#define INTERPOLATOR_INTERPOLATINGPOLYNOMIAL_H
+
+
+
+/**
+* @struct Interpolator::InterpolatingPolynomial
+* Auxiliary class ( a struct ) to be used as an argument for various class
+* constructors.
+* Filled by reading data from the database or from other data sources by means of factory classes
+*/
+
+
+
+namespace Interpolator{
+
+struct InterpolatingPolynomial
+{
+ double xMin; ///< minimum X value
+ double xMax; ///< maximum X value
+ doubleVector coefficient; ///< vector of polynomial cofficients. cofficient[i] is for polynomial coefficient of degree i
+};
+
+/**
+ * std::vector of InperpolatingPolynomial
+ */
+
+typedef std::vector<InterpolatingPolynomial> InterpolatingPolynoamialVector;
+
+} // close namespace Interpolator
+
+#endif // INTERPOLATOR_INTERPOLATINGPOLYNOMIAL_H
diff --git a/interpolator.cpp b/interpolator.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..dee53a6634f1746f8c5029854f0aae5e3e89b30e
--- /dev/null
+++ b/interpolator.cpp
@@ -0,0 +1,132 @@
+/*
+ * class Interpolator - implementation
+ *
+ * $Author: claudio $
+ *
+ * $Revision: 1.2 $
+ *
+ * $Log: interpolator.cpp,v $
+ * Revision 1.2 2017-03-03 15:32:26 claudio
+ * correctetd spline::inverse_calculate if returned value is at range limit
+ *
+ * Revision 1.1.1.1 2012-09-05 08:17:22 claudio
+ * Interpolator utility classes
+ *
+ *
+ * $Name: release_01_01_01 $
+ */
+
+#include "interpolator.h"
+#include <cmath>
+#include <limits>
+
+
+//mark revision name for easier control of installations
+#define CVSVERSION "$Name: release_01_01_01 $"
+const char *cvs_version = CVSVERSION;
+
+
+namespace Interpolator{
+
+Interpolator::Interpolator():_min(NAN),_max(NAN),_tolbracket(1.0e-4),_tolsolve(1.0e-8),_maxitersolve(100),_ntry(50){}
+
+Interpolator::~Interpolator () {}
+
+void Interpolator::get_range(double& min, double& max )
+{
+ min=_min;
+ max=_max;
+ return;
+}
+
+
+double Interpolator::inverse_evaluate(double goal,double initialEstimate)
+{
+ double xmin=0,xmax=0;
+ if(bracket(goal,initialEstimate,_tolbracket,xmin,xmax))
+ return solve(goal,xmin,xmax,_tolsolve,_maxitersolve);
+ else return initialEstimate;
+}
+
+bool Interpolator::bracket(double goal,double initial,double tol, double& xlow, double& xhigh)
+{
+ if(tol<=std::numeric_limits<double>::epsilon()) throw std::range_error("bracket: tol too small");
+ const double fact=1.618;
+ const unsigned int ntry=_ntry;
+
+ double x1,x2;
+ x1=initial-tol; if(x1 < _min) x1=_min;
+ x2=initial+tol; if(x2 > _max) x2=_max;
+ double f1,f2;
+ f1=goal - evaluate(x1);
+ f2=goal - evaluate(x2);
+ for(unsigned int j=0; j < ntry; j++){
+ if(f1*f2 <=0.0){
+ xlow=x1;
+ xhigh=x2;
+ return true;
+ }
+ if(fabs(f1) < fabs(f2)){
+ x1 += fact*(x1-x2);
+ if(x1 < _min) x1=_min;
+ if(x1 > _max) x1=_max;
+ f1=goal - evaluate(x1);
+ }
+ else{
+ x2 += fact *(x2-x1);
+ if(x2 < _min) x2=_min;
+ if(x2 > _max) x2=_max;
+ f2=goal - evaluate(x2);
+ }
+ }
+ return false;
+}
+
+
+double Interpolator::solve(double goal,double x1,double x2,double tol,unsigned int maxiter) throw (std::range_error)
+{
+ double fl,fh,xl,xh,del,f,rtf,dx;
+ fl = goal - evaluate(x1);
+ fh = goal - evaluate(x2);
+ if( fl*fh > 0.0) throw std::range_error("solve: initial interval does not braket solution");
+ if(fl < 0.0 ){
+ xl=x1;
+ xh=x2;
+ }
+ else{
+ xl=x2;
+ xh=x1;
+ double swap=fl;
+ fl=fh;
+ fh=swap;
+ }
+ dx=xh-xl;
+ for(unsigned int j=0; j < maxiter; j++){
+ rtf=xl+dx*fl/(fl-fh);
+ f= goal - evaluate(rtf);
+ if(f < 0.0){
+ del=xl-rtf;
+ xl=rtf;
+ fl=f;
+ }
+ else{
+ del=xh-rtf;
+ xh=rtf;
+ fh=f;
+ }
+ dx=xh-xl;
+ if(fabs(del) < tol || f == 0.0) return rtf;
+ }
+ return NAN; //maxiter exceeded
+}
+
+void Interpolator::set_solve_params(double tolbracket,double tolsolve, unsigned int maxitersolve, unsigned int ntry)
+{
+ _tolbracket=tolbracket;
+ _tolsolve=tolsolve;
+ _maxitersolve=maxitersolve;
+ _ntry=ntry;
+ return;
+}
+
+} // close namespace Interpolator
diff --git a/interpolator.h b/interpolator.h
new file mode 100644
index 0000000000000000000000000000000000000000..d2dbae40c0c69a28421d33fe46505deb782cb026
--- /dev/null
+++ b/interpolator.h
@@ -0,0 +1,135 @@
+#ifndef INTERPOLATOR_INTERPOLATOR_H
+#define INTERPOLATOR_INTERPOLATOR_H
+
+#include <stdexcept>
+#include <vector>
+#include <cmath>
+
+
+/**
+* @class Interpolator::Interpolator
+* Base class for interpolation, typically used for calibration tables of fields,currents etc.
+* It is used to evalute y=f(x) with x and y both real [ R(1) -> R(1) ].
+* It is implemented by means of different interpolation strategies
+*/
+
+/* $Author: claudio $
+*
+* $Revision: 1.1.1.1 $
+*
+* $Log: interpolator.h,v $
+* Revision 1.1.1.1 2012-09-05 08:17:22 claudio
+* Interpolator utility classes
+*
+*
+*/
+
+namespace Interpolator{
+ /**
+ * typedef for easier declaration of some virtual methods
+ * in application programs std::vector<double> may be used wherever doubleVector is required
+ */
+typedef std::vector<double> doubleVector;
+class Interpolator
+{
+public:
+ /**
+ * default Constructor
+ * min and max are set to a "legal" but un-initialized values
+ */
+ Interpolator ( );
+
+ virtual ~Interpolator ( );
+
+ /**
+ * calculate the interpolated value v=f(in_val) y=f(x)
+ * @return interpolated value
+ * @param in_val input value (x)
+ * @return interpolated value y=f(x)
+ * @exception std::range_error
+ */
+ virtual double evaluate (double in_val ) throw (std::range_error)=0;
+
+
+ /**
+ * estimate the value of the inverse iterpolator x=g(y)
+ * We make the hypothesys that the interpolated function is invertible!!
+ * @param goal the recquired output (y)
+ * @param initialEstimate initial estimate for x
+ * @return x, NAN if the algorithm does not converge
+ */
+
+ virtual double inverse_evaluate(double goal,double initialEstimate);
+
+ /**
+ * calculate the interpolated value v[ ]=f(in_val[ ]). Vector version
+ * @return doubleVector
+ * @param in_val
+ * @exception std::range_error
+ */
+ virtual doubleVector evaluate (doubleVector in_val ) throw(std::range_error)=0;
+
+ /**
+ * returns the existence filed of y=f(x). It is assumed the f(x) has values for all
+ * the range [min,max).
+ * @param min
+ * @param max
+ */
+ void get_range(double& min, double& max );
+
+
+ /**
+ * auxiliary function used in inverse_evaluate
+ * finds an interval in x where the interpolated functions brackets the goal values
+ * (See "Numerical Recipes in C++", zbrak algorithm )
+ * @param goal the value to be braketed
+ * @param initial x value around which the interaval is searched
+ * @param tol half-width of initial x interval
+ * @param xlow [out] lower limit of interval
+ * @param xhigh [out] upper limit of interval
+ * @return true if successful
+ */
+
+ bool bracket(double goal,double initial,double tol, double& xlow, double& xhigh);
+
+
+ /**
+ * auxiliary function used in inverse_evaluate
+ * finds an finds the requested goal value by means of the "regula falsi" root finding
+ * algorithm (See "Numerical Recipes in C++", rtflsp algorithm )
+ * @param goal the y value to be reached
+ * @param xlow lower bound for x
+ * @param xhigh upper bound for x
+ * @param tol precision of x;
+ * @param maxiter maximum number of iterations
+ * @return x stasfying relation goal == evaluate(x)
+ * @exception std::range_error initial range does not bracket a solution
+ */
+
+ double solve(double goal,double xlow,double xhigh,double tol,unsigned int maxiter=100) throw (std::range_error);
+
+
+ /**
+ * special auxiliary function to modify the default convergence parameters of bracket() and solve()
+ * @param tolbracket tolerance for bracketing algorithm
+ * @param tolsolve tolerance for solve algorithm
+ * @param maxitersolve maximum number of iterations for the solve algorithm
+ * @param ntry maximum number of iterations for the bracket algorithm
+ */
+
+ void set_solve_params(double tolbracket,double tolsolve, unsigned int maxitersolve, unsigned int ntry);
+
+protected:
+ double _min; /**< minimum value of x */
+ double _max; /**< maximum value of x */
+ double _tolbracket; /**< tolerance for bracketing algorithm */
+ double _tolsolve; /**< tolerance for solve algorithm */
+ unsigned int _maxitersolve; /**< maximum number of iterations for the solve algorithm */
+ unsigned int _ntry; /**< maximum number of iterations for the bracket algorithm */
+
+};
+
+
+} //close namespace Interpolator;
+
+#endif // INTERPOLATOR_INTERPOLATOR_H
diff --git a/multipolynomial.cpp b/multipolynomial.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a6432bf2d2cd36dd732e4b5e471f2220b57f416e
--- /dev/null
+++ b/multipolynomial.cpp
@@ -0,0 +1,124 @@
+/*
+ * class multiPolynomial - implementation
+ *
+ * $Author: claudio $
+ *
+ * $Revision: 1.2 $
+ *
+ * $Log: multipolynomial.cpp,v $
+ * Revision 1.2 2017-03-03 15:32:26 claudio
+ * correctetd spline::inverse_calculate if returned value is at range limit
+ *
+ * Revision 1.1.1.1 2012-09-05 08:17:22 claudio
+ * Interpolator utility classes
+ *
+ */
+
+#include "multipolynomial.h"
+#include <stdexcept>
+
+namespace Interpolator{
+
+multiPolynomial::multiPolynomial ( InterpolatingPolynoamialVector& poly){
+ n_poly=poly.size();
+ if( ! n_poly ) throw std::length_error("multiPolynomial(): at least one polynomial is recquired");
+
+ //scan limits and check that they are in order and continguos
+ {
+ double low,high;
+ unsigned int i=0;
+ low=poly[i].xMin;
+ high=poly[i].xMax;
+ while(++i < n_poly){
+ if(low >= high) throw std::domain_error("multiPolynomial: xMin >= xMax");
+ if(high != poly[i].xMin) throw std::domain_error("multiPolynomial: not contiguos");
+ low=poly[i].xMin;
+ high=poly[i].xMax;
+ }
+ }
+ //good - create structure to store the data
+ range = new double[n_poly+1];
+ range[0]=poly[0].xMin; //first is special
+ for (unsigned int i=0; i<n_poly ; i++){
+ range[i+1]=poly[i].xMax;
+ _poly.push_back(poly[i].coefficient);
+ }
+ Interpolator::_min=range[0];
+ Interpolator::_max=range[n_poly];
+}
+
+
+multiPolynomial::multiPolynomial(multiPolynomial& src):Interpolator()
+{
+ _min=src._min;
+ _max=src._max;
+ n_poly=src.n_poly;
+ range = new double[n_poly+1];
+ range[0]=src.range[0];
+ for (unsigned int i=0; i<n_poly ; i++){
+ range[i+1]=src.range[i+1];
+ }
+ _poly=src._poly; //exploit assignment operator from stl
+}
+
+multiPolynomial& multiPolynomial::operator=(const multiPolynomial& src)
+{
+ _min=src._min;
+ _max=src._max;
+ n_poly=src.n_poly;
+ range = new double[n_poly+1];
+ range[0]=src.range[0];
+ for (unsigned int i=0; i<n_poly ; i++){
+ range[i+1]=src.range[i+1];
+ }
+ _poly=src._poly; //exploit assignment operator from stl
+ return *this;
+}
+
+multiPolynomial::~multiPolynomial ( )
+{
+ delete [] range;
+}
+
+double multiPolynomial::evaluate(double v) throw (std::range_error)
+{
+ if ( v < _min || v > _max ) throw std::range_error("input value out of range");
+ //find the correct polynomial to use
+
+ unsigned int klo,khi,k;
+ klo=0;
+ khi=1;
+ while( khi < n_poly ){
+ if( v< range[khi] && v >= range[klo]) break;
+ klo++;khi++;
+ }
+ k=klo;
+
+ // Horner rule - can also evaluate the derivative
+ unsigned int degree=_poly[k].size()-1;
+ doubleVector* poly=&_poly[k];
+ double p=poly->at(degree);;
+ //double dp=0; - derivative , disabled
+ for( int i= degree-1; i >= 0; i-- ){
+ p = v*p + poly->at(i);
+ //dp = v*dp + p; //derivative evaluation
+ }
+ return p;
+}
+
+doubleVector multiPolynomial::evaluate (doubleVector in_val ) throw(std::range_error)
+{
+ //simplicistic implementation. Should make better use of the indexing
+ doubleVector outval;
+ for(unsigned int i=0; i < in_val.size(); i++) outval.push_back(evaluate(in_val[i]));
+ return outval;
+}
+
+double multiPolynomial::inverse_evaluate(double goal,double initialEstimate)
+{
+ //dirty work-around untill the algorithm is tested also for multy-polinomial!!!
+ //TODO : implement inverse_evaluate
+ return initialEstimate;
+}
+
+} // close namespace Interpolator
diff --git a/multipolynomial.h b/multipolynomial.h
new file mode 100644
index 0000000000000000000000000000000000000000..7c509c64cc628779f8b3d9607ea68ef7f7581105
--- /dev/null
+++ b/multipolynomial.h
@@ -0,0 +1,102 @@
+/*
+* $Author: claudio $
+*
+* $Revision: 1.1.1.1 $
+*
+* $Log: multipolynomial.h,v $
+* Revision 1.1.1.1 2012-09-05 08:17:22 claudio
+* Interpolator utility classes
+*
+*
+*/
+#ifndef INTERPOLATOR_MULTIPOLYNOMIAL_H
+#define INTERPOLATOR_MULTIPOLYNOMIAL_H
+
+#include <interpolator.h>
+#include <interpolatingpolynomial.h>
+
+/**
+* @class Interpolator::multiPolynomial
+* Class for interpolating fields,currents etc.
+* It is used to evalute y=f(x) with x,y belonging to R(1).
+* It is implemented by means of a set of interpolating polynomials,
+* Each polynomial is valid for the input range [xmini,xmaxi). The polinomial ranges must be contiguos
+* and ordered: thath is xmaxi==xmini+1. The range of validity thus becomes [xmin0,xmaxn).
+* The ranges are checked by the constructor.
+*/
+
+
+
+namespace Interpolator{
+
+class multiPolynomial : public Interpolator
+{
+public:
+
+ // Constructors/Destructors
+ //
+
+
+ /**
+ * Constructor with data
+ * @param polynomials vector of interpolatingPolynomial. must be built outside
+ * @exception std::domain_error: range of polynomial not contiguos or not in order
+ */
+ multiPolynomial ( InterpolatingPolynoamialVector& polynomials );
+
+ /**
+ * Copy constructor
+ * @param src the multiPolynomial to copied
+ */
+ multiPolynomial( multiPolynomial& src);
+
+ /**
+ * Assignment operator
+ * @param src the multiPolynomial to copied
+ */
+ multiPolynomial& operator=(const multiPolynomial& src);
+
+ /**
+ * Destructor
+ */
+ ~multiPolynomial ( );
+
+ /**
+ * Calculate the inetrpolated value y=f(x)
+ * @return interpolated value y
+ * @param in_val input value (x)
+ * @exception lenght_error wrong number of polynomial (>=1)
+ * @exception range_error at least one input data is outside the [min,max] range
+ */
+ double evaluate(double in_val) throw (std::range_error);
+
+ /**
+ * estimate the value of the inverse iterpolator x=g(y)
+ * We make the hypothesys that the interpolated function is invertible!!
+ * @param goal the recquired output (y)
+ * @param initialEstimate initial estimate for x
+ * @return x, NAN if the algorithm did not converge
+ */
+
+ double inverse_evaluate(double goal,double initialEstimate);
+
+ /**
+ * calculate the interpolated value v[]=f(in_val[]). Vector version
+ * @return doubleVector
+ * @param in_val doubleVector of input values
+ * @exception range_error: at least one input data is outside the [min,max] range
+ */
+ doubleVector evaluate (doubleVector in_val ) throw(std::range_error);
+
+protected:
+ double* range; /**< auxiliary struct to search the polynomial to use */
+ unsigned int n_poly; /**< number ofinterpolating polynomials polynomials */
+ std::vector<doubleVector> _poly; /**< the polynomials, stored as vectors of coefficients */
+private:
+
+
+};
+
+} //close namespace Interpolator
+
+#endif // INTERPOLATOR_MULTIPOLYNOMIAL_H
diff --git a/periodicspline.cpp b/periodicspline.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4765ae9417063eab2d38014f1110b517e7319523
--- /dev/null
+++ b/periodicspline.cpp
@@ -0,0 +1,231 @@
+/*
+ * class periodicSpline - implementation
+ *
+ * Algorithm taken from:
+ * Qwt Widget Library
+ * Copyright (C) 1997 Josef Wilgen
+ * Copyright (C) 2002 Uwe Rathmann
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the Qwt License, Version 1.0
+ *
+ * $Author: claudio $
+ *
+ * $Revision: 1.2 $
+ *
+ * $Log: periodicspline.cpp,v $
+ * Revision 1.2 2017-03-06 13:05:54 claudio
+ * use unsigned for comparison
+ *
+ * Revision 1.1.1.1 2012-09-05 08:17:22 claudio
+ * Interpolator utility classes
+ *
+ */
+
+#include "periodicspline.h"
+#include <cmath>
+
+namespace Interpolator{
+
+periodicSpline::periodicSpline (doubleVector& xvalues, doubleVector& yvalues ):Interpolator(){
+ if(xvalues.size() != yvalues.size()) throw std::length_error("periodicSpline(): data vectors must be same length");
+ if(xvalues.size() < 2) throw std::length_error("periodicSpline(): data vector length must be >=2");
+ n_base=xvalues.size();
+ x_base = new double [n_base];
+ y_base= new double [n_base];
+ //from qwt_spline
+ a.resize(n_base-1);
+ b.resize(n_base-1);
+ c.resize(n_base-1);
+
+ //copy the base values into the just allocated arrays
+ //and check that xvalues are sorted correctly
+ double lowbound=-HUGE_VAL;
+ for(unsigned int i=0;i < n_base;i++){
+ if(xvalues[i]>lowbound){
+ lowbound=xvalues[i];
+ }
+ else{
+ delete [] y_base;
+ delete [] x_base;
+ throw std::domain_error("periodicSpline(): x_values must be in ascending order");
+ }
+
+ x_base[i]=xvalues[i];
+ y_base[i]=yvalues[i];
+ }
+ //record range limits
+ _min=x_base[0];
+ _max=x_base[n_base-1];
+
+ {
+ //from qwt_spline build cofficents for periodic spline
+ std::vector<double> d(n_base-1);
+ std::vector<double> h(n_base-1);
+ std::vector<double> s(n_base);
+ //
+ // setup equation system; use coefficient
+ // vectors as temporary buffers
+ //
+ for ( unsigned int i = 0; i < n_base - 1; i++)
+ {
+ h[i] = x_base[i+1] - x_base[i];
+ if (h[i] <= 0.0)
+ throw std::domain_error("periodicSpline(): x_values must be in ascending order");
+ }
+
+ const int imax = n_base - 2;
+ double htmp = h[imax];
+ double dy1 = (y_base[0] - y_base[imax]) / htmp;
+ for (int i = 0; i <= imax; i++)
+ {
+ b[i] = c[i] = h[i];
+ a[i] = 2.0 * (htmp + h[i]);
+ const double dy2 = (y_base[i+1] - y_base[i]) / h[i];
+ d[i] = 6.0 * ( dy1 - dy2);
+ dy1 = dy2;
+ htmp = h[i];
+ }
+
+ //
+ // solve it
+ //
+
+ // L-U Factorization
+ a[0] = sqrt(a[0]);
+ c[0] = h[imax] / a[0];
+ double sum = 0;
+
+ for( int i = 0; i < imax - 1; i++)
+ {
+ b[i] /= a[i];
+ if (i > 0)
+ c[i] = - c[i-1] * b[i-1] / a[i];
+ a[i+1] = sqrt( a[i+1] - (b[i] * b[i]));
+ sum += (c[i] * c[i]);
+ }
+ b[imax-1] = (b[imax-1] - c[imax-2] * b[imax-2]) / a[imax-1];
+ a[imax] = sqrt(a[imax] - (b[imax-1] * b[imax-1]) - sum);
+
+
+ // forward elimination
+ s[0] = d[0] / a[0];
+ sum = 0;
+ for( int i = 1; i < imax; i++)
+ {
+ s[i] = (d[i] - b[i-1] * s[i-1]) / a[i];
+ sum += c[i-1] * s[i-1];
+ }
+ s[imax] = (d[imax] - b[imax-1] * s[imax-1] - sum) / a[imax];
+
+
+ // backward elimination
+ s[imax] = - s[imax] / a[imax];
+ s[imax-1] = -(s[imax-1] + b[imax-1] * s[imax]) / a[imax-1];
+ for (int i= imax - 2; i >= 0; i--)
+ s[i] = - (s[i] + b[i] * s[i+1] + c[i] * s[imax]) / a[i];
+
+ //
+ // Finally, determine the spline coefficients
+ //
+ s[n_base-1] = s[0];
+ for ( unsigned int i=0; i < n_base-1; i++)
+ {
+ a[i] = ( s[i+1] - s[i] ) / ( 6.0 * h[i]);
+ b[i] = 0.5 * s[i];
+ c[i] = ( y_base[i+1] - y_base[i] ) / h[i] - (s[i+1] + 2.0 * s[i] ) * h[i] / 6.0;
+ }
+ }
+
+}
+
+periodicSpline::periodicSpline(periodicSpline& src):Interpolator()
+{
+ if( src.n_base < 2 ) throw std::length_error("periodicSpline(): data vector length must be >=2");
+ n_base=src.n_base;
+ x_base = new double [n_base];
+ y_base= new double [n_base];
+ _min=src._min;
+ _max=src._max;
+ for(unsigned int i=0; i<n_base;i++){
+ x_base[i]=src.x_base[i];
+ y_base[i]=src.y_base[i];
+ }
+ a=src.a;
+ b=src.b;
+ c=src.c;
+}
+
+periodicSpline& periodicSpline::operator=(const periodicSpline& src)
+{
+ if( src.n_base < 2 ) throw std::length_error("periodicSpline(): data vector length must be >=2");
+ if(n_base){
+ delete [] y_base;
+ delete [] x_base;
+ }
+ n_base=src.n_base;
+ x_base = new double [n_base];
+ y_base= new double [n_base];
+ _min=src._min;
+ _max=src._max;
+ for(unsigned int i=0; i<n_base;i++){
+ x_base[i]=src.x_base[i];
+ y_base[i]=src.y_base[i];
+ }
+ a=src.a;
+ b=src.b;
+ c=src.c;
+ return *this;
+}
+
+periodicSpline::~periodicSpline ( )
+{
+ if(n_base){
+ delete [] y_base;
+ delete [] x_base;
+ }
+}
+
+/**
+ * @param xvalues returns the abscissae of the base point of the periodicspline
+ * @param yvalues returns the ordinates of the base point of the periodicspline
+ */
+void periodicSpline::get_base_points (doubleVector& xvalues, doubleVector& yvalues ) const
+{
+ xvalues.clear();
+ yvalues.clear();
+ for(unsigned int i=0;i<n_base;i++){
+ xvalues.push_back(x_base[i]);
+ yvalues.push_back(y_base[i]);
+ }
+ return;
+}
+
+double periodicSpline::evaluate (double v ) throw(std::range_error)
+{
+
+ int klo,khi,k;
+
+ if (v < _min || v > _max ) throw std::range_error("input value out of range");
+
+ klo=0;
+ khi=n_base-1; //era n_base !!!!!
+ while (khi-klo > 1) {
+ k=(khi+klo) >> 1;//k=(khi+klo)/2;
+ if (x_base[k] > v) khi=k;
+ else klo=k;
+ }
+
+ const double delta = v - x_base[klo];
+ return( ( ( ( a[klo] * delta) + b[klo] ) * delta + c[klo] ) * delta + y_base[klo]);
+}
+
+doubleVector periodicSpline::evaluate (doubleVector in_val ) throw(std::range_error)
+{
+ //simplicistic implementation. Should make better use of the indexing
+ doubleVector outval;
+ for(unsigned int i=0; i < in_val.size(); i++) outval.push_back(evaluate(in_val[i]));
+ return outval;
+}
+
+} //close namespace Interpolator
diff --git a/periodicspline.h b/periodicspline.h
new file mode 100644
index 0000000000000000000000000000000000000000..8218419d6928b63540c7ceac67c7a79d20895fd2
--- /dev/null
+++ b/periodicspline.h
@@ -0,0 +1,116 @@
+/* $Author: claudio $
+ *
+ * $Revision: 1.1.1.1 $
+ *
+ * $Log: periodicspline.h,v $
+ * Revision 1.1.1.1 2012-09-05 08:17:22 claudio
+ * Interpolator utility classes
+ *
+ *
+ */
+#ifndef INTERPOLATOR_PERIODICSPLINE_H
+#define INTERPOLATOR_PERIODICSPLINE_H
+
+#include <interpolator.h>
+#include <vector>
+#include <stdexcept>
+
+/**
+ * @class Interpolator::periodicSpline
+ * Class for interpolating tabulated data (fields,currents etc.).
+ * It is used to evalute y=f(x) with x,y belonging to R(1).
+ * It is implemented by means of the cubic periodic spline algortihm.If y values at the exterems are equal, then the
+ * interpolated function will have equal vallues for the 1st and 2nd derivatives at the extremes.
+ * Algorithm taken from:
+ * Qwt Widget Library
+ * Copyright (C) 1997 Josef Wilgen
+ * Copyright (C) 2002 Uwe Rathmann
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the Qwt License, Version 1.0
+ * \see Numerical Recipes, Stoer
+ */
+namespace Interpolator{
+
+class periodicSpline : public Interpolator
+{
+public:
+
+ /**
+ * Constructor with starting data
+ * @param xvalues abscissae of the base point of the periodicspline, must be in ascending order : X[j] > X[i] for j>i .
+ * @param yvalues ordinates of the base point of the periodicspline
+ * @exception std::length_error input data vectors too small or of different sizes
+ */
+ periodicSpline ( doubleVector& xvalues, doubleVector& yvalues);
+
+ /**
+ * Copy constructor
+ * @param src the periodicSpline to be copied
+ */
+
+ periodicSpline( periodicSpline& src);
+
+ /**
+ * Destructor
+ */
+ ~periodicSpline ();
+
+ /**
+ * Oveload copy operator
+ * @param src the periodicSpline to be copied
+ */
+
+ periodicSpline& operator=(const periodicSpline&);
+
+ /**
+ * @param xvalues returns the abscissae of the base point of the periodicSpline
+ * @param yvalues returns the ordinates of the base point of the periodicSpline
+ */
+ void get_base_points (doubleVector& xvalues, doubleVector& yvalues ) const;
+
+ /**
+ * Calculate the ineterpolated value v=f(in_val)
+ * @return interpolated value
+ * @param in_val input value
+ * @exception range_error input value is outside the [min,max] range
+ */
+ double evaluate(double in_val) throw (std::range_error);
+
+ /**
+ * calculate the interpolated value v[]=f(in_val[]). Vector version
+ * @return doubleVector
+ * @param in_val doubleVector of input values
+ * @exception range_error at least one input data is outside the [min,max] range
+ */
+ doubleVector evaluate (doubleVector in_val ) throw(std::range_error);
+
+protected:
+
+
+ // plain double arrays for holding data - More efficient than std::vectors
+ double* x_base; /**< x values of the base points */
+ double* y_base; /**< y values of the base points */
+
+ //number of points of the tabulated data;
+ unsigned int n_base;/**< number of base points */
+
+ /** coefficient vectors */
+ std::vector<double> a;
+ std::vector<double> b;
+ std::vector<double> c;
+
+
+private:
+
+
+ // Private attribute accessor methods
+ //
+
+
+
+};
+
+} // close namespce Interpolator
+
+#endif // INTERPOLATOR_PERIODICSPLINE_H
diff --git a/spline.cpp b/spline.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..624818ff17905461a5e2b6cba4a914b5d3b7028f
--- /dev/null
+++ b/spline.cpp
@@ -0,0 +1,259 @@
+/*
+ * class Spline - implementation
+ *
+ * $Author: claudio $
+ *
+ * $Revision: 1.2 $
+ *
+ * $Log: spline.cpp,v $
+ * Revision 1.2 2019-12-06 13:03:02 claudio
+ * gcc 7.4.0
+ *
+ * Revision 1.1.1.1 2012-09-05 08:17:22 claudio
+ * Interpolator utility classes
+ *
+ */
+
+#include "spline.h"
+#include <cmath>
+#include <cstdlib>
+#include <limits>
+
+namespace Interpolator{
+
+Spline::Spline (doubleVector& xvalues, doubleVector& yvalues, double yp0 , double ypN ):Interpolator(){
+ if(xvalues.size() != yvalues.size()) throw std::length_error("Spline(): data vectors must be same length");
+ if(xvalues.size() < 2) throw std::length_error("Spline(): data vector length must be >=2");
+ n_base=xvalues.size();
+ x_base= static_cast<double*>(malloc(n_base*sizeof(double))); //get one element more the needed - just for testing bug
+ y_base= static_cast<double*>(malloc(n_base*sizeof(double)));
+ y2 = static_cast<double*>(malloc((n_base+1)*sizeof(double)));
+ u = static_cast<double*>(malloc((n_base-1)*sizeof(double)));
+
+ //copy the base values into the just allocated arrays
+ //and check that xvalues are sorted correctly
+ double lowbound=-HUGE_VAL;
+ for(unsigned int i=0;i < n_base;i++){
+ if(xvalues[i]>lowbound){
+ lowbound=xvalues[i];
+ }
+ else{
+ free(u);
+ free(y2);
+ free(y_base);
+ free(x_base);
+ throw std::domain_error("Spline(): x_values must be in ascending order");
+ }
+
+ x_base[i]=xvalues[i];
+ y_base[i]=yvalues[i];
+ }
+ //record range limits
+ _min=x_base[0];
+ _max=x_base[n_base-1];
+ // do the initial calulations for the spline algorithm
+ // got from an old Pascal based text about Numerical Algorithms
+ {
+ double p,sig,qn,un;
+ unsigned int i;
+ int k;
+ //set boundary conditions for intial point
+ if(std::isnan(yp0)){
+ u[0]=0.0; y2[0]=0.0;//natural spline condition
+ }
+ else{
+ u[0]=(3.0/(x_base[1]-x_base[0]))*((y_base[1]-y_base[0])/(x_base[1]-x_base[0])-yp0);
+ y2[0]=-0.5;
+ }
+
+ for (i=1;i<n_base-1;i++) {
+ sig=(x_base[i]-x_base[i-1])/(x_base[i+1]-x_base[i-1]);
+ p=sig*y2[i-1]+2.0;
+ y2[i]=(sig-1.0)/p;
+ u[i]=(y_base[i+1]-y_base[i])/(x_base[i+1]-x_base[i]) - (y_base[i]-y_base[i-1])/(x_base[i]-x_base[i-1]);
+ u[i]=(6.0*u[i]/(x_base[i+1]-x_base[i-1])-sig*u[i-1])/p;
+ }
+ //set boundary conditions for intial point
+ if(std::isnan(ypN)){
+ qn=un=0.0; //natural spline condition
+ }
+ else{
+ un=(3.0/(x_base[n_base-1]-x_base[n_base-2]))* (ypN-(y_base[n_base-1]-y_base[n_base-2])/(x_base[n_base-1]-x_base[n_base-2]));
+ qn=0.5;
+ }
+ y2[n_base-1]=(un-qn*u[n_base-2])/(qn*y2[n_base-2]+1.0);
+ for (k=n_base-2;k>=0;k--)
+ y2[k]=y2[k]*y2[k+1]+u[k];
+ }
+}
+
+Spline::Spline(Spline& src):Interpolator()
+{
+ if( src.n_base < 2 ) throw std::length_error("Spline(): data vector length must be >=2");
+ n_base=src.n_base;
+ x_base= static_cast<double*>(malloc(n_base*sizeof(double)));
+ y_base= static_cast<double*>(malloc(n_base*sizeof(double)));
+ y2 = static_cast<double*>(malloc((n_base+1)*sizeof(double)));
+ u = static_cast<double*>(malloc((n_base-1)*sizeof(double)));
+ _min=src._min;
+ _max=src._max;
+ for(unsigned int i=0; i<n_base;i++){
+ x_base[i]=src.x_base[i];
+ y_base[i]=src.y_base[i];
+ y2[i]=src.y2[i];
+ }
+ y2[n_base]=src.y2[n_base];
+ for(unsigned int i=0; i<n_base-1;i++){
+ u[i]=src.u[i];
+ }
+}
+
+Spline& Spline::operator=(const Spline& src)
+{
+
+ n_base=src.n_base;
+ x_base= static_cast<double*>(malloc(n_base*sizeof(double)));
+ y_base= static_cast<double*>(malloc(n_base*sizeof(double)));
+ y2 = static_cast<double*>(malloc((n_base+1)*sizeof(double)));
+ u = static_cast<double*>(malloc((n_base+1)*sizeof(double)));
+ _min=src._min;
+ _max=src._max;
+ for(unsigned int i=0; i<n_base;i++){
+ x_base[i]=src.x_base[i];
+ y_base[i]=src.y_base[i];
+ y2[i]=src.y2[i];
+ }
+ y2[n_base]=src.y2[n_base];
+ for(unsigned int i=0; i<n_base-1;i++){
+ u[i]=src.u[i];
+ }
+ return *this;
+}
+
+Spline::~Spline ( )
+{
+ if(n_base){
+ free(u);
+ //free(pp);
+ free(y2);
+ free(y_base);
+ free(x_base);
+ }
+}
+
+/**
+ * @param xvalues returns the abscissae of the base point of the spline
+ * @param yvalues returns the ordinates of the base point of the spline
+ */
+void Spline::get_base_points (doubleVector& xvalues, doubleVector& yvalues ) const
+{
+ xvalues.clear();
+ yvalues.clear();
+ for(unsigned int i=0;i<n_base;i++){
+ xvalues.push_back(x_base[i]);
+ yvalues.push_back(y_base[i]);
+ }
+ return;
+}
+
+double Spline::evaluate (double v ) throw(std::range_error)
+{
+ int klo,khi,k;
+ double h,b,a;
+ if (v < _min || v > _max ) throw std::range_error("input value out of range");
+
+ klo=0;
+ khi=n_base-1; //era n_base !!!!!
+ while (khi-klo > 1) {
+ k=(khi+klo) >> 1;//k=(khi+klo)/2;
+ if (x_base[k] > v) khi=k;
+ else klo=k;
+ }
+ h=x_base[khi]-x_base[klo];
+ if (h == 0.0) throw std::range_error("input value out of range");
+ a=(x_base[khi]-v)/h;
+ b=(v-x_base[klo])/h;
+ return a*y_base[klo]+b*y_base[khi]+((a*a*a-a)*y2[klo]+(b*b*b-b)*y2[khi])*(h*h)/6.0;
+}
+
+doubleVector Spline::evaluate (doubleVector in_val ) throw(std::range_error)
+{
+ //simplicistic implementation. Should make better use of the indexing
+ doubleVector outval;
+ for(unsigned int i=0; i < in_val.size(); i++) outval.push_back(evaluate(in_val[i]));
+ return outval;
+}
+
+
+// bool Spline::bracket(double goal,double initial,double tol, double& xlow, double& xhigh)
+// {
+// if(tol<=std::numeric_limits<double>::epsilon()) throw std::range_error("bracket: tol too small");
+// const double fact=1.618;
+// const unsigned int ntry=50;
+// double x1,x2;
+// x1=initial-tol; if(x1 < _min) x1=_min;
+// x2=initial+tol; if(x2 > _max) x2=_max;
+// double f1,f2;
+// f1=goal - evaluate(x1);
+// f2=goal - evaluate(x2);
+// for(unsigned int j=0; j < ntry; j++){
+// if(f1*f2 <0.0){
+// xlow=x1;
+// xhigh=x2;
+// return true;
+// }
+// if(fabs(f1) < fabs(f2)){
+// x1 += fact*(x1-x2);
+// if(x1 < _min) x1=_min;
+// if(x1 > _max) x1=_max;
+// f1=goal - evaluate(x1);
+// }
+// else{
+// x2 += fact *(x2-x1);
+// if(x2 < _min) x2=_min;
+// if(x2 > _max) x2=_max;
+// f2=goal - evaluate(x2);
+// }
+// }
+// return false;
+// }
+//
+//
+// double Spline::solve(double goal,double x1,double x2,double tol,unsigned int maxiter)
+// {
+// double fl,fh,xl,xh,del,f,rtf,dx;
+// fl = goal - evaluate(x1);
+// fh = goal - evaluate(x2);
+// if( fl*fh > 0.0) throw std::range_error("solve: initial interval does not braket solution");
+// if(fl < 0.0 ){
+// xl=x1;
+// xh=x2;
+// }
+// else{
+// xl=x2;
+// xh=x1;
+// double swap=fl;
+// fl=fh;
+// fh=swap;
+// }
+// dx=xh-xl;
+// for(unsigned int j=0; j < maxiter; j++){
+// rtf=xl+dx*fl/(fl-fh);
+// f= goal - evaluate(rtf);
+// if(f < 0.0){
+// del=xl-rtf;
+// xl=rtf;
+// fl=f;
+// }
+// else{
+// del=xh-rtf;
+// xh=rtf;
+// fh=f;
+// }
+// dx=xh-xl;
+// if(fabs(del) < tol || f == 0.0) return rtf;
+// }
+// return NAN; //maxiter exceeded
+// }
+
+} //close namespace Interpolator
diff --git a/spline.h b/spline.h
new file mode 100644
index 0000000000000000000000000000000000000000..34e28dfd0be83434ed379d1aef2f880371eb7c40
--- /dev/null
+++ b/spline.h
@@ -0,0 +1,110 @@
+/* $Author: claudio $
+ *
+ * $Revision: 1.1.1.1 $
+ *
+ * $Log: spline.h,v $
+ * Revision 1.1.1.1 2012-09-05 08:17:22 claudio
+ * Interpolator utility classes
+ *
+ *
+ *
+ */
+#ifndef INTERPOLATOR_SPLINE_H
+#define INTERPOLATOR_SPLINE_H
+
+#include <interpolator.h>
+
+/**
+ * @class Interpolator::Spline
+ * Class for interpolating tabulated data (fields,currents etc.).
+ * It is used to evalute y=f(x) with x,y belonging to R(1).
+ * It is implemented by means of the cubic spline algortihm.
+ * By default it builds a "natural spline" Interpolator with special conditions at the extreme: zero 2nd derivative.
+ * It is also possible to specify the values of the 1st derivative at both extremes.
+ * \see Numerical Recipes, Stoer
+ */
+
+
+namespace Interpolator{
+
+class Spline : public Interpolator
+{
+public:
+
+ /**
+ * Constructor with starting data
+ * @param xvalues abscissae of the base point of the spline, must be in ascending order : X[j] > Y[i] for j>i .
+ * @param yvalues ordinates of the base point of the spline
+ * @param yp0 1st derivate value at the start of interval. Use NAN for free spline condition
+ * @param ypN 1st derivate value at the end of interval. Use NAN for free spline condition
+ * @exception std::length_error input data vectors too small or of different sizes
+ */
+ Spline ( doubleVector& xvalues, doubleVector& yvalues, double yp0=NAN , double ypN=NAN);
+
+ /**
+ * Copy constructor
+ * @param src the Spline to be copied
+ */
+
+ Spline( Spline& src);
+
+ /**
+ * Destructor
+ */
+ ~Spline ();
+
+ /**
+ * Oveload copy operator
+ * @param src the Spline to be copied
+ */
+
+ Spline& operator=(const Spline&);
+
+ /**
+ * @param xvalues returns the abscissae of the base point of the Spline
+ * @param yvalues returns the ordinates of the base point of the Spline
+ */
+ void get_base_points (doubleVector& xvalues, doubleVector& yvalues ) const;
+
+ /**
+ * Calculate the ineterpolated value v=f(in_val)
+ * @return interpolated value
+ * @param in_val input value
+ * @exception range_error input value is outside the [min,max] range
+ */
+ double evaluate(double in_val) throw (std::range_error);
+
+ /**
+ * calculate the interpolated value v[]=f(in_val[]). Vector version
+ * @return doubleVector
+ * @param in_val doubleVector of input values
+ * @exception range_error at least one input data is outside the [min,max] range
+ */
+ doubleVector evaluate (doubleVector in_val ) throw(std::range_error);
+
+protected:
+
+
+ // plain double arrays for holding data - More efficient than std::vectors
+ double* x_base; /**< x values of the base points */
+ double* y_base; /**< y values of the base points */
+
+ //array for storing precalculated data spline parameters. Actually it contains the 2nd derivative of the tabulated function
+ double* y2; /**< array for storing precalculated data spline parameters */
+ double* u; /**< array for storing precalculated data spline parameters */
+ //number of points of the tabulated data;
+ unsigned int n_base;/**< number of base points */
+
+private:
+
+
+ // Private attribute accessor methods
+ //
+
+
+
+};
+
+} // close namespce Interpolator
+
+#endif // INTERPOLATOR_SPLINE_H
diff --git a/spltest.cpp b/spltest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3bd43a21d7b03442e7a865455c7cb2303f9c6b03
--- /dev/null
+++ b/spltest.cpp
@@ -0,0 +1,138 @@
+
+
+/* $Author: claudio $
+ *
+ * $Revision: 1.1.1.1 $
+ *
+ * $Log: spltest.cpp,v $
+ * Revision 1.1.1.1 2012-09-05 08:17:22 claudio
+ * Interpolator utility classes
+ *
+ */
+
+#include <cstdlib>
+#include <iostream>
+#include <iomanip>
+#include "spline.h"
+
+
+using namespace std;
+using namespace Interpolator;
+
+int main(){
+
+ try{
+ double i_qd[]={0.0,14.8056,74.028,148.056,222.084,296.112,333.126,370.14};
+ double g_qd[]={0.0,0.7745,3.9519,7.9231,11.8187,15.7069,17.5649,19.2667};
+// doubleVector iqd;
+// doubleVector fieldqd;
+ vector<double> iqd;
+ vector<double> fieldqd;
+
+ //controle printout of table
+ int PRECISION=16;
+ int WIDTH=20;
+
+ for (unsigned int i=0;i<8;i++){
+ iqd.push_back(i_qd[i]);
+ fieldqd.push_back(g_qd[i]);
+
+ cout << setprecision(PRECISION) << setw(WIDTH) << right;
+ cout << iqd[i];
+
+ cout << setprecision(PRECISION) << setw(WIDTH) << right;
+ cout << fieldqd[i] << endl;
+ }
+ Spline* dir= new Spline(iqd,fieldqd);
+ Spline* inv= new Spline(fieldqd,iqd);
+
+
+ double min=0.0,max=i_qd[7];
+ unsigned int steps=200;
+ double mincur=min;
+ double inccur=(max-min)/steps;
+ double I,g,inv_I;
+ cout << setw(WIDTH) << right;
+ cout << "I";
+
+ cout << setw(WIDTH) << right;
+ cout << "g" ;
+
+ cout << setw(WIDTH) << right;
+ cout << "inv_I" ;
+
+
+
+ cout << endl;
+ for(unsigned int i=0;i<=steps;i++){
+ I=mincur+ i*inccur;
+ g=dir->evaluate(I);
+ inv_I=inv->evaluate(g);
+
+
+
+ cout << setprecision(PRECISION) << setw(WIDTH) << right;
+ cout << I;
+
+ cout << setprecision(PRECISION) << setw(WIDTH) << right;
+ cout << g;
+
+ cout << setprecision(PRECISION) << setw(WIDTH) << right;
+ cout << inv_I;
+
+
+ if(I != 0.0) {
+ cout << setprecision(PRECISION) << setw(WIDTH) << right;
+ cout << 10000.0*((I-inv_I)/I);
+ }
+ else{
+ cout << setprecision(PRECISION) << setw(WIDTH) << right;
+ cout << 0.0;
+ }
+
+// if(dir->bracket(g,I,1.0e-4,xmin,xmax)){
+// cout << setprecision(PRECISION) << setw(WIDTH) << right;
+// cout << xmin;
+//
+// cout << setprecision(PRECISION) << setw(WIDTH) << right;
+// cout << xmax;
+//
+// double I2=dir->solve(g,xmin,xmax,1.0e-8,100);
+//
+// cout << setprecision(PRECISION) << setw(WIDTH) << right;
+// cout << I2;
+//
+// }
+// else{
+// cout << setprecision(PRECISION) << setw(WIDTH) << right;
+// cout << "------";
+//
+// cout << setprecision(PRECISION) << setw(WIDTH) << right;
+// cout << "++++++";
+//
+// cout << setprecision(PRECISION) << setw(WIDTH) << right;
+// cout << "//////";
+// }
+ {
+
+
+ double I2=dir->inverse_evaluate(g,I);
+
+ cout << setprecision(PRECISION) << setw(WIDTH) << right;
+ cout << I2;
+
+ }
+
+
+ cout << endl;
+ }
+
+ delete inv;
+ delete dir;
+ }
+ catch(std::invalid_argument& ex)
+ {
+ cout<<ex.what()<<endl;
+ }
+ return 0;
+}
diff --git a/test_freespline.cpp b/test_freespline.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..185684042fa571342cf6806b1a2bc5033f8f294b
--- /dev/null
+++ b/test_freespline.cpp
@@ -0,0 +1,263 @@
+/*
+ * program to test spline class
+ *
+ * $Author: claudio $
+ *
+ * $Revision: 1.1.1.1 $
+ *
+ * $Log: test_freespline.cpp,v $
+ * Revision 1.1.1.1 2012-09-05 08:17:22 claudio
+ * Interpolator utility classes
+ *
+ *
+ */
+
+#include <iostream>
+#include <fstream>
+#include "spline.h"
+#include <cmath>
+
+
+using namespace std;
+using namespace Interpolator;
+
+int main(){
+ cout<< "test of Spline interpolator class" << endl;
+ doubleVector x0,y0;
+ const int NP=30;
+ for(int i=0; i < NP; i++){
+ x0.push_back(i*i*.02);
+ y0.push_back(sin(i*i*.02));
+ }
+ Spline* test_spline = new Spline(x0,y0,NAN,0.0);
+ {
+ doubleVector a;
+ vector<double> b;
+ test_spline->get_base_points(a,b);
+ cout<<"test get_base_points()"<<endl;
+ for(int i=0; i < NP; i++ )
+ if((x0[i] != a[i]) ||(y0[i] != b[i])){
+ cout << i << " " << x0[i] << " " << a[i] << " | "<< y0[i] << " " << b[i] << endl;
+ return -1;
+ }
+
+ }
+ {
+ cout << "test get_range()" << endl;
+ double minx,maxx;
+ test_spline->get_range(minx,maxx);
+ if( (minx != x0[0]) || (maxx != x0[NP-1]) ){
+ cout << minx << " " << x0[0] << " | " << maxx << " " << x0[NP-1] << endl;
+ return -1;
+ }
+ }
+ {
+ cout << "test evaluate()" << endl;
+ double xn=0.02;
+ double yn=test_spline->evaluate(xn);
+ cout << xn << " " << yn << endl;
+ yn=test_spline->evaluate(x0[5]);
+ cout << xn << " " << yn << " " << y0[5] << endl;
+ }
+ {
+ cout << "test evaluate()" << endl;
+ double xn=0.02;
+ double yn=test_spline->evaluate(xn);
+ cout << xn << " " << yn << endl;
+ yn=test_spline->evaluate(x0[5]);
+ cout << xn << " " << yn << " " << y0[5] << endl;
+ }
+ { cout << "test evaluate()for plot" << endl;
+ doubleVector X;
+ doubleVector Y;
+ double range=x0[NP-1];
+ double nsteps=300;
+ for(int i=0;i<=nsteps;i++){
+ X.push_back((range*i)/nsteps);
+ Y.push_back(test_spline->evaluate(X[i]));
+ }
+ //print out data
+ cout<<"---------------------------"<<endl;
+ for(int i=0; i < NP; i++) cout<<x0[i]<<"\t"<<y0[i]<<endl;
+ cout<<"---------------------------"<<endl;
+ for(int i=0; i < nsteps; i++) cout<<X[i]<<"\t"<<Y[i]<<endl;
+ cout<<"+++++++++++++++++++++++++++"<<endl;
+ }
+ {
+ cout << "test evaluate()for doubleVector" << endl;
+ doubleVector X;
+ doubleVector Y;
+ double range=x0[NP-1];
+ double nsteps=300;
+ for(int i=0;i<=nsteps;i++)
+ X.push_back((range*i)/nsteps);
+ Y=test_spline->evaluate(X);
+ //print out data
+ cout<<"###########################"<<endl;
+ for(int i=0; i < nsteps; i++) cout<<X[i]<<"\t"<<Y[i]<<endl;
+ cout<<"***************************"<<endl;
+ }
+ {
+ cout << "test range exception" << endl;
+ try{
+ test_spline->evaluate(-10.0);
+ cerr << "range_error not thrown!" << endl;
+ }
+ catch(std::range_error& ex){
+ cout << "OK, catched expected exception: "<< ex.what() << endl;
+ }
+ catch(...){
+ cerr << "unexpected exception thrown!" << endl;
+ }
+ }
+ {
+ cout << "test constructor with unequal lengths" << endl;
+ try{
+ doubleVector x,y;
+ x.push_back(1.0);
+ x.push_back(2.0);
+ x.push_back(3.0);
+
+ y.push_back(4.0);
+ y.push_back(5.0);
+ Spline spl(x,y,NAN,0.0);
+ cerr << "length_error not thrown!" << endl;
+ }
+ catch(std::length_error& ex){
+ cout << "OK, catched expected exception: "<< ex.what() << endl;
+ }
+ catch(...){
+ cerr << "unexpected exception thrown!" << endl;
+ }
+ }
+ {
+ cout << "test constructor with out of order data" << endl;
+ try{
+ doubleVector x,y;
+ x.push_back(1.0);
+ x.push_back(2.0);
+ x.push_back(1.5);
+
+ y.push_back(4.0);
+ y.push_back(5.0);
+ y.push_back(6.0);
+ Spline spl(x,y,NAN,0.0);
+ cerr << "domain_error not thrown!" << endl;
+ }
+ catch(std::domain_error& ex){
+ cout << "OK, catched expected exception: "<< ex.what() << endl;
+ }
+ catch(...){
+ cerr << "unexpected exception thrown!" << endl;
+ }
+ }
+ {
+ cout << "Test copy constructor" << endl;
+ Spline copy_spline(*test_spline);
+ double minx,maxx,copyminx,copymaxx;
+ test_spline->get_range(minx,maxx);
+ copy_spline.get_range(copyminx,copymaxx);
+ if( (minx != copyminx) || (maxx != copymaxx) ){
+ cerr << "COPY FAILED!" << endl;
+ return -1;
+ }
+ double range=x0[NP-1];
+ double nsteps=100;
+ for(int i=0;i<nsteps;i++){
+ double val=(range*i)/nsteps;
+ double y,ycopy;
+ ycopy=copy_spline.evaluate(val);
+ y=test_spline->evaluate(val);
+ if(y != ycopy){
+ cerr << "COPY FAILED " << val << " " << y << " " << ycopy << endl;
+ }
+ }
+ }
+ {
+ cout << "Test copy assignment operator" << endl;
+ Spline oper_spline=*test_spline;
+ double minx,maxx,operminx,opermaxx;
+ test_spline->get_range(minx,maxx);
+ oper_spline.get_range(operminx,opermaxx);
+ if( (minx != operminx) || (maxx != opermaxx) ){
+ cerr << "ASSIGN COPY FAILED!" << endl;
+ return -1;
+ }
+ double range=x0[NP-1];
+ double nsteps=100;
+ for(int i=0;i<nsteps;i++){
+ double val=(range*i)/nsteps;
+ double y,yoper;
+ yoper=oper_spline.evaluate(val);
+ y=test_spline->evaluate(val);
+ if(y != yoper){
+ cerr << "COPY FAILED " << val << " " << y << " " << yoper << endl;
+ }
+ }
+ }
+ cout << "test delete" << endl;
+ delete test_spline;
+ {
+ cout << "test leaks - new/delete" << endl;
+ Spline *dynspl;
+ doubleVector x,y;
+ x.push_back(1.0);
+ x.push_back(2.0);
+ x.push_back(3.5);
+
+ y.push_back(4.0);
+ y.push_back(5.0);
+ y.push_back(6.0);
+ for(int i=0; i < 100 ; i++){
+ dynspl = new Spline(x,y,NAN,0.0);
+ delete dynspl;
+ }
+ }
+ {
+ cout << "test leaks - stack" << endl;
+ doubleVector x,y;
+ x.push_back(1.0);
+ x.push_back(2.0);
+ x.push_back(3.5);
+
+ y.push_back(4.0);
+ y.push_back(5.0);
+ y.push_back(6.0);
+ for(int i=0; i < 100 ; i++){
+ Spline stspl(x,y,NAN,0.0);
+ }
+ }
+
+ {
+ cout <<"test end conditions" << endl;
+ std::ofstream outfile;
+ outfile.open("splineval.csv",ios::out); //CHECK success
+ if ( (outfile.rdstate() & ifstream::failbit ) != 0 ) return -1;
+ doubleVector x,y;
+ x.push_back(0.0);y.push_back(0.0);
+ x.push_back(1.0);y.push_back(1.0);
+ x.push_back(1.5);y.push_back(1.0);
+ x.push_back(2.0);y.push_back(1.0);
+ x.push_back(3.0);y.push_back(0.0);
+ Spline nsp(x,y); //natural spline
+ Spline zsp(x,y,0.0,0.0); // spline with 0 derivate start and stop
+ Spline nzsp(x,y,NAN,0.0); // spline with 0 derivate at stop
+ Spline znsp(x,y,0.0,NAN); // spline with 0 derivate at start
+ Spline ddsp(x,y,-0.5,0.5); // spline with derivatives at start and stop
+ double X;
+ double Yn,Yz,Ynz,Yzn,Ydd;
+ double range=3.0;
+ int nsteps=300;
+ for(int i=0;i<=nsteps;i++){
+ X=((range*i)/nsteps);
+ Yn=nsp.evaluate(X);
+ Yz=zsp.evaluate(X);
+ Ynz=nzsp.evaluate(X);
+ Yzn=znsp.evaluate(X);
+ Ydd=ddsp.evaluate(X);
+ outfile<< X << "\t\t" << Yn << "\t\t" << Yz << "\t\t" << Ynz << "\t\t" << Yzn << "\t\t" << Ydd << std::endl;
+ }
+ outfile.close();
+ }
+ return 0;
+}
diff --git a/test_inverse.cpp b/test_inverse.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0fa48ce0b5f087194e339875280f590d7587b1b5
--- /dev/null
+++ b/test_inverse.cpp
@@ -0,0 +1,67 @@
+/*
+ * program to test spline classes
+ *
+ * $Author: claudio $
+ *
+ * $Revision: 1.1 $
+ *
+ * $Log: test_inverse.cpp,v $
+ * Revision 1.1 2017-03-03 15:32:40 claudio
+ * correctetd spline::inverse_calculate if returned value is at range limit
+ *
+ * Revision 1.1.1.1 2012-09-05 08:17:22 claudio
+ * Interpolator utility classes
+ *
+ *
+ */
+
+#include <iostream>
+#include "spline.h"
+#include <cmath>
+//#include <stdexcept>
+
+using namespace std;
+using namespace Interpolator;
+int main(){
+ cout << "test spline::inverse_evaluate" << endl;
+ {
+
+ doubleVector x,y;
+ x.push_back(10.0); y.push_back(83.6);
+ x.push_back(12); y.push_back(68.3);
+ x.push_back( 14); y.push_back( 55.7);
+ x.push_back( 16); y.push_back(45.5);
+ x.push_back( 18); y.push_back(37.3);
+ x.push_back( 20); y.push_back(30.7);
+ x.push_back( 22); y.push_back( 25.4);
+ x.push_back( 24); y.push_back( 21.2);
+ x.push_back( 26); y.push_back( 17.9);
+ x.push_back( 28); y.push_back( 15.2);
+ x.push_back( 30); y.push_back( 13.1);
+ x.push_back( 32); y.push_back( 11.4);
+ x.push_back( 34); y.push_back( 10.1);
+ x.push_back( 36); y.push_back( 9);
+ x.push_back( 40); y.push_back( 7.5);
+ x.push_back( 50); y.push_back( 5.77);
+ Spline *lambdasplie = new Spline(x,y);
+ double gmin,gmax;
+
+ lambdasplie->get_range(gmin,gmax);
+ double max = lambdasplie->evaluate(gmin);
+ double min = lambdasplie->evaluate(gmax);
+ double initial=min +((max-min)/2);
+
+ double l=min;
+ double gap, linv;
+ cout << "l\tgap\tlinverse" << endl;
+ cout.precision(8);
+ while (l<max){
+ gap = lambdasplie->inverse_evaluate(l,initial);
+ linv = lambdasplie->evaluate(gap);
+ cout << l <<"\t" << gap << "\t" << linv <<endl;
+ l+=0.05;
+ }
+ delete lambdasplie;
+ }
+ return 0;
+}
diff --git a/test_multipoly.cpp b/test_multipoly.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d67022b86af83f73e0b22dbc1311184d72a77f45
--- /dev/null
+++ b/test_multipoly.cpp
@@ -0,0 +1,282 @@
+/*
+ * program to test multipolynomial class
+ *
+ * $Author: claudio $
+ *
+ * $Revision: 1.1.1.1 $
+ *
+ * $Log: test_multipoly.cpp,v $
+ * Revision 1.1.1.1 2012-09-05 08:17:22 claudio
+ * Interpolator utility classes
+ *
+ */
+
+#include <iostream>
+#include "multipolynomial.h"
+
+using namespace std;
+using namespace Interpolator;
+
+int main(){
+ cout << "test of multiPolynomial interpolator class" << endl;
+ try{
+ //create interpolating polynomials
+
+ InterpolatingPolynomial p0,p1,p2;
+ p0.xMin=0.5;
+ p0.xMax=1.5;
+ p0.coefficient.push_back(1.1);
+
+ p1.xMin=1.5;
+ p1.xMax=2.5;
+ p1.coefficient.push_back(1.1);
+ p1.coefficient.push_back(2.2);
+
+ p2.xMin=2.5;
+ p2.xMax=3.5;
+ p2.coefficient.push_back(1.1);
+ p2.coefficient.push_back(2.2);
+ p2.coefficient.push_back(3.3);
+
+ //stuff polynomial into vector
+
+ InterpolatingPolynoamialVector pvector;
+ pvector.push_back(p0);
+ pvector.push_back(p1);
+ pvector.push_back(p2);
+
+ multiPolynomial *multipoly = new multiPolynomial(pvector);
+ double minX,maxX;
+ multipoly->get_range(minX,maxX);
+ if( (minX != p0.xMin) || ( maxX != p2.xMax) ){
+ cout << "multiPoly.get_range(): error " << minX << " " << maxX << " " << p0.xMin << " " << p2.xMax << endl;
+ return -1; //stop further testing
+ }
+ delete multipoly;
+ }
+ catch(...){
+ cerr << " unexpected excpetion!" << endl;
+ return -1;
+ }
+
+ {
+ cout << "test a single 2nd degree polynomial " << endl;
+ try{
+ InterpolatingPolynomial p0;
+ p0.xMin=0.0;
+ p0.xMax=4.0;
+ p0.coefficient.push_back(0.0); //a0
+ p0.coefficient.push_back(0.0); //a1
+ p0.coefficient.push_back(1.0); //a2 --> Y=X*X
+ InterpolatingPolynoamialVector pvector;
+ pvector.push_back(p0);
+ multiPolynomial *multipoly = new multiPolynomial(pvector);
+ double v=0.0,y;
+ y=multipoly->evaluate(v); if( y != 0.0 ){ cerr << "value not correct (0.0): " << y << endl; return -1; }
+ v=1.0;
+ y=multipoly->evaluate(v); if( y != 1.0 ){ cerr << "value not correct (1.0): " << y << endl; return -1; }
+ v=4.0;
+ y=multipoly->evaluate(v); if( y != 16.0 ){ cerr << "value not correct (16.0): " << y << endl; return -1; }
+ doubleVector X;
+ for (int i=0; i<=400;i++) X.push_back(i/100.0); //create range for testing
+ doubleVector Y=multipoly->evaluate(X);
+ cout << "----------------------------------" << endl;
+ for (unsigned int i=0; i < X.size(); i++) cout << X[i] << "\t" << Y[i] << endl;
+ cout << "----------------------------------" << endl;
+ delete multipoly;
+ }
+ catch(std::range_error& ex){
+ std::cerr << ex.what() << std::endl;
+ return -1;
+ }
+ catch(std::length_error& ex){
+ std::cerr << ex.what() << std::endl;
+ return -1;
+ }
+ try{
+ //create interpolating polynomials
+ cout << "test non contiguos" << endl;
+ InterpolatingPolynomial p0,p1,p2;
+ p0.xMin=0.5;
+ p0.xMax=1.5;
+ p0.coefficient.push_back(1.1);
+
+ p1.xMin=1.6; // <-------
+ p1.xMax=2.5;
+ p1.coefficient.push_back(1.1);
+ p1.coefficient.push_back(2.2);
+
+ p2.xMin=2.5;
+ p2.xMax=3.5;
+ p2.coefficient.push_back(1.1);
+ p2.coefficient.push_back(2.2);
+ p2.coefficient.push_back(3.3);
+
+ //stuff polynomial into vector
+
+ InterpolatingPolynoamialVector pvector;
+ pvector.push_back(p0);
+ pvector.push_back(p1);
+ pvector.push_back(p2);
+
+ multiPolynomial *multipoly = new multiPolynomial(pvector);
+ delete multipoly;
+ }
+ catch(std::domain_error& ex){
+ cout << "OK expected exception: " << ex.what() << endl;
+
+ }
+ catch(...){
+ cerr << " unexpected excpetion!" << endl;
+ return -1;
+ }
+
+ try{
+ //create interpolating polynomials
+ cout << "test 3 polynomials" << endl;
+ InterpolatingPolynomial p0,p1,p2;
+ p0.xMin=0.0;
+ p0.xMax=1.5;
+ p0.coefficient.push_back(1.0); //c[0]
+
+ p1.xMin=1.5;
+ p1.xMax=2.5;
+ p1.coefficient.push_back(0.0); //c[0]
+ p1.coefficient.push_back(0.0); //c[1]
+ p1.coefficient.push_back(1.0); //c[2]
+
+ p2.xMin=2.5;
+ p2.xMax=3.5;
+ p2.coefficient.push_back(5.0); //c[0]
+ p2.coefficient.push_back(-2.0); //c[1]
+ //stuff polynomial into vector
+
+ InterpolatingPolynoamialVector pvector;
+ pvector.push_back(p0);
+ pvector.push_back(p1);
+ pvector.push_back(p2);
+
+ multiPolynomial *multipoly = new multiPolynomial(pvector);
+ doubleVector X;
+ for (int i=0; i<=350;i++) X.push_back(i/100.0); //create range for testing
+ doubleVector Y=multipoly->evaluate(X);
+ cout << "-+++++++++++++++++++++++++++++++++" << endl;
+ for (unsigned int i=0; i < X.size(); i++) cout << X[i] << "\t" << Y[i] << endl;
+ cout << "++++++++++++++++++++++++++++++++++" << endl;
+ delete multipoly;
+ }
+ catch(std::domain_error& ex){
+ cout << "unexpected exception: " << ex.what() << endl;
+ return -1;
+ }
+ catch(std::range_error& ex){
+ std::cerr << ex.what() << std::endl;
+ return -1;
+ }
+ catch(std::length_error& ex){
+ std::cerr << ex.what() << std::endl;
+ return -1;
+ }
+ catch(...){
+ cerr << " unexpected excpetion!" << endl;
+ return -1;
+ }
+
+ try{
+ //create interpolating polynomials
+ cout << "test 3 identical polynomials" << endl;
+ InterpolatingPolynomial p0,p1,p2;
+ p0.xMin=0.0;
+ p0.xMax=1.5;
+ p0.coefficient.push_back(1.0); //c[0]
+
+
+ InterpolatingPolynoamialVector pvector;
+ pvector.push_back(p0);
+ pvector.push_back(p0);
+ pvector.push_back(p0);
+
+ multiPolynomial *multipoly = new multiPolynomial(pvector);
+ doubleVector X;
+ for (int i=0; i<=150;i++) X.push_back(i/100.0); //create range for testing
+ doubleVector Y=multipoly->evaluate(X);
+ cout << "**********************************" << endl;
+ for (unsigned int i=0; i < X.size(); i++) cout << X[i] << "\t" << Y[i] << endl;
+ cout << "**********************************" << endl;
+ delete multipoly;
+ }
+ catch(std::domain_error& ex){
+ cout << "OK, exception caught: " << ex.what() << endl;
+ }
+ catch(std::range_error& ex){
+ std::cerr << ex.what() << std::endl;
+ return -1;
+ }
+ catch(std::length_error& ex){
+ std::cerr << ex.what() << std::endl;
+ return -1;
+ }
+ catch(...){
+ cerr << " unexpected excpetion!" << endl;
+ return -1;
+ }
+ try{
+ //create interpolating polynomials
+ cout << "test copy constructor" << endl;
+ InterpolatingPolynomial p0,p1,p2;
+ p0.xMin=0.0;
+ p0.xMax=1.5;
+ p0.coefficient.push_back(1.0); //c[0]
+
+ p1.xMin=1.5;
+ p1.xMax=2.5;
+ p1.coefficient.push_back(0.0); //c[0]
+ p1.coefficient.push_back(0.0); //c[1]
+ p1.coefficient.push_back(1.0); //c[2]
+
+ p2.xMin=2.5;
+ p2.xMax=3.5;
+ p2.coefficient.push_back(5.0); //c[0]
+ p2.coefficient.push_back(-2.0); //c[1]
+ //stuff polynomial into vector
+
+ InterpolatingPolynoamialVector pvector;
+ pvector.push_back(p0);
+ pvector.push_back(p1);
+ pvector.push_back(p2);
+
+ multiPolynomial *multipoly = new multiPolynomial(pvector);
+ multiPolynomial copypoly=*multipoly;
+ doubleVector X;
+ for (int i=0; i<=350;i++) X.push_back(i/100.0); //create range for testing
+ doubleVector Y=multipoly->evaluate(X);
+ doubleVector Yc=copypoly.evaluate(X);
+
+ for (unsigned int i=0; i < X.size(); i++){
+ if (Y[i] != Yc[i])
+ cerr << "COPY ERROR: " << X[i] << " " << Y[i] << " " << Yc[i] <<endl;
+ }
+ delete multipoly;
+ }
+ catch(...){
+ cerr << " test copy constructor: unexpected excpetion!" << endl;
+ return -1;
+ }
+ }
+ {
+ cout << "test interpolatore for constat value on simmetric range" << endl;
+ InterpolatingPolynomial p0;
+ InterpolatingPolynoamialVector pvector;
+ p0.xMin=-5.0;
+ p0.xMax=5.0;
+ p0.coefficient.push_back(0.52); //c[0]
+ pvector.push_back(p0);
+ multiPolynomial *multipoly = new multiPolynomial(pvector);
+ double a0=-4.99;
+ cout << "for "<< a0 << " "<< multipoly->evaluate(a0);
+ a0=4.99;
+ cout << "for "<< a0 << " "<< multipoly->evaluate(a0);
+ delete multipoly;
+ }
+ return 0;
+}
diff --git a/test_periodicspline.cpp b/test_periodicspline.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2d8af4e338944262aaeae7d968faba30e40947b1
--- /dev/null
+++ b/test_periodicspline.cpp
@@ -0,0 +1,242 @@
+/*
+ * program to test peridic spline class
+ *
+ * $Author: claudio $
+ *
+ * $Revision: 1.2 $
+ *
+ * $Log: test_periodicspline.cpp,v $
+ * Revision 1.2 2017-03-06 13:13:20 claudio
+ * use unsigned for comparison
+ *
+ * Revision 1.1.1.1 2012-09-05 08:17:22 claudio
+ * Interpolator utility classes
+ *
+ *
+ */
+
+#include <iostream>
+#include "periodicspline.h"
+#include <cmath>
+
+using namespace std;
+using namespace Interpolator;
+
+int main(){
+ cout<< "test of periodicSpline interpolator class" << endl;
+ doubleVector x0,y0;
+ const int NP=15;
+ const double k = (2*M_PI)/NP;
+ for(int i=0; i <= NP; i++){
+ x0.push_back(i*k);
+ y0.push_back(sin(i*k));
+ }
+ periodicSpline* test_spline = new periodicSpline(x0,y0);
+ {
+ doubleVector a;
+ doubleVector b;
+ test_spline->get_base_points(a,b);
+ cout<<"test get_base_points()"<<endl;
+ for(unsigned int i=0; i < a.size()-1; i++ )
+ if((x0[i] != a[i]) ||(y0[i] != b[i])){
+ cout << i << " " << x0[i] << " " << a[i] << " | "<< y0[i] << " " << b[i] << endl;
+ return -1;
+ }
+
+ }
+ {
+ cout << "test get_range()" << endl;
+ double minx,maxx;
+ test_spline->get_range(minx,maxx);
+ int I=x0.size()-1;
+ if( (minx != x0[0]) || (maxx != x0[I]) ){
+ cout << minx << " " << x0[0] << " | " << maxx << " " << x0[I] << endl;
+ return -1;
+ }
+ }
+ {
+ cout << "test evaluate()" << endl;
+ double xn=0.02;
+ double yn=test_spline->evaluate(xn);
+ cout << xn << " " << yn << endl;
+ yn=test_spline->evaluate(x0[5]);
+ cout << xn << " " << yn << " " << y0[5] << endl;
+ }
+ {
+ cout << "test evaluate()" << endl;
+ double xn=0.02;
+ double yn=test_spline->evaluate(xn);
+ cout << xn << " " << yn << endl;
+ yn=test_spline->evaluate(x0[5]);
+ cout << xn << " " << yn << " " << y0[5] << endl;
+ }
+ { cout << "test evaluate()for plot" << endl;
+ doubleVector X;
+ doubleVector Y;
+ double minx;
+ double range;
+ test_spline->get_range(minx,range);
+ double nsteps=100;
+ for(int i=0;i<=nsteps;i++){
+ X.push_back((range*i)/nsteps);
+ Y.push_back(test_spline->evaluate(X[i]));
+ }
+ //print out data
+ cout<<"---------------------------"<<endl;
+ for(unsigned int i=0; i < x0.size(); i++) cout<<x0[i]<<"\t"<<y0[i]<<endl;
+ cout<<"---------------------------"<<endl;
+ for(int i=0; i <= nsteps; i++) cout<<X[i]<<"\t"<<Y[i]<<endl;
+ cout<<"+++++++++++++++++++++++++++"<<endl;
+ }
+ {
+ cout << "test evaluate()for doubleVector" << endl;
+ doubleVector X;
+ doubleVector Y;
+ double minx;
+ double range;
+ test_spline->get_range(minx,range);
+ double nsteps=100;
+ for(int i=0;i<=nsteps;i++)
+ X.push_back((range*i)/nsteps);
+ Y=test_spline->evaluate(X);
+ //print out data
+ cout<<"###########################"<<endl;
+ for(int i=0; i <= nsteps; i++) cout<<X[i]<<"\t"<<Y[i]<<endl;
+ cout<<"***************************"<<endl;
+ }
+ {
+ cout << "test range exception" << endl;
+ try{
+ test_spline->evaluate(-10.0);
+ cerr << "range_error not thrown!" << endl;
+ }
+ catch(std::range_error& ex){
+ cout << "OK, catched expected exception: "<< ex.what() << endl;
+ }
+ catch(...){
+ cerr << "unexpected exception thrown!" << endl;
+ }
+ }
+ {
+ cout << "test constructor with unequal lengths" << endl;
+ try{
+ doubleVector x,y;
+ x.push_back(1.0);
+ x.push_back(2.0);
+ x.push_back(3.0);
+
+ y.push_back(4.0);
+ y.push_back(5.0);
+ periodicSpline spl(x,y);
+ cerr << "length_error not thrown!" << endl;
+ }
+ catch(std::length_error& ex){
+ cout << "OK, catched expected exception: "<< ex.what() << endl;
+ }
+ catch(...){
+ cerr << "unexpected exception thrown!" << endl;
+ }
+ }
+ {
+ cout << "test constructor with out of order data" << endl;
+ try{
+ doubleVector x,y;
+ x.push_back(1.0);
+ x.push_back(2.0);
+ x.push_back(1.5);
+
+ y.push_back(4.0);
+ y.push_back(5.0);
+ y.push_back(6.0);
+ periodicSpline spl(x,y);
+ cerr << "domain_error not thrown!" << endl;
+ }
+ catch(std::domain_error& ex){
+ cout << "OK, catched expected exception: "<< ex.what() << endl;
+ }
+ catch(...){
+ cerr << "unexpected exception thrown!" << endl;
+ }
+ }
+ {
+ cout << "Test copy constructor" << endl;
+ periodicSpline copy_spline(*test_spline);
+ double minx,maxx,copyminx,copymaxx;
+ test_spline->get_range(minx,maxx);
+ copy_spline.get_range(copyminx,copymaxx);
+ if( (minx != copyminx) || (maxx != copymaxx) ){
+ cerr << "COPY FAILED!" << endl;
+ return -1;
+ }
+ double range;
+ test_spline->get_range(minx,range);
+ double nsteps=100;
+ for(int i=0;i<nsteps;i++){
+ double val=(range*i)/nsteps;
+ double y,ycopy;
+ ycopy=copy_spline.evaluate(val);
+ y=test_spline->evaluate(val);
+ if(y != ycopy){
+ cerr << "COPY FAILED " << val << " " << y << " " << ycopy << endl;
+ }
+ }
+ }
+ {
+ cout << "Test copy operator" << endl;
+ periodicSpline copy_spline(*test_spline);
+ copy_spline=*test_spline;
+ double minx,maxx,copyminx,copymaxx;
+ test_spline->get_range(minx,maxx);
+ copy_spline.get_range(copyminx,copymaxx);
+ if( (minx != copyminx) || (maxx != copymaxx) ){
+ cerr << "COPY FAILED!" << endl;
+ return -1;
+ }
+ double range;
+ test_spline->get_range(minx,range);
+ double nsteps=100;
+ for(int i=0;i<nsteps;i++){
+ double val=(range*i)/nsteps;
+ double y,ycopy;
+ ycopy=copy_spline.evaluate(val);
+ y=test_spline->evaluate(val);
+ if(y != ycopy){
+ cerr << "COPY FAILED " << val << " " << y << " " << ycopy << endl;
+ }
+ }
+ }
+ cout << "test delete" << endl;
+ delete test_spline;
+ {
+ cout << "test leaks - new/delete" << endl;
+ periodicSpline *dynspl;
+ doubleVector x,y;
+ x.push_back(1.0);
+ x.push_back(2.0);
+ x.push_back(3.5);
+
+ y.push_back(4.0);
+ y.push_back(5.0);
+ y.push_back(6.0);
+ for(int i=0; i < 100 ; i++){
+ dynspl = new periodicSpline(x,y);
+ delete dynspl;
+ }
+ }
+ {
+ cout << "test leaks - stack" << endl;
+ doubleVector x,y;
+ x.push_back(1.0);
+ x.push_back(2.0);
+ x.push_back(3.5);
+
+ y.push_back(4.0);
+ y.push_back(5.0);
+ y.push_back(6.0);
+ for(int i=0; i < 100 ; i++){
+ periodicSpline stspl(x,y);
+ }
+ }
+
+ return 0;
+}
diff --git a/test_spline.cpp b/test_spline.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b2f0a97fccf92507e2f263af3fa33fc4bb569fad
--- /dev/null
+++ b/test_spline.cpp
@@ -0,0 +1,209 @@
+/*
+ * program to test spline classes
+ *
+ * $Author: claudio $
+ *
+ * $Revision: 1.1.1.1 $
+ *
+ * $Log: test_spline.cpp,v $
+ * Revision 1.1.1.1 2012-09-05 08:17:22 claudio
+ * Interpolator utility classes
+ *
+ *
+ */
+
+#include <iostream>
+#include "spline.h"
+#include <cmath>
+//#include <stdexcept>
+
+using namespace std;
+using namespace Interpolator;
+int main(){
+ cout<< "test of Spline interpolator class" << endl;
+ doubleVector x0,y0;
+ const int NP=30;
+ for(int i=0; i < NP; i++){
+ x0.push_back(i*i*.02);
+ y0.push_back(sin(i*i*.02));
+ }
+ Spline* test_spline = new Spline(x0,y0);
+ {
+ doubleVector a;
+ doubleVector b;
+ test_spline->get_base_points(a,b);
+ cout<<"test get_base_points()"<<endl;
+ for(int i=0; i < NP; i++ )
+ if((x0[i] != a[i]) ||(y0[i] != b[i])){
+ cout << i << " " << x0[i] << " " << a[i] << " | "<< y0[i] << " " << b[i] << endl;
+ return -1;
+ }
+
+ }
+ {
+ cout << "test get_range()" << endl;
+ double minx,maxx;
+ test_spline->get_range(minx,maxx);
+ if( (minx != x0[0]) || (maxx != x0[NP-1]) ){
+ cout << minx << " " << x0[0] << " | " << maxx << " " << x0[NP-1] << endl;
+ return -1;
+ }
+ }
+ {
+ cout << "test evaluate()" << endl;
+ double xn=0.02;
+ double yn=test_spline->evaluate(xn);
+ cout << xn << " " << yn << endl;
+ yn=test_spline->evaluate(x0[5]);
+ cout << xn << " " << yn << " " << y0[5] << endl;
+ }
+ {
+ cout << "test evaluate()" << endl;
+ double xn=0.02;
+ double yn=test_spline->evaluate(xn);
+ cout << xn << " " << yn << endl;
+ yn=test_spline->evaluate(x0[5]);
+ cout << xn << " " << yn << " " << y0[5] << endl;
+ }
+ { cout << "test evaluate()for plot" << endl;
+ doubleVector X;
+ doubleVector Y;
+ double range=x0[NP-1];
+ double nsteps=100;
+ for(int i=0;i<nsteps;i++){
+ X.push_back((range*i)/nsteps);
+ Y.push_back(test_spline->evaluate(X[i]));
+ }
+ //print out data
+ cout<<"---------------------------"<<endl;
+ for(int i=0; i < NP; i++) cout<<x0[i]<<"\t"<<y0[i]<<endl;
+ cout<<"---------------------------"<<endl;
+ for(int i=0; i < nsteps; i++) cout<<X[i]<<"\t"<<Y[i]<<endl;
+ cout<<"+++++++++++++++++++++++++++"<<endl;
+ }
+ {
+ cout << "test evaluate()for doubleVector" << endl;
+ doubleVector X;
+ doubleVector Y;
+ double range=x0[NP-1];
+ double nsteps=100;
+ for(int i=0;i<nsteps;i++)
+ X.push_back((range*i)/nsteps);
+ Y=test_spline->evaluate(X);
+ //print out data
+ cout<<"###########################"<<endl;
+ for(int i=0; i < nsteps; i++) cout<<X[i]<<"\t"<<Y[i]<<endl;
+ cout<<"***************************"<<endl;
+ }
+ {
+ cout << "test range exception" << endl;
+ try{
+ test_spline->evaluate(-10.0);
+ cerr << "range_error not thrown!" << endl;
+ }
+ catch(std::range_error& ex){
+ cout << "OK, catched expected exception: "<< ex.what() << endl;
+ }
+ catch(...){
+ cerr << "unexpected exception thrown!" << endl;
+ }
+ }
+ {
+ cout << "test constructor with unequal lengths" << endl;
+ try{
+ doubleVector x,y;
+ x.push_back(1.0);
+ x.push_back(2.0);
+ x.push_back(3.0);
+
+ y.push_back(4.0);
+ y.push_back(5.0);
+ Spline spl(x,y);
+ cerr << "length_error not thrown!" << endl;
+ }
+ catch(std::length_error& ex){
+ cout << "OK, catched expected exception: "<< ex.what() << endl;
+ }
+ catch(...){
+ cerr << "unexpected exception thrown!" << endl;
+ }
+ }
+ {
+ cout << "test constructor with out of order data" << endl;
+ try{
+ doubleVector x,y;
+ x.push_back(1.0);
+ x.push_back(2.0);
+ x.push_back(1.5);
+
+ y.push_back(4.0);
+ y.push_back(5.0);
+ y.push_back(6.0);
+ Spline spl(x,y);
+ cerr << "domain_error not thrown!" << endl;
+ }
+ catch(std::domain_error& ex){
+ cout << "OK, catched expected exception: "<< ex.what() << endl;
+ }
+ catch(...){
+ cerr << "unexpected exception thrown!" << endl;
+ }
+ }
+ {
+ cout << "Test copy constructor" << endl;
+ //Spline copy_spline=*test_spline;
+ Spline copy_spline(*test_spline);
+ double minx,maxx,copyminx,copymaxx;
+ test_spline->get_range(minx,maxx);
+ copy_spline.get_range(copyminx,copymaxx);
+ if( (minx != copyminx) || (maxx != copymaxx) ){
+ cerr << "COPY FAILED!" << endl;
+ return -1;
+ }
+ double range=x0[NP-1];
+ double nsteps=100;
+ for(int i=0;i<nsteps;i++){
+ double val=(range*i)/nsteps;
+ double y,ycopy;
+ ycopy=copy_spline.evaluate(val);
+ y=test_spline->evaluate(val);
+ if(y != ycopy){
+ cerr << "COPY FAILED " << val << " " << y << " " << ycopy << endl;
+ }
+ }
+ }
+ cout << "test delete" << endl;
+ delete test_spline;
+ {
+ cout << "test leaks - new/delete" << endl;
+ Spline *dynspl;
+ doubleVector x,y;
+ x.push_back(1.0);
+ x.push_back(2.0);
+ x.push_back(3.5);
+
+ y.push_back(4.0);
+ y.push_back(5.0);
+ y.push_back(6.0);
+ for(int i=0; i < 100 ; i++){
+ dynspl = new Spline(x,y);
+ delete dynspl;
+ }
+ }
+ {
+ cout << "test leaks - stack" << endl;
+ doubleVector x,y;
+ x.push_back(1.0);
+ x.push_back(2.0);
+ x.push_back(3.5);
+
+ y.push_back(4.0);
+ y.push_back(5.0);
+ y.push_back(6.0);
+ for(int i=0; i < 100 ; i++){
+ Spline stspl(x,y);
+ }
+ }
+
+ return 0;
+}