diff --git a/diffractsim/colour_functions.py b/diffractsim/colour_functions.py
index de4c6ef..b9c9622 100644
--- a/diffractsim/colour_functions.py
+++ b/diffractsim/colour_functions.py
@@ -3,7 +3,14 @@
 from scipy.interpolate import CubicSpline
 from .util.backend_functions import backend as bd
 
-illuminant_d65 = np.loadtxt(Path(__file__).parent / "./data/illuminant_d65.txt", usecols=(1))
+illuminant_d65 = np.loadtxt(Path(__file__).parent / "./data/illuminant_d65.txt", usecols = 1)
+high_pressure_sodium = np.loadtxt(Path(__file__).parent / "./data/high_pressure_sodium.txt", usecols = 1)
+incandescent_tugsten = np.loadtxt(Path(__file__).parent / "./data/incandescent_tugsten.txt", usecols = 1)
+compact_fluorescent_lamp = np.loadtxt(Path(__file__).parent / "./data/compact_fluorescent_lamp.txt", usecols = 1)
+mercury_vapor = np.loadtxt(Path(__file__).parent / "./data/mercury_vapor.txt", usecols = 1)
+LED_6770K = np.loadtxt(Path(__file__).parent / "./data/LED_6770K.txt", usecols = 1)
+ceramic_metal_halide = np.loadtxt(Path(__file__).parent / "./data/ceramic_metal_halide.txt", usecols = 1)
+cie_cmf = np.loadtxt(Path(__file__).parent / "./data/cie-cmf.txt", usecols = 1)
 
 """
 MPL 2.0 Clause License 
@@ -12,28 +19,29 @@
 All rights reserved.
 """
 
+
 class ColourSystem:
-    def __init__(self, spectrum_size = 400, spec_divisions = 40, clip_method = 1):
+    def __init__(self, spectrum_size=400, spec_divisions=40, clip_method=1):
         global bd
         from .util.backend_functions import backend as bd
 
         self.spectrum_size = spectrum_size
         # import CIE XYZ standard observer color matching functions
 
-        cmf = np.loadtxt(Path(__file__).parent / "./data/cie-cmf.txt", usecols=(1, 2, 3))
-        
-        self.Δλ = (779-380)/spectrum_size
-        self.λ_list = np.linspace(380,779, spectrum_size)
+        cmf = np.loadtxt(Path(__file__).parent / "./data/cie-cmf.txt", usecols = (1, 2, 3))
+
+        self.Δλ = (779 - 380) / spectrum_size
+        self.λ_list = np.linspace(380, 779, spectrum_size)
+
+        if spectrum_size == 400:
 
-        if spectrum_size == 400: 
-            
             # CIE XYZ standard observer color matching functions
-            self.cie_x =  cmf.T[0]
-            self.cie_y =  cmf.T[1]
-            self.cie_z =  cmf.T[2]
+            self.cie_x = cmf.T[0]
+            self.cie_y = cmf.T[1]
+            self.cie_z = cmf.T[2]
 
-        else: #by default spectrum has a size of 400. If new size, we interpolate
-            λ_list_old = np.linspace(380,779, 400)
+        else:  # by default spectrum has a size of 400. If new size, we interpolate
+            λ_list_old = np.linspace(380, 779, 400)
             self.cie_x = np.interp(self.λ_list, λ_list_old, cmf.T[0])
             self.cie_y = np.interp(self.λ_list, λ_list_old, cmf.T[1])
             self.cie_z = np.interp(self.λ_list, λ_list_old, cmf.T[2])
@@ -54,8 +62,8 @@ def __init__(self, spectrum_size = 400, spec_divisions = 40, clip_method = 1):
         self.T = bd.vstack(
             [[3.2406, -1.5372, -0.4986], [-0.9689, 1.8758, 0.0415], [0.0557, -0.2040, 1.0570]]
         )
-        
-        #clip methods for negative sRGB values:
+
+        # clip methods for negative sRGB values:
 
         self.CLIP_CLAMP_TO_ZERO = 0
         self.CLIP_ADD_WHITE = 1
@@ -73,7 +81,7 @@ def XYZ_to_sRGB_linear(self, XYZ):
                                                                          [Z1,Z2,Z3...]])
         """
 
-        rgb = bd.tensordot(self.T, XYZ, axes=([1, 0]))
+        rgb = bd.tensordot(self.T, XYZ, axes = ([1, 0]))
 
         if self.clip_method == self.CLIP_CLAMP_TO_ZERO:
             # set negative rgb values to zero
@@ -83,9 +91,9 @@ def XYZ_to_sRGB_linear(self, XYZ):
         if self.clip_method == self.CLIP_ADD_WHITE:
             # add enough white to make all rgb values nonnegative
             # find max negative rgb (or 0.0 if all non-negative), we need that much white
-            rgb_min = bd.amin(rgb, axis=0)
+            rgb_min = bd.amin(rgb, axis = 0)
             # get max positive component
-            rgb_max = bd.amax(rgb, axis=0)
+            rgb_max = bd.amax(rgb, axis = 0)
 
             # get scaling factor to maintain max rgb after adding white
             scaling = bd.where(rgb_max > 0.0, rgb_max / (rgb_max - rgb_min + 0.00001), bd.ones(rgb.shape))
@@ -94,8 +102,7 @@ def XYZ_to_sRGB_linear(self, XYZ):
             rgb = bd.where(rgb_min < 0.0, scaling * (rgb - rgb_min), rgb)
             return rgb
 
-
-    def sRGB_linear_to_sRGB(self,rgb_linear):
+    def sRGB_linear_to_sRGB(self, rgb_linear):
 
         """
         Convert a linear RGB color to a non linear RGB color (gamma correction).
@@ -114,14 +121,13 @@ def sRGB_linear_to_sRGB(self,rgb_linear):
         )
 
         # clip intensity if needed (rgb values > 1.0) by scaling
-        rgb_max = bd.amax(rgb, axis=0) + 0.00001  # avoid division by zero
+        rgb_max = bd.amax(rgb, axis = 0) + 0.00001  # avoid division by zero
         intensity_cutoff = 1.0
         rgb = bd.where(rgb_max > intensity_cutoff, rgb * intensity_cutoff / (rgb_max), rgb)
 
         return rgb
 
-
-    def sRGB_to_sRGB_linear(self,rgb):
+    def sRGB_to_sRGB_linear(self, rgb):
         """
         Convert a RGB color to a linear RGB color.
 
@@ -132,7 +138,6 @@ def sRGB_to_sRGB_linear(self,rgb):
         """
         return bd.where(rgb <= 0.03928, rgb / 12.92, bd.power((rgb + 0.055) / 1.055, 2.4))
 
-
     def XYZ_to_sRGB(self, XYZ):
         """
         Convert a XYZ to an RGB color.
@@ -148,7 +153,6 @@ def XYZ_to_sRGB(self, XYZ):
 
         return rgb
 
-
     def spec_to_XYZ(self, spec):
         """
         Convert a spectrum to an XYZ color.
@@ -158,8 +162,6 @@ def spec_to_XYZ(self, spec):
         Number of samples of each spectral intensity list doesn't matter, but they must be equally spaced.
         """
 
-        
-        
         if spec.ndim == 1:
 
             X = bd.dot(spec, self.cie_x) * self.Δλ * 0.003975 * 683.002
@@ -168,10 +170,9 @@ def spec_to_XYZ(self, spec):
             return bd.array([X, Y, Z])
 
         else:
-            return bd.tensordot(spec, self.cie_xyz, axes=([1, 1])).T * self.Δλ * 0.003975 * 683.002
-
+            return bd.tensordot(spec, self.cie_xyz, axes = ([1, 1])).T * self.Δλ * 0.003975 * 683.002
 
-    def spec_partition_to_XYZ(self, spec_partition, index = 0):
+    def spec_partition_to_XYZ(self, spec_partition, index=0):
         """
         Convert a spectrum to an XYZ color.
 
@@ -187,9 +188,7 @@ def spec_partition_to_XYZ(self, spec_partition, index = 0):
             return bd.array([X, Y, Z])
 
         else:
-            return bd.tensordot(spec_partition, self.cie_xyz_partitions[index], axes=([1, 1])).T * self.Δλ * 0.003975 * 683.002
-
-
+            return bd.tensordot(spec_partition, self.cie_xyz_partitions[index], axes = ([1, 1])).T * self.Δλ * 0.003975 * 683.002
 
     def spec_to_sRGB(self, spec):
         """
@@ -204,12 +203,10 @@ def spec_to_sRGB(self, spec):
         XYZ = self.spec_to_XYZ(spec)
         return self.XYZ_to_sRGB(XYZ)
 
-
-    def wavelength_to_XYZ(self,wavelength, intensity):
-
+    def wavelength_to_XYZ(self, wavelength, intensity):
 
         if (wavelength > 380) and (wavelength < 780):
-            index = int(wavelength-380)
+            index = int(wavelength - 380)
             X = intensity * self.cie_x[index] * self.Δλ * 0.003975 * 683.002
             Y = intensity * self.cie_y[index] * self.Δλ * 0.003975 * 683.002
             Z = intensity * self.cie_z[index] * self.Δλ * 0.003975 * 683.002
@@ -220,7 +217,6 @@ def wavelength_to_XYZ(self,wavelength, intensity):
 
         return bd.array([X, Y, Z])
 
-
     def wavelength_to_sRGB(self, wavelength, intensity):
 
         XYZ = self.wavelength_to_XYZ(wavelength, intensity)
@@ -229,4 +225,4 @@ def wavelength_to_sRGB(self, wavelength, intensity):
     def wavelength_to_sRGB_linear(self, wavelength, intensity):
 
         XYZ = self.wavelength_to_XYZ(wavelength, intensity)
-        return self.XYZ_to_sRGB_linear(XYZ)
\ No newline at end of file
+        return self.XYZ_to_sRGB_linear(XYZ)
diff --git a/diffractsim/monochromatic_simulator.py b/diffractsim/monochromatic_simulator.py
index 0589514..d237141 100644
--- a/diffractsim/monochromatic_simulator.py
+++ b/diffractsim/monochromatic_simulator.py
@@ -32,9 +32,7 @@ def __init__(self,  wavelength, extent_x, extent_y, Nx, Ny, intensity = 0.1 * W
         intensity: intensity of the field
         """
         global bd
-        global backend_name
         from .util.backend_functions import backend as bd
-        from .util.backend_functions import backend_name
         
         self.extent_x = extent_x
         self.extent_y = extent_y
@@ -109,8 +107,6 @@ def zoom_propagate(self, z, x_interval, y_interval):
         self.E = bluestein_method(self, self.E, z, self.λ, x_interval, y_interval)
 
 
-
-
     def propagate_to_image_plane(self, pupil, M, zi, z0, scale_factor = 1):
         """
         Assuming an optical system with linear response and assuming the system is only diffraction-limited by
@@ -237,21 +233,23 @@ def compute_colors_at(self, z):
 
     def interpolate(self, Nx, Ny):
         """Interpolate the field to the new shape (Nx,Ny)"""
-        from scipy.interpolate import RectBivariateSpline
+        from scipy.interpolate import interp2d
 
 
-        if backend_name == 'cupy':
+        if bd != np:
             self.E = self.E.get()
 
-        fun_real = RectBivariateSpline(
+        fun_real = interp2d(
                     self.dx*(np.arange(self.Nx)-self.Nx//2),
                     self.dy*(np.arange(self.Ny)-self.Ny//2),
-                    np.real(self.E))
+                    np.real(self.E),
+                    kind="cubic",)
 
-        fun_imag = RectBivariateSpline(
+        fun_imag = interp2d(
                     self.dx*(np.arange(self.Nx)-self.Nx//2),
                     self.dy*(np.arange(self.Ny)-self.Ny//2),
-                    np.imag(self.E))
+                    np.imag(self.E),
+                    kind="cubic",)
 
 
         self.Nx = Nx
@@ -306,12 +304,8 @@ def get_longitudinal_profile(self, start_distance, end_distance, steps, scale_fa
                 self.yy/=scale_factor
 
             rgb = self.get_colors()
-            if backend_name == 'jax':
-                longitudinal_profile_rgb = longitudinal_profile_rgb.at[i,:,:].set( rgb[self.Ny//2,:,:])
-                longitudinal_profile_E = longitudinal_profile_E.at[i,:].set(self.E[self.Ny//2,:])
-            else:
-                longitudinal_profile_rgb[i,:,:]  = rgb[self.Ny//2,:,:]
-                longitudinal_profile_E[i,:] = self.E[self.Ny//2,:]
+            longitudinal_profile_rgb[i,:,:]  = rgb[self.Ny//2,:,:]
+            longitudinal_profile_E[i,:] = self.E[self.Ny//2,:]
             self.E = np.copy(self.E0)
 
 
@@ -340,4 +334,4 @@ def __add__(self, Field):
             "The wavelength, dimensions and sampling of the interfering fields must be identical")
 
 
-    from .visualization import plot_colors, plot_phase, plot_intensity, plot_longitudinal_profile_colors, plot_longitudinal_profile_intensity
+    from .visualization import save_plot, plot_colors, plot_phase, plot_intensity, plot_longitudinal_profile_colors, plot_longitudinal_profile_intensity
diff --git a/diffractsim/polychromatic_simulator.py b/diffractsim/polychromatic_simulator.py
index 6d855eb..24158ec 100644
--- a/diffractsim/polychromatic_simulator.py
+++ b/diffractsim/polychromatic_simulator.py
@@ -23,9 +23,7 @@
 class PolychromaticField:
     def __init__(self, spectrum, extent_x, extent_y, Nx, Ny, spectrum_size = 180, spectrum_divisions = 30):
         global bd
-        global backend_name
         from .util.backend_functions import backend as bd
-        from .util.backend_functions import backend_name
 
         self.extent_x = extent_x
         self.extent_y = extent_y
@@ -97,7 +95,7 @@ def get_colors(self):
 
         bar = progressbar.ProgressBar()
 
-        # We compute the pattern of each wavelength separately, and associate it to small spectrum interval dλ = (780- 380)/spectrum_divisions . We approximately the final colour
+        # We compute the pattern of each wavelength separately, and associate it to small spectrum interval dλ = (780- 380)/spectrum_divisions . We approximate the final colour
         # by summing the contribution of each small spectrum interval converting its intensity distribution to a RGB space.
         
 
@@ -126,8 +124,7 @@ def get_colors(self):
                         sRGB_linear += self.cs.XYZ_to_sRGB_linear(XYZ)
 
 
-
-        if backend_name == 'cupy':
+        if bd != np:
             bd.cuda.Stream.null.synchronize()
         rgb = self.cs.sRGB_linear_to_sRGB(sRGB_linear)
         rgb = (rgb.T).reshape((self.Ny, self.Nx, 3))
@@ -135,6 +132,12 @@ def get_colors(self):
         return rgb
 
 
+    def get_field(self):
+        """get field of the cross-section profile at the current distance"""
+
+        return self.E
+
+
     def scale_propagate(self, z, scale_factor):
         """
         #raise NotImplementedError(self.__class__.__name__ + '.scale_propagate')
@@ -176,7 +179,6 @@ def get_colors_at_image_plane(self, pupil, M, zi, z0, scale_factor = 1):
         """
 
 
-
         for j in range(len(self.optical_elements)):
             self.E = self.E * self.optical_elements[j].get_transmittance(self.xx, self.yy, 0)
 
@@ -215,7 +217,7 @@ def get_colors_at_image_plane(self, pupil, M, zi, z0, scale_factor = 1):
             XYZ = self.cs.spec_partition_to_XYZ(bd.outer(Iλ, self.spec_partitions[i]),i)
             sRGB_linear += self.cs.XYZ_to_sRGB_linear(XYZ)
 
-        if backend_name == 'cupy':
+        if bd != np:
             bd.cuda.Stream.null.synchronize()
 
         self.xx = M_abs * self.xx
@@ -233,4 +235,4 @@ def get_colors_at_image_plane(self, pupil, M, zi, z0, scale_factor = 1):
 
 
 
-    from .visualization import plot_colors
+    from .visualization import save_plot, plot_colors
\ No newline at end of file
diff --git a/diffractsim/visualization/__init__.py b/diffractsim/visualization/__init__.py
index 7b97397..21a7711 100644
--- a/diffractsim/visualization/__init__.py
+++ b/diffractsim/visualization/__init__.py
@@ -1,4 +1,4 @@
-from .plot_colors import plot_colors
+from .plot_colors import save_plot, plot_colors
 from .plot_intensity import plot_intensity
 from .plot_phase import plot_phase
 from .plot_longitudinal_profile import plot_longitudinal_profile_colors, plot_longitudinal_profile_intensity
\ No newline at end of file
diff --git a/diffractsim/visualization/plot_colors.py b/diffractsim/visualization/plot_colors.py
index 131e15c..4ce181a 100644
--- a/diffractsim/visualization/plot_colors.py
+++ b/diffractsim/visualization/plot_colors.py
@@ -2,7 +2,6 @@
 import numpy as np
 from ..util.constants import *
 
-
 """
 
 MPL 2.0 License 
@@ -12,27 +11,69 @@
 
 """
 
-def plot_colors(self, rgb, figsize=(6, 6), xlim=None, ylim=None, text = None, units = mm, dark_background = False):
+def save_plot(self, rgb, figsize=(16, 16), xlim=None, ylim=None, text=None, path="pic.png", dark_background=True, tight=True):
+    from ..util.backend_functions import backend as bd
+    if bd != np:
+        rgb = rgb.get()
+
+    if dark_background:
+        plt.style.use("dark_background")
+    else:
+        plt.style.use("default")
+
+    fig = plt.figure(figsize = figsize)
+    ax = fig.add_subplot(1, 1, 1)
+
+    if xlim is not None:
+        ax.set_xlim(np.array(xlim) / mm)
+
+    if ylim is not None:
+        ax.set_ylim(np.array(ylim) / mm)
+
+    ax.set_xlabel("[mm]")
+    ax.set_ylabel("[mm]")
+
+    if text is None:
+        ax.set_title("Screen distance = " + str(self.z * 100) + " cm")
+    else:
+        ax.set_title(text)
+
+    im = ax.imshow(
+        rgb,
+        extent = [
+            float(self.x[0] - self.dx / 2) / mm,
+            float(self.x[-1] + self.dx / 2) / mm,
+            float(self.y[0] - self.dy / 2) / mm,
+            float(self.y[-1] + self.dy / 2) / mm,
+        ],
+        interpolation = "spline36", origin = "lower"
+    )
+    if tight:
+        fig.savefig(path, bbox_inches = 'tight')
+    else:
+        fig.savefig(path)  # save the figure to file
+    plt.close(fig)  # close the figure window
+
+
+def plot_colors(self, rgb, figsize=(6, 6), xlim=None, ylim=None, text=None, units=mm, dark_background=True):
     """visualize the diffraction pattern colors with matplotlib"""
 
     from ..util.backend_functions import backend as bd
-    from ..util.backend_functions import backend_name
-    
-    if dark_background == True:
+    if dark_background:
         plt.style.use("dark_background")
     else:
         plt.style.use("default")
-    if backend_name == 'cupy':
+    if bd != np:
         rgb = rgb.get()
 
-    fig = plt.figure(figsize=figsize)
+    fig = plt.figure(figsize = figsize)
     ax = fig.add_subplot(1, 1, 1)
 
-    if xlim != None:
-        ax.set_xlim(np.array(xlim)/units)
+    if xlim is not None:
+        ax.set_xlim(np.array(xlim) / units)
 
-    if ylim != None:
-        ax.set_ylim(np.array(ylim)/units)
+    if ylim is not None:
+        ax.set_ylim(np.array(ylim) / units)
 
     # we use mm by default
     if units == mm:
@@ -51,19 +92,19 @@ def plot_colors(self, rgb, figsize=(6, 6), xlim=None, ylim=None, text = None, un
         ax.set_xlabel("[m]")
         ax.set_ylabel("[m]")
 
-    if text == None:
+    if text is None:
         ax.set_title("Screen distance = " + str(self.z * 100) + " cm")
-    else: 
+    else:
         ax.set_title(text)
 
     im = ax.imshow(
-        (rgb),
-        extent=[
-            float(self.x[0] - self.dx/2) / units,
-            float(self.x[-1] + self.dx/2) / units,
-            float(self.y[0] - self.dy/2)/ units,
-            float(self.y[-1] + self.dy/2) / units,
+        rgb,
+        extent = [
+            float(self.x[0] - self.dx / 2) / units,
+            float(self.x[-1] + self.dx / 2) / units,
+            float(self.y[0] - self.dy / 2) / units,
+            float(self.y[-1] + self.dy / 2) / units,
         ],
-        interpolation="spline36", origin = "lower"
+        interpolation = "spline36", origin = "lower"
     )
     plt.show()
diff --git a/examples/animation/dummy b/examples/animation/dummy
new file mode 100644
index 0000000..8b13789
--- /dev/null
+++ b/examples/animation/dummy
@@ -0,0 +1 @@
+
diff --git a/examples/create_animation.py b/examples/create_animation.py
new file mode 100644
index 0000000..296a2cd
--- /dev/null
+++ b/examples/create_animation.py
@@ -0,0 +1,85 @@
+import subprocess
+from concurrent.futures import ThreadPoolExecutor
+import ffmpeg
+
+
+def run_simulation_for_frame(slave, frame_number, is_monochromatic, wave_length, spectrum, phase_mask_name, spectrum_sz, spectrum_div,
+                                img_size_x, img_size_y, out_filename_pattern):
+    # Build the command to run diffractsim_simulation.py with all configurable arguments
+    cmd = [
+        "python", slave,
+        "--frame", str(frame_number),
+        "--wavelength", str(wave_length),
+        "--light_spectrum", spectrum,
+        "--phase_mask", phase_mask_name,
+        "--spectrum_size", str(spectrum_sz),
+        "--spectrum_divisions", str(spectrum_div),
+        "--image_size", str(img_size_x), str(img_size_y),
+        "--output_filename_pattern", out_filename_pattern
+    ]
+
+    # Add the --monochromatic flag if necessary
+    if is_monochromatic:
+        cmd.append("--monochromatic")
+
+    # Run the command using subprocess
+    print(f"Computing frame {frame_number}")
+    subprocess.run(cmd)
+
+
+if __name__ == "__main__":
+    # Number of concurrent instances to run
+    max_workers = 11
+
+    # Animation start and end frame as position on z axis (ie distance to lens/sensor)
+    start_frame = 1
+    end_frame = 100
+    encode_to_video = True
+    framerate = 20
+
+    # Simulation slave to use
+    sim_slave = "phase_hologram_reconstruction_slave.py"
+
+    # Configuration for the simulation (outer scope variables)
+    monochromatic = False
+    wavelength = 591.0
+    light_spectrum = "illuminant_d65"
+    phase_mask = "github_logo.png"
+    spectrum_size = 400  # by default spectrum has a size of 400. If new size, we interpolate
+    spectrum_divisions = 400
+    image_size_x = 23.0  # in mm
+    image_size_y = 23.0  # in mm
+    # Options for light_spectrum are:
+    # illuminant_d65, high_pressure_sodium, incandescent_tugsten,
+    # compact_fluorescent_lamp, LED_6770K, ceramic_metal_halide,
+    # mercury_vapor, cie_cmf
+
+    path = "./animation/"
+    filename = "frame_"
+    filetype = ".png"  # Can save as jpg or png
+    output_filename_pattern = path + filename + "{frame:06d}" + filetype
+
+    # ThreadPoolExecutor to handle parallel execution
+    with ThreadPoolExecutor(max_workers = max_workers) as executor:
+        # Submitting tasks for frames
+        futures = []
+        for frame in range(start_frame, end_frame + 1):
+            futures.append(executor.submit(
+                run_simulation_for_frame, sim_slave, str(frame), monochromatic, str(wavelength), light_spectrum, phase_mask,
+                str(spectrum_size), str(spectrum_divisions), str(image_size_x), str(image_size_y), output_filename_pattern))
+
+        # Waiting for all tasks to complete
+        for future in futures:
+            try:
+                future.result()
+            except Exception as e:
+                print(f"An error occurred: {e}")
+
+    if encode_to_video:
+        (
+            ffmpeg
+            .input(path + filename + '%06d' + filetype, start_number = start_frame, framerate = framerate)  # Note: pattern_type='glob' and "*" wildcard only work in linux
+            .output(path + 'animation.mp4', vcodec = 'png', movflags = 'faststart')  # crf=1, **{'qscale:v': 1}
+            .run(overwrite_output = True)
+            # for Android use: vcodec='libvpx-vp9', pix_fmt='yuv420p', video_bitrate=200000000, movflags='faststart'
+        )
diff --git a/examples/phase_hologram_reconstruction_slave.py b/examples/phase_hologram_reconstruction_slave.py
new file mode 100644
index 0000000..2a09d9e
--- /dev/null
+++ b/examples/phase_hologram_reconstruction_slave.py
@@ -0,0 +1,96 @@
+import argparse
+import diffractsim
+from diffractsim import PolychromaticField, MonochromaticField, ApertureFromImage, Lens, mm, nm, cm, colour_functions as cf
+
+# Set the backend
+diffractsim.set_backend("CPU")
+
+
+def run_simulation(frame: int, monochromatic: bool, wavelength: float, spectrum: str, phase_mask: str, spectrum_size: int,
+                    spectrum_divisions: int, image_size_x: float, image_size_y: float,
+                    output_filename_pattern: str):
+
+    match spectrum:
+        case "illuminant_d65":
+            light_spectrum = cf.illuminant_d65
+        case "high_pressure_sodium":
+            light_spectrum = cf.high_pressure_sodium
+        case "incandescent_tugsten":
+            light_spectrum = cf.incandescent_tugsten
+        case "compact_fluorescent_lamp":
+            light_spectrum = cf.compact_fluorescent_lamp
+        case "mercury_vapor":
+            light_spectrum = cf.mercury_vapor
+        case "LED_6770K":
+            light_spectrum = cf.LED_6770K
+        case "ceramic_metal_halide":
+            light_spectrum = cf.ceramic_metal_halide
+        case "cie_cmf":
+            light_spectrum = cf.cie_cmf
+
+    # Handle monochromatic or polychromatic fields
+    if monochromatic:
+        F = MonochromaticField(
+            wavelength = wavelength * nm, extent_x = 30 * mm, extent_y = 30 * mm, Nx = 2400, Ny = 2400, intensity = 0.0075)
+    else:
+        F = PolychromaticField(
+            spectrum = light_spectrum, extent_x = 30 * mm, extent_y = 30 * mm, Nx = 2400, Ny = 2400,
+            spectrum_size = spectrum_size, spectrum_divisions = spectrum_divisions)
+
+    # Load the hologram as a phase mask aperture
+    F.add(ApertureFromImage(
+        amplitude_mask_path = "./apertures/white_background.png",
+        phase_mask_path = phase_mask,
+        image_size = (image_size_x * mm, image_size_y * mm), simulation = F))
+
+    # Add a lens with a focal length
+    F.add(Lens(f = 80 * cm))
+
+    # Propagate the field based on the frame argument
+    F.propagate(frame * cm)
+
+    # Get the reconstructed image (colors) at the specified z position
+    rgb = F.get_colors()
+
+    # Format the filename using the provided pattern
+    imagefile = output_filename_pattern.format(frame = frame)
+
+    print(imagefile + " done")
+
+    # Save the plot to the specified path
+    F.save_plot(rgb, figsize = (16, 16), path = imagefile, tight = True)
+
+
+if __name__ == "__main__":
+    # Set up argument parser to accept all configurable arguments
+    parser = argparse.ArgumentParser(description = "Run a phase hologram reconstruction with customizable options.")
+
+    # Adding arguments for frame, monochromatic field, and the new parameters
+    parser.add_argument('--frame', required=True, type = int, default = 1, help = "Set the frame number (default is 1).")
+    parser.add_argument('--monochromatic', action=argparse.BooleanOptionalAction, default = False, help = "Use monochromatic field instead of polychromatic.")
+    parser.add_argument('--wavelength', required=True, type = float, default = 591.0, help = "Wavelength in nm for the monochromatic field.")
+    parser.add_argument('--light_spectrum', required=True, type = str, default = "illuminant_d65", help = "Light spectrum to use. Options are:"
+                                                                                                "illuminant_d65, high_pressure_sodium, incandescent_tugsten,"
+                                                                                                "compact_fluorescent_lamp, LED_6770K, ceramic_metal_halide,"
+                                                                                                "mercury_vapor, cie_cmfSee ./diffractsim/data/*.")
+    parser.add_argument('--phase_mask', required=True, type = str, default = "github_logo.png", help = "Path to the phase mask image.")
+    parser.add_argument('--spectrum_size', required=True, type = int, default = 400, help = "Spectrum size for the polychromatic field.")
+    parser.add_argument('--spectrum_divisions', required=True, type = int, default = 400, help = "Number of spectrum divisions for the polychromatic field.")
+    parser.add_argument('--image_size', required=True, nargs = 2, type = float, default = [23.0, 23.0], help = "Size of the image in mm (x, y).")
+    parser.add_argument('--output_filename_pattern', required=True, type = str, default = "./animation/frame_{frame:06d}.png",
+                        help = "Pattern for naming the output files. Use '{frame}' to insert the frame number.")
+
+    # Parse the arguments from the command line
+    args = parser.parse_args()
+
+    # Run the simulation with the provided arguments
+    run_simulation(frame = args.frame,
+                    monochromatic = args.monochromatic,
+                    wavelength = args.wavelength,
+                    spectrum = args.light_spectrum,
+                    phase_mask = args.phase_mask,
+                    spectrum_size = args.spectrum_size,
+                    spectrum_divisions = args.spectrum_divisions,
+                    image_size_x = args.image_size[0],
+                    image_size_y = args.image_size[1],
+                    output_filename_pattern = args.output_filename_pattern)
diff --git a/examples/run_ffmpeg_manually.py b/examples/run_ffmpeg_manually.py
new file mode 100644
index 0000000..9426e60
--- /dev/null
+++ b/examples/run_ffmpeg_manually.py
@@ -0,0 +1,16 @@
+import ffmpeg
+
+framerate = 25
+start_frame = 1
+path = "./animation/"
+filename = "frame_"
+filetype = ".png"  # Can save as jpg or png
+
+# pip install ffmpeg-python // Add ffmpeg.exe to system path variable (Or just put in same folder)
+(
+    ffmpeg
+    .input(path + filename + '%06d' + filetype, start_number = start_frame, framerate = framerate)  # Note: pattern_type='glob' and "*" wildcard only work in linux
+    .output(path + 'animation.mp4', vcodec='png', movflags='faststart')  # crf=1, **{'qscale:v': 1}
+    .run(overwrite_output = True)
+    # for Android use: vcodec='libvpx-vp9', pix_fmt='yuv420p', video_bitrate=200000000, movflags='faststart'
+)