diff --git a/docs/example_test_script.py b/docs/example_test_script.py
new file mode 100644
index 0000000..2badd5f
--- /dev/null
+++ b/docs/example_test_script.py
@@ -0,0 +1,162 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+## example_test_script.py
+## Created by Aurélien STCHERBININE
+## Last modified by Aurélien STCHERBININE : 29/05/2024
+
+##----------------------------------------------------------------------------------------
+"""Example test for OMEGA-Py: import, apply corrections, compute band depth and display
+figures for OMEGA observation 0979_2.
+"""
+##----------------------------------------------------------------------------------------
+##----------------------------------------------------------------------------------------
+## Package importation
+import omegapy.omega_data as od
+import omegapy.omega_plots as op
+import omegapy.useful_functions as uf
+
+import numpy as np
+import matplotlib.pyplot as plt
+plt.ion()   # Activation of Matplotlib interactive mode
+
+##----------------------------------------------------------------------------------------
+## Parameters to adjust
+# Number of simultaneous processes uses to perform the thermal correction
+npool = 10
+
+# Colormap for reflectance maps
+cmap_refl = 'Greys_r'
+
+# Colormap for BD maps
+cmap_bd = 'Blues'
+
+########
+## Alternative with cmocean (if installed)
+# import cmocean.cm as cmo
+# cmap_bd = cmo.ice
+########
+
+## If needed, set the paths with the functions below
+## (see https://astcherbinine.github.io/omegapy/configuration/#windows-or-if-you-have-troubles-using-the-environment-variables)
+# od.set_omega_bin_path("/path/to/binary/files/folder/")
+# od.set_omega_py_path("/path/to/omegapy-made/files/folder/")
+##
+
+##----------------------------------------------------------------------------------------
+## Data importation & correction
+# Load the data cube
+omega0 = od.OMEGAdata('0979_2')
+
+# Apply thermal and atmospheric corrections
+omega = od.corr_therm_atm(
+    omega0, 
+    npool = npool,    # Adjust npool according to your system
+    )
+
+##----------------------------------------------------------------------------------------
+## OMEGA mask to hide bad pixels / calibration lines
+mask = od.omega_mask(
+    omega,
+    hide_128 = True,
+    emer_lim = 10,
+    inci_lim = 70,
+    tempc_lim = -194,
+    limsat_c = 500
+    )
+
+##----------------------------------------------------------------------------------------
+## Compute band depth at 1.5μm
+bd_15 = od.BD_omega(
+    omega, 
+    [1.5, 1.51], 
+    1.30, 
+    1.71, 
+    norm = True
+    )
+
+# Mask on band depth values for overplotting
+mask_bd15 = np.ones(bd_15.shape)    # Initialisation with array of 1
+mask_bd15[bd_15 < 0.1] = np.nan     # NaN for the pixels with no ice that we want to hide
+
+##----------------------------------------------------------------------------------------
+## Display figures
+# Surface reflectance
+op.show_omega_v2(
+    omega, 
+    lam = 1.085, 
+    polar = True, 
+    vmin = 0, 
+    vmax = 0.6,
+    Nfig = "reflectance",
+    cmap = cmap_refl,
+    )
+
+# 1.5μm Band Depth
+op.show_data_v2(
+    omega, 
+    data = bd_15, 
+    cb_title = "1.5μm BD", 
+    polar = True, 
+    cmap = cmap_bd, 
+    vmin = 0, 
+    vmax = 0.75,
+    Nfig = "BD15",
+    )
+
+# Overplotting 1.5μm BD over surface reflectance
+op.show_omega_v2(
+    omega,
+    lam = 1.085,
+    polar = True,
+    vmin = 0,
+    vmax = 0.5,
+    cbar = False,
+    Nfig = "overplot",
+    cmap = cmap_refl,
+    )
+
+op.show_data_v2(
+    omega,
+    data = bd_15,
+    mask = mask_bd15,
+    polar = True,
+    cmap = cmap_bd,
+    cb_title = "1.5 μm BD",
+    vmin = 0,
+    vmax = 0.75,
+    Nfig = "overplot",
+    title = 'Overplot map for ORB0979_2',
+    )
+
+plt.figure('reflectance').savefig('omegapy_ORB0979_2_reflectance.png')
+plt.figure('BD15').savefig('omegapy_ORB0979_2_BD15.png')
+plt.figure('overplot').savefig('omegapy_ORB0979_2_overplot_reflectance_BD15.png')
+
+# Interactive display of the observation (@ λ = 1.085 µm)
+op.show_omega_interactif_v2(
+    omega, 
+    lam = 1.085, 
+    cmap = cmap_refl, 
+    vmin = 0, 
+    vmax = 0.5, 
+    polar = True,
+    mask = mask,
+    title = 'Interactive map',
+    autoyscale = False,
+    ylim_sp = (0, 0.6),
+    )
+
+plt.show()
+
+##----------------------------------------------------------------------------------------
+# Search for the index of λ = 1.085 µm in the wavelength array
+i_lam = uf.where_closer(1.085, omega.lam)
+
+##----------------------------------------------------------------------------------------
+## Print completed
+print("\n\033[32;1m[OMEGA-Py test example completed]\033[0m")
+
+##----------------------------------------------------------------------------------------
+## End of code
+##----------------------------------------------------------------------------------------
diff --git a/docs/tests.md b/docs/tests.md
index d689f1e..3d26197 100644
--- a/docs/tests.md
+++ b/docs/tests.md
@@ -74,13 +74,39 @@ You can check the values stored in them by displaying the output of
 
 The final test would be to actually download and process an OMEGA observation.
 
-For instance, download the data from cube 3 of orbit 0979:
+For instance, download the data from cube 2 of orbit 0979:
 
- * [ORB0979.QUB](https://archives.esac.esa.int/psa/ftp/MARS-EXPRESS/OMEGA/MEX-M-OMEGA-2-EDR-FLIGHT-V1.0/DATA/ORB09/ORB0979_3.QUB)
- * [ORB0979.NAV](https://archives.esac.esa.int/psa/ftp/MARS-EXPRESS/OMEGA/MEX-M-OMEGA-2-EDR-FLIGHT-V1.0/DATA/GEM09/ORB0979_3.NAV)
+ * [ORB0979_2.QUB](https://archives.esac.esa.int/psa/ftp/MARS-EXPRESS/OMEGA/MEX-M-OMEGA-2-EDR-FLIGHT-V1.0/DATA/ORB09/ORB0979_2.QUB)
+ * [ORB0979_2.NAV](https://archives.esac.esa.int/psa/ftp/MARS-EXPRESS/OMEGA/MEX-M-OMEGA-2-EDR-FLIGHT-V1.0/DATA/GEM09/ORB0979_2.NAV)
+
+Then download the [`example_test_script.py`](example_test_script.py) Python file and run it
+to load, process, and display this OMEGA observation.
+
+??? abstract "Source code for [`example_test_script.py`](example_test_script.py)"
+    ~~~python title="example_test_script.py"
+    --8<-- "docs/example_test_script.py"
+    ~~~
+
+This example script will:
+
+ * load the binary files;
+ * apply the atmospheric and thermal corrections to the OMEGA observation;
+ * compute the 1.5μm band depth criteria (see the [band depth example](../example1_band_depth));
+ * display and save the reflectance, 1.5μm BD, and overplotted maps 
+   (see the [overplotting example](../example2_overplot));
+ * display the interactive visualization tool (if running the code in an IPython shell).
+
+You can then check that you have reproduced the figures shown in the examples.
+
+!!! tip "Notes"
+    If running the code directly from a bash terminal (with `python3 example_test_script.py`),
+    the data will be processed, and the figures saved in your current directory.
+    But you will not have access to the interactive visualization of the cube.
+
+    To use the interactive aspect of *OMEGA-Py*, open an IPython terminal to run the code
+    with `#!python run example_test_script.py` (or execute it from your favorite IDE with
+    a IPython console).
 
-Then run the code of the [basic usage](../basic_usage/) example. 
-You should be able to load, process, and display this OMEGA observation.
 
 !!! failure
     * If encountering an error in the data loading/processing, check your installation with the above steps.