diff --git a/.gitignore b/.gitignore
index 581a345..e925824 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,17 +1,13 @@
-Airy_2d/Airy2d-out/*
-Laguerre_Gauss_3d/LaguerreGauss3d-out/*
-Gauss_2d/planar/*
-Gauss_2d/concave/*
-Gauss_2d/convex/*
+scripts/Airy_2d/Airy2d-out/*
+scripts/Laguerre_Gauss_3d/LaguerreGauss3d-out/*
+scripts/Gauss_2d/planar/*
+scripts/Gauss_2d/concave/*
+scripts/Gauss_2d/convex/*
 
 ### Jupyter ###
 .ipynb_checkpoints
 */.ipynb_checkpoints/*
 
-### Cython ###
-Laguerre_Gauss_3d/optbeam.c
-Laguerre_Gauss_3d/optbeam.html
-
 ### Python ###
 # Byte-compiled / optimized / DLL files
 __pycache__/
@@ -124,3 +120,6 @@ dmypy.json
 
 # Pyre type checker
 .pyre/
+
+# Virtual environment created via python -m venv optbeam-env
+optbeam-env
diff --git a/Airy_2d/optbeam.pxd b/Airy_2d/optbeam.pxd
deleted file mode 100644
index 95570f9..0000000
--- a/Airy_2d/optbeam.pxd
+++ /dev/null
@@ -1,58 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-file:    optbeam.pxd
-brief:   ...
-author:  Daniel Kotik
-version: 1.5-beta
-release date: xx.xx.2020
-creation date: 03.04.2020
-"""
-cimport cython
-from cpython.pycapsule cimport PyCapsule_New
-from cpython cimport bool
-
-# -----------------------------------------------------------------------------
-# declare C functions as "cpdef" to export them to the module
-# -----------------------------------------------------------------------------
-cdef extern from "math.h":
-    cpdef double _exp "exp" (double x) nogil
-    cpdef double _sqrt "sqrt" (double x) nogil
-
-cdef extern from "complex.h":
-    cpdef double complex _cexp "cexp" (double complex z) nogil
-
-# -----------------------------------------------------------------------------
-# function declarations
-# -----------------------------------------------------------------------------
-cdef double _imag_1d_func_c(int n, double *arr, void *func_ptr)
-cdef double _real_1d_func_c(int n, double *arr, void *func_ptr)
-
-@cython.locals(real=cython.double, imag=cython.double, real_tol=cython.double,
-               imag_tol=cython.double)
-cdef (double complex, double, double) _complex_quad(func, double a, double b, dict kwargs=*)
-
-# -----------------------------------------------------------------------------
-# class declarations
-# -----------------------------------------------------------------------------
-cdef class Beam2dCartesian:
-    cdef:
-        dict __dict__
-        double x, _k
-        public bool called
-        double _ry, _rz
-        double _a, _b
-
-    cdef double complex spectrum(self, double k_y) nogil
-    cdef double _phase(self, double k_y, double x, double y) nogil
-    cdef double complex _integrand(self, double k_y) nogil
-
-cdef class IncAiry2d(Beam2dCartesian):
-    cdef:
-        double _W_y
-        double _M
-        double _W
-        double xi
-
-    cdef double _heaviside(self, double x) nogil
-    cdef double complex _f_Airy(self, double k_y, double W_y, double M, double W) nogil
-    cdef double complex spectrum(self, double k_y) nogil
diff --git a/Airy_2d/optbeam.py b/Airy_2d/optbeam.py
deleted file mode 100644
index 72aeca0..0000000
--- a/Airy_2d/optbeam.py
+++ /dev/null
@@ -1,229 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-file:    optbeam.py
-brief:   ...
-author:  Daniel Kotik
-version: 1.5-beta
-release date: xx.xx.2020
-creation date: 12.04.2020
-"""
-import cython
-import math
-import sys
-
-from scipy.integrate import quad
-from types import MappingProxyType
-
-
-if not cython.compiled:
-    from math import sqrt as _sqrt
-    from cmath import exp as _cexp
-    print("\nPlease consider compiling `%s.py` via Cython: "
-          "`$ cythonize -3 -i %s.py`\n" % (__name__, __name__))
-else:
-    from scipy import LowLevelCallable
-
-
-def _real_1d_func(x, func):
-    """Return real part of a 1d function."""
-    return func(x).real
-
-
-def _imag_1d_func(x, func):
-    """Return imag part of a 1d function."""
-    return func(x).imag
-
-
-def _imag_1d_func_c(n, arr, func_ptr):
-    """Return imag part of a 1d function.
-
-    Cython implementation.
-    """
-    # pure python formulation of:
-    # return (<Beam2dCartesian>func_ptr)(arr[0]).imag
-    return cython.cast(Beam2dCartesian, func_ptr)._integrand(arr[0]).imag
-
-
-def _real_1d_func_c(n, arr, func_ptr):
-    """Return real part of a 1d function.
-
-    Cython implementation.
-    """
-    # pure python formulation of:
-    # return (<Beam2dCartesian>func_ptr)(arr[0]).real
-    return cython.cast(Beam2dCartesian, func_ptr)._integrand(arr[0]).real
-
-
-def _complex_quad(func, a, b, kwargs={}):
-    """Integrate real and imaginary part of the given function."""
-    if cython.compiled:
-        # pure python formulation of: cdef void *f_ptr = <void*>func
-        f_ptr = cython.declare(cython.p_void, cython.cast(cython.p_void, func))
-
-        func_capsule = PyCapsule_New(f_ptr, cython.NULL, cython.NULL)
-
-        current_module = sys.modules[__name__]
-
-        ll_real_1d_func_c = LowLevelCallable.from_cython(current_module,
-                                                         '_real_1d_func_c',
-                                                         func_capsule)
-        ll_imag_1d_func_c = LowLevelCallable.from_cython(current_module,
-                                                         '_imag_1d_func_c',
-                                                         func_capsule)
-        real, real_tol = quad(ll_real_1d_func_c, a, b, **kwargs)
-        imag, imag_tol = quad(ll_imag_1d_func_c, a, b, **kwargs)
-    else:
-        real, real_tol = quad(_real_1d_func, a, b, (func,), **kwargs)
-        imag, imag_tol = quad(_imag_1d_func, a, b, (func,), **kwargs)
-
-    return real + 1j*imag, real_tol, imag_tol
-
-
-def critical(n1, n2):
-    """Calculate critical angle in degrees."""
-    assert n1 > n2, "\nWarning: Critical angle is not defined, since n1 <= n2!"
-    return math.degrees(math.asin(n2/n1))
-
-
-def brewster(n1, n2):
-    """Calculate Brewster angle in degrees."""
-    return math.degrees(math.atan(n2/n1))
-
-
-class Beam2dCartesian:
-    """..."""
-
-    def __init__(self, x, params, called=False):
-        """..."""
-        self.x = x
-        self._k = params['k']
-        self._params = MappingProxyType(params)  # read-only view of a dict
-        self.called = called
-
-        # integration boundaries
-        self._a = -self._k
-        self._b = self._k
-
-    @property
-    def params(self):
-        """Beam specific parameters.
-
-        This is a read-only property.
-        """
-        return self._params
-
-    def profile(self, r):
-        """Field amplitude function.
-
-        Plane wave decomposition: calculate field amplitude at light source
-        position if not coinciding with beam waist.
-        """
-        self._ry = r.y
-        self._rz = r.z
-
-        if not self.called:
-            print("Calculating inital field configuration. "
-                  "This will take some time...")
-            self.called = True
-
-        try:
-            (result,
-             real_tol,
-             imag_tol) = _complex_quad(self if cython.compiled else self._integrand,
-                                       self._a, self._b, dict(limit=100))
-        except Exception as e:
-            print(type(e).__name__ + ":", e)
-            sys.exit()
-
-        return result
-
-    def spectrum(self, k_y):
-        """Spectrum amplitude distribution function, f."""
-        raise NotImplementedError
-
-    def _phase(self, k_y, x, y):
-        """Phase function."""
-        return x*_sqrt(self._k**2 - k_y**2) + k_y*y
-
-    def _integrand(self, k_y):
-        """Integrand function."""
-        return self.spectrum(k_y) * _cexp(1j*self._phase(k_y, self.x, self._ry))
-
-
-# -----------------------------------------------------------------------------
-# specific beam implementations based on Beam2dCartesian base class
-# -----------------------------------------------------------------------------
-class IncAiry2d(Beam2dCartesian):
-    """2d incomplete Airy beam."""
-
-    def __init__(self, x, params, called=False):
-        """..."""
-        self._W_y = params['W_y']
-        self._M = params['M']
-        self._W = params['W']
-        super().__init__(x, params, called)
-
-    def profile(self, r):
-        """..."""
-        if self.x == 0:
-            # adjust integration boundaries
-            self._a = self._M-self._W
-            self._b = self._M+self._W
-
-        return super().profile(r)
-
-    def spectrum(self, k_y):
-        """..."""
-        return self._f_Airy(k_y, self._W_y, self._M, self._W)
-
-    def _f_Airy(self, k_y, W_y, M, W):
-        """Airy spectrum amplitude."""
-        return W_y * _cexp(1.0j*(-1/3)*(k_y*W_y)**3) \
-            * self._heaviside(W_y * k_y - (M - W)) \
-            * self._heaviside((M + W) - W_y * k_y)
-
-    def _heaviside(self, x):
-        """Theta (Heaviside step) function."""
-        return 0 if x < 0 else 1
-
-    def _integrand(self, k_y):
-        """..."""
-        if self.x == 0:
-            xi = k_y
-            return _cexp(1.0j*(-(xi**3)/3 + (xi * self._ry/self._W_y)))
-        else:
-            # next line needs _integrand declared cpdef without nogil attribute,
-            # and will execute slower than repeating the super class integration
-            # routine function here
-            #return super(IncAiry2d, self)._integrand(k_y)
-            return self.spectrum(k_y) * _cexp(1j*self._phase(k_y, self.x, self._ry))
-
-
-def main():
-    class Vector3:
-        """Simple 3d vector class."""
-
-        def __init__(self, x, y, z):
-            self.x = x
-            self.y = y
-            self.z = z
-
-    x, y, z = -2.15, 0.3, 0.5
-
-    r = Vector3(0, y, z)
-    # alternative:
-    # import meep; r = meep.Vector3(0, y, z)
-
-    k1 = 75.39822368615503
-    w_0 = 0.09284038347027228
-    M = 0
-    W = 4
-    params = dict(W_y=w_0, k=k1, M=M, W=W)
-
-    beam = IncAiry2d(x=x, params=params)
-
-    return (beam, r)
-
-
-if __name__ == '__main__':
-    main()
diff --git a/Gauss_2d/optbeam.pxd b/Gauss_2d/optbeam.pxd
deleted file mode 100644
index 0a0bcef..0000000
--- a/Gauss_2d/optbeam.pxd
+++ /dev/null
@@ -1,55 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-file:    optbeam.pxd
-brief:   ...
-author:  Daniel Kotik
-version: 1.5-beta
-release date: xx.xx.2020
-creation date: 03.04.2020
-"""
-cimport cython
-from cpython.pycapsule cimport PyCapsule_New
-from cpython cimport bool
-
-# -----------------------------------------------------------------------------
-# declare C functions as "cpdef" to export them to the module
-# -----------------------------------------------------------------------------
-cdef extern from "math.h":
-    cpdef double _exp "exp" (double x) nogil
-    cpdef double _sqrt "sqrt" (double x) nogil
-
-cdef extern from "complex.h":
-    cpdef double complex _cexp "cexp" (double complex z) nogil
-
-# -----------------------------------------------------------------------------
-# function declarations
-# -----------------------------------------------------------------------------
-cdef double _imag_1d_func_c(int n, double *arr, void *func_ptr)
-cdef double _real_1d_func_c(int n, double *arr, void *func_ptr)
-
-@cython.locals(real=cython.double, imag=cython.double, real_tol=cython.double,
-               imag_tol=cython.double)
-cdef (double complex, double, double) _complex_quad(func, double a, double b, dict kwargs=*)
-
-# -----------------------------------------------------------------------------
-# class declarations
-# -----------------------------------------------------------------------------
-cdef class Beam2dCartesian:
-    cdef:
-        dict __dict__
-        double x, _k
-        public bool called
-
-        double _ry, _rz
-
-    cdef double spectrum(self, double k_y) nogil
-    cdef double _phase(self, double k_y, double x, double y) nogil
-    cdef double complex _integrand(self, double k_y) nogil
-
-cdef class Gauss2d(Beam2dCartesian):
-    cdef:
-        double _W_y
-        double _norm
-
-    cdef double _f_Gauss(self, double k_y, double W_y) nogil
-    cdef double spectrum(self, double k_y) nogil
diff --git a/Gauss_2d/optbeam.py b/Gauss_2d/optbeam.py
deleted file mode 100644
index 2bcf4f2..0000000
--- a/Gauss_2d/optbeam.py
+++ /dev/null
@@ -1,204 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-file:    optbeam.py
-brief:   ...
-author:  Daniel Kotik
-version: 1.5-beta
-release date: xx.xx.2020
-creation date: 03.04.2020
-"""
-import cython
-import math
-import sys
-
-from scipy.integrate import quad
-from types import MappingProxyType
-
-
-if not cython.compiled:
-    from math import (exp as _exp,
-                      sqrt as _sqrt)
-    from cmath import exp as _cexp
-    print("\nPlease consider compiling `%s.py` via Cython: "
-          "`$ cythonize -3 -i %s.py`\n" % (__name__, __name__))
-else:
-    from scipy import LowLevelCallable
-
-
-def _real_1d_func(x, func):
-    """Return real part of a 1d function."""
-    return func(x).real
-
-
-def _imag_1d_func(x, func):
-    """Return imag part of a 1d function."""
-    return func(x).imag
-
-
-def _imag_1d_func_c(n, arr, func_ptr):
-    """Return imag part of a 1d function.
-
-    Cython implementation.
-    """
-    # pure python formulation of:
-    # return (<Beam2dCartesian>func_ptr)(arr[0]).imag
-    return cython.cast(Beam2dCartesian, func_ptr)._integrand(arr[0]).imag
-
-
-def _real_1d_func_c(n, arr, func_ptr):
-    """Return real part of a 1d function.
-
-    Cython implementation.
-    """
-    # pure python formulation of:
-    # return (<Beam2dCartesian>func_ptr)(arr[0]).real
-    return cython.cast(Beam2dCartesian, func_ptr)._integrand(arr[0]).real
-
-
-def _complex_quad(func, a, b, kwargs={}):
-    """Integrate real and imaginary part of the given function."""
-    if cython.compiled:
-        # pure python formulation of: cdef void *f_ptr = <void*>func
-        f_ptr = cython.declare(cython.p_void, cython.cast(cython.p_void, func))
-
-        func_capsule = PyCapsule_New(f_ptr, cython.NULL, cython.NULL)
-
-        current_module = sys.modules[__name__]
-
-        ll_real_1d_func_c = LowLevelCallable.from_cython(current_module,
-                                                         '_real_1d_func_c',
-                                                         func_capsule)
-        ll_imag_1d_func_c = LowLevelCallable.from_cython(current_module,
-                                                         '_imag_1d_func_c',
-                                                         func_capsule)
-        real, real_tol = quad(ll_real_1d_func_c, a, b, **kwargs)
-        imag, imag_tol = quad(ll_imag_1d_func_c, a, b, **kwargs)
-    else:
-        real, real_tol = quad(_real_1d_func, a, b, (func,), **kwargs)
-        imag, imag_tol = quad(_imag_1d_func, a, b, (func,), **kwargs)
-
-    return real + 1j*imag, real_tol, imag_tol
-
-
-def critical(n1, n2):
-    """Calculate critical angle in degrees."""
-    assert n1 > n2, "\nWarning: Critical angle is not defined, since n1 <= n2!"
-    return math.degrees(math.asin(n2/n1))
-
-
-def brewster(n1, n2):
-    """Calculate Brewster angle in degrees."""
-    return math.degrees(math.atan(n2/n1))
-
-
-class Beam2dCartesian:
-    """..."""
-
-    def __init__(self, x, params, called=False):
-        """..."""
-        self.x = x
-        self._k = params['k']
-        self._params = MappingProxyType(params)  # read-only view of a dict
-        self.called = called
-
-    @property
-    def params(self):
-        """Beam specific parameters.
-
-        This is a read-only property.
-        """
-        return self._params
-
-    def profile(self, r):
-        """Field amplitude function.
-
-        Plane wave decomposition: calculate field amplitude at light source
-        position if not coinciding with beam waist.
-        """
-        self._ry = r.y
-        self._rz = r.z
-
-        if not self.called:
-            print("Calculating inital field configuration. "
-                  "This will take some time...")
-            self.called = True
-
-        try:
-            (result,
-             real_tol,
-             imag_tol) = _complex_quad(self if cython.compiled else self._integrand,
-                                       -self._k, self._k)
-        except Exception as e:
-            print(type(e).__name__ + ":", e)
-            sys.exit()
-
-        return result
-
-    def spectrum(self, k_y):
-        """Spectrum amplitude distribution function, f."""
-        raise NotImplementedError
-
-    def _phase(self, k_y, x, y):
-        """Phase function."""
-        return x*_sqrt(self._k**2 - k_y**2) + k_y*y
-
-    def _integrand(self, k_y):
-        """Integrand function."""
-        return self.spectrum(k_y) * _cexp(1j*self._phase(k_y, self.x, self._ry))
-
-
-# -----------------------------------------------------------------------------
-# specific beam implementations based on Beam2dCartesian base class
-# -----------------------------------------------------------------------------
-class Gauss2d(Beam2dCartesian):
-    """2d Gauss beam."""
-
-    def __init__(self, x, params, called=False):
-        """..."""
-        self._W_y = params['W_y']
-        self._norm = 2 * _sqrt(math.pi) / self._W_y
-        super().__init__(x, params, called)
-
-    def profile(self, r):
-        """..."""
-        # beam profile distribution (field amplitude) at the waist of the beam
-        if self.x == 0:
-            return self._norm * _exp(-(r.y / self._W_y)**2)
-        else:
-            return super().profile(r)
-
-    def spectrum(self, k_y):
-        """Spectrum amplitude function, f."""
-        return self._f_Gauss(k_y, self._W_y)
-
-    def _f_Gauss(self, k_y, W_y):
-        """Gaussian spectrum amplitude."""
-        return _exp(-(k_y*W_y/2)**2)
-
-
-def main():
-    class Vector3:
-        """Simple 3d vector class."""
-
-        def __init__(self, x, y, z):
-            self.x = x
-            self.y = y
-            self.z = z
-
-    x, y, z = -2.15, 0.3, 0.5
-
-    r = Vector3(0, y, z)
-    # alternative:
-    # import meep; r = meep.Vector3(0, y, z)
-
-    k1 = 116.11326447667875
-    w_0 = 0.1061032953945969
-    params = dict(W_y=w_0, k=k1)
-
-    beam = Gauss2d(x=x, params=params)
-
-    return (beam, r)
-
-
-if __name__ == '__main__':
-    main()
diff --git a/Laguerre_Gauss_3d/optbeam.pxd b/Laguerre_Gauss_3d/optbeam.pxd
deleted file mode 100644
index 2d3aa1c..0000000
--- a/Laguerre_Gauss_3d/optbeam.pxd
+++ /dev/null
@@ -1,94 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-file:    optbeam.pxd
-brief:   ...
-author:  Daniel Kotik
-version: 1.5-beta
-release date: xx.xx.2020
-creation date: 22.02.2020
-"""
-cimport cython
-from cpython.pycapsule cimport PyCapsule_New
-from cpython cimport bool
-
-# -----------------------------------------------------------------------------
-# declare C functions as "cpdef" to export them to the module
-# -----------------------------------------------------------------------------
-cdef extern from "stdlib.h":
-    cpdef int _abs "abs" (int n) nogil
-
-cdef extern from "math.h":
-    cpdef double _sin "sin" (double x) nogil
-    cpdef double _cos "cos" (double x) nogil
-    cpdef double _exp "exp" (double x) nogil
-    cpdef double _acos "acos" (double x) nogil
-    cpdef double _atan2 "atan2" (double y, double x) nogil
-
-cdef extern from "complex.h":
-    cpdef double complex _cexp "cexp" (double complex z) nogil
-    cpdef double complex _csqrt "csqrt" (double complex z) nogil
-
-# -----------------------------------------------------------------------------
-# function declarations
-# -----------------------------------------------------------------------------
-cdef double _imag_2d_func_c(int n, double *arr, void *func_ptr)
-cdef double _real_2d_func_c(int n, double *arr, void *func_ptr)
-
-@cython.locals(real=cython.double, imag=cython.double, real_tol=cython.double,
-               imag_tol=cython.double)
-cdef (double complex, double, double) _complex_dblquad(Beam3d func,
-                                                       double a, double b,
-                                                       double gfun, double hfun,
-                                                       dict kwargs=*)
-
-# -----------------------------------------------------------------------------
-# class declarations
-# -----------------------------------------------------------------------------
-cdef class Beam3d:
-    cdef:
-        dict __dict__
-        double x, _k
-        public bool called
-
-    cdef double complex _integrand(self, double x, double y) nogil
-
-cdef class Beam3dSpherical(Beam3d):
-    cdef:
-        double _ry, _rz
-
-    cdef double _phase(self, double sin_theta, double cos_theta, double phi,
-                       double x, double y, double z) nogil
-    cdef double complex spectrum(self, double sin_theta, double theta, double phi) nogil
-    cdef double complex _integrand(self, double theta, double phi) nogil
-
-cdef class Beam3dCartesian(Beam3d):
-    cdef:
-        double _ry, _rz
-
-    cdef double _phase(self, double k_y, double k_z, double x, double y, double z) nogil
-    cdef double complex spectrum(self, double k_y, double k_z) nogil
-    cdef double complex _integrand(self, double k_y, double k_z) nogil
-
-cdef class LaguerreGauss3d(Beam3dSpherical):
-    cdef:
-      int _m
-      double _W_y
-      double _norm
-
-    cdef double complex _f_Gauss_spherical(self, double sin_theta, double _W_y, double k) nogil
-    cdef double complex _f_Laguerre_Gauss_spherical(self, double sin_theta, double theta, double phi,
-                                                    double _W_y, double k, int _m) nogil
-    cdef double complex spectrum(self, double sin_theta, double theta, double phi) nogil
-
-cdef class LaguerreGauss3dCartesian(Beam3dCartesian):
-    cdef:
-      int _m
-      double _W_y
-      double _norm
-
-    cdef double _phi(self, double k_y, double k_z) nogil
-    cdef double _theta(self, double k_y, double k_z, double k) nogil
-    cdef double complex _f_Gauss_cartesian(self, double k_y, double k_z, double _W_y) nogil
-    cdef double complex _f_Laguerre_Gauss_cartesian(self, double k_y, double k_z,
-                                                    double _W_y, double k, int _m) nogil
-    cdef double complex spectrum(self, double k_y, double k_z) nogil
diff --git a/README.md b/README.md
index bce3cf2..5efb24b 100644
--- a/README.md
+++ b/README.md
@@ -1,12 +1,12 @@
-![concave](Gauss_2d/img/concave_intensity_cropped_rotated_resized.png)
-![planar](Gauss_2d/img/planar_intensity_cropped_rotated_resized.png)
-![convex](Gauss_2d/img/convex_intensity_cropped_rotated_resized.png)
-![Airy](Airy_2d/img/Airy_beam_M_0_W_4_scattering.png)
+![concave](scripts/Gauss_2d/img/concave_intensity_cropped_rotated_resized.png)
+![planar](scripts/Gauss_2d/img/planar_intensity_cropped_rotated_resized.png)
+![convex](scripts/Gauss_2d/img/convex_intensity_cropped_rotated_resized.png)
+![Airy](scripts/Airy_2d/img/Airy_beam_M_0_W_4_scattering.png)
 
-![snap](Laguerre_Gauss_3d/img/vortex_beam_m_2_longitudinal_resized.png)
-![snap](Laguerre_Gauss_3d/img/vortex_beam_m_2_transverse_resized.png)
-![snap](Laguerre_Gauss_3d/img/vortex_beam_m_2_3d_half_resized.png)
-![Airy](Airy_2d/img/Airy_beam_M_0_W_4_free_space.png)
+![snap](scripts/Laguerre_Gauss_3d/img/vortex_beam_m_2_longitudinal_resized.png)
+![snap](scripts/Laguerre_Gauss_3d/img/vortex_beam_m_2_transverse_resized.png)
+![snap](scripts/Laguerre_Gauss_3d/img/vortex_beam_m_2_3d_half_resized.png)
+![Airy](scripts/Airy_2d/img/Airy_beam_M_0_W_4_free_space.png)
 
 # Optical-beams-MEEP
 [![DOI](https://zenodo.org/badge/91711821.svg)](https://zenodo.org/badge/latestdoi/91711821)
@@ -31,8 +31,19 @@ Originally, these files have been used in studying optical beam shifts providing
 We highly recommend to install the parallel version of PyMeep via Conda:
 
 ```shell
-$ conda create -n pmp -c conda-forge pymeep=*=mpi_mpich_*
+# create conda virtual environment "pmp"
+$ conda create -n pmp -c conda-forge pymeep=*=mpi_mpich_* scipy
+
+# next command is optional, but recommended to enforce environment isolation (no local
+# user site packages in conda environment that may shadow conda installed dependencies)
+$ conda env config vars set PYTHONNOUSERSITE=True -n pmp
+
+# activate environment
 $ conda activate pmp
+
+# install optbeam package inside environment (-e flag is optional; it makes an
+# editable install for developers)
+$ python -m pip install [-e] .
 ```
 
 For detailed installation instructions, see the [Meep documentation](https://meep.readthedocs.io/en/latest/Installation/#conda-packages).
diff --git a/optbeam/_2d/__init__.py b/optbeam/_2d/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/optbeam/_2d/beams.py b/optbeam/_2d/beams.py
new file mode 100644
index 0000000..5e39280
--- /dev/null
+++ b/optbeam/_2d/beams.py
@@ -0,0 +1,136 @@
+
+import math
+import sys
+
+try:
+    import cython
+except ModuleNotFoundError:
+    cython = None
+
+if not cython or not cython.compiled:
+    from math import (exp as _exp,
+                      sqrt as _sqrt)
+    from cmath import exp as _cexp
+
+if cython and not cython.compiled:
+    print("\nPlease consider compiling `%s.py` via Cython: "
+          "`$ cythonize -3 -i %s.py`\n" % (__name__, __name__))
+
+from .helpers import _complex_quad, Beam2d
+
+
+class Beam2dCartesian(Beam2d):
+    """..."""
+
+    def profile(self, r):
+        """Field amplitude function.
+
+        Plane wave decomposition: calculate field amplitude at light source
+        position if not coinciding with beam waist.
+        """
+        self._ry = r.y
+        self._rz = r.z
+
+        if not self.called:
+            print("Calculating inital field configuration. "
+                  "This will take some time...")
+            self.called = True
+
+        try:
+            (result,
+             real_tol,
+             imag_tol) = _complex_quad(self if cython
+                                       and cython.compiled else self._integrand,
+                                       -self._k, self._k)
+        except Exception as e:
+            print(type(e).__name__ + ":", e)
+            sys.exit()
+
+        return result
+
+    def spectrum(self, k_y):
+        """Spectrum amplitude distribution function, f."""
+        raise NotImplementedError
+
+    def _phase(self, k_y, x, y):
+        """Phase function."""
+        return x*_sqrt(self._k**2 - k_y**2) + k_y*y
+
+    def _integrand(self, k_y):
+        """Integrand function."""
+        return self.spectrum(k_y) * _cexp(1j*self._phase(k_y, self.x, self._ry))
+
+
+# -----------------------------------------------------------------------------
+# specific beam implementations based on Beam2dCartesian base class
+# -----------------------------------------------------------------------------
+class Gauss2d(Beam2dCartesian):
+    """2d Gauss beam."""
+
+    def __init__(self, x, params, called=False):
+        """..."""
+        self._W_y = params['W_y']
+        self._norm = 2 * _sqrt(math.pi) / self._W_y
+        super().__init__(x, params, called)
+
+    def profile(self, r):
+        """..."""
+        # beam profile distribution (field amplitude) at the waist of the beam
+        if self.x == 0:
+            return self._norm * _exp(-(r.y / self._W_y)**2)
+        else:
+            return super().profile(r)
+
+    def spectrum(self, k_y):
+        """Spectrum amplitude function, f."""
+        return self._f_Gauss(k_y, self._W_y)
+
+    def _f_Gauss(self, k_y, W_y):
+        """Gaussian spectrum amplitude."""
+        return _exp(-(k_y*W_y/2)**2)
+
+
+class IncAiry2d(Beam2dCartesian):
+    """2d incomplete Airy beam."""
+
+    def __init__(self, x, params, called=False):
+        """..."""
+        self._W_y = params['W_y']
+        self._M = params['M']
+        self._W = params['W']
+        super().__init__(x, params, called)
+
+    def profile(self, r):
+        """..."""
+        if self.x == 0:
+            # adjust integration boundaries
+            self._a = self._M-self._W
+            self._b = self._M+self._W
+
+        return super().profile(r)
+
+    def spectrum(self, k_y):
+        """..."""
+        return self._f_Airy(k_y, self._W_y, self._M, self._W)
+
+    def _f_Airy(self, k_y, W_y, M, W):
+        """Airy spectrum amplitude."""
+        return W_y * _cexp(1.0j*(-1/3)*(k_y*W_y)**3) \
+            * self._heaviside(W_y * k_y - (M - W)) \
+            * self._heaviside((M + W) - W_y * k_y)
+
+    def _heaviside(self, x):
+        """Theta (Heaviside step) function."""
+        return 0 if x < 0 else 1
+
+    def _integrand(self, k_y):
+        """..."""
+        if self.x == 0:
+            xi = k_y
+            return _cexp(1.0j*(-(xi**3)/3 + (xi * self._ry/self._W_y)))
+        else:
+            # next line needs _integrand declared cpdef without nogil attribute,
+            # and will execute slower than repeating the super class integration
+            # routine function here
+            #return super(IncAiry2d, self)._integrand(k_y)
+            return self.spectrum(k_y) * _cexp(1j*self._phase(k_y, self.x, self._ry))
diff --git a/optbeam/_2d/helpers.py b/optbeam/_2d/helpers.py
new file mode 100644
index 0000000..0ce9909
--- /dev/null
+++ b/optbeam/_2d/helpers.py
@@ -0,0 +1,95 @@
+
+import sys
+
+try:
+    import cython
+except ModuleNotFoundError:
+    cython = None
+
+if cython:
+    if cython.compiled:
+        from scipy import LowLevelCallable
+    else:
+        print("\nPlease consider compiling `%s.py` via Cython: "
+              "`$ cythonize -3 -i %s.py`\n" % (__name__, __name__))
+
+from scipy.integrate import quad
+from types import MappingProxyType
+
+
+def _real_1d_func(x, func):
+    """Return real part of a 1d function."""
+    return func(x).real
+
+
+def _imag_1d_func(x, func):
+    """Return imag part of a 1d function."""
+    return func(x).imag
+
+
+def _imag_1d_func_c(n, arr, func_ptr):
+    """Return imag part of a 1d function.
+
+    Cython implementation.
+    """
+    # pure python formulation of:
+    # return (<Beam2dCartesian>func_ptr)(arr[0]).imag
+    return cython.cast(Beam2d, func_ptr)._integrand(arr[0]).imag
+
+
+def _real_1d_func_c(n, arr, func_ptr):
+    """Return real part of a 1d function.
+
+    Cython implementation.
+    """
+    # pure python formulation of:
+    # return (<Beam2dCartesian>func_ptr)(arr[0]).real
+    return cython.cast(Beam2d, func_ptr)._integrand(arr[0]).real
+
+
+def _complex_quad(func, a, b, kwargs={}):
+    """Integrate real and imaginary part of the given function."""
+    if cython and cython.compiled:
+        # pure python formulation of: cdef void *f_ptr = <void*>func
+        f_ptr = cython.declare(cython.p_void, cython.cast(cython.p_void, func))
+
+        func_capsule = PyCapsule_New(f_ptr, cython.NULL, cython.NULL)
+
+        current_module = sys.modules[__name__]
+
+        ll_real_1d_func_c = LowLevelCallable.from_cython(current_module,
+                                                         '_real_1d_func_c',
+                                                         func_capsule)
+        ll_imag_1d_func_c = LowLevelCallable.from_cython(current_module,
+                                                         '_imag_1d_func_c',
+                                                         func_capsule)
+        real, real_tol = quad(ll_real_1d_func_c, a, b, **kwargs)
+        imag, imag_tol = quad(ll_imag_1d_func_c, a, b, **kwargs)
+    else:
+        real, real_tol = quad(_real_1d_func, a, b, (func,), **kwargs)
+        imag, imag_tol = quad(_imag_1d_func, a, b, (func,), **kwargs)
+
+    return real + 1j*imag, real_tol, imag_tol
+
+
+class Beam2d:
+    """Abstract base class."""
+
+    def __init__(self, x, params, called=False):
+        """..."""
+        self.x = x
+        self._k = params['k']
+        self._params = MappingProxyType(params)  # read-only view of a dict
+        self.called = called
+
+    @property
+    def params(self):
+        """Beam specific parameters.
+
+        This is a read-only property.
+        """
+        return self._params
+
+    def _integrand(self, x):
+        """Integrand function over one coordinate x."""
+        raise NotImplementedError
diff --git a/optbeam/_3d/__init__.py b/optbeam/_3d/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/Laguerre_Gauss_3d/optbeam.py b/optbeam/_3d/beams.py
similarity index 61%
rename from Laguerre_Gauss_3d/optbeam.py
rename to optbeam/_3d/beams.py
index 9756e7c..7dcc1ef 100644
--- a/Laguerre_Gauss_3d/optbeam.py
+++ b/optbeam/_3d/beams.py
@@ -1,21 +1,13 @@
-# -*- coding: utf-8 -*-
-"""
-file:    optbeam.py
-brief:   ...
-author:  Daniel Kotik
-version: 1.5-beta
-release date: xx.xx.2020
-creation date: 22.02.2020
-"""
-import cython
+
 import math
 import sys
 
-from scipy.integrate import dblquad
-from types import MappingProxyType
-
+try:
+    import cython
+except ModuleNotFoundError:
+    cython = None
 
-if not cython.compiled:
+if not cython or not cython.compiled:
     from math import (sin as _sin,
                       cos as _cos,
                       exp as _exp,
@@ -24,99 +16,12 @@
     from cmath import (exp as _cexp,
                        sqrt as _csqrt)
     from builtins import abs as _abs
+
+if cython and not cython.compiled:
     print("\nPlease consider compiling `%s.py` via Cython: "
           "`$ cythonize -3 -i %s.py`\n" % (__name__, __name__))
-else:
-    from scipy import LowLevelCallable
-
-
-def _real_2d_func(x, y, func):
-    """Return real part of a 2d function."""
-    return func(x, y).real
-
-
-def _imag_2d_func(x, y, func):
-    """Return imag part of a 2d function."""
-    return func(x, y).imag
-
-
-def _imag_2d_func_c(n, arr, func_ptr):
-    """Return imag part of a 2d function.
-
-    Cython implementation.
-    """
-    # pure python formulation of:
-    # return (<Beam3d>func_ptr)(arr[0], arr[1]).imag
-    return cython.cast(Beam3d, func_ptr)._integrand(arr[0], arr[1]).imag
-
-
-def _real_2d_func_c(n, arr, func_ptr):
-    """Return real part of a 2d function.
-
-    Cython implementation.
-    """
-    # pure python formulation of:
-    # return (<Beam3d>func_ptr)(arr[0], arr[1]).real
-    return cython.cast(Beam3d, func_ptr)._integrand(arr[0], arr[1]).real
-
-
-def _complex_dblquad(func, a, b, gfun, hfun, kwargs={}):
-    """Integrate real and imaginary part of the given function."""
-    if cython.compiled:
-        # pure python formulation of: cdef void *f_ptr = <void*>func
-        f_ptr = cython.declare(cython.p_void, cython.cast(cython.p_void, func))
-
-        func_capsule = PyCapsule_New(f_ptr, cython.NULL, cython.NULL)
-
-        current_module = sys.modules[__name__]
-
-        ll_real_2d_func_c = LowLevelCallable.from_cython(current_module,
-                                                         '_real_2d_func_c',
-                                                         func_capsule)
-        ll_imag_2d_func_c = LowLevelCallable.from_cython(current_module,
-                                                         '_imag_2d_func_c',
-                                                         func_capsule)
-        real, real_tol = dblquad(ll_real_2d_func_c, a, b, gfun, hfun, **kwargs)
-        imag, imag_tol = dblquad(ll_imag_2d_func_c, a, b, gfun, hfun, **kwargs)
-    else:
-        real, real_tol = dblquad(_real_2d_func, a, b, gfun, hfun, (func,), **kwargs)
-        imag, imag_tol = dblquad(_imag_2d_func, a, b, gfun, hfun, (func,), **kwargs)
-
-    return real + 1j*imag, real_tol, imag_tol
-
-
-def critical(n1, n2):
-    """Calculate critical angle in degrees."""
-    assert n1 > n2, "\nWarning: Critical angle is not defined, since n1 <= n2!"
-    return math.degrees(math.asin(n2/n1))
-
 
-def brewster(n1, n2):
-    """Calculate Brewster angle in degrees."""
-    return math.degrees(math.atan(n2/n1))
-
-
-class Beam3d:
-    """Abstract base class."""
-
-    def __init__(self, x, params, called=False):
-        """..."""
-        self.x = x   # TODO: rename x to x_shift
-        self._k = params['k']
-        self._params = MappingProxyType(params)  # read-only view of a dict
-        self.called = called
-
-    @property
-    def params(self):
-        """Beam specific parameters.
-
-        This is a read-only property.
-        """
-        return self._params
-
-    def _integrand(self, x, y):
-        """Integrand function over two coordinates x and y."""
-        raise NotImplementedError
+from .helpers import _complex_dblquad, Beam3d
 
 
 class Beam3dSpherical(Beam3d):
@@ -135,7 +40,8 @@ def profile(self, r):
         try:
             (result,
              real_tol,
-             imag_tol) = _complex_dblquad(self if cython.compiled else self._integrand,
+             imag_tol) = _complex_dblquad(self if cython and cython.compiled
+                                          else self._integrand,
                                           0, 2*math.pi, 0, math.pi/2)
         except Exception as e:
             print(type(e).__name__ + ":", e)
@@ -160,7 +66,8 @@ def _integrand(self, theta, phi):
         cos_theta = _cos(theta)
 
         return sin_theta * cos_theta * self.spectrum(sin_theta, theta, phi) * \
-            _cexp(1j*self._phase(sin_theta, cos_theta, phi, self.x, self._ry, self._rz))
+            _cexp(1j*self._phase(sin_theta, cos_theta,
+                                 phi, self.x, self._ry, self._rz))
 
 
 class Beam3dCartesian(Beam3d):
@@ -179,7 +86,8 @@ def profile(self, r):
         try:
             (result,
              real_tol,
-             imag_tol) = _complex_dblquad(self if cython.compiled else self._integrand,
+             imag_tol) = _complex_dblquad(self if cython
+                                          and cython.compiled else self._integrand,
                                           -self._k, self._k, -self._k, self._k)
         except Exception as e:
             print(type(e).__name__ + ":", e)
@@ -309,34 +217,5 @@ def _f_Laguerre_Gauss_spherical(self, sin_theta, theta, phi, W_y, k, m):
         """
         return self._f_Gauss_spherical(sin_theta, W_y, k) * theta**_abs(m) * \
             _cexp(1j*m*phi)
-
-
-def main():
-    class Vector3:
-        """Simple vector class."""
-
-        def __init__(self, x, y, z):
-            self.x = x
-            self.y = y
-            self.z = z
-
-    x, y, z = -2.15, 0.3, 0.5
-
-    r = Vector3(0, y, z)
-    # alternative:
-    # import meep; r = meep.Vector3(0, y, z)
-
-    k1 = 31.41592653589793
-    w_0 = 0.25464790894703254
-    m_charge = 2
-
-    params = dict(W_y=w_0, m=m_charge, k=k1)
-
-    LGbeam = LaguerreGauss3d(x=x, params=params)
-    LGbeamCartesian = LaguerreGauss3dCartesian(x=x, params=params)
-
-    return (LGbeam, LGbeamCartesian, r)
-
-
-if __name__ == '__main__':
-    main()
+        return self._f_Gauss_spherical(sin_theta, W_y, k) * theta**_abs(m) * \
+            _cexp(1j*m*phi)
diff --git a/optbeam/_3d/helpers.py b/optbeam/_3d/helpers.py
new file mode 100644
index 0000000..d8a4f86
--- /dev/null
+++ b/optbeam/_3d/helpers.py
@@ -0,0 +1,97 @@
+
+import sys
+
+try:
+    import cython
+except ModuleNotFoundError:
+    cython = None
+
+if cython:
+    if cython.compiled:
+        from scipy import LowLevelCallable
+    else:
+        print("\nPlease consider compiling `%s.py` via Cython: "
+              "`$ cythonize -3 -i %s.py`\n" % (__name__, __name__))
+
+from scipy.integrate import dblquad
+from types import MappingProxyType
+
+
+def _real_2d_func(x, y, func):
+    """Return real part of a 2d function."""
+    return func(x, y).real
+
+
+def _imag_2d_func(x, y, func):
+    """Return imag part of a 2d function."""
+    return func(x, y).imag
+
+
+def _imag_2d_func_c(n, arr, func_ptr):
+    """Return imag part of a 2d function.
+
+    Cython implementation.
+    """
+    # pure python formulation of:
+    # return (<Beam3d>func_ptr)(arr[0], arr[1]).imag
+    return cython.cast(Beam3d, func_ptr)._integrand(arr[0], arr[1]).imag
+
+
+def _real_2d_func_c(n, arr, func_ptr):
+    """Return real part of a 2d function.
+
+    Cython implementation.
+    """
+    # pure python formulation of:
+    # return (<Beam3d>func_ptr)(arr[0], arr[1]).real
+    return cython.cast(Beam3d, func_ptr)._integrand(arr[0], arr[1]).real
+
+
+def _complex_dblquad(func, a, b, gfun, hfun, kwargs={}):
+    """Integrate real and imaginary part of the given function."""
+    if cython and cython.compiled:
+        # pure python formulation of: cdef void *f_ptr = <void*>func
+        f_ptr = cython.declare(cython.p_void, cython.cast(cython.p_void, func))
+
+        func_capsule = PyCapsule_New(f_ptr, cython.NULL, cython.NULL)
+
+        current_module = sys.modules[__name__]
+
+        ll_real_2d_func_c = LowLevelCallable.from_cython(current_module,
+                                                         '_real_2d_func_c',
+                                                         func_capsule)
+        ll_imag_2d_func_c = LowLevelCallable.from_cython(current_module,
+                                                         '_imag_2d_func_c',
+                                                         func_capsule)
+        real, real_tol = dblquad(ll_real_2d_func_c, a, b, gfun, hfun, **kwargs)
+        imag, imag_tol = dblquad(ll_imag_2d_func_c, a, b, gfun, hfun, **kwargs)
+    else:
+        real, real_tol = dblquad(
+            _real_2d_func, a, b, gfun, hfun, (func,), **kwargs)
+        imag, imag_tol = dblquad(
+            _imag_2d_func, a, b, gfun, hfun, (func,), **kwargs)
+
+    return real + 1j*imag, real_tol, imag_tol
+
+
+class Beam3d:
+    """Abstract base class."""
+
+    def __init__(self, x, params, called=False):
+        """..."""
+        self.x = x   # TODO: rename x to x_shift
+        self._k = params['k']
+        self._params = MappingProxyType(params)  # read-only view of a dict
+        self.called = called
+
+    @property
+    def params(self):
+        """Beam specific parameters.
+
+        This is a read-only property.
+        """
+        return self._params
+
+    def _integrand(self, x, y):
+        """Integrand function over two coordinates x and y."""
+        raise NotImplementedError
diff --git a/optbeam/__init__.py b/optbeam/__init__.py
new file mode 100644
index 0000000..cbb9e7a
--- /dev/null
+++ b/optbeam/__init__.py
@@ -0,0 +1,20 @@
+
+__version__ = "1.4.2"
+
+from math import (degrees as _degrees,
+                  asin as _asin,
+                  atan as _atan)
+
+from ._2d.beams import Gauss2d, IncAiry2d
+from ._3d.beams import LaguerreGauss3d
+
+
+def critical(n1, n2):
+    """Calculate critical angle in degrees."""
+    assert n1 > n2, "\nWarning: Critical angle is not defined, since n1 <= n2!"
+    return _degrees(_asin(n2/n1))
+
+
+def brewster(n1, n2):
+    """Calculate Brewster angle in degrees."""
+    return _degrees(_atan(n2/n1))
diff --git a/Airy_2d/Airy2d.ctl b/scripts/Airy_2d/Airy2d.ctl
similarity index 100%
rename from Airy_2d/Airy2d.ctl
rename to scripts/Airy_2d/Airy2d.ctl
diff --git a/Airy_2d/Airy2d.py b/scripts/Airy_2d/Airy2d.py
similarity index 100%
rename from Airy_2d/Airy2d.py
rename to scripts/Airy_2d/Airy2d.py
diff --git a/Airy_2d/img/Airy_beam_M_0_W_4_free_space.png b/scripts/Airy_2d/img/Airy_beam_M_0_W_4_free_space.png
similarity index 100%
rename from Airy_2d/img/Airy_beam_M_0_W_4_free_space.png
rename to scripts/Airy_2d/img/Airy_beam_M_0_W_4_free_space.png
diff --git a/Airy_2d/img/Airy_beam_M_0_W_4_scattering.png b/scripts/Airy_2d/img/Airy_beam_M_0_W_4_scattering.png
similarity index 100%
rename from Airy_2d/img/Airy_beam_M_0_W_4_scattering.png
rename to scripts/Airy_2d/img/Airy_beam_M_0_W_4_scattering.png
diff --git a/Gauss_2d/Gauss2d.ctl b/scripts/Gauss_2d/Gauss2d.ctl
similarity index 100%
rename from Gauss_2d/Gauss2d.ctl
rename to scripts/Gauss_2d/Gauss2d.ctl
diff --git a/Gauss_2d/Gauss2d.py b/scripts/Gauss_2d/Gauss2d.py
similarity index 100%
rename from Gauss_2d/Gauss2d.py
rename to scripts/Gauss_2d/Gauss2d.py
diff --git a/Gauss_2d/img/concave_cropped_rotated_resized.png b/scripts/Gauss_2d/img/concave_cropped_rotated_resized.png
similarity index 100%
rename from Gauss_2d/img/concave_cropped_rotated_resized.png
rename to scripts/Gauss_2d/img/concave_cropped_rotated_resized.png
diff --git a/Gauss_2d/img/concave_intensity_cropped_rotated_resized.png b/scripts/Gauss_2d/img/concave_intensity_cropped_rotated_resized.png
similarity index 100%
rename from Gauss_2d/img/concave_intensity_cropped_rotated_resized.png
rename to scripts/Gauss_2d/img/concave_intensity_cropped_rotated_resized.png
diff --git a/Gauss_2d/img/convex_cropped_rotated_resized.png b/scripts/Gauss_2d/img/convex_cropped_rotated_resized.png
similarity index 100%
rename from Gauss_2d/img/convex_cropped_rotated_resized.png
rename to scripts/Gauss_2d/img/convex_cropped_rotated_resized.png
diff --git a/Gauss_2d/img/convex_intensity_cropped_rotated_resized.png b/scripts/Gauss_2d/img/convex_intensity_cropped_rotated_resized.png
similarity index 100%
rename from Gauss_2d/img/convex_intensity_cropped_rotated_resized.png
rename to scripts/Gauss_2d/img/convex_intensity_cropped_rotated_resized.png
diff --git a/Gauss_2d/img/planar_cropped_rotated_resized.png b/scripts/Gauss_2d/img/planar_cropped_rotated_resized.png
similarity index 100%
rename from Gauss_2d/img/planar_cropped_rotated_resized.png
rename to scripts/Gauss_2d/img/planar_cropped_rotated_resized.png
diff --git a/Gauss_2d/img/planar_intensity_cropped_rotated_resized.png b/scripts/Gauss_2d/img/planar_intensity_cropped_rotated_resized.png
similarity index 100%
rename from Gauss_2d/img/planar_intensity_cropped_rotated_resized.png
rename to scripts/Gauss_2d/img/planar_intensity_cropped_rotated_resized.png
diff --git a/scripts/Gauss_2d/planar/Gauss2d-e2_s-000003696.h5 b/scripts/Gauss_2d/planar/Gauss2d-e2_s-000003696.h5
new file mode 100644
index 0000000..ee3b7f1
Binary files /dev/null and b/scripts/Gauss_2d/planar/Gauss2d-e2_s-000003696.h5 differ
diff --git a/scripts/Gauss_2d/planar/Gauss2d-e2_s-000003696.png b/scripts/Gauss_2d/planar/Gauss2d-e2_s-000003696.png
new file mode 100644
index 0000000..ca8deb0
Binary files /dev/null and b/scripts/Gauss_2d/planar/Gauss2d-e2_s-000003696.png differ
diff --git a/scripts/Gauss_2d/planar/Gauss2d-eps-000000000.h5 b/scripts/Gauss_2d/planar/Gauss2d-eps-000000000.h5
new file mode 100644
index 0000000..9ee04d3
Binary files /dev/null and b/scripts/Gauss_2d/planar/Gauss2d-eps-000000000.h5 differ
diff --git a/scripts/Gauss_2d/planar/Gauss2d-ez-000003696.h5 b/scripts/Gauss_2d/planar/Gauss2d-ez-000003696.h5
new file mode 100644
index 0000000..fda31b2
Binary files /dev/null and b/scripts/Gauss_2d/planar/Gauss2d-ez-000003696.h5 differ
diff --git a/scripts/Gauss_2d/planar/Gauss2d-ez-000003696.png b/scripts/Gauss_2d/planar/Gauss2d-ez-000003696.png
new file mode 100644
index 0000000..a7f2033
Binary files /dev/null and b/scripts/Gauss_2d/planar/Gauss2d-ez-000003696.png differ
diff --git a/Laguerre_Gauss_3d/LaguerreGauss3d.ctl b/scripts/Laguerre_Gauss_3d/LaguerreGauss3d.ctl
similarity index 100%
rename from Laguerre_Gauss_3d/LaguerreGauss3d.ctl
rename to scripts/Laguerre_Gauss_3d/LaguerreGauss3d.ctl
diff --git a/Laguerre_Gauss_3d/LaguerreGauss3d.py b/scripts/Laguerre_Gauss_3d/LaguerreGauss3d.py
similarity index 100%
rename from Laguerre_Gauss_3d/LaguerreGauss3d.py
rename to scripts/Laguerre_Gauss_3d/LaguerreGauss3d.py
diff --git a/Laguerre_Gauss_3d/img/vortex_beam_m_2_3d_full_resized.png b/scripts/Laguerre_Gauss_3d/img/vortex_beam_m_2_3d_full_resized.png
similarity index 100%
rename from Laguerre_Gauss_3d/img/vortex_beam_m_2_3d_full_resized.png
rename to scripts/Laguerre_Gauss_3d/img/vortex_beam_m_2_3d_full_resized.png
diff --git a/Laguerre_Gauss_3d/img/vortex_beam_m_2_3d_half_resized.png b/scripts/Laguerre_Gauss_3d/img/vortex_beam_m_2_3d_half_resized.png
similarity index 100%
rename from Laguerre_Gauss_3d/img/vortex_beam_m_2_3d_half_resized.png
rename to scripts/Laguerre_Gauss_3d/img/vortex_beam_m_2_3d_half_resized.png
diff --git a/Laguerre_Gauss_3d/img/vortex_beam_m_2_longitudinal_resized.png b/scripts/Laguerre_Gauss_3d/img/vortex_beam_m_2_longitudinal_resized.png
similarity index 100%
rename from Laguerre_Gauss_3d/img/vortex_beam_m_2_longitudinal_resized.png
rename to scripts/Laguerre_Gauss_3d/img/vortex_beam_m_2_longitudinal_resized.png
diff --git a/Laguerre_Gauss_3d/img/vortex_beam_m_2_transverse_resized.png b/scripts/Laguerre_Gauss_3d/img/vortex_beam_m_2_transverse_resized.png
similarity index 100%
rename from Laguerre_Gauss_3d/img/vortex_beam_m_2_transverse_resized.png
rename to scripts/Laguerre_Gauss_3d/img/vortex_beam_m_2_transverse_resized.png
diff --git a/Laguerre_Gauss_3d/integration_test.ipynb b/scripts/Laguerre_Gauss_3d/integration_test.ipynb
similarity index 100%
rename from Laguerre_Gauss_3d/integration_test.ipynb
rename to scripts/Laguerre_Gauss_3d/integration_test.ipynb
diff --git a/Laguerre_Gauss_3d/plot_2d_matplotlib.py b/scripts/Laguerre_Gauss_3d/plot_2d_matplotlib.py
similarity index 100%
rename from Laguerre_Gauss_3d/plot_2d_matplotlib.py
rename to scripts/Laguerre_Gauss_3d/plot_2d_matplotlib.py
diff --git a/Laguerre_Gauss_3d/plot_3d_mayavi.py b/scripts/Laguerre_Gauss_3d/plot_3d_mayavi.py
similarity index 100%
rename from Laguerre_Gauss_3d/plot_3d_mayavi.py
rename to scripts/Laguerre_Gauss_3d/plot_3d_mayavi.py
diff --git a/setup.py b/setup.py
new file mode 100644
index 0000000..a61060c
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,24 @@
+from setuptools import setup, find_packages
+
+
+with open("README.md") as f:
+    readme = f.read()
+
+with open("LICENSE") as f:
+    license = f.read()
+
+setup(
+    name="optbeam",
+    version="1.4.2",
+    description=("Simulation of reflection and refraction of polarized "
+                 "opticial beams at plane and curved dielectric interfaces"),
+    long_description=readme,
+    long_description_content_type='text/markdown',
+    url="https://github.com/DanielKotik/Optical-beams-MEEP",
+    author="Daniel Kotik et al.",
+    author_email="kotik@physics.org",
+    license=license,
+    packages=find_packages(exclude=("scripts")),
+    include_package_data=True,
+    install_requires=["scipy"],
+)