diff --git a/.coveragerc b/.coveragerc
new file mode 100644
index 0000000..6392437
--- /dev/null
+++ b/.coveragerc
@@ -0,0 +1,33 @@
+[run]
+
+include = numjuggler/*
+
+omit = 
+    *_tab.py
+
+# branch = True
+
+
+
+[report]
+
+# Regexes for lines to exclude from consideration
+exclude_lines =
+    # Have to re-enable the standard pragma
+    pragma: no cover
+
+    # Don't complain about missing debug-only code:
+    def __repr__
+    if self\.debug
+
+    # Don't complain if tests don't hit defensive assertion code:
+    raise AssertionError
+    raise NotImplementedError
+
+    # Don't complain if non-runnable code isn't run:
+    if 0:
+    if __name__ == .__main__.:
+
+ignore_errors = True
+sort = Cover
+
diff --git a/.gitattributes b/.gitattributes
new file mode 100644
index 0000000..0910289
--- /dev/null
+++ b/.gitattributes
@@ -0,0 +1,15 @@
+#
+# Special handling for binary files
+# See https://git-scm.com/book/tr/v2/Customizing-Git-Git-Attributes
+# Install doc2txt https://sourceforge.net/projects/docx2txt/?source=navbar
+#
+*.npz binary
+*.xz binary
+*.zip binary
+*.bz2 binary
+*.gz binary
+*.npy binary
+*.docx binary
+#*.docx binary diff=word
+*.png binary
+#*.png diff=exif
\ No newline at end of file
diff --git a/.gitignore b/.gitignore
index 198e48d..fd9267d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,7 @@
+# Python setup folders
+build/
+*.egg-info
+
 # Byte-compiled
 *.pyc
 
@@ -5,3 +9,18 @@
 .ropeproject/
 
 /.gtm/
+
+# Environments
+.cache/
+.idea/
+
+# Dev.scripts results
+numjuggler.dot
+numjuggler.pdf
+numjuggler.sfood1
+numjuggler.sfood2
+
+# Pytest and Coverage
+*.log
+.coverage
+htmlcov/
diff --git a/.pylintrc b/.pylintrc
new file mode 100644
index 0000000..c7bc4d7
--- /dev/null
+++ b/.pylintrc
@@ -0,0 +1,407 @@
+[MASTER]
+
+# Specify a configuration file.
+#rcfile=
+
+# Python code to execute, usually for sys.path manipulation such as
+# pygtk.require().
+#init-hook=
+
+# Add files or directories to the blacklist. They should be base names, not
+# paths.
+ignore=CVS
+
+# Add files or directories matching the regex patterns to the blacklist. The
+# regex matches against base names, not paths.
+ignore-patterns=
+
+# Pickle collected data for later comparisons.
+persistent=yes
+
+# List of plugins (as comma separated values of python modules names) to load,
+# usually to register additional checkers.
+load-plugins=
+
+# Use multiple processes to speed up Pylint.
+jobs=1
+
+# Allow loading of arbitrary C extensions. Extensions are imported into the
+# active Python interpreter and may run arbitrary code.
+unsafe-load-any-extension=no
+
+# A comma-separated list of package or module names from where C extensions may
+# be loaded. Extensions are loading into the active Python interpreter and may
+# run arbitrary code
+extension-pkg-whitelist=
+
+# Allow optimization of some AST trees. This will activate a peephole AST
+# optimizer, which will apply various small optimizations. For instance, it can
+# be used to obtain the result of joining multiple strings with the addition
+# operator. Joining a lot of strings can lead to a maximum recursion error in
+# Pylint and this flag can prevent that. It has one side effect, the resulting
+# AST will be different than the one from reality. This option is deprecated
+# and it will be removed in Pylint 2.0.
+optimize-ast=no
+
+
+[MESSAGES CONTROL]
+
+# Only show warnings with the listed confidence levels. Leave empty to show
+# all. Valid levels: HIGH, INFERENCE, INFERENCE_FAILURE, UNDEFINED
+confidence=
+
+# Enable the message, report, category or checker with the given id(s). You can
+# either give multiple identifier separated by comma (,) or put this option
+# multiple time (only on the command line, not in the configuration file where
+# it should appear only once). See also the "--disable" option for examples.
+#enable=
+
+# Disable the message, report, category or checker with the given id(s). You
+# can either give multiple identifiers separated by comma (,) or put this
+# option multiple times (only on the command line, not in the configuration
+# file where it should appear only once).You can also use "--disable=all" to
+# disable everything first and then reenable specific checks. For example, if
+# you want to run only the similarities checker, you can use "--disable=all
+# --enable=similarities". If you want to run only the classes checker, but have
+# no Warning level messages displayed, use"--disable=all --enable=classes
+# --disable=W"
+disable=missing-docstring,print-statement,parameter-unpacking,unpacking-in-except,old-raise-syntax,backtick,import-star-module-level,apply-builtin,basestring-builtin,buffer-builtin,cmp-builtin,coerce-builtin,execfile-builtin,file-builtin,long-builtin,raw_input-builtin,reduce-builtin,standarderror-builtin,unicode-builtin,xrange-builtin,coerce-method,delslice-method,getslice-method,setslice-method,no-absolute-import,old-division,dict-iter-method,dict-view-method,next-method-called,metaclass-assignment,indexing-exception,raising-string,reload-builtin,oct-method,hex-method,nonzero-method,cmp-method,input-builtin,round-builtin,intern-builtin,unichr-builtin,map-builtin-not-iterating,zip-builtin-not-iterating,range-builtin-not-iterating,filter-builtin-not-iterating,using-cmp-argument,long-suffix,old-ne-operator,old-octal-literal,suppressed-message,useless-suppression
+
+
+[REPORTS]
+
+# Set the output format. Available formats are text, parseable, colorized, msvs
+# (visual studio) and html. You can also give a reporter class, eg
+# mypackage.mymodule.MyReporterClass.
+output-format=text
+
+# Put messages in a separate file for each module / package specified on the
+# command line instead of printing them on stdout. Reports (if any) will be
+# written in a file name "pylint_global.[txt|html]". This option is deprecated
+# and it will be removed in Pylint 2.0.
+files-output=no
+
+# Tells whether to display a full report or only the messages
+reports=yes
+
+# Python expression which should return a note less than 10 (10 is the highest
+# note). You have access to the variables errors warning, statement which
+# respectively contain the number of errors / warnings messages and the total
+# number of statements analyzed. This is used by the global evaluation report
+# (RP0004).
+evaluation=10.0 - ((float(5 * error + warning + refactor + convention) / statement) * 10)
+
+# Template used to display messages. This is a python new-style format string
+# used to format the message information. See doc for all details
+#msg-template=
+
+
+[VARIABLES]
+
+# Tells whether we should check for unused import in __init__ files.
+init-import=no
+
+# A regular expression matching the name of dummy variables (i.e. expectedly
+# not used).
+dummy-variables-rgx=(_+[a-zA-Z0-9]*?$)|dummy
+
+# List of additional names supposed to be defined in builtins. Remember that
+# you should avoid to define new builtins when possible.
+additional-builtins=
+
+# List of strings which can identify a callback function by name. A callback
+# name must start or end with one of those strings.
+callbacks=cb_,_cb
+
+# List of qualified module names which can have objects that can redefine
+# builtins.
+redefining-builtins-modules=six.moves,future.builtins
+
+
+[LOGGING]
+
+# Logging modules to check that the string format arguments are in logging
+# function parameter format
+logging-modules=logging
+
+
+[TYPECHECK]
+
+# Tells whether missing members accessed in mixin class should be ignored. A
+# mixin class is detected if its name ends with "mixin" (case insensitive).
+ignore-mixin-members=yes
+
+# List of module names for which member attributes should not be checked
+# (useful for modules/projects where namespaces are manipulated during runtime
+# and thus existing member attributes cannot be deduced by static analysis. It
+# supports qualified module names, as well as Unix pattern matching.
+ignored-modules=
+
+# List of class names for which member attributes should not be checked (useful
+# for classes with dynamically set attributes). This supports the use of
+# qualified names.
+ignored-classes=optparse.Values,thread._local,_thread._local
+
+# List of members which are set dynamically and missed by pylint inference
+# system, and so shouldn't trigger E1101 when accessed. Python regular
+# expressions are accepted.
+generated-members=
+
+# List of decorators that produce context managers, such as
+# contextlib.contextmanager. Add to this list to register other decorators that
+# produce valid context managers.
+contextmanager-decorators=contextlib.contextmanager
+
+
+[SIMILARITIES]
+
+# Minimum lines number of a similarity.
+min-similarity-lines=4
+
+# Ignore comments when computing similarities.
+ignore-comments=yes
+
+# Ignore docstrings when computing similarities.
+ignore-docstrings=yes
+
+# Ignore imports when computing similarities.
+ignore-imports=no
+
+
+[SPELLING]
+
+# Spelling dictionary name. Available dictionaries: none. To make it working
+# install python-enchant package.
+spelling-dict=
+
+# List of comma separated words that should not be checked.
+spelling-ignore-words=
+
+# A path to a file that contains private dictionary; one word per line.
+spelling-private-dict-file=
+
+# Tells whether to store unknown words to indicated private dictionary in
+# --spelling-private-dict-file option instead of raising a message.
+spelling-store-unknown-words=no
+
+
+[FORMAT]
+
+# Maximum number of characters on a single line.
+max-line-length=132
+
+# Regexp for a line that is allowed to be longer than the limit.
+ignore-long-lines=^\s*(# )?<?https?://\S+>?$
+
+# Allow the body of an if to be on the same line as the test if there is no
+# else.
+single-line-if-stmt=no
+
+# List of optional constructs for which whitespace checking is disabled. `dict-
+# separator` is used to allow tabulation in dicts, etc.: {1  : 1,\n222: 2}.
+# `trailing-comma` allows a space between comma and closing bracket: (a, ).
+# `empty-line` allows space-only lines.
+no-space-check=trailing-comma,dict-separator
+
+# Maximum number of lines in a module
+max-module-lines=1000
+
+# String used as indentation unit. This is usually "    " (4 spaces) or "\t" (1
+# tab).
+indent-string='    '
+
+# Number of spaces of indent required inside a hanging  or continued line.
+indent-after-paren=4
+
+# Expected format of line ending, e.g. empty (any line ending), LF or CRLF.
+expected-line-ending-format=
+
+
+[MISCELLANEOUS]
+
+# List of note tags to take in consideration, separated by a comma.
+notes=FIXME,XXX,TODO
+
+
+[BASIC]
+
+# Good variable names which should always be accepted, separated by a comma
+good-names=i,j,k,ex,Run,_
+
+# Bad variable names which should always be refused, separated by a comma
+bad-names=foo,bar,baz,toto,tutu,tata
+
+# Colon-delimited sets of names that determine each other's naming style when
+# the name regexes allow several styles.
+name-group=
+
+# Include a hint for the correct naming format with invalid-name
+include-naming-hint=no
+
+# List of decorators that produce properties, such as abc.abstractproperty. Add
+# to this list to register other decorators that produce valid properties.
+property-classes=abc.abstractproperty
+
+# Regular expression matching correct module names
+module-rgx=(([a-z_][a-z0-9_]*)|([A-Z][a-zA-Z0-9]+))$
+
+# Naming hint for module names
+module-name-hint=(([a-z_][a-z0-9_]*)|([A-Z][a-zA-Z0-9]+))$
+
+# Regular expression matching correct constant names
+const-rgx=(([A-Z_][A-Z0-9_]*)|(__.*__))$
+
+# Naming hint for constant names
+const-name-hint=(([A-Z_][A-Z0-9_]*)|(__.*__))$
+
+# Regular expression matching correct class names
+class-rgx=[A-Z_][a-zA-Z0-9]+$
+
+# Naming hint for class names
+class-name-hint=[A-Z_][a-zA-Z0-9]+$
+
+# Regular expression matching correct function names
+function-rgx=[a-z_][a-z0-9_]{2,30}$
+
+# Naming hint for function names
+function-name-hint=[a-z_][a-z0-9_]{2,30}$
+
+# Regular expression matching correct method names
+method-rgx=[a-z_][a-z0-9_]{2,30}$
+
+# Naming hint for method names
+method-name-hint=[a-z_][a-z0-9_]{2,30}$
+
+# Regular expression matching correct attribute names
+attr-rgx=[a-z_][a-z0-9_]{2,30}$
+
+# Naming hint for attribute names
+attr-name-hint=[a-z_][a-z0-9_]{2,30}$
+
+# Regular expression matching correct argument names
+argument-rgx=[a-z_][a-z0-9_]{2,30}$
+
+# Naming hint for argument names
+argument-name-hint=[a-z_][a-z0-9_]{2,30}$
+
+# Regular expression matching correct variable names
+variable-rgx=[a-z_][a-z0-9_]{2,30}$
+
+# Naming hint for variable names
+variable-name-hint=[a-z_][a-z0-9_]{2,30}$
+
+# Regular expression matching correct class attribute names
+class-attribute-rgx=([A-Za-z_][A-Za-z0-9_]{2,30}|(__.*__))$
+
+# Naming hint for class attribute names
+class-attribute-name-hint=([A-Za-z_][A-Za-z0-9_]{2,30}|(__.*__))$
+
+# Regular expression matching correct inline iteration names
+inlinevar-rgx=[A-Za-z_][A-Za-z0-9_]*$
+
+# Naming hint for inline iteration names
+inlinevar-name-hint=[A-Za-z_][A-Za-z0-9_]*$
+
+# Regular expression which should only match function or class names that do
+# not require a docstring.
+no-docstring-rgx=^_
+
+# Minimum line length for functions/classes that require docstrings, shorter
+# ones are exempt.
+docstring-min-length=-1
+
+
+[ELIF]
+
+# Maximum number of nested blocks for function / method body
+max-nested-blocks=5
+
+
+[CLASSES]
+
+# List of method names used to declare (i.e. assign) instance attributes.
+defining-attr-methods=__init__,__new__,setUp
+
+# List of valid names for the first argument in a class method.
+valid-classmethod-first-arg=cls
+
+# List of valid names for the first argument in a metaclass class method.
+valid-metaclass-classmethod-first-arg=mcs
+
+# List of member names, which should be excluded from the protected access
+# warning.
+exclude-protected=_asdict,_fields,_replace,_source,_make
+
+
+[DESIGN]
+
+# Maximum number of arguments for function / method
+max-args=5
+
+# Argument names that match this expression will be ignored. Default to name
+# with leading underscore
+ignored-argument-names=_.*
+
+# Maximum number of locals for function / method body
+max-locals=15
+
+# Maximum number of return / yield for function / method body
+max-returns=6
+
+# Maximum number of branch for function / method body
+max-branches=12
+
+# Maximum number of statements in function / method body
+max-statements=50
+
+# Maximum number of parents for a class (see R0901).
+max-parents=7
+
+# Maximum number of attributes for a class (see R0902).
+max-attributes=7
+
+# Minimum number of public methods for a class (see R0903).
+min-public-methods=2
+
+# Maximum number of public methods for a class (see R0904).
+max-public-methods=20
+
+# Maximum number of boolean expressions in a if statement
+max-bool-expr=5
+
+
+[IMPORTS]
+
+# Deprecated modules which should not be used, separated by a comma
+deprecated-modules=optparse
+
+# Create a graph of every (i.e. internal and external) dependencies in the
+# given file (report RP0402 must not be disabled)
+import-graph=
+
+# Create a graph of external dependencies in the given file (report RP0402 must
+# not be disabled)
+ext-import-graph=
+
+# Create a graph of internal dependencies in the given file (report RP0402 must
+# not be disabled)
+int-import-graph=
+
+# Force import order to recognize a module as part of the standard
+# compatibility libraries.
+known-standard-library=
+
+# Force import order to recognize a module as part of a third party library.
+known-third-party=enchant
+
+# Analyse import fallback blocks. This can be used to support both Python 2 and
+# 3 compatible code, which means that the block might have code that exists
+# only in one or another interpreter, leading to false positives when analysed.
+analyse-fallback-blocks=no
+
+
+[EXCEPTIONS]
+
+# Exceptions that will emit a warning when being caught. Defaults to
+# "Exception"
+overgeneral-exceptions=Exception
diff --git a/numjuggler/__init__.py b/numjuggler/__init__.py
index 8f73a57..d4660be 100644
--- a/numjuggler/__init__.py
+++ b/numjuggler/__init__.py
@@ -4,5 +4,5 @@
 Package provides script mcnp.juggler. See its help for description.
 
 """
-
+from numjuggler.utils.PartialFormatter import PartialFormatter
 from .__version__ import __version__ as version
diff --git a/numjuggler/likefunc.py b/numjuggler/likefunc.py
index 167ddb0..7943dcf 100644
--- a/numjuggler/likefunc.py
+++ b/numjuggler/likefunc.py
@@ -7,7 +7,8 @@
 
 from __future__ import print_function
 from collections import OrderedDict
-import sys
+
+from numjuggler.utils.io import resolve_fname_or_stream
 
 
 def trivial(x):
@@ -59,7 +60,7 @@ def __init__(self, log=False):
         self.ld = {}
 
         # Flag to store log
-        self.log = log
+        self.__lf = log
 
         # Default mapping, applied to elements not in ranges
         self.default = trivial
@@ -68,30 +69,19 @@ def __init__(self, log=False):
         self.doc = ""
         return
 
-    def __call1(self, x):
-        res = self.get_value(x)
-        return res
-
-    def __call2(self, x):
+    def __call__(self, x):
         res = self.get_value(x)
-        self.ld[x] = res
+        if self.log:
+            self.ld[x] = res
         return res
 
-    def __call__(self, x):
-        return self.__call1(x)
-
     @property
     def log(self):
         return self.__lf
 
     @log.setter
-    def log(self, v):
-        v = bool(v)
-        if v:
-            self.__call__ = self.__call2
-        else:
-            self.__call__ = self.__call1
-        self.__lf = v
+    def log(self, _log):
+        self.__lf = _log
 
     def __str__(self):
         res = []
@@ -102,25 +92,12 @@ def __str__(self):
         res.append('other -> {}'.format(self.default._mydoc))
         return '\n'.join(res)
 
-    def _get_log_file(self, fname):
-        if self.log is None:
+    def write_log_as_map(self, t, fname_or_stream=None):
+        if not self.log:
             raise ValueError("Cannon write log for unlogged mapping.")
-
-        # Define where to print
-        if fname is None:
-            fout = sys.stdout
-        else:
-            try:
-                fout = open(fname, 'w')
-            except:
-                print("Cannot open {}. Print to stdout".format(repr(fname)))
-                fout = sys.stdout
-        return fout
-
-    def write_log_as_map(self, t, fname=None):
-        fout = self._get_log_file(fname)
-        for nold, nnew in self.ld.items():
-            print('{} {}: {}'.format(t, nnew, nold), file=fout)
+        with resolve_fname_or_stream(fname_or_stream, "w") as fout:
+            for nold, nnew in self.ld.items():
+                print('{} {}: {}'.format(t, nnew, nold), file=fout)
 
 
 class LikeFunction(LikeFunctionBase):
@@ -262,7 +239,7 @@ def read_map_file(fname, log=False):
     # Dictionary type -> LikeFunction
     maps = {}
 
-    with open(fname, 'r') as mapfile:
+    with resolve_fname_or_stream(fname, 'r') as mapfile:
         for l in mapfile:
             t, ranges, f = _parse_map_line(l)
             if t is None:
diff --git a/numjuggler/main.py b/numjuggler/main.py
index e59ff45..176059a 100644
--- a/numjuggler/main.py
+++ b/numjuggler/main.py
@@ -3,14 +3,23 @@
 from __future__ import print_function
 
 import argparse as ap
+import sys
 from math import pi as Pi
 import os.path
 from numjuggler import numbering as mn
 from numjuggler import parser as mp
 from numjuggler import ri_notation as rin
 from numjuggler import string_cells as stc
+from numjuggler import likefunc as lf
 from numjuggler import version
 
+try:
+    import pirs.mcnp.mctal.Mctal as Mctal, Vector3, Material
+except:
+    Material = None
+    Mctal = None
+    Vector3 = None
+
 
 def multiline(lines, prefix=''):
     return prefix + ('\n' + prefix).join(lines)
@@ -77,7 +86,7 @@ def tr2str(pl, fmt1='{:12.9f}', fmte='{:16.8e}'):
          'annotate', 'getc', 'mnew', 'combinec', 'cdens')
 
 
-def main():
+def main(args=sys.argv[1:]):
     p = ap.ArgumentParser(prog='numjuggler', description=descr, epilog=epilog)
     p.add_argument('--version', action='version',
                    version='%(prog)s {}'.format(version))
@@ -122,7 +131,7 @@ def main():
                     nargs='?',
                     default='',
                     const='gen')
-    harg, clo = ph.parse_known_args()
+    harg, clo = ph.parse_known_args(args)
     if harg.h:
         if harg.h == 'gen':
             p.print_help()
@@ -260,7 +269,7 @@ def main():
             # as well.
 
             import re
-            r = re.compile('(u)(\d+)')
+            r = re.compile(r'(u)(\d+)')
 
             csets = {}
             ulst = []
@@ -530,7 +539,7 @@ def cfunc(n):
                     rms.update(rml)
 
             # read reference materials and create Materials
-            from pirs.mcnp import Material
+            assert Material is not None,  "pirs Material class is not imported"
             rmd = {}
             for c in cards:
                 if c.ctype == mp.CID.data:
@@ -747,9 +756,9 @@ def cfunc(n):
 
         elif args.mode == 'zrotate':
 
-            from pirs.core.trageom import Vector3, pi
+            assert  Vector3 is not None, "pirs Vector3 class is not imported"
             ag = float(args.c)    # in grad
-            ar = ag * pi / 180.   # in radians
+            ar = ag * Pi / 180.   # in radians
 
             # new transformation number:
             trn = args.t
@@ -1234,7 +1243,7 @@ def check(v, s):
             # -m argument is the mctal name followed by tally number of the
             # tally containing cell volumes.
             if args.m != '0':
-                from pirs.mcnp.mctal import Mctal
+                assert Mctal is not None, "pirs Mctal class is not imported"
                 mctal = Mctal()
                 fname, tn = args.m.split()
                 tn = int(tn)
@@ -1604,7 +1613,6 @@ def print_spherical(s, r):
                 print(c.card(), end='')
 
         elif args.mode == 'renum':
-            import likefunc as lf
             if args.map:
                 maps = lf.read_map_file(args.map, log=args.log != '')
             else:
diff --git a/numjuggler/nogq2.py b/numjuggler/nogq2.py
index 75db636..8a47193 100644
--- a/numjuggler/nogq2.py
+++ b/numjuggler/nogq2.py
@@ -1,444 +1,451 @@
-from math import copysign, isnan
-from pirs.core.trageom import vector
-
-
-def areclose(l, rtol=1e-4, atol=1e-7, cmnt=None, name='', detailed=True):
-    """
-    Check if all elements in the list l are close to each other.
-
-    Implementation originally based on numpy.isclose formula
-    """
-    f1 = '{}'
-    f2 = '{:15.8e}'
-
-    A = min(l)
-    B = max(l)
-    if atol is None:
-        ares = False
-        astr = f1.format(atol)
-    else:
-        ares = B - A <= atol
-        astr = f2.format(atol)
-    if rtol is None:
-        rres = False
-        rstr = f1.format(rtol)
-    else:
-        b = max(map(abs, l))
-        rres = B - A <= rtol*b
-        rstr = f2.format(rtol*b)
-
-    result = ares or rres
-
-    if cmnt is not None:
-        # assume it is a list of comments. Add here information
-        c = cmnt.append
-        c('Are close check: ' + name + ': {}'.format(result))
-        if detailed:
-            c('  values: ' + ' '.join('{:15.8e}'.format(v) for v in l))
-            c('  B:      {:15.8e}'.format(B))
-            c('  A:      {:15.8e}'.format(A))
-            c('  B - A:  {:15.8e}'.format(B - A))
-            c('  atol:   ' + astr)
-            c('  rtol*b: ' + rstr)
-    return result
-
-
-def get_params(card):
-    """
-    Return parameters of the GQ card.
-    """
-    p1, p2 = card.lower().split('gq')
-    # Check if there is transformation:
-    tr = len(p1.split()) > 1
-    return tr, map(float, p2.split())
-
-
-def get_cone_or_cyl(pl):
-    """
-    See notes from 14.03.2018
-    """
-    cmnt = []
-
-    # Define normalization coeff gamma:
-    A, B, C, D, E, F, G, H, J, K = pl
-    ABC = sum((A, B, C))
-    cmnt.append(' A, B, C:' + ' '.join('{:15.8e}'.format(v) for v in (A, B, C)))
-    cmnt.append(' D, E, F:' + ' '.join('{:15.8e}'.format(v) for v in (D, E, F)))
-    cmnt.append(' G, H, J:' + ' '.join('{:15.8e}'.format(v) for v in (G, H, J)))
-    cmnt.append('       K:' + ' {:15.8e}'.format(K))
-
-    # List of normalization coefficients. 1 is always assumed
-    gammas = [1.0, -1.0]
-    tDEF = 5e-5
-    tABC = 1e-9
-    if areclose((0, E), atol=tDEF, rtol=None):
-        if areclose((1, B), atol=tABC, rtol=None):
-            g = 1.0/B
-            gammas.append(g)
-            cmnt.append('gamma_B = {:15.8e}'.format(g))
-        if areclose((1, C), atol=tABC, rtol=None):
-            g = 1.0/C
-            gammas.append(g)
-            cmnt.append('gamma_C = {:15.8e}'.format(g))
-    else:
-        g = 2*E/(2*E*A - D*F)
-        gammas.append(g)
-        cmnt.append('gamma_E = {:15.8e}'.format(g))
-
-    if areclose((0, F), atol=tDEF, rtol=None):
-        if areclose((1, A), atol=tABC, rtol=None):
-            g = 1.0/A
-            gammas.append(g)
-            cmnt.append('gamma_A = {:15.8e}'.format(g))
-        if areclose((1, C), atol=tABC, rtol=None):
-            g = 1.0/C
-            gammas.append(g)
-            cmnt.append('gamma_C = {:15.8e}'.format(g))
-    else:
-        g = 2*F/(2*F*B - D*E)
-        gammas.append(g)
-        cmnt.append('gamma_F = {:15.8e}'.format(g))
-
-    if areclose((0, D), atol=tDEF, rtol=None):
-        if areclose((1, B), atol=tABC, rtol=None):
-            g = 1.0/B
-            gammas.append(g)
-            cmnt.append('gamma_B = {:15.8e}'.format(g))
-        if areclose((1, A), atol=tABC, rtol=None):
-            g = 1.0/A
+from __future__ import print_function, division, nested_scopes
+import six
+try:
+    import pirs.core.trageom.vector as vector
+except:
+    vector = None
+
+if six.PY2 and vector is not None:
+    from math import copysign, isnan
+
+
+    def areclose(l, rtol=1e-4, atol=1e-7, cmnt=None, name='', detailed=True):
+        """
+        Check if all elements in the list l are close to each other.
+
+        Implementation originally based on numpy.isclose formula
+        """
+        f1 = '{}'
+        f2 = '{:15.8e}'
+
+        A = min(l)
+        B = max(l)
+        if atol is None:
+            ares = False
+            astr = f1.format(atol)
+        else:
+            ares = B - A <= atol
+            astr = f2.format(atol)
+        if rtol is None:
+            rres = False
+            rstr = f1.format(rtol)
+        else:
+            b = max(map(abs, l))
+            rres = B - A <= rtol*b
+            rstr = f2.format(rtol*b)
+
+        result = ares or rres
+
+        if cmnt is not None:
+            # assume it is a list of comments. Add here information
+            c = cmnt.append
+            c('Are close check: ' + name + ': {}'.format(result))
+            if detailed:
+                c('  values: ' + ' '.join('{:15.8e}'.format(v) for v in l))
+                c('  B:      {:15.8e}'.format(B))
+                c('  A:      {:15.8e}'.format(A))
+                c('  B - A:  {:15.8e}'.format(B - A))
+                c('  atol:   ' + astr)
+                c('  rtol*b: ' + rstr)
+        return result
+
+
+    def get_params(card):
+        """
+        Return parameters of the GQ card.
+        """
+        p1, p2 = card.lower().split('gq')
+        # Check if there is transformation:
+        tr = len(p1.split()) > 1
+        return tr, map(float, p2.split())
+
+
+    def get_cone_or_cyl(pl):
+        """
+        See notes from 14.03.2018
+        """
+        cmnt = []
+
+        # Define normalization coeff gamma:
+        A, B, C, D, E, F, G, H, J, K = pl
+        ABC = sum((A, B, C))
+        cmnt.append(' A, B, C:' + ' '.join('{:15.8e}'.format(v) for v in (A, B, C)))
+        cmnt.append(' D, E, F:' + ' '.join('{:15.8e}'.format(v) for v in (D, E, F)))
+        cmnt.append(' G, H, J:' + ' '.join('{:15.8e}'.format(v) for v in (G, H, J)))
+        cmnt.append('       K:' + ' {:15.8e}'.format(K))
+
+        # List of normalization coefficients. 1 is always assumed
+        gammas = [1.0, -1.0]
+        tDEF = 5e-5
+        tABC = 1e-9
+        if areclose((0, E), atol=tDEF, rtol=None):
+            if areclose((1, B), atol=tABC, rtol=None):
+                g = 1.0/B
+                gammas.append(g)
+                cmnt.append('gamma_B = {:15.8e}'.format(g))
+            if areclose((1, C), atol=tABC, rtol=None):
+                g = 1.0/C
+                gammas.append(g)
+                cmnt.append('gamma_C = {:15.8e}'.format(g))
+        else:
+            g = 2*E/(2*E*A - D*F)
             gammas.append(g)
-            cmnt.append('gamma_A = {:15.8e}'.format(g))
-    else:
-        g = 2*D/(2*D*C - E*F)
-        gammas.append(g)
-        cmnt.append('gamma_d = {:15.8e}'.format(g))
-
-    # Ensure that gamma=1 is considered first
-    gammas = set(gammas)
-    gammas.remove(1)
-    gammas.remove(-1)
-    gammas = (1.0, -1.0) + tuple(gammas)
-    for gamma in gammas:
-        t2, n, R0c, r2c, R0k, r2k, c1, c2 = get_surface_parameters(gamma, pl)
-        n2 = scalar_product(n, n)
-        R0kn = scalar_product(n, R0k)
-        R0k2 = scalar_product(R0k, R0k)
-        R0cn = scalar_product(n, R0c)
-        R0c2 = scalar_product(R0c, R0c)
-        tt = t2**0.5 if t2 >= 0 else float('nan')
-        rc = r2c**0.5 if r2c >= 0 else float('nan')
-        rk = r2k**0.5 if r2k >= 0 else float('nan')
-        cmnt.append('gamma: 1 + {:15.8e}'.format(gamma - 1.0))
-        cmnt.append('     t^2: {:15.8e}'.format(t2))
-        cmnt.append('      t : {:15.8e}'.format(tt))
-        cmnt.append('      n : ' + ' '.join('{:15.8e}'.format(v) for v in n))
-        cmnt.append('   (n,n): {:15.8e}'.format(n2))
-        cmnt.append('     r^2: {:15.8e}  {:15.8e}'.format(r2c, r2k))
-        cmnt.append('       r: {:15.8e}  {:15.8e}'.format(rc, rk))
-        cmnt.append('     R0c: ' + ' '.join('{:15.8e}'.format(v) for v in R0c))
-        cmnt.append('     R0k: ' + ' '.join('{:15.8e}'.format(v) for v in R0k))
-        cmnt.append('  (n,R0): {:15.8e}  {:15.8e}'.format(R0cn, R0kn))
-        cmnt.append(' (R0,R0): {:15.8e}  {:15.8e}'.format(R0c2, R0k2))
-        cmnt.append('      c1: {:15.8e}'.format(c1))
-        cmnt.append('      c2: {:15.8e}'.format(c2))
-
-        # Check parameters common for cone and cylinder
-        if isnan(n2) or areclose((0, n2), atol=1e-4, rtol=None):
-            cmnt.append('        gamma sorted out due to n')
-            continue
-
-        # Evaluate surface at some points
-        distances = (-100, -10, -1, 0, 1, 10, 100)
-        rsdc = {}  # residuals for cylinder
-        rsdk = {}  # residuals for cone
-        ni, nj, nk = vector.Vector3(car=n).basis()
-        if areclose((0, r2c), atol=1e-6, rtol=None) or isnan(rc) or isnan(R0c2):
-            can_be_cylinder = False
-            cmnt.append('        cannot be cylinder due to r2 or R0')
+            cmnt.append('gamma_E = {:15.8e}'.format(g))
+
+        if areclose((0, F), atol=tDEF, rtol=None):
+            if areclose((1, A), atol=tABC, rtol=None):
+                g = 1.0/A
+                gammas.append(g)
+                cmnt.append('gamma_A = {:15.8e}'.format(g))
+            if areclose((1, C), atol=tABC, rtol=None):
+                g = 1.0/C
+                gammas.append(g)
+                cmnt.append('gamma_C = {:15.8e}'.format(g))
         else:
-            can_be_cylinder = True
-            r0 = vector.Vector3(car=R0c)
-            nj.R = rc
-            nk.R = rc
-            for d in distances:
-                r00 = r0 + d*ni
-                d1 = evaluate_gq(pl, (r00 + nj).car)
-                d2 = evaluate_gq(pl, (r00 - nj).car)
-                d3 = evaluate_gq(pl, (r00 + nk).car)
-                d4 = evaluate_gq(pl, (r00 - nk).car)
-                rsdc[d] = (d1, d2, d3, d4)
-        if areclose((0, t2), atol=1e-6, rtol=None) or isnan(tt) or isnan(R0k2):
-            can_be_cone = False
-            cmnt.append('        cannot be cone due to t2 or R0')
+            g = 2*F/(2*F*B - D*E)
+            gammas.append(g)
+            cmnt.append('gamma_F = {:15.8e}'.format(g))
+
+        if areclose((0, D), atol=tDEF, rtol=None):
+            if areclose((1, B), atol=tABC, rtol=None):
+                g = 1.0/B
+                gammas.append(g)
+                cmnt.append('gamma_B = {:15.8e}'.format(g))
+            if areclose((1, A), atol=tABC, rtol=None):
+                g = 1.0/A
+                gammas.append(g)
+                cmnt.append('gamma_A = {:15.8e}'.format(g))
         else:
-            can_be_cone = True
-            r0 = vector.Vector3(car=R0k)
-            for d in distances:
-                rkd = tt * d
-                nj.R = rkd
-                nk.R = rkd
-                r00 = r0 + d*ni
-                d1 = evaluate_gq(pl, (r00 + nj).car)
-                d2 = evaluate_gq(pl, (r00 - nj).car)
-                d3 = evaluate_gq(pl, (r00 + nk).car)
-                d4 = evaluate_gq(pl, (r00 - nk).car)
-                rsdk[d] = (d1, d2, d3, d4)
-
-        # Compare residuals for cone and cylinder and choose one
-        if can_be_cone and can_be_cylinder:
-            typ = 0
-            for d in distances:
-                dc = scalar_product(rsdc[d], rsdc[d])
-                dk = scalar_product(rsdk[d], rsdk[d])
-                if dc <= dk:
-                    typ += 1
+            g = 2*D/(2*D*C - E*F)
+            gammas.append(g)
+            cmnt.append('gamma_d = {:15.8e}'.format(g))
+
+        # Ensure that gamma=1 is considered first
+        gammas = set(gammas)
+        gammas.remove(1)
+        gammas.remove(-1)
+        gammas = (1.0, -1.0) + tuple(gammas)
+        for gamma in gammas:
+            t2, n, R0c, r2c, R0k, r2k, c1, c2 = get_surface_parameters(gamma, pl)
+            n2 = scalar_product(n, n)
+            R0kn = scalar_product(n, R0k)
+            R0k2 = scalar_product(R0k, R0k)
+            R0cn = scalar_product(n, R0c)
+            R0c2 = scalar_product(R0c, R0c)
+            tt = t2**0.5 if t2 >= 0 else float('nan')
+            rc = r2c**0.5 if r2c >= 0 else float('nan')
+            rk = r2k**0.5 if r2k >= 0 else float('nan')
+            cmnt.append('gamma: 1 + {:15.8e}'.format(gamma - 1.0))
+            cmnt.append('     t^2: {:15.8e}'.format(t2))
+            cmnt.append('      t : {:15.8e}'.format(tt))
+            cmnt.append('      n : ' + ' '.join('{:15.8e}'.format(v) for v in n))
+            cmnt.append('   (n,n): {:15.8e}'.format(n2))
+            cmnt.append('     r^2: {:15.8e}  {:15.8e}'.format(r2c, r2k))
+            cmnt.append('       r: {:15.8e}  {:15.8e}'.format(rc, rk))
+            cmnt.append('     R0c: ' + ' '.join('{:15.8e}'.format(v) for v in R0c))
+            cmnt.append('     R0k: ' + ' '.join('{:15.8e}'.format(v) for v in R0k))
+            cmnt.append('  (n,R0): {:15.8e}  {:15.8e}'.format(R0cn, R0kn))
+            cmnt.append(' (R0,R0): {:15.8e}  {:15.8e}'.format(R0c2, R0k2))
+            cmnt.append('      c1: {:15.8e}'.format(c1))
+            cmnt.append('      c2: {:15.8e}'.format(c2))
+
+            # Check parameters common for cone and cylinder
+            if isnan(n2) or areclose((0, n2), atol=1e-4, rtol=None):
+                cmnt.append('        gamma sorted out due to n')
+                continue
+
+            # Evaluate surface at some points
+            distances = (-100, -10, -1, 0, 1, 10, 100)
+            rsdc = {}  # residuals for cylinder
+            rsdk = {}  # residuals for cone
+            ni, nj, nk = vector.Vector3(car=n).basis()
+            if areclose((0, r2c), atol=1e-6, rtol=None) or isnan(rc) or isnan(R0c2):
+                can_be_cylinder = False
+                cmnt.append('        cannot be cylinder due to r2 or R0')
+            else:
+                can_be_cylinder = True
+                r0 = vector.Vector3(car=R0c)
+                nj.R = rc
+                nk.R = rc
+                for d in distances:
+                    r00 = r0 + d*ni
+                    d1 = evaluate_gq(pl, (r00 + nj).car)
+                    d2 = evaluate_gq(pl, (r00 - nj).car)
+                    d3 = evaluate_gq(pl, (r00 + nk).car)
+                    d4 = evaluate_gq(pl, (r00 - nk).car)
+                    rsdc[d] = (d1, d2, d3, d4)
+            if areclose((0, t2), atol=1e-6, rtol=None) or isnan(tt) or isnan(R0k2):
+                can_be_cone = False
+                cmnt.append('        cannot be cone due to t2 or R0')
+            else:
+                can_be_cone = True
+                r0 = vector.Vector3(car=R0k)
+                for d in distances:
+                    rkd = tt * d
+                    nj.R = rkd
+                    nk.R = rkd
+                    r00 = r0 + d*ni
+                    d1 = evaluate_gq(pl, (r00 + nj).car)
+                    d2 = evaluate_gq(pl, (r00 - nj).car)
+                    d3 = evaluate_gq(pl, (r00 + nk).car)
+                    d4 = evaluate_gq(pl, (r00 - nk).car)
+                    rsdk[d] = (d1, d2, d3, d4)
+
+            # Compare residuals for cone and cylinder and choose one
+            if can_be_cone and can_be_cylinder:
+                typ = 0
+                for d in distances:
+                    dc = scalar_product(rsdc[d], rsdc[d])
+                    dk = scalar_product(rsdk[d], rsdk[d])
+                    if dc <= dk:
+                        typ += 1
+                    else:
+                        typ -= 1
+                if typ > 0:
+                    typ = 'c'
                 else:
-                    typ -= 1
-            if typ > 0:
+                    typ = 'k'
+            elif can_be_cone:
+                typ = 'k'
+            elif can_be_cylinder:
                 typ = 'c'
             else:
-                typ = 'k'
-        elif can_be_cone:
-            typ = 'k'
-        elif can_be_cylinder:
-            typ = 'c'
-        else:
-            typ = 'o'
-            continue
-
-        # Prepare output
-        if typ == 'k':
-            org = R0k
-            r2 = r2k
-            rsd = rsdk
-        elif typ == 'c':
-            org = R0c
-            r2 = r2c
-            rsd = rsdc
-
-        rsdmax = max(map(abs, sum(rsd.values(), ())))
-        cmnt.append(' Residuals for {}, {:15.8e}'.format(typ, rsdmax))
-        for d in distances:
-            cmnt.append(' at d={:10.3e}:'.format(d) + ' '.join('{:15.8e}'.format(v) for v in rsd[d]))
-        if rsdmax > 1e-1:
-            typ = 'o'
-            continue
-        else:
-            cmnt.append(' Final max. residual for {}, {:15.8e}'.format(typ, rsdmax))
+                typ = 'o'
+                continue
+
+            # Prepare output
+            if typ == 'k':
+                org = R0k
+                r2 = r2k
+                rsd = rsdk
+            elif typ == 'c':
+                org = R0c
+                r2 = r2c
+                rsd = rsdc
+
+            rsdmax = max(map(abs, sum(rsd.values(), ())))
+            cmnt.append(' Residuals for {}, {:15.8e}'.format(typ, rsdmax))
+            for d in distances:
+                cmnt.append(' at d={:10.3e}:'.format(d) + ' '.join('{:15.8e}'.format(v) for v in rsd[d]))
+            if rsdmax > 1e-1:
+                typ = 'o'
+                continue
+            else:
+                cmnt.append(' Final max. residual for {}, {:15.8e}'.format(typ, rsdmax))
 
+            cmnt = ['c ' + c for c in cmnt]
+            return typ, n, org, t2, r2, cmnt
+        typ = 'o'
         cmnt = ['c ' + c for c in cmnt]
-        return typ, n, org, t2, r2, cmnt
-    typ = 'o'
-    cmnt = ['c ' + c for c in cmnt]
-    return typ, None, None, None, None, cmnt
-
-
-def get_surface_parameters(gamma, pl):
-    """
-    pl -- list of original GQ parameters,
-    gamma -- normlaization coefficient
-    """
-    # normliaze GQ parameters
-    pl = (gamma*v for v in pl)
-    A, B, C, D, E, F, G, H, J, K = pl
-
-    # parameter t^2
-    t2 = sum((2.0, -A, -B, -C))
-
-    # parameters xn, yn, zn
-    t21 = sum((3.0, -A, -B, -C))
-    xn2 = (1 - A)/t21
-    yn2 = (1 - B)/t21
-    zn2 = (1 - C)/t21
-    xn = float('nan')
-    yn = float('nan')
-    zn = float('nan')
-    if xn2 >= 0:
-        xn = xn2 ** 0.5
-    if yn2 >= 0:
-        yn = yn2 ** 0.5
-    if zn2 >= 0:
-        zn = zn2 ** 0.5
-    mn2 = max((xn2, yn2, zn2))
-    if xn2 == mn2:
-        yn = copysign(yn, -D)
-        zn = copysign(zn, -F)
-    elif yn2 == mn2:
-        zn = copysign(zn, -E)
-        xn = copysign(xn, -D)
-    else:
-        xn = copysign(xn, -F)
-        yn = copysign(yn, -E)
-
-    # Parameters x0, y0, z0
-    # and r2 (cylinder radius, or error in cone focus position)
-
-    # Assuming that t2 is 0 (i.e. formulae for cylinder)
-    x0c = -G * 0.5
-    y0c = -H * 0.5
-    z0c = -J * 0.5
-    r2c = x0c**2 + y0c**2 + z0c**2 - K
-
-    # Formulae for cone
-    cR0n = (G*xn + H*yn + J*zn)*t21/2.0
-    if t2 > 0:
-        cR0n = cR0n / t2
-    else:
-        cR0n = copysign(float('inf'), cR0n*t2)
-    x0k = xn*cR0n - G*0.5
-    y0k = yn*cR0n - H*0.5
-    z0k = zn*cR0n - J*0.5
-    r2k = -(G*x0k + H*y0k + J*z0k)/2 - K
-
-    # Consistency checks
-    rgh1 = D*E*F
-    lft1 = -8*(1 - A)*(1 - B)*(1 - C)
-    rgh2 = (G*xn + H*yn + J*zn)**2 * t21
-    lft2 = G**2*(1 - A) + H**2*(1 - B) + J**2*(1 - C) - G*H*D - J*G*F - J*H*E
-    c1 = rgh1 - lft1
-    c2 = rgh2 - lft2
-    return (t2, (xn, yn, zn),
-            (x0c, y0c, z0c), r2c,    # cylinder formulae
-            (x0k, y0k, z0k), r2k,    # cone formulae
-            c1, c2)
-
-
-def evaluate_gq(pl, p):
-    """
-    Evaluate GQ equation at point p.
-    """
-    x, y, z = p
-    A, B, C, D, E, F, G, H, J, K = pl
-    d = (A*x**2 + B*y**2 + C*z**2 +
-         D*x*y  + E*y*z  + F*z*x  +
-         G*x    + H*y    + J*z    + K)
-    return d
-
-
-def basis_on_axis(axis):
-    """
-    Return basis.
-
-    One of the basis vectors has components xn, yn, zn. The others are chosen
-    to be as close as possible to the original basis vectors.
-
-    Components of `axis` are normalized to 1, and the component that has the
-    maximal absolute value is positive.
-    """
-    m = max(axis)
-    n = 0      # number of shifts
-    o = 'xyz'  # axis names
-    t = tuple(axis)
-    while t[0] != m:
-        t = (t[2], t[0], t[1])
-        o = o[2] + o[0] + o[1]
-        n += 1
-    xn, yn, zn = t
-    # After i shifts, xn component is the largest one. Now the original basis
-    # vectors are named so, that the largest component of the axis is xn, and
-    # the 1-st basis vector, i', is `axis`. The second basis vector, j', is most
-    # closest to the 2-nd basis vector of the original CS, j. This means, that
-    # j' is perpendicular to i' and on the plane built on i' and j.
-    i = (xn, yn, zn)
-    xn2 = xn**2
-    yn2 = yn**2
-    zn2 = zn**2
-    b = 1.0/((1.0 - yn2)**2 + yn2*(xn2 + zn2))**0.5
-    a = -b*yn
-    j = (xn*a, yn*a + b, zn*a)
-    k = cross_product(i, j)
-
-    # shift vector names back n times
-    b = (i, j, k)
-    for kkk in range(n):
-        i = b[1][1], b[1][2], b[1][0]
-        j = b[2][1], b[2][2], b[2][0]
-        k = b[0][1], b[0][2], b[0][0]
+        return typ, None, None, None, None, cmnt
+
+
+    def get_surface_parameters(gamma, pl):
+        """
+        pl -- list of original GQ parameters,
+        gamma -- normlaization coefficient
+        """
+        # normliaze GQ parameters
+        pl = (gamma*v for v in pl)
+        A, B, C, D, E, F, G, H, J, K = pl
+
+        # parameter t^2
+        t2 = sum((2.0, -A, -B, -C))
+
+        # parameters xn, yn, zn
+        t21 = sum((3.0, -A, -B, -C))
+        xn2 = (1 - A)/t21
+        yn2 = (1 - B)/t21
+        zn2 = (1 - C)/t21
+        xn = float('nan')
+        yn = float('nan')
+        zn = float('nan')
+        if xn2 >= 0:
+            xn = xn2 ** 0.5
+        if yn2 >= 0:
+            yn = yn2 ** 0.5
+        if zn2 >= 0:
+            zn = zn2 ** 0.5
+        mn2 = max((xn2, yn2, zn2))
+        if xn2 == mn2:
+            yn = copysign(yn, -D)
+            zn = copysign(zn, -F)
+        elif yn2 == mn2:
+            zn = copysign(zn, -E)
+            xn = copysign(xn, -D)
+        else:
+            xn = copysign(xn, -F)
+            yn = copysign(yn, -E)
+
+        # Parameters x0, y0, z0
+        # and r2 (cylinder radius, or error in cone focus position)
+
+        # Assuming that t2 is 0 (i.e. formulae for cylinder)
+        x0c = -G * 0.5
+        y0c = -H * 0.5
+        z0c = -J * 0.5
+        r2c = x0c**2 + y0c**2 + z0c**2 - K
+
+        # Formulae for cone
+        cR0n = (G*xn + H*yn + J*zn)*t21/2.0
+        if t2 > 0:
+            cR0n = cR0n / t2
+        else:
+            cR0n = copysign(float('inf'), cR0n*t2)
+        x0k = xn*cR0n - G*0.5
+        y0k = yn*cR0n - H*0.5
+        z0k = zn*cR0n - J*0.5
+        r2k = -(G*x0k + H*y0k + J*z0k)/2 - K
+
+        # Consistency checks
+        rgh1 = D*E*F
+        lft1 = -8*(1 - A)*(1 - B)*(1 - C)
+        rgh2 = (G*xn + H*yn + J*zn)**2 * t21
+        lft2 = G**2*(1 - A) + H**2*(1 - B) + J**2*(1 - C) - G*H*D - J*G*F - J*H*E
+        c1 = rgh1 - lft1
+        c2 = rgh2 - lft2
+        return (t2, (xn, yn, zn),
+                (x0c, y0c, z0c), r2c,    # cylinder formulae
+                (x0k, y0k, z0k), r2k,    # cone formulae
+                c1, c2)
+
+
+    def evaluate_gq(pl, p):
+        """
+        Evaluate GQ equation at point p.
+        """
+        x, y, z = p
+        A, B, C, D, E, F, G, H, J, K = pl
+        d = (A*x**2 + B*y**2 + C*z**2 +
+             D*x*y  + E*y*z  + F*z*x  +
+             G*x    + H*y    + J*z    + K)
+        return d
+
+
+    def basis_on_axis(axis):
+        """
+        Return basis.
+
+        One of the basis vectors has components xn, yn, zn. The others are chosen
+        to be as close as possible to the original basis vectors.
+
+        Components of `axis` are normalized to 1, and the component that has the
+        maximal absolute value is positive.
+        """
+        m = max(axis)
+        n = 0      # number of shifts
+        o = 'xyz'  # axis names
+        t = tuple(axis)
+        while t[0] != m:
+            t = (t[2], t[0], t[1])
+            o = o[2] + o[0] + o[1]
+            n += 1
+        xn, yn, zn = t
+        # After i shifts, xn component is the largest one. Now the original basis
+        # vectors are named so, that the largest component of the axis is xn, and
+        # the 1-st basis vector, i', is `axis`. The second basis vector, j', is most
+        # closest to the 2-nd basis vector of the original CS, j. This means, that
+        # j' is perpendicular to i' and on the plane built on i' and j.
+        i = (xn, yn, zn)
+        xn2 = xn**2
+        yn2 = yn**2
+        zn2 = zn**2
+        b = 1.0/((1.0 - yn2)**2 + yn2*(xn2 + zn2))**0.5
+        a = -b*yn
+        j = (xn*a, yn*a + b, zn*a)
+        k = cross_product(i, j)
+
+        # shift vector names back n times
         b = (i, j, k)
-    check_basis(*b)
-    # print 'Basis on axis'
-    # print axis
-    # print b, o
-    return b, o[0]
-
-
-def cross_product(a, b):
-    """
-    Return cross-product [a, b]
-    """
-    x = a[1]*b[2] - a[2]*b[1]
-    y = b[0]*a[2] - b[2]*a[0]
-    z = a[0]*b[1] - a[1]*b[0]
-    return x, y, z
-
-
-def scalar_product(a, b):
-    """
-    Return scalar product (a, b)
-    """
-    return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]
-
-
-def check_basis(i, j, k):
-    """
-    Check that i, j, k are orthogonal, unit and right-hand.
-    """
-    f = '{:16.8e}'
-    vf = f*3
-    il = scalar_product(i, i)
-    jl = scalar_product(j, j)
-    kl = scalar_product(k, k)
-    cmnt = []
-    if not areclose((il, jl, kl, 1.0), atol=1e-7, rtol=None, cmnt=cmnt, name='Basis vector normalization'):
-        for comment in cmnt:
-            print(comment)
-        print('Basis vectors not normal')
-        # raise ValueError('Basis vectors not normal')
-
-    ij = scalar_product(i, j)
-    ik = scalar_product(i, k)
-    jk = scalar_product(j, k)
-    if not areclose((ij, ik, jk, 0.0), atol=1e-7, rtol=None):
-        print('i', vf.format(*i))
-        print('j', vf.format(*j))
-        print('k', vf.format(*k))
-        print('products', vf.format(ij, ik, jk))
-        print('Basis vectors not orthogonal')
-        # raise ValueError('Basis vectors not orthogonal')
-
-    ii = scalar_product(i, cross_product(j, k))
-    jj = scalar_product(j, cross_product(k, i))
-    kk = scalar_product(k, cross_product(i, j))
-    if not areclose((ii, jj, kk, 1.0)):
-        print('(i, [j, k])', f.format(ii))
-        print('(j, [k, i])', f.format(jj))
-        print('(k, [i, j])', f.format(kk))
-        print('Basis vectors not right-hand')
-        # raise ValueError('Basis vectors not right-hand')
-    return i, j, k
-
-
-def transform(p, b, o):
-    """
-    Return coordinates of point `p` in the basis `b`, which center has
-    coordinates `o`.
-
-    Vectors of basis `b` must be orthonormal.
-    """
-    x, y, z = p
-    i, j, k = b
-    X, Y, Z = o
-    dv = (x-X, y-Y, z-Z)
-    xp = scalar_product(dv, i)
-    yp = scalar_product(dv, j)
-    zp = scalar_product(dv, k)
-    return xp, yp, zp
-
-
-if __name__ == '__main__':
-    a = 1.0
-    b = a + 1e-7
-    c = a - 1e-7
-    l = (a, b, c, 0)
-    print l
-    print areclose(l, atol=1e-6)
+        for kkk in range(n):
+            i = b[1][1], b[1][2], b[1][0]
+            j = b[2][1], b[2][2], b[2][0]
+            k = b[0][1], b[0][2], b[0][0]
+            b = (i, j, k)
+        check_basis(*b)
+        # print 'Basis on axis'
+        # print axis
+        # print b, o
+        return b, o[0]
+
+
+    def cross_product(a, b):
+        """
+        Return cross-product [a, b]
+        """
+        x = a[1]*b[2] - a[2]*b[1]
+        y = b[0]*a[2] - b[2]*a[0]
+        z = a[0]*b[1] - a[1]*b[0]
+        return x, y, z
+
+
+    def scalar_product(a, b):
+        """
+        Return scalar product (a, b)
+        """
+        return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]
+
+
+    def check_basis(i, j, k):
+        """
+        Check that i, j, k are orthogonal, unit and right-hand.
+        """
+        f = '{:16.8e}'
+        vf = f*3
+        il = scalar_product(i, i)
+        jl = scalar_product(j, j)
+        kl = scalar_product(k, k)
+        cmnt = []
+        if not areclose((il, jl, kl, 1.0), atol=1e-7, rtol=None, cmnt=cmnt, name='Basis vector normalization'):
+            for comment in cmnt:
+                print(comment)
+            print('Basis vectors not normal')
+            # raise ValueError('Basis vectors not normal')
+
+        ij = scalar_product(i, j)
+        ik = scalar_product(i, k)
+        jk = scalar_product(j, k)
+        if not areclose((ij, ik, jk, 0.0), atol=1e-7, rtol=None):
+            print('i', vf.format(*i))
+            print('j', vf.format(*j))
+            print('k', vf.format(*k))
+            print('products', vf.format(ij, ik, jk))
+            print('Basis vectors not orthogonal')
+            # raise ValueError('Basis vectors not orthogonal')
+
+        ii = scalar_product(i, cross_product(j, k))
+        jj = scalar_product(j, cross_product(k, i))
+        kk = scalar_product(k, cross_product(i, j))
+        if not areclose((ii, jj, kk, 1.0)):
+            print('(i, [j, k])', f.format(ii))
+            print('(j, [k, i])', f.format(jj))
+            print('(k, [i, j])', f.format(kk))
+            print('Basis vectors not right-hand')
+            # raise ValueError('Basis vectors not right-hand')
+        return i, j, k
+
+
+    def transform(p, b, o):
+        """
+        Return coordinates of point `p` in the basis `b`, which center has
+        coordinates `o`.
+
+        Vectors of basis `b` must be orthonormal.
+        """
+        x, y, z = p
+        i, j, k = b
+        X, Y, Z = o
+        dv = (x-X, y-Y, z-Z)
+        xp = scalar_product(dv, i)
+        yp = scalar_product(dv, j)
+        zp = scalar_product(dv, k)
+        return xp, yp, zp
+
+
+    if __name__ == '__main__':
+        a = 1.0
+        b = a + 1e-7
+        c = a - 1e-7
+        l = (a, b, c, 0)
+        print(l)
+        print(areclose(l, atol=1e-6))
diff --git a/numjuggler/parser.py b/numjuggler/parser.py
index 704e8c8..c8ea353 100644
--- a/numjuggler/parser.py
+++ b/numjuggler/parser.py
@@ -8,6 +8,9 @@
 
 import re
 import warnings
+import six
+import os
+from numjuggler import PartialFormatter
 
 try:
     # This clause define the fallback for cPickle, which is an accelerated
@@ -19,21 +22,21 @@
     import pickle as cPickle
 
 # integer with one prefix character
-re_int = re.compile('\D{0,1}\d+')
+re_int = re.compile(r'\D{0,1}\d+')
 
 # interior of square brackets for tally in lattices
-re_ind = re.compile('\[.+\]', flags=re.DOTALL)
+re_ind = re.compile(r'\[.+\]', flags=re.DOTALL)
 
 # repitition syntax of MCNP input file
-re_rpt = re.compile('\d+[ri]', flags=re.IGNORECASE)
+re_rpt = re.compile(r'\d+[ri]', flags=re.IGNORECASE)
 
 # imp or tmp parameters in cell card
-re_prm = re.compile('((imp:n|imp:p|tmp)\s+\S+)')
-re_prm = re.compile('[it]mp:*[npe]*[=\s]+\S+', flags=re.IGNORECASE)
-re_prm = re.compile('([it]mp:*[npe]*[=\s]+)(\S+)', flags=re.IGNORECASE)
+re_prm = re.compile(r'((imp:n|imp:p|tmp)\s+\S+)')
+re_prm = re.compile(r'[it]mp:*[npe]*[=\s]+\S+', flags=re.IGNORECASE)
+re_prm = re.compile(r'([it]mp:*[npe]*[=\s]+)(\S+)', flags=re.IGNORECASE)
 
 # fill keyword
-re_fll = re.compile('\*{0,1}fill[=\s]+', flags=re.IGNORECASE)
+re_fll = re.compile(r'\*{0,1}fill[=\s]+', flags=re.IGNORECASE)
 
 
 # If type specifier not given, any data type can be formatted:
@@ -43,6 +46,7 @@ def fmt_gen(s):
 fmt_g = fmt_gen
 fmt_s = fmt_gen
 
+partial_formmatter = PartialFormatter()
 
 class __CIDClass(object):
     """
@@ -599,7 +603,7 @@ def card(self, wrap=False, comment=True):
                 indent = ' '*5
                 if self.ctype == CID.title:
                     indent = 'c' + indent
-                tparts = re.split('\{.*?\}', self.template)[1:]
+                tparts = re.split(r'\{.*?\}', self.template)[1:]
                 # print 'wrapped inp', repr(self.template)
                 # print 'wrapped spl', repr(tparts)
                 newt = ['']  # new template parts
@@ -645,7 +649,8 @@ def card(self, wrap=False, comment=True):
             else:
                 tmpl = self.template
 
-            card = tmpl.format(*inpt)
+            card = partial_formmatter.format(tmpl, *inpt)
+            # card = tmpl.format(*inpt)
         else:
             card = self.template
         return card
@@ -1150,44 +1155,54 @@ def is_blankline(l):
     """
     return l.strip() == ''
 
+if six.PY2:
+    def get_cards(inp, debug=None):
+        """
+        Check first existence of a dump file
 
-def get_cards(inp, debug=None):
-    """
-    Check first existence of a dump file
+        If dump exists and it is newwer than the input file, read the dump file
+        """
+        from os import stat
+        iname = inp
+        dname = '.{}.~'.format(os.path.basename(inp))
+        try:
+            it = stat(iname).st_mtime
+        except OSError as e:
+            raise e
 
-    If dump exists and it is newwer than the input file, read the dump file
-    """
-    from os import stat
-    iname = inp
-    dname = '.{}.~'.format(inp)
-    try:
-        it = stat(iname).st_mtime
-    except OSError as e:
-        raise e
-
-    try:
-        dt = stat(dname).st_mtime
-    except OSError:
-        # print('No dump file exists')
-        dt = it - 1.0
-    if it < dt and debug is None:
-        # print('Reading from dump')
-        # dump is youger
-        dfile = open(dname, 'r')
-        cl = cPickle.load(dfile)
-        for c in cl:
-            yield c
-    else:
-        # print('Reading from input')
-        cl = []
-        for c in get_cards_from_input(inp, debug=debug):
-            yield c
-            cl.append(c)
+        try:
+            dt = stat(dname).st_mtime
+        except OSError:
+            # print('No dump file exists')
+            dt = it - 1.0
+        if it < dt and debug is None:
+            # print('Reading from dump')
+            # dump is youger
+            dfile = open(dname, 'r')
+            cl = cPickle.load(dfile)
+            for c in cl:
+                yield c
+        else:
+            # print('Reading from input')
+            cl = []
+            for c in get_cards_from_input(inp, debug=debug):
+                yield c
+                cl.append(c)
         if debug is None:
             # otherwise the instances of c contain the file object, which
             # cannot be dumped.
             dfile = open(dname, 'w')
             cPickle.dump(cl, dfile)
+else:
+    def get_cards(inp, debug=None):
+        """
+        Check first existence of a dump file
+
+        If dump exists and it is newwer than the input file, read the dump file
+        """
+        iname = inp
+        for c in get_cards_from_input(inp, debug=debug):
+            yield c
 
 
 def get_cards_from_input(inp, debug=None):
diff --git a/numjuggler/shortener.py b/numjuggler/shortener.py
index 11c6e12..94669cc 100644
--- a/numjuggler/shortener.py
+++ b/numjuggler/shortener.py
@@ -42,6 +42,7 @@ def f(l):
                 iD = e - es
 
 if __name__ == '__main__':
+    pass
 
 
 
diff --git a/numjuggler/string_cells.py b/numjuggler/string_cells.py
index 3df11e1..99b44ba 100644
--- a/numjuggler/string_cells.py
+++ b/numjuggler/string_cells.py
@@ -142,7 +142,7 @@ def complementary(ccell) :
     ccell.remove_comments()
 
     # insert external parenthesis
-    ccell.str= re.sub(leftp,"(\g<left>",ccell.str,count=1)
+    ccell.str= re.sub(leftp,r"(\g<left>",ccell.str,count=1)
     ccell.str= re.sub(rightp,r")\g<right>",ccell.str,count=1)
 
     # insert '&' as intersection operator
@@ -161,7 +161,7 @@ def complementary(ccell) :
     ccell.str=re.sub(number,chgsign,ccell.str)
 
     # insert external parenthesis
-    ccell.str= re.sub(leftp,"(\g<left>",ccell.str,count=1)
+    ccell.str= re.sub(leftp,r"(\g<left>",ccell.str,count=1)
     ccell.str= re.sub(rightp,r")\g<right>",ccell.str,count=1)
 
 # restore original comments
@@ -253,13 +253,13 @@ def remove_redundant(self,remove_com=True,remopt='nochg'):
       geom=re.sub(blnkline,'',geom)
 
       # ensure 5 blanks continous line
-      geom=re.sub(contline,'\n     \g<start>',geom)
+      geom=re.sub(contline, r'\n     \g<start>',geom)
 
       if remopt != 'all' :
         # add parenthesis to set geom as MCNP complex cell
         if geom.find(':') == -1 and geom.find('#') == -1 :
-           geom=re.sub(startgeom,'\g<previous>(\g<start>',geom)
-           geom=re.sub(endgeom,'\g<last>)\g<next>',geom)
+           geom=re.sub(startgeom,r'\g<previous>(\g<start>',geom)
+           geom=re.sub(endgeom,r'\g<last>)\g<next>',geom)
 
       # restore original comments
       self.str = geom
@@ -267,7 +267,7 @@ def remove_redundant(self,remove_com=True,remopt='nochg'):
       if remove_com: self.restore_comments()
 
       # subtitute comment $ with  blank line
-      self.str=re.sub(comdollar,'\nC\g<blnk>',self.str)
+      self.str=re.sub(comdollar,r'\nC\g<blnk>',self.str)
       pdiff = [x-y for x,y in zip(pmod,porg)]
       self.removedp=pdiff
       return
@@ -385,10 +385,10 @@ def get_lines(self):
       card=re.sub(blnkline,'',card)
 
       # subtitute comment $ with  blank line
-      card=re.sub(comdollar,'\nC\g<blnk>',card)
+      card=re.sub(comdollar, r'\nC\g<blnk>',card)
 
       # ensure 5 blanks continous line
-      card=re.sub(contline,'\n     \g<start>',card)
+      card=re.sub(contline, r'\n     \g<start>',card)
 
       if (card[-1] == '\n') : card=card[:-1]
       return  map(lambda x: x+'\n' ,card.split('\n'))
diff --git a/numjuggler/utils/PartialFormatter.py b/numjuggler/utils/PartialFormatter.py
new file mode 100644
index 0000000..5f977de
--- /dev/null
+++ b/numjuggler/utils/PartialFormatter.py
@@ -0,0 +1,102 @@
+from __future__ import print_function, absolute_import, division
+
+import six
+import string
+
+
+def make_label(item):
+    if 0 < len(item):
+        return '{' + item + '}'
+    else:
+        return item
+
+
+class SafeDict(dict):
+
+    def __getitem__(self, item):
+        return super(SafeDict, self).__getitem__(item) or ''
+
+    def __missing__(self, key):
+        return make_label(key)
+
+
+# ' 1 #'
+class PartialFormatter(string.Formatter):
+
+    def format(*args, **kwargs):
+        if not args:
+            raise TypeError("descriptor 'format' of 'Formatter' object "
+                            "needs an argument")
+        self, args = args[0], args[1:]  # allow the "self" keyword be passed
+        try:
+            format_string, args = args[0], args[1:] # allow the "format_string" keyword be passed
+        except IndexError:
+            if 'format_string' in kwargs:
+                format_string = kwargs.pop('format_string')
+            else:
+                raise TypeError("format() missing 1 required positional "
+                                "argument: 'format_string'")
+        return string.Formatter.vformat(self, format_string, args, SafeDict(kwargs))
+
+    def _vformat(self, format_string, args, kwargs, used_args, recursion_depth, auto_arg_index=0):
+        """Clone from Python3 version to fix Python2 mess with unnamed {} format specifiers"""
+        if recursion_depth < 0:
+            raise ValueError('Max string recursion exceeded')
+        result = []
+        for literal_text, field_name, format_spec, conversion in \
+                self.parse(format_string):
+
+            # output the literal text
+            if literal_text:
+                result.append(literal_text)
+
+            # if there's a field, output it
+            if field_name is not None:
+                # this is some markup, find the object and do
+                #  the formatting
+
+                # handle arg indexing when empty field_names are given.
+                if field_name == '':
+                    if auto_arg_index is False:
+                        raise ValueError('cannot switch from manual field '
+                                         'specification to automatic field '
+                                         'numbering')
+                    field_name = str(auto_arg_index)
+                    auto_arg_index += 1
+                elif field_name.isdigit():
+                    if auto_arg_index:
+                        raise ValueError('cannot switch from manual field '
+                                         'specification to automatic field '
+                                         'numbering')
+                    # disable auto arg incrementing, if it gets
+                    # used later on, then an exception will be raised
+                    auto_arg_index = False
+
+                # given the field_name, find the object it references
+                #  and the argument it came from
+                obj, arg_used = self.get_field(field_name, args, kwargs)
+                used_args.add(arg_used)
+
+                # do any conversion on the resulting object
+                obj = self.convert_field(obj, conversion)
+
+                # expand the format spec, if needed
+                if six.PY2:
+                    format_spec = self._vformat(
+                        format_spec, args, kwargs,
+                        used_args, recursion_depth-1,
+                        auto_arg_index=auto_arg_index)
+                else:
+                    format_spec, auto_arg_index = self._vformat(
+                        format_spec, args, kwargs,
+                        used_args, recursion_depth-1,
+                        auto_arg_index=auto_arg_index)
+                # format the object and append to the result
+                result.append(self.format_field(obj, format_spec))
+
+        result = ''.join(result)
+
+        if six.PY2:
+            return result
+        else:
+            return result, auto_arg_index
diff --git a/numjuggler/utils/__init__.py b/numjuggler/utils/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/numjuggler/utils/io.py b/numjuggler/utils/io.py
new file mode 100644
index 0000000..20e6d13
--- /dev/null
+++ b/numjuggler/utils/io.py
@@ -0,0 +1,30 @@
+import os
+import sys
+
+from .resource import Path
+from contextlib import contextmanager
+
+
+@contextmanager
+def cd_temporarily(cd_to):
+    cur_dir = str(Path.cwd())
+    try:
+        os.chdir(str(cd_to))
+        yield
+    finally:
+        os.chdir(cur_dir)
+
+
+@contextmanager
+def resolve_fname_or_stream(fname_or_stream, mode="r"):
+    is_input = mode == 'r'
+    if fname_or_stream is None:
+        if is_input:
+            yield sys.stdin
+        else:
+            yield sys.stdout
+    elif is_input and hasattr(fname_or_stream, "read") or not is_input and hasattr(fname_or_stream, "write"):
+        yield fname_or_stream
+    else:
+        with open(fname_or_stream, mode=mode) as fid:
+            yield fid
diff --git a/numjuggler/utils/resource.py b/numjuggler/utils/resource.py
new file mode 100644
index 0000000..09a1178
--- /dev/null
+++ b/numjuggler/utils/resource.py
@@ -0,0 +1,41 @@
+import pkg_resources as pkg
+import inspect
+
+# noinspection PyCompatibility
+try:
+    from pathlib import Path
+except:
+    from pathlib2 import Path
+
+
+def filename_resolver(package=None):
+    if package is None:
+        caller_package = inspect.getmodule(inspect.stack()[1][0]).__name__
+        package = caller_package
+
+    resource_manager = pkg.ResourceManager()
+
+    def func(resource):
+        return resource_manager.resource_filename(package, resource)
+
+    # noinspection PyCompatibility
+    func.__doc__ = "Computes file names for resources located in {}".format(package)
+
+    return func
+
+
+def path_resolver(package=None):
+
+    resolver = filename_resolver(package)
+
+    def func(resource):
+        filename = resolver(resource)
+        return Path(filename)
+
+    if package is None:
+        package = "caller package"
+
+    # noinspection PyCompatibility
+    func.__doc__ = "Computes Path for resources located in {}".format(package)
+
+    return func
diff --git a/pytest.ini b/pytest.ini
new file mode 100644
index 0000000..ae660f4
--- /dev/null
+++ b/pytest.ini
@@ -0,0 +1,14 @@
+[pytest]
+#markers (linelist)       markers for test functions
+norecursedirs = .* build dist {arch} *.egg adhoc examples notebooks experiment dev_tests data out wrk
+testpaths=tests
+#usefixtures (args)       list of default fixtures to be used with this project
+python_files=test*.py    # (args)      glob-style file patterns for Python test module discovery
+python_classes=Test*     # (args)    prefixes or glob names for Python test class discovery
+python_functions=test*   # (args)  prefixes or glob names for Python test function and method discovery
+#xfail_strict (bool)      default for the strict parameter of xfail markers when not given explicitly (default: Fals
+doctest_optionflags = NORMALIZE_WHITESPACE IGNORE_EXCEPTION_DETAIL ALLOW_UNICODE ALLOW_BYTES
+addopts = --ignore setup.py --ignore *_tab.py --doctest-modules --color=yes
+# coverage doesn't allow to work with PyCharm debugger, run test_coverage.sh script to update coverage
+#addopts = --ignore setup.py --doctest-modules --doctest-glob='*.rst' --cov-report term-missing --cov m2t
+#minversion (string)      minimally required pytest version
diff --git a/pytest_coverage b/pytest_coverage
new file mode 100755
index 0000000..2454e9e
--- /dev/null
+++ b/pytest_coverage
@@ -0,0 +1,4 @@
+#!/usr/bin/env bash
+
+pytest --cov-report=term-missing:skip-covered --cov=numjuggler $*
+
diff --git a/run_coverage b/run_coverage
new file mode 100755
index 0000000..d03db57
--- /dev/null
+++ b/run_coverage
@@ -0,0 +1,6 @@
+#!/usr/bin/env bash
+
+coverage run -m pytest tests
+coverage report --show-missing --skip-covered
+coverage html
+
diff --git a/setup.cfg b/setup.cfg
new file mode 100644
index 0000000..a662eae
--- /dev/null
+++ b/setup.cfg
@@ -0,0 +1,2 @@
+[aliaces]
+test=pytest
diff --git a/setup.py b/setup.py
index d662bd7..65d33bd 100644
--- a/setup.py
+++ b/setup.py
@@ -1,19 +1,46 @@
-# from distutils.core import setup
-from setuptools import setup
-
-# Get the version from numjuggler/__init__.py
-fd = {}
-with open('./numjuggler/__version__.py', 'r') as f:
-    exec(f.read(), fd)
-    version = fd['__version__']
-
-
-setup(name='numjuggler',
-      version=version,  # __version__,  # '2.41a.27',
-      description='MCNP input file renumbering tool',
-      author='A. Travleev',
-      author_email='anton.travleev@gmail.com',
-      packages=['numjuggler', ],
-      # scripts = ['numjuggler/numjuggler'],
-      entry_points={'console_scripts': ['numjuggler = numjuggler.main:main']},
-      )
+import sys
+from setuptools import setup, find_packages
+from setuptools.command.test import test as test_command
+
+
+# noinspection PyAttributeOutsideInit
+class PyTest(test_command):
+    """
+    See recomendations at https://docs.pytest.org/en/latest/goodpractices.html
+    """
+    user_options = [('pytest-args=', 'a', "Arguments to pass to pytest")]
+
+    def initialize_options(self):
+        test_command.initialize_options(self)
+        self.pytest_args = ""
+
+    def run_tests(self):
+        import shlex
+        # import here, cause outside the eggs aren't loaded
+        import pytest
+        errno = pytest.main(shlex.split(self.pytest_args))
+        sys.exit(errno)
+
+
+def load_version():
+    fd = {}
+    with open('./numjuggler/__version__.py', 'r') as f:
+        exec(f.read(), fd)
+        return fd['__version__']
+
+
+packages = find_packages(
+    include=("numjuggler", "numjuggler.*",),
+)
+
+setup(
+    name='numjuggler',
+    version=load_version(),
+    description='MCNP input file renumbering tool',
+    author='A.Travleev',
+    author_email='anton.travleev@gmail.com',
+    packages=packages,
+    tests_require=['pytest', 'pytest-cov>=2.3.1'],
+    cmdclass={'test': PyTest},
+    entry_points={'console_scripts': ['numjuggler = numjuggler.main:main']},
+)
diff --git a/tests/__init__.py b/tests/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/data/simple_cubes.mcnp b/tests/data/simple_cubes.mcnp
new file mode 100644
index 0000000..d5c39af
--- /dev/null
+++ b/tests/data/simple_cubes.mcnp
@@ -0,0 +1,49 @@
+testing integration scripts on simple cubic cells
+c 
+c given two envelopes and fillers insert new envelop and filler intersecting
+c the original fillers
+c
+1 0 -1 : 2 : -3 : 4 : -5 : 6 imp:n=0
+        $ outer space
+2 0  1 -2  3 -7  5 -6 imp:n=1 fill=1
+        $ envelop #1
+3 0  1 -2  7 -4  5 -6 imp:n=1 fill=2
+        $ envelop #2
+4 1 1.0  20 -21  22 -23  24 -25 imp:n=1 u=1
+        $ filler #1 body
+5 0     -20 : 21 : -22 : 23 : -24 : 25 imp:n=1 u=1
+        $ filler #1 outer space
+6 2 1.0  30 -31  32 -33  34 -35 imp:n=1 u=2
+        $ filler #2 body
+7 0     -30 : 31 : -32 : 33 :-34 : 35 imp:n=1 u=2
+        $ filler #2 outer space
+
+c envelopes
+1  px -50
+2  px  50
+3  py -50
+4  py  50
+5  pz -50
+6  pz  50
+7  py   0
+c filler #1
+20 px -25
+21 px  25
+22 py -25
+23 py   0
+24 pz -25
+25 pz  25
+c filler #2
+30 px -35
+31 px  35
+32 py   0
+33 py  35
+34 pz -35
+35 pz  35
+
+m1 1001.31c 1.0
+m2 2004.31c 1.0
+c
+c Comment may contain entries with braces {31c}
+c
+mode n
diff --git a/tests/test_likefunc.py b/tests/test_likefunc.py
new file mode 100644
index 0000000..a5a0a66
--- /dev/null
+++ b/tests/test_likefunc.py
@@ -0,0 +1,54 @@
+import pytest
+from numjuggler.utils.io import cd_temporarily
+from six import StringIO
+import numjuggler.likefunc as lf
+
+
+
+@pytest.mark.parametrize("data, log, present, absent, expected_present_value, expected_absent_value, expected_text", [
+    (
+        """
+            c 1: 12
+            c 2: 14
+        """,
+        False,
+        1, 3,
+        12, 3,
+        "",
+    ),
+    (
+        """
+            c 1: 12
+            c 2: 14
+        """,
+        True,
+        1, 3,
+        12, 3,
+        "cel 12: 1\ncel 3: 3\n",
+    ),
+])
+def test_LikeFunction(
+    data,
+    log,
+    present,
+    absent,
+    expected_present_value,
+    expected_absent_value,
+    expected_text
+):
+    input = StringIO(data)
+    maps = lf.read_map_file(input, log)
+    actual = StringIO()
+    for k in maps:
+        like_function = maps[k]
+        assert like_function.log == log
+        present_value = like_function(present)
+        assert present_value == expected_present_value
+        absent_value = like_function(absent)
+        assert absent_value == expected_absent_value
+        if log:
+            like_function.write_log_as_map(k, actual)
+        else:
+            with pytest.raises(ValueError):
+                like_function.write_log_as_map(k, actual)
+    assert actual.getvalue() == expected_text
diff --git a/tests/test_main.py b/tests/test_main.py
new file mode 100644
index 0000000..9c6b93b
--- /dev/null
+++ b/tests/test_main.py
@@ -0,0 +1,50 @@
+# -*- coding: utf-8 -*-
+from __future__ import print_function, division, nested_scopes
+
+import pytest
+import six
+from numjuggler.utils.resource import path_resolver
+from numjuggler.utils.io import cd_temporarily
+from numjuggler.main import main
+
+test_data_path = path_resolver('tests')('data')
+assert test_data_path.exists(), "Cannot access test data files"
+
+
+def load_line_heading_numbers(lines):
+    res = []
+    for line in lines:
+        line = six.ensure_str(line, encoding='utf8')
+        if line and str.isdigit(line[0]):
+            if six.PY2:
+                card_no = int(line.split()[0])
+            else:
+                card_no = int(line.split(maxsplit=2)[0])
+            res.append(card_no)
+    return res
+
+
+@pytest.mark.parametrize("inp,command,expected", [
+    (
+        'simple_cubes.mcnp',
+        "-c 10",
+        "11 12 13 14 15 16 17 1 2 3 4 5 6 7 20 21 22 23 24 25 30 31 32 33 34 35"
+    ),
+    (
+        'simple_cubes.mcnp',
+        "-c 20 -s 10",
+        "21 22 23 24 25 26 27 11 12 13 14 15 16 17 30 31 32 33 34 35 40 41 42 43 44 45"
+    ),
+])
+def test_test_main(tmpdir, capsys, inp, command, expected):
+    source = test_data_path / inp
+    command = command.split()
+    command.append(str(source))
+    with cd_temporarily(tmpdir):
+        main(command)
+    out, err = capsys.readouterr()
+    actual_numbers = load_line_heading_numbers(out.split('\n'))
+    expected_numbers = list(f for f in map(int, expected.split()))
+    assert expected_numbers == actual_numbers, "Output of numjuggler is wrong"
+
+