diff --git a/bfps/fluid_converter.py b/bfps/FluidConvert.py
similarity index 77%
rename from bfps/fluid_converter.py
rename to bfps/FluidConvert.py
index 04f113c4ed44d99fb020d4a782bb7d28ae09adbb..cbb9d9c3f478463426afbb63c18bb26abba49285 100644
--- a/bfps/fluid_converter.py
+++ b/bfps/FluidConvert.py
@@ -24,22 +24,26 @@
 
 
 
-import bfps
-import bfps.fluid_base
-import bfps.tools
 import numpy as np
 import pickle
 import os
+from ._fluid_base import _fluid_particle_base
 
-class fluid_converter(bfps.fluid_base.fluid_particle_base):
+class FluidConvert(_fluid_particle_base):
+    """This class is meant to be used for conversion of native DNS field
+    representations to real-space representations of velocity/vorticity
+    fields.
+    It may be superseeded by streamlined functionality in the future...
+    """
     def __init__(
             self,
-            name = 'fluid_converter',
+            name = 'FluidConvert',
             work_dir = './',
             simname = 'test',
             fluid_precision = 'single',
             use_fftw_wisdom = True):
-        super(fluid_converter, self).__init__(
+        _fluid_particle_base.__init__(
+                self,
                 name = name,
                 work_dir = work_dir,
                 simname = simname,
@@ -86,4 +90,24 @@ class fluid_converter(bfps.fluid_base.fluid_particle_base):
                 """
         self.fluid_end += 'delete fs;\n'
         return None
+    def specific_parser_arguments(
+            self,
+            parser):
+        _fluid_particle_base.specific_parser_arguments(self, parser)
+        parser.add_argument(
+                '--src-wd',
+                type = str,
+                dest = 'src_work_dir',
+                default = './')
+        parser.add_argument(
+                '--src-simname',
+                type = str,
+                dest = 'src_simname',
+                default = '')
+        parser.add_argument(
+                '--src-iteration',
+                type = int,
+                dest = 'src_iteration',
+                default = 0)
+        return None
 
diff --git a/bfps/fluid_resize.py b/bfps/FluidResize.py
similarity index 68%
rename from bfps/fluid_resize.py
rename to bfps/FluidResize.py
index 5c4e6c5ed001c9ca698d19824677c335a6a0ff53..0cb0f561163ed3a7f080da5157562b90ab2c26dd 100644
--- a/bfps/fluid_resize.py
+++ b/bfps/FluidResize.py
@@ -24,20 +24,27 @@
 
 
 
-import bfps
-import bfps.fluid_base
-
+import os
+import argparse
 import numpy as np
 
-class fluid_resize(bfps.fluid_base.fluid_particle_base):
+import bfps
+from ._fluid_base import _fluid_particle_base
+
+class FluidResize(_fluid_particle_base):
+    """This class is meant to resize snapshots of DNS states to new grids.
+    Typical stuff for DNS of turbulence.
+    It will become superfluous when HDF5 is used for field I/O.
+    """
     def __init__(
             self,
-            name = 'fluid_resize',
+            name = 'FluidResize',
             work_dir = './',
             simname = 'test',
             dtype = np.float32,
             use_fftw_wisdom = False):
-        super(fluid_resize, self).__init__(
+        _fluid_particle_base.__init__(
+                self,
                 name = name,
                 work_dir = work_dir,
                 simname = simname,
@@ -100,4 +107,50 @@ class fluid_resize(bfps.fluid_base.fluid_particle_base):
                 //endcpp
                 """
         return None
+    def specific_parser_arguments(
+            self,
+            parser):
+        _fluid_particle_base.specific_parser_arguments(self, parser)
+        parser.add_argument(
+                '-m',
+                type = int,
+                dest = 'm',
+                default = 32,
+                metavar = 'M',
+                help = 'resize from N to M')
+        parser.add_argument(
+                '--src_wd',
+                type = str,
+                dest = 'src_work_dir',
+                required = True)
+        parser.add_argument(
+                '--src_iteration',
+                type = int,
+                dest = 'src_iteration',
+                required = True)
+        return None
+    def launch(
+            self,
+            args = [],
+            **kwargs):
+        parser = argparse.ArgumentParser('bfps ' + type(self).__name__)
+        self.add_parser_arguments(parser)
+        opt = parser.parse_args(args)
+        cmd_line_pars = vars(opt)
+        for k in ['dst_nx', 'dst_ny', 'dst_nz']:
+            if type(cmd_line_pars[k]) == type(None):
+                cmd_line_pars[k] = opt.m
+        self.pars_from_namespace(opt)
+        src_file = os.path.join(
+                opt.src_work_dir,
+                opt.src_simname + '_cvorticity_i{0:0>5x}'.format(opt.src_iteration))
+        read_file = os.path.join(
+                self.work_dir,
+                opt.src_simname + '_cvorticity_i{0:0>5x}'.format(opt.src_iteration))
+        if not os.path.exists(read_file):
+            os.symlink(src_file, read_file)
+        self.set_host_info(bfps.host_info)
+        self.write_par(iter0 = opt.src_iteration)
+        self.run(ncpu = opt.ncpu)
+        return None
 
diff --git a/bfps/Launcher.py b/bfps/Launcher.py
deleted file mode 100644
index ad09a1f22188426b99c94b799c6a82e600fdedf4..0000000000000000000000000000000000000000
--- a/bfps/Launcher.py
+++ /dev/null
@@ -1,156 +0,0 @@
-import os
-
-import argparse
-import bfps
-
-class Launcher:
-    def __init__(self):
-        self.base_class = bfps.NavierStokes
-        self.parser = argparse.ArgumentParser(prog = 'bfps')
-        self.parser.add_argument(
-                '-v', '--version',
-                action = 'version',
-                version = '%(prog)s ' + bfps.__version__)
-        self.parser.add_argument(
-                '-n',
-                type = int,
-                dest = 'n',
-                default = 32,
-                metavar = 'N',
-                help = 'code is run by default in a grid of NxNxN')
-        self.parser.add_argument(
-                '--run',
-                dest = 'run',
-                action = 'store_true')
-        self.parser.add_argument(
-                '--ncpu',
-                type = int, dest = 'ncpu',
-                default = 2)
-        self.parser.add_argument(
-                '--precision',
-                type = str, dest = 'precision',
-                default = 'single')
-        self.parser.add_argument(
-                '--simname',
-                type = str, dest = 'simname',
-                default = 'test')
-        self.parser.add_argument(
-                '--wd',
-                type = str, dest = 'work_dir',
-                default = './')
-        self.parser.add_argument(
-                '--njobs',
-                type = int, dest = 'njobs',
-                default = 1)
-        self.parser.add_argument(
-                '--QR-stats',
-                action = 'store_true',
-                dest = 'QR_stats',
-                help = 'add this option if you want to compute velocity gradient and QR stats')
-        self.parser.add_argument(
-                '--kMeta',
-                type = float,
-                dest = 'kMeta',
-                default = 2.0)
-        self.parser.add_argument(
-                '--dtfactor',
-                type = float,
-                dest = 'dtfactor',
-                default = 0.5,
-                help = 'dt is computed as DTFACTOR / N')
-        self.parser.add_argument(
-                '--environment',
-                type = str,
-                dest = 'environment',
-                default = '')
-        self.parser.add_argument(
-                '--src-simname',
-                type = str,
-                dest = 'src_simname',
-                default = '')
-        self.parser.add_argument(
-                '--src-iteration',
-                type = int,
-                dest = 'src_iteration',
-                default = 0)
-        self.parser.add_argument(
-                '--particle-rand-seed',
-                type = int,
-                dest = 'particle_rand_seed',
-                default = None)
-        c = self.base_class()
-        for k in sorted(c.parameters.keys()):
-            self.parser.add_argument(
-                    '--{0}'.format(k),
-                    type = type(c.parameters[k]),
-                    dest = k,
-                    default = None)
-        return None
-    def __call__(
-            self,
-            args = None):
-        opt = self.parser.parse_args(args)
-        if opt.environment != '':
-            bfps.host_info['environment'] = opt.environment
-        opt.work_dir = os.path.join(
-                os.path.realpath(opt.work_dir),
-                'N{0:0>4}'.format(opt.n))
-        c = self.base_class(
-                work_dir = opt.work_dir,
-                fluid_precision = opt.precision,
-                simname = opt.simname,
-                QR_stats_on = opt.QR_stats)
-        # with the default Lundgren forcing, I can estimate the dissipation
-        # with nondefault forcing, figure out the amplitude for this viscosity
-        # yourself
-        c.parameters['nu'] = (opt.kMeta * 2 / opt.n)**(4./3)
-        c.parameters['dt'] = (opt.dtfactor / opt.n)
-        if ((c.parameters['niter_todo'] % c.parameters['niter_out']) != 0):
-            c.parameters['niter_out'] = c.parameters['niter_todo']
-        if c.QR_stats_on:
-            # max_Q_estimate and max_R_estimate are just used for the 2D pdf
-            # therefore I just want them to be small multiples of mean trS2
-            # I'm already estimating the dissipation with kMeta...
-            meantrS2 = (opt.n//2 / opt.kMeta)**4 * c.parameters['nu']**2
-            c.parameters['max_Q_estimate'] = meantrS2
-            c.parameters['max_R_estimate'] = .4*meantrS2**1.5
-
-        # command line parameters will overwrite any defaults
-        cmd_line_pars = vars(opt)
-        for k in ['nx', 'ny', 'nz']:
-            if type(cmd_line_pars[k]) == type(None):
-                cmd_line_pars[k] = opt.n
-        for k in c.parameters.keys():
-            if k in cmd_line_pars.keys():
-                if not type(cmd_line_pars[k]) == type(None):
-                    c.parameters[k] = cmd_line_pars[k]
-        c.fill_up_fluid_code()
-        c.finalize_code()
-        c.write_src()
-        c.set_host_info(bfps.host_info)
-        if opt.run:
-            if not os.path.exists(os.path.join(c.work_dir, c.simname + '.h5')):
-                c.write_par()
-                if c.parameters['nparticles'] > 0:
-                    data = c.generate_tracer_state(species = 0, rseed = opt.particle_rand_seed)
-                    for s in range(1, c.particle_species):
-                        c.generate_tracer_state(species = s, data = data)
-                init_condition_file = os.path.join(
-                        c.work_dir,
-                        c.simname + '_cvorticity_i{0:0>5x}'.format(0))
-                if not os.path.exists(init_condition_file):
-                    if len(opt.src_simname) > 0:
-                        src_file = os.path.join(
-                                c.work_dir,
-                                opt.src_simname + '_cvorticity_i{0:0>5x}'.format(opt.src_iteration))
-                        os.symlink(src_file, init_condition_file)
-                    else:
-                       c.generate_vector_field(
-                               write_to_file = True,
-                               spectra_slope = 2.0,
-                               amplitude = 0.25)
-            c.run(ncpu = opt.ncpu,
-                  njobs = opt.njobs)
-        return c
-
-
diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index aa11892d21ecb1e2af658a3d7885c81e0d84d44b..7ccc7ed2b6dbf8ca706c6e40a54ccfd10ae6462d 100644
--- a/bfps/NavierStokes.py
+++ b/bfps/NavierStokes.py
@@ -27,12 +27,17 @@
 import os
 import numpy as np
 import h5py
+import argparse
 
 import bfps
-import bfps.fluid_base
-import bfps.tools
+from ._fluid_base import _fluid_particle_base
 
-class NavierStokes(bfps.fluid_base.fluid_particle_base):
+class NavierStokes(_fluid_particle_base):
+    """Objects of this class can be used to generate production DNS codes.
+    Any functionality that users require should be available through this class,
+    in the sense that they can implement whatever they need by simply inheriting
+    this class.
+    """
     def __init__(
             self,
             name = 'NavierStokes',
@@ -46,12 +51,26 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
         self.QR_stats_on = QR_stats_on
         self.frozen_fields = frozen_fields
         self.fftw_plan_rigor = fftw_plan_rigor
-        super(NavierStokes, self).__init__(
+        _fluid_particle_base.__init__(
+                self,
                 name = name,
                 work_dir = work_dir,
                 simname = simname,
                 dtype = fluid_precision,
                 use_fftw_wisdom = use_fftw_wisdom)
+        self.parameters['nu'] = 0.1
+        self.parameters['fmode'] = 1
+        self.parameters['famplitude'] = 0.5
+        self.parameters['fk0'] = 1.5
+        self.parameters['fk1'] = 3.0
+        self.parameters['forcing_type'] = 'linear'
+        self.parameters['histogram_bins'] = 256
+        self.parameters['max_velocity_estimate'] = 1.0
+        self.parameters['max_vorticity_estimate'] = 1.0
+        self.parameters['QR2D_histogram_bins'] = 64
+        self.parameters['max_trS2_estimate'] = 1.0
+        self.parameters['max_Q_estimate'] = 1.0
+        self.parameters['max_R_estimate'] = 1.0
         self.file_datasets_grow = """
                 //begincpp
                 std::string temp_string;
@@ -715,10 +734,87 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
                          self.parameters['nx']//2+1,
                          3))
     def write_par(self, iter0 = 0):
-        super(NavierStokes, self).write_par(iter0 = iter0)
+        _fluid_particle_base.write_par(self, iter0 = iter0)
         with h5py.File(os.path.join(self.work_dir, self.simname + '.h5'), 'r+') as ofile:
             kspace = self.get_kspace()
             nshells = kspace['nshell'].shape[0]
+            for k in ['velocity', 'vorticity']:
+                time_chunk = 2**20//(8*3*self.parameters['nx'])
+                time_chunk = max(time_chunk, 1)
+                ofile.create_dataset('statistics/xlines/' + k,
+                                     (1, self.parameters['nx'], 3),
+                                     chunks = (time_chunk, self.parameters['nx'], 3),
+                                     maxshape = (None, self.parameters['nx'], 3),
+                                     dtype = self.dtype,
+                                     compression = 'gzip')
+                time_chunk = 2**20//(8*3*3*nshells)
+                time_chunk = max(time_chunk, 1)
+                ofile.create_dataset('statistics/spectra/' + k + '_' + k,
+                                     (1, nshells, 3, 3),
+                                     chunks = (time_chunk, nshells, 3, 3),
+                                     maxshape = (None, nshells, 3, 3),
+                                     dtype = np.float64,
+                                     compression = 'gzip')
+                time_chunk = 2**20//(8*4*10)
+                time_chunk = max(time_chunk, 1)
+                a = ofile.create_dataset('statistics/moments/' + k,
+                                     (1, 10, 4),
+                                     chunks = (time_chunk, 10, 4),
+                                     maxshape = (None, 10, 4),
+                                     dtype = np.float64,
+                                     compression = 'gzip')
+                time_chunk = 2**20//(8*4*self.parameters['histogram_bins'])
+                time_chunk = max(time_chunk, 1)
+                ofile.create_dataset('statistics/histograms/' + k,
+                                     (1,
+                                      self.parameters['histogram_bins'],
+                                      4),
+                                     chunks = (time_chunk,
+                                               self.parameters['histogram_bins'],
+                                               4),
+                                     maxshape = (None,
+                                                 self.parameters['histogram_bins'],
+                                                 4),
+                                     dtype = np.int64,
+                                     compression = 'gzip')
+            for s in range(self.particle_species):
+                time_chunk = 2**20 // (8*3*
+                                       self.parameters['nparticles']*
+                                       self.parameters['tracers{0}_integration_steps'.format(s)])
+                time_chunk = max(time_chunk, 1)
+                ofile.create_dataset('particles/tracers{0}/rhs'.format(s),
+                                     (1,
+                                      self.parameters['tracers{0}_integration_steps'.format(s)],
+                                      self.parameters['nparticles'],
+                                      3),
+                                     maxshape = (None,
+                                                 self.parameters['tracers{0}_integration_steps'.format(s)],
+                                                 self.parameters['nparticles'],
+                                                 3),
+                                     chunks =  (time_chunk,
+                                                self.parameters['tracers{0}_integration_steps'.format(s)],
+                                                self.parameters['nparticles'],
+                                                3),
+                                     dtype = np.float64)
+                time_chunk = 2**20 // (8*3*self.parameters['nparticles'])
+                time_chunk = max(time_chunk, 1)
+                ofile.create_dataset(
+                    '/particles/tracers{0}/velocity'.format(s),
+                    (1,
+                     self.parameters['nparticles'],
+                     3),
+                    chunks = (time_chunk, self.parameters['nparticles'], 3),
+                    maxshape = (None, self.parameters['nparticles'], 3),
+                    dtype = np.float64)
+                if self.parameters['tracers{0}_acc_on'.format(s)]:
+                    ofile.create_dataset(
+                        '/particles/tracers{0}/acceleration'.format(s),
+                        (1,
+                         self.parameters['nparticles'],
+                         3),
+                        chunks = (time_chunk, self.parameters['nparticles'], 3),
+                        maxshape = (None, self.parameters['nparticles'], 3),
+                        dtype = np.float64)
             if self.QR_stats_on:
                 time_chunk = 2**20//(8*3*self.parameters['histogram_bins'])
                 time_chunk = max(time_chunk, 1)
@@ -827,33 +923,107 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
         self.fluid_start += update_fields
         self.fluid_loop += update_fields
         return None
+    def specific_parser_arguments(
+            self,
+            parser):
+        _fluid_particle_base.specific_parser_arguments(self, parser)
+        parser.add_argument(
+                '--src-wd',
+                type = str,
+                dest = 'src_work_dir',
+                default = './')
+        parser.add_argument(
+                '--src-simname',
+                type = str,
+                dest = 'src_simname',
+                default = '')
+        parser.add_argument(
+                '--src-iteration',
+                type = int,
+                dest = 'src_iteration',
+                default = 0)
+        parser.add_argument(
+               '--precision',
+               type = str, dest = 'precision',
+               default = 'single')
+        parser.add_argument(
+               '--njobs',
+               type = int, dest = 'njobs',
+               default = 1)
+        parser.add_argument(
+               '--QR-stats',
+               action = 'store_true',
+               dest = 'QR_stats',
+               help = 'add this option if you want to compute velocity gradient and QR stats')
+        parser.add_argument(
+               '--kMeta',
+               type = float,
+               dest = 'kMeta',
+               default = 2.0)
+        parser.add_argument(
+               '--dtfactor',
+               type = float,
+               dest = 'dtfactor',
+               default = 0.5,
+               help = 'dt is computed as DTFACTOR / N')
+        parser.add_argument(
+               '--particle-rand-seed',
+               type = int,
+               dest = 'particle_rand_seed',
+               default = None)
+        return None
+    def launch(
+            self,
+            args = [],
+            **kwargs):
+        # with the default Lundgren forcing, I can estimate the dissipation
+        # with nondefault forcing, figure out the amplitude for this viscosity
+        # yourself
+        parser = argparse.ArgumentParser('bfps ' + type(self).__name__)
+        self.add_parser_arguments(parser)
+        opt = parser.parse_args(args)
+        self.QR_stats_on = opt.QR_stats
+        self.parameters['nu'] = (opt.kMeta * 2 / opt.n)**(4./3)
+        self.parameters['dt'] = (opt.dtfactor / opt.n)
+        if ((self.parameters['niter_todo'] % self.parameters['niter_out']) != 0):
+            self.parameters['niter_out'] = self.parameters['niter_todo']
+        if self.QR_stats_on:
+            # max_Q_estimate and max_R_estimate are just used for the 2D pdf
+            # therefore I just want them to be small multiples of mean trS2
+            # I'm already estimating the dissipation with kMeta...
+            meantrS2 = (opt.n//2 / opt.kMeta)**4 * self.parameters['nu']**2
+            self.parameters['max_Q_estimate'] = meantrS2
+            self.parameters['max_R_estimate'] = .4*meantrS2**1.5
 
-def launch(
-        opt,
-        nu = None):
-    c = NavierStokes(work_dir = opt.work_dir)
-    assert((opt.nsteps % 4) == 0)
-    c.parameters['nx'] = opt.n
-    c.parameters['ny'] = opt.n
-    c.parameters['nz'] = opt.n
-    if type(nu) == type(None):
-        c.parameters['nu'] = 5.5*opt.n**(-4./3)
-    else:
-        c.parameters['nu'] = nu
-    c.parameters['dt'] = 5e-3 * (64. / opt.n)
-    c.parameters['niter_todo'] = opt.nsteps
-    c.parameters['niter_part'] = 2
-    c.parameters['famplitude'] = 0.2
-    c.parameters['nparticles'] = 32
-    if opt.particles:
-        c.add_particles()
-        c.add_particles(kcut = 'fs->kM/2')
-    c.finalize_code()
-    c.write_src()
-    c.write_par()
-    if opt.run:
-        if opt.iteration == 0 and opt.initialize:
-            c.generate_initial_condition()
-        c.run(ncpu = opt.ncpu, njobs = opt.njobs)
-    return c
+        self.pars_from_namespace(opt)
+        self.fill_up_fluid_code()
+        self.finalize_code()
+        self.write_src()
+        self.set_host_info(bfps.host_info)
+        if not os.path.exists(os.path.join(self.work_dir, self.simname + '.h5')):
+            self.write_par()
+            if self.parameters['nparticles'] > 0:
+                data = self.generate_tracer_state(
+                        species = 0,
+                        rseed = opt.particle_rand_seed)
+                for s in range(1, self.particle_species):
+                    self.generate_tracer_state(species = s, data = data)
+            init_condition_file = os.path.join(
+                    self.work_dir,
+                    self.simname + '_cvorticity_i{0:0>5x}'.format(0))
+            if not os.path.exists(init_condition_file):
+                if len(opt.src_simname) > 0:
+                    src_file = os.path.join(
+                            self.work_dir,
+                            opt.src_simname + '_cvorticity_i{0:0>5x}'.format(opt.src_iteration))
+                    os.symlink(src_file, init_condition_file)
+                else:
+                   self.generate_vector_field(
+                           write_to_file = True,
+                           spectra_slope = 2.0,
+                           amplitude = 0.25)
+        self.run(
+                ncpu = opt.ncpu,
+                njobs = opt.njobs)
+        return None
 
diff --git a/bfps/__init__.py b/bfps/__init__.py
index b9bc3a35ffafd7a181aa3521f2789aaa9f470244..bc841ed576fcd171fa3a2e6af08b1770678e4235 100644
--- a/bfps/__init__.py
+++ b/bfps/__init__.py
@@ -48,8 +48,7 @@ bfpsfolder = os.path.join(homefolder, '.config/', 'bfps')
 sys.path.append(bfpsfolder)
 from host_information import host_info
 
-from .code import code
-from .fluid_converter import fluid_converter
-from .fluid_resize import fluid_resize
+from .FluidConvert import FluidConvert
+from .FluidResize import FluidResize
 from .NavierStokes import NavierStokes
 
diff --git a/bfps/__main__.py b/bfps/__main__.py
index e9b883fb48e35685ecca95ce7f809ce7136f58b3..9e89964bd5830a65891e6f0905702bef3ddd858a 100644
--- a/bfps/__main__.py
+++ b/bfps/__main__.py
@@ -1,9 +1,55 @@
+#######################################################################
+#                                                                     #
+#  Copyright 2015 Max Planck Institute                                #
+#                 for Dynamics and Self-Organization                  #
+#                                                                     #
+#  This file is part of bfps.                                         #
+#                                                                     #
+#  bfps is free software: you can redistribute it and/or modify       #
+#  it under the terms of the GNU General Public License as published  #
+#  by the Free Software Foundation, either version 3 of the License,  #
+#  or (at your option) any later version.                             #
+#                                                                     #
+#  bfps is distributed in the hope that it will be useful,            #
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of     #
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the      #
+#  GNU General Public License for more details.                       #
+#                                                                     #
+#  You should have received a copy of the GNU General Public License  #
+#  along with bfps.  If not, see <http://www.gnu.org/licenses/>       #
+#                                                                     #
+# Contact: Cristian.Lalescu@ds.mpg.de                                 #
+#                                                                     #
+#######################################################################
+
+
+
 import sys
-from .Launcher import Launcher
+import argparse
+
+import bfps
+from .NavierStokes import NavierStokes
+from .FluidConvert import FluidConvert
+from .FluidResize import FluidResize
 
 def main():
-    l = Launcher()
-    l(sys.argv[1:] + ['--run'])
+    parser = argparse.ArgumentParser(prog = 'bfps')
+    parser.add_argument(
+            '-v', '--version',
+            action = 'version',
+            version = '%(prog)s ' + bfps.__version__)
+    parser.add_argument(
+            'base_class',
+            choices = ['NavierStokes',
+                       'FluidResize'],
+            type = str)
+    # first option is the choice of base class or -h or -v
+    # all other options are passed on to the base_class instance
+    opt = parser.parse_args(sys.argv[1:2])
+    # error is thrown if first option is not a base class, so Launcher
+    # cannot be executed by mistake.
+    c = eval('{0}()'.format(opt.base_class))
+    c.launch(args = sys.argv[2:])
     return None
 
 if __name__ == '__main__':
diff --git a/bfps/base.py b/bfps/_base.py
similarity index 73%
rename from bfps/base.py
rename to bfps/_base.py
index 1301fe65d465066134755c1e10d6b72d707c3236..78667424a6e455f20158e317f0e82868f5b21904 100644
--- a/bfps/base.py
+++ b/bfps/_base.py
@@ -28,12 +28,12 @@ import os
 import sys
 import numpy as np
 import h5py
-import bfps
+from bfps import install_info
+from bfps import __version__
 
-class base(object):
-    """
-        This class contains simulation parameters, and handles parameter related
-        functionalities of both python objects and C++ codes.
+class _base(object):
+    """This class contains simulation parameters, and handles parameter related
+    functionalities of both python objects and C++ codes.
     """
     def __init__(
             self,
@@ -109,8 +109,8 @@ class base(object):
             else:
                 ofile['parameters/' + k] = self.parameters[k]
         ofile['iteration'] = int(iter0)
-        for k in bfps.install_info.keys():
-            ofile['install_info/' + k] = str(bfps.install_info[k])
+        for k in install_info.keys():
+            ofile['install_info/' + k] = str(install_info[k])
         ofile.close()
         return None
     def read_parameters(self):
@@ -120,13 +120,66 @@ class base(object):
                     self.parameters[k] = type(self.parameters[k])(data_file['parameters/' + k].value)
         return None
     def pars_from_namespace(self, opt):
-        new_pars = vars(opt)
+        cmd_line_pars = vars(opt)
+        for k in ['nx', 'ny', 'nz']:
+            if type(cmd_line_pars[k]) == type(None):
+                cmd_line_pars[k] = opt.n
+        for k in self.parameters.keys():
+            if k in cmd_line_pars.keys():
+                if not type(cmd_line_pars[k]) == type(None):
+                    self.parameters[k] = cmd_line_pars[k]
         self.simname = opt.simname
         self.work_dir = opt.work_dir
-        for k in self.parameters.keys():
-            self.parameters[k] = new_pars[k]
         return None
     def get_coord(self, direction):
         assert(direction == 'x' or direction == 'y' or direction == 'z')
         return np.arange(.0, self.parameters['n' + direction])*2*np.pi / self.parameters['n' + direction]
+    def add_parser_arguments(
+            self,
+            parser):
+        self.specific_parser_arguments(parser)
+        self.parameters_to_parser_arguments(parser)
+        return None
+    def specific_parser_arguments(
+            self,
+            parser):
+        parser.add_argument(
+                '-v', '--version',
+                action = 'version',
+                version = '%(prog)s ' + __version__)
+        parser.add_argument(
+               '-n', '--cube-size',
+               type = int,
+               dest = 'n',
+               default = 32,
+               metavar = 'N',
+               help = 'code is run by default in a grid of NxNxN')
+        parser.add_argument(
+                '--ncpu',
+                type = int, dest = 'ncpu',
+                default = 2)
+        parser.add_argument(
+                '--simname',
+                type = str, dest = 'simname',
+                default = 'test')
+        parser.add_argument(
+                '--environment',
+                type = str,
+                dest = 'environment',
+                default = '')
+        parser.add_argument(
+                '--wd',
+                type = str, dest = 'work_dir',
+                default = './')
+        return None
+    def parameters_to_parser_arguments(
+            self,
+            parser):
+        for k in sorted(self.parameters.keys()):
+            parser.add_argument(
+                    '--{0}'.format(k),
+                    type = type(self.parameters[k]),
+                    dest = k,
+                    default = None)
+        return None
 
diff --git a/bfps/code.py b/bfps/_code.py
similarity index 97%
rename from bfps/code.py
rename to bfps/_code.py
index e8baced47d29812354974f3b8f44cda896b5ff67..c0f019fc9dbe1903a921a4c319c85818469bc5e4 100644
--- a/bfps/code.py
+++ b/bfps/_code.py
@@ -31,19 +31,19 @@ import subprocess
 import h5py
 from datetime import datetime
 import math
+
 import bfps
-from bfps.base import base
+from ._base import _base
 
-class code(base):
-    """
-        This class is meant to stitch together the C++ code into a final source file,
-        compile it, and handle all job launching.
+class _code(_base):
+    """This class is meant to stitch together the C++ code into a final source file,
+    compile it, and handle all job launching.
     """
     def __init__(
             self,
             work_dir = './',
             simname = 'test'):
-        super(code, self).__init__(work_dir = work_dir, simname = simname)
+        _base.__init__(self, work_dir = work_dir, simname = simname)
         self.version_message = ('/***********************************************************************\n' +
                                 '* this code automatically generated by bfps\n' +
                                 '* version {0}\n'.format(bfps.__version__) +
diff --git a/bfps/fluid_base.py b/bfps/_fluid_base.py
similarity index 76%
rename from bfps/fluid_base.py
rename to bfps/_fluid_base.py
index 06c3098be516f440ee0649155d26fd30ae97f08b..8c313b05f57f44d2f555a4ae53c440e958f5da29 100644
--- a/bfps/fluid_base.py
+++ b/bfps/_fluid_base.py
@@ -24,15 +24,18 @@
 
 
 
-import bfps
-import bfps.code
-import bfps.tools
+from ._code import _code
+from bfps import tools
 
 import os
 import numpy as np
 import h5py
 
-class fluid_particle_base(bfps.code):
+class _fluid_particle_base(_code):
+    """This class is meant to put together all common code between the
+    different C++ solvers/postprocessing tools, so that development of
+    specific functionalities is not overwhelming.
+    """
     def __init__(
             self,
             name = 'solver',
@@ -40,7 +43,8 @@ class fluid_particle_base(bfps.code):
             simname = 'test',
             dtype = np.float32,
             use_fftw_wisdom = True):
-        super(fluid_particle_base, self).__init__(
+        _code.__init__(
+                self,
                 work_dir = work_dir,
                 simname = simname)
         self.use_fftw_wisdom = use_fftw_wisdom
@@ -70,19 +74,6 @@ class fluid_particle_base(bfps.code):
         self.parameters['niter_out'] = 1024
         self.parameters['nparticles'] = 0
         self.parameters['dt'] = 0.01
-        self.parameters['nu'] = 0.1
-        self.parameters['fmode'] = 1
-        self.parameters['famplitude'] = 0.5
-        self.parameters['fk0'] = 1.5
-        self.parameters['fk1'] = 3.0
-        self.parameters['forcing_type'] = 'linear'
-        self.parameters['histogram_bins'] = 256
-        self.parameters['max_velocity_estimate'] = 1.0
-        self.parameters['max_vorticity_estimate'] = 1.0
-        self.parameters['QR2D_histogram_bins'] = 64
-        self.parameters['max_trS2_estimate'] = 1.0
-        self.parameters['max_Q_estimate'] = 1.0
-        self.parameters['max_R_estimate'] = 1.0
         self.fluid_includes = '#include "fluid_solver.hpp"\n'
         self.fluid_variables = ''
         self.fluid_definitions = ''
@@ -284,19 +275,19 @@ class fluid_particle_base(bfps.code):
             field_name = 'vorticity',
             write_to_file = False):
         np.random.seed(rseed)
-        Kdata00 = bfps.tools.generate_data_3D(
+        Kdata00 = tools.generate_data_3D(
                 self.parameters['nz']//2,
                 self.parameters['ny']//2,
                 self.parameters['nx']//2,
                 p = spectra_slope,
                 amplitude = amplitude).astype(self.ctype)
-        Kdata01 = bfps.tools.generate_data_3D(
+        Kdata01 = tools.generate_data_3D(
                 self.parameters['nz']//2,
                 self.parameters['ny']//2,
                 self.parameters['nx']//2,
                 p = spectra_slope,
                 amplitude = amplitude).astype(self.ctype)
-        Kdata02 = bfps.tools.generate_data_3D(
+        Kdata02 = tools.generate_data_3D(
                 self.parameters['nz']//2,
                 self.parameters['ny']//2,
                 self.parameters['nx']//2,
@@ -308,7 +299,7 @@ class fluid_particle_base(bfps.code):
         Kdata0[..., 0] = Kdata00
         Kdata0[..., 1] = Kdata01
         Kdata0[..., 2] = Kdata02
-        Kdata1 = bfps.tools.padd_with_zeros(
+        Kdata1 = tools.padd_with_zeros(
                 Kdata0,
                 self.parameters['nz'],
                 self.parameters['ny'],
@@ -396,90 +387,18 @@ class fluid_particle_base(bfps.code):
         assert (self.parameters['niter_todo'] % self.parameters['niter_part'] == 0)
         assert (self.parameters['niter_out']  % self.parameters['niter_stat'] == 0)
         assert (self.parameters['niter_out']  % self.parameters['niter_part'] == 0)
-        super(fluid_particle_base, self).write_par(iter0 = iter0)
+        _code.write_par(self, iter0 = iter0)
         with h5py.File(os.path.join(self.work_dir, self.simname + '.h5'), 'r+') as ofile:
             ofile['field_dtype'] = np.dtype(self.dtype).str
             kspace = self.get_kspace()
             for k in kspace.keys():
                 ofile['kspace/' + k] = kspace[k]
             nshells = kspace['nshell'].shape[0]
-            for k in ['velocity', 'vorticity']:
-                time_chunk = 2**20//(8*3*self.parameters['nx'])
-                time_chunk = max(time_chunk, 1)
-                ofile.create_dataset('statistics/xlines/' + k,
-                                     (1, self.parameters['nx'], 3),
-                                     chunks = (time_chunk, self.parameters['nx'], 3),
-                                     maxshape = (None, self.parameters['nx'], 3),
-                                     dtype = self.dtype,
-                                     compression = 'gzip')
-                time_chunk = 2**20//(8*3*3*nshells)
-                time_chunk = max(time_chunk, 1)
-                ofile.create_dataset('statistics/spectra/' + k + '_' + k,
-                                     (1, nshells, 3, 3),
-                                     chunks = (time_chunk, nshells, 3, 3),
-                                     maxshape = (None, nshells, 3, 3),
-                                     dtype = np.float64,
-                                     compression = 'gzip')
-                time_chunk = 2**20//(8*4*10)
-                time_chunk = max(time_chunk, 1)
-                a = ofile.create_dataset('statistics/moments/' + k,
-                                     (1, 10, 4),
-                                     chunks = (time_chunk, 10, 4),
-                                     maxshape = (None, 10, 4),
-                                     dtype = np.float64,
-                                     compression = 'gzip')
-                time_chunk = 2**20//(8*4*self.parameters['histogram_bins'])
-                time_chunk = max(time_chunk, 1)
-                ofile.create_dataset('statistics/histograms/' + k,
-                                     (1,
-                                      self.parameters['histogram_bins'],
-                                      4),
-                                     chunks = (time_chunk,
-                                               self.parameters['histogram_bins'],
-                                               4),
-                                     maxshape = (None,
-                                                 self.parameters['histogram_bins'],
-                                                 4),
-                                     dtype = np.int64,
-                                     compression = 'gzip')
-            for s in range(self.particle_species):
-                time_chunk = 2**20 // (8*3*
-                                       self.parameters['nparticles']*
-                                       self.parameters['tracers{0}_integration_steps'.format(s)])
-                time_chunk = max(time_chunk, 1)
-                ofile.create_dataset('particles/tracers{0}/rhs'.format(s),
-                                     (1,
-                                      self.parameters['tracers{0}_integration_steps'.format(s)],
-                                      self.parameters['nparticles'],
-                                      3),
-                                     maxshape = (None,
-                                                 self.parameters['tracers{0}_integration_steps'.format(s)],
-                                                 self.parameters['nparticles'],
-                                                 3),
-                                     chunks =  (time_chunk,
-                                                self.parameters['tracers{0}_integration_steps'.format(s)],
-                                                self.parameters['nparticles'],
-                                                3),
-                                     dtype = np.float64)
-                time_chunk = 2**20 // (8*3*self.parameters['nparticles'])
-                time_chunk = max(time_chunk, 1)
-                ofile.create_dataset(
-                    '/particles/tracers{0}/velocity'.format(s),
-                    (1,
-                     self.parameters['nparticles'],
-                     3),
-                    chunks = (time_chunk, self.parameters['nparticles'], 3),
-                    maxshape = (None, self.parameters['nparticles'], 3),
-                    dtype = np.float64)
-                if self.parameters['tracers{0}_acc_on'.format(s)]:
-                    ofile.create_dataset(
-                        '/particles/tracers{0}/acceleration'.format(s),
-                        (1,
-                         self.parameters['nparticles'],
-                         3),
-                        chunks = (time_chunk, self.parameters['nparticles'], 3),
-                        maxshape = (None, self.parameters['nparticles'], 3),
-                        dtype = np.float64)
             ofile.close()
         return None
+    def launch(
+            self,
+            args = [],
+            **kwargs):
+        return None
 
diff --git a/documentation/_static/api.rst b/documentation/_static/api.rst
index 68278ad80242f410c4e205337f1cfec0f11059de..a1288bad8a875cff458e68ce13f7bcb85debed9c 100644
--- a/documentation/_static/api.rst
+++ b/documentation/_static/api.rst
@@ -2,35 +2,43 @@
 API
 ===
 
+
 ----------
 bfps.tools
 ----------
 
-.. automodule:: tools
+.. automodule:: bfps.tools
     :members:
     :undoc-members:
-    :inherited-members:
     :show-inheritance:
 
 
--------------
-bfps.Launcher
--------------
+-----------------
+bfps.NavierStokes
+-----------------
 
-.. autoclass:: Launcher.Launcher
+.. automodule:: bfps.NavierStokes
     :members:
     :undoc-members:
-    :inherited-members:
     :show-inheritance:
 
 
 -----------------
-bfps.NavierStokes
+bfps.FluidConvert
 -----------------
 
-.. autoclass:: NavierStokes.NavierStokes
+.. automodule:: bfps.FluidConvert
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+
+----------------
+bfps.FluidResize
+----------------
+
+.. automodule:: bfps.FluidResize
     :members:
     :undoc-members:
-    :inherited-members:
     :show-inheritance:
 
diff --git a/documentation/conf.py b/documentation/conf.py
index 6e1cb48af7aeac96bcf022498c80c4e042c7caaa..5fc2f86699b6f446001efb8dda00fab733c434a4 100644
--- a/documentation/conf.py
+++ b/documentation/conf.py
@@ -21,7 +21,7 @@ import shlex
 # add these directories to sys.path here. If the directory is relative to the
 # documentation root, use os.path.abspath to make it absolute, like shown here.
 #sys.path.insert(0, os.path.abspath('.'))
-sys.path.append(os.path.abspath('../bfps'))
+sys.path.insert(0, os.path.abspath('..'))
 
 # -- General configuration ------------------------------------------------
 
diff --git a/done.txt b/done.txt
index 14afa7cd4dfa184999bb3084776733bf2968dc75..4bb2b052fc594db523bc5d35961f2a786cef596c 100644
--- a/done.txt
+++ b/done.txt
@@ -12,3 +12,6 @@ x 2016-01-08 simplify tracer/field addition mechanism
 x 2016-01-08 add stat choice parameter to add_particles                              @design +v1.0 +particle_api
 x 2016-01-15 particle output is broken when niter_part != 1                          @bugfix
 x 2016-01-19 clean up machine_settings mess                                          @design @documentation +v2.0
+x 2016-01-24 clear delimitation of public API                                        @documentation +v1.0
+x 2016-01-24 document coordinate conventions                                         @documentation +v1.0
+x 2016-01-24 move parameters from _fluid_particle_base to NavierStokes etc           @design
diff --git a/tests/base.py b/tests/base.py
index 948a1f8a0cadce1f85d8c5aa36c1d887e2a8f49e..2264de23b013fd2f497561fed0cab2b40cd54e56 100644
--- a/tests/base.py
+++ b/tests/base.py
@@ -33,10 +33,48 @@ import numpy as np
 import matplotlib.pyplot as plt
 
 import bfps
-from bfps import fluid_resize
+from bfps import FluidResize
 from bfps.tools import particle_finite_diff_test as acceleration_test
 
-parser = bfps.get_parser()
+import argparse
+
+def get_parser(base_class = bfps.NavierStokes,
+               n = 32,
+               ncpu = 2,
+               precision = 'single',
+               simname = 'test',
+               work_dir = './',
+               njobs = 1):
+    parser = argparse.ArgumentParser()
+    parser.add_argument('--run', dest = 'run', action = 'store_true')
+    parser.add_argument('-n',
+            type = int, dest = 'n',
+            default = n)
+    parser.add_argument('--ncpu',
+            type = int, dest = 'ncpu',
+            default = ncpu)
+    parser.add_argument('--precision',
+            type = str, dest = 'precision',
+            default = precision)
+    parser.add_argument('--simname',
+            type = str, dest = 'simname',
+            default = simname)
+    parser.add_argument('--wd',
+            type = str, dest = 'work_dir',
+            default = work_dir)
+    parser.add_argument('--njobs',
+            type = int, dest = 'njobs',
+            default = njobs)
+    c = base_class(simname = simname)
+    for k in sorted(c.parameters.keys()):
+        parser.add_argument(
+                '--{0}'.format(k),
+                type = type(c.parameters[k]),
+                dest = k,
+                default = c.parameters[k])
+    return parser
+
+parser = get_parser()
 parser.add_argument('--initialize', dest = 'initialize', action = 'store_true')
 parser.add_argument('--frozen', dest = 'frozen', action = 'store_true')
 parser.add_argument('--iteration',
@@ -54,25 +92,21 @@ parser.add_argument(
 def double(opt):
     old_simname = 'N{0:0>3x}'.format(opt.n)
     new_simname = 'N{0:0>3x}'.format(opt.n*2)
-    c = fluid_resize(
-            work_dir = opt.work_dir,
-            simname = old_simname + '_double',
-            dtype = opt.precision)
-    c.parameters['nx'] = opt.n
-    c.parameters['ny'] = opt.n
-    c.parameters['nz'] = opt.n
-    c.parameters['dst_nx'] = 2*opt.n
-    c.parameters['dst_ny'] = 2*opt.n
-    c.parameters['dst_nz'] = 2*opt.n
-    c.parameters['dst_simname'] = new_simname
-    c.parameters['src_simname'] = old_simname
-    c.parameters['niter_todo'] = 0
-    c.write_src()
-    c.set_host_info({'type' : 'pc'})
-    c.write_par()
-    c.run(ncpu = opt.ncpu,
-          err_file = 'err_',
-          out_file = 'out_')
+    c = FluidResize(dtype = opt.precision)
+    c.launch(
+            args = ['--simname', old_simname + '_double',
+                    '--wd', opt.work_dir,
+                    '--nx', '{0}'.format(opt.n),
+                    '--ny', '{0}'.format(opt.n),
+                    '--nz', '{0}'.format(opt.n),
+                    '--dst_nx', '{0}'.format(2*opt.n),
+                    '--dst_ny', '{0}'.format(2*opt.n),
+                    '--dst_nz', '{0}'.format(2*opt.n),
+                    '--dst_simname', new_simname,
+                    '--src_simname', old_simname,
+                    '--src_iteration', '0',
+                    '--src_wd', './',
+                    '--niter_todo', '0'])
     return None
 
 def launch(
diff --git a/tests/test_convergence.py b/tests/test_convergence.py
index cd9d52828f69617c303d9b87818c5fd652721a8b..25fa796edd6e55b3d7978fc318ca3b40bf4e2d06 100644
--- a/tests/test_convergence.py
+++ b/tests/test_convergence.py
@@ -72,7 +72,7 @@ def convergence_test(
             code_class = code_class,
             tracer_state_file = h5py.File(os.path.join(c0.work_dir, c0.simname + '.h5'), 'r'))
     # get real space fields
-    converter = bfps.fluid_converter(
+    converter = bfps.FluidConvert(
             fluid_precision = opt.precision,
             use_fftw_wisdom = False)
     converter.write_src()
diff --git a/tests/test_io.py b/tests/test_io.py
index 447cb405406a61573cc7d7b2b5942d012e6f55f6..327889ec96885f93eb43500952cb3e85d832c24d 100644
--- a/tests/test_io.py
+++ b/tests/test_io.py
@@ -26,7 +26,10 @@
 
 from base import *
 
-class test_io(bfps.code):
+import bfps
+from bfps._code import _code
+
+class test_io(_code):
     def __init__(
             self,
             name = 'test_io',
@@ -44,8 +47,6 @@ class test_io(bfps.code):
 if __name__ == '__main__':
     opt = parser.parse_args(
             ['-n', '32',
-             '--run',
-             '--initialize',
              '--ncpu', '2'] +
             sys.argv[1:])
     c = test_io(work_dir = opt.work_dir + '/io')
diff --git a/tests/test_main.py b/tests/test_main.py
new file mode 100644
index 0000000000000000000000000000000000000000..d748b157c01a44ae15bf169355f53e36fe5b1bc9
--- /dev/null
+++ b/tests/test_main.py
@@ -0,0 +1,32 @@
+#######################################################################
+#                                                                     #
+#  Copyright 2015 Max Planck Institute                                #
+#                 for Dynamics and Self-Organization                  #
+#                                                                     #
+#  This file is part of bfps.                                         #
+#                                                                     #
+#  bfps is free software: you can redistribute it and/or modify       #
+#  it under the terms of the GNU General Public License as published  #
+#  by the Free Software Foundation, either version 3 of the License,  #
+#  or (at your option) any later version.                             #
+#                                                                     #
+#  bfps is distributed in the hope that it will be useful,            #
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of     #
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the      #
+#  GNU General Public License for more details.                       #
+#                                                                     #
+#  You should have received a copy of the GNU General Public License  #
+#  along with bfps.  If not, see <http://www.gnu.org/licenses/>       #
+#                                                                     #
+# Contact: Cristian.Lalescu@ds.mpg.de                                 #
+#                                                                     #
+#######################################################################
+
+
+
+import bfps
+from bfps.__main__ import main
+
+if __name__ == '__main__':
+    main()
+
diff --git a/todo.txt b/todo.txt
index 75e6f29952cbd26b18df7ebe8d650079a5cf1a0c..04aaf5dbf519e6681b33a087a070115865250bfe 100644
--- a/todo.txt
+++ b/todo.txt
@@ -1,5 +1,4 @@
-(B) compute z polynomials only when needed                                  @optimization +v1.0
-(B) document coordinate conventions                                         @documentation +v1.0
+(B) compute z polynomials only when needed                                  @optimization
 (B) read https://www.xsede.org/documents/271087/369161/ExtScale-Koziol.pdf  @optimization @HDF5 +I/O
 (B) set up mechanism for adding in new PDEs                                 @design +v2.0 +alternate_algorithms
 (B) tweak HDF5 settings                                                     @optimization @HDF5 +I/O
@@ -7,8 +6,8 @@
 (C) code overview                                                           @documentation
 (C) move stat I/O to cpp lib                                                @design @HDF5
 (C) test involving hydrodynamic similarity                                  @tests
+(C) tests should use launch instead of get_parser                           @design @tests
 (C) use HDF5 io for fields                                                  @design @HDF5 +I/O
-(D) clear delimitation of public API                                        @documentation +v1.0
 (D) generalize interpolation comparison test                                @tests
 (D) generate separate lib(s) with extra classes                             @tests +alternate_algorithms
 (D) test anisotropic grids                                                  @tests