diff --git a/docs/example1_band_depth.md b/docs/example1_band_depth.md
index be76e6f..833d0e9 100644
--- a/docs/example1_band_depth.md
+++ b/docs/example1_band_depth.md
@@ -1,2 +1,161 @@
-**Coming soon.**
+Here is an example on how to compute and display a band depth on an OMEGA observation using
+OMEGA-Py.
+
+Let's assume you have downloaded the OMEGA data cube *ORB0979_2* from the PSA
+(or from these links:
+[QUB](https://archives.esac.esa.int/psa/ftp/MARS-EXPRESS/OMEGA/MEX-M-OMEGA-2-EDR-FLIGHT-V1.0/DATA/ORB09/ORB0979_2.QUB),
+[NAV](https://archives.esac.esa.int/psa/ftp/MARS-EXPRESS/OMEGA/MEX-M-OMEGA-2-EDR-FLIGHT-V1.0/DATA/GEM09/ORB0979_2.NAV))
+and loaded the module and data as follows:
+~~~python
+import omegapy.omega_data as od
+import omegapy.omega_plots as op
+~~~
+
+## Step 1 - Loading the data
+Option 1: Load the binary files and apply the thermal and atmospheric corrections if needed.
+~~~python
+# Load the data cube
+omega0 = od.OMEGAdata('0979_2')
+# Apply thermal and atmospheric corrections
+omega = od.corr_therm_atm(omega0, npool=15) # Adjust npool according to your system
+~~~
+
+Option 2: Load directly a previously saved [`OMEGAdata`](../reference/omega_data/#omega_data.OMEGAdata) object
+(with or without the corrections already applied).
+~~~python
+# Load directly the already corrected data cube
+omega = od.autoload_omega('0979_2', therm_corr=True, atm_corr=True)
+~~~
+
+## Step 2 - Generating the data mask (optional)
+
+~~~python
+mask = od.omega_mask(
+ omega,
+ hide_128=True,
+ emer_lim=10,
+ inci_lim=70,
+ tempc_lim=-194,
+ limsat_c=500
+ )
+~~~
+
+## Step 3 - Computing band depth at 1.5 μm
+
+Let's compute the water ice 1.5 μm band depth as defined in Poulet et al. (2007)[^1]:
+
+[^1]: F. Poulet, C. Gomez, J.-P. Bibring, et al. (2007).
+Martian surface mineralogy from Observatoire pour la Minéralogie, l'Eau, les Glaces et l'Activité on board the Mars Express
+spacecraft (OMEGA/MEx): Global mineral maps.
+*JGR, 112*, E08S02.
+[doi:10.1029/2006JE002840](https://doi.org/10.1029/2006JE002840)
+
+~~~python
+nbd_15 = BD_omega(omega, [1.5, 1.51], 1.30, 1.71, norm=True)
+~~~
+
+??? note "Source code for function `#!python BD_omega(omega, lam0, lamc1, lamc2, norm=True)`"
+ ~~~python
+ import omegapy.useful_functions as uf
+ import numpy as np
+
+ def BD_omega(omega, lam0, lamc1, lamc2, norm=True):
+ """Compute the band depth on an OMEGA observation cube.
+ Continuum linear between lamc1 and lamc2.
+
+ If an array is passed as argument for a wavelength value, the average is used.
+
+ Parameters
+ ==========
+ omega : OMEGAdata
+ The OMEGA/MEx observation.
+ lam0 : float or array-like
+ The wavelength of the center of the band.
+ lamc1 : float or array-like
+ The wavelength of the bluer point for the continuum determination.
+ lamc2 : float or array-like
+ The wavelength of the redder point for the continuum determination.
+ norm : bool, optional (default True)
+ | True -> band_depth output is the normalized BD values.
+ | False -> band_depth output is the BD values.
+
+ Returns
+ =======
+ band_depth : 2D array
+ The array of the band depth values for the observation
+ (normalized or not depending on norm).
+ rf_c : 2D array
+ The value of the continuum used to measure the band depth.
+ """
+ if not omega.therm_corr:
+ print("\033[01;33mWarning: No thermal correction applied.\033[0m")
+ # Initialisation
+ refl_cube = omega.cube_rf
+ nx, ny, nlam = refl_cube.shape
+ # Conversion floats -> list
+ if isinstance(lam0, (int, float)):
+ lam0 = [lam0]
+ if isinstance(lamc1, (int, float)):
+ lamc1 = [lamc1]
+ if isinstance(lamc2, (int, float)):
+ lamc2 = [lamc2]
+ # Search for wavelength indexes
+ i_lam0 = uf.where_closer_array(lam0, omega.lam)
+ i_lamc1 = uf.where_closer_array(lamc1, omega.lam)
+ i_lamc2 = uf.where_closer_array(lamc2, omega.lam)
+ # Average wavelengths
+ lam0 = np.mean(omega.lam[i_lam0])
+ lamc1 = np.mean(omega.lam[i_lamc1])
+ lamc2 = np.mean(omega.lam[i_lamc2])
+ # Average reflectances
+ rf_band = np.mean(refl_cube[:, :, i_lam0], axis=2)
+ rf_c1 = np.mean(refl_cube[:, :, i_lamc1], axis=2)
+ rf_c2 = np.mean(refl_cube[:, :, i_lamc2], axis=2)
+ # Average continuum
+ rf_c = rf_c1 + (rf_c2 - rf_c1) * (lam0 - lamc1) / (lamc2 - lamc1)
+ # Compute BD over the OMEGA cube
+ if norm:
+ band_depth = (rf_c - rf_band) / rf_c
+ else:
+ band_depth = rf_c - rf_band
+ # Output
+ return band_depth
+ ~~~
+
+## Step 4 - Displaying the band depth map
+
+Now that we have computed the band depth map, we can display it with
+the [`show_data_v2`](../reference/omega_plots/#omega_plots.show_data_v2) function,
+as described in the [data visualization](../data_visualization/#reflectance-vs-previously-computed-data) page.
+
+!!! tip
+ Install and import the [cmocean](https://matplotlib.org/cmocean/) module to access some more nice colormaps:
+
+ ~~~python
+ import cmocean.cm as cmo
+ ~~~
+
+=== "1.5μm Band Depth"
+ ~~~python
+ op.show_data_v2(omega, data=nbd_15, cb_title="1.5μm BD", polar=True, cmap=cmo.ice, vmin=0, vmax=0.75)
+ ~~~
+
+
+ 
+
+ ORB0979_2 1.5μm BD – Polar projection
+
+
+
+=== "Surface reflectance"
+ ~~~python
+ op.show_omega_v2(omega, lam=1.085, polar=True, vmin=0, vmax=0.6)
+ ~~~
+
+
+ 
+
+ ORB0979_2 reflectance – Polar projection
+
+
diff --git a/docs/img/ORB0979_2__show_data_v2_nbd15_polar.png b/docs/img/ORB0979_2__show_data_v2_nbd15_polar.png
new file mode 100644
index 0000000..8de8b82
Binary files /dev/null and b/docs/img/ORB0979_2__show_data_v2_nbd15_polar.png differ
diff --git a/docs/img/ORB0979_2__show_omega_v2_alb108_polar.png b/docs/img/ORB0979_2__show_omega_v2_alb108_polar.png
new file mode 100644
index 0000000..39cd023
Binary files /dev/null and b/docs/img/ORB0979_2__show_omega_v2_alb108_polar.png differ