diff --git a/opengrid/datasets/gas_012017_30s.pkl b/opengrid/datasets/gas_012017_30s.pkl
new file mode 100644
index 0000000..06a38f2
Binary files /dev/null and b/opengrid/datasets/gas_012017_30s.pkl differ
diff --git a/opengrid/library/analysis.py b/opengrid/library/analysis.py
index 3fa9a4a..c5e4cbb 100644
--- a/opengrid/library/analysis.py
+++ b/opengrid/library/analysis.py
@@ -11,6 +11,7 @@
 import numpy as np
 import numbers
 from opengrid.library.exceptions import EmptyDataFrame
+from collections import namedtuple
 
 
 class Analysis(object):
@@ -197,3 +198,65 @@ def load_factor(ts, resolution=None, norm=None):
     lf = ts / norm
 
     return lf
+
+
+def load_duration(df, trim_zeros=False):
+    """
+    Create descending load duration series
+    (mainly for use in a load duration curve)
+
+    Parameters
+    ----------
+    df : pd.DataFrame or pd.Series
+    trim_zeros : bool
+        trim trailing zero's
+
+    Returns
+    -------
+    pd.DataFrame or pd.Series
+    """
+    df = pd.DataFrame(df)  # in case a series is passed, wrap it in a dataframe
+    load_durations = (df[column].reset_index(drop=True).sort_values(ascending=False).reset_index(drop=True) for column in df)
+    if trim_zeros:
+        load_durations = (np.trim_zeros(s, trim='b') for s in load_durations)
+    df = pd.concat(load_durations, axis=1)
+    result = df.squeeze()
+    return result
+
+
+def modulation_detection(ts, min_level=0.1):
+    """
+    Detect the modulation levels of a gas boiler
+
+    Parameters
+    ----------
+    ts : pd.Series
+    min_level : float
+        Physically, a gas boiler cannot modulate under a certain percentage of its maximum power
+        So we use this percentage to cut off any noise
+
+    Returns
+    -------
+    namedtuple(median, minimum)
+    """
+    # drop all values below the minimum level
+    ts = ts[ts >= (ts.max() * min_level)]
+
+    # load duration curve
+    ld = load_duration(ts, trim_zeros=True)
+
+    # find the part in the load duration curve with the highest number of consecutive identical values
+    # a.k.a. find the longest 'flat part' in the curve
+    median_modul = ld.round().groupby(ld.round()).size().sort_values(ascending=False).index[0]
+
+    # take the second derivative of the whatever happens after the flat part
+    dif2 = ld[ld < median_modul].diff().diff().shift(-2)
+    # find the maximum in this second derivative
+    # this is where the curve 'drops off'
+    min_modul_ix = dif2.sort_values(ascending=False).index[0]
+    min_modul = ld[min_modul_ix]
+
+    # return as a namedtuple
+    ModulationLevel = namedtuple('ModulationLevel', ['median', 'minimum'])
+    ml = ModulationLevel(median_modul, min_modul)
+    return ml
diff --git a/opengrid/library/plotting.py b/opengrid/library/plotting.py
index 4136e06..f659c50 100644
--- a/opengrid/library/plotting.py
+++ b/opengrid/library/plotting.py
@@ -1,10 +1,7 @@
 import os
-import os
 import numpy as np
 import pandas as pd
 import matplotlib
-import pandas as pd
-import numpy as np
 import matplotlib.cm as cm
 import matplotlib.pyplot as plt
 from matplotlib.dates import date2num, num2date, HourLocator, DayLocator, AutoDateLocator, DateFormatter
@@ -124,7 +121,6 @@ def carpet(timeseries, **kwargs):
 
     return im
 
-
 def boxplot(df, plot_mean=False, plot_ids=None, title=None, xlabel=None, ylabel=None):
     """
     Plot boxplots
@@ -175,4 +171,3 @@ def boxplot(df, plot_mean=False, plot_ids=None, title=None, xlabel=None, ylabel=
         plt.ylabel(ylabel)
 
     return plt.gcf()
-
diff --git a/opengrid/tests/test_analyses.py b/opengrid/tests/test_analyses.py
index b18633c..3786f9a 100644
--- a/opengrid/tests/test_analyses.py
+++ b/opengrid/tests/test_analyses.py
@@ -64,6 +64,13 @@ def test_load_factor(self):
         self.assertIsInstance(ts, pd.Series)
         self.assertAlmostEqual(175.0345212009457, (lf2 * 800).iloc[0])
 
+    def test_load_duration(self):
+        ts = pd.Series([1, 5, 7, 0, 6, 0, 3, 2])
+        ld = og.analysis.load_duration(ts)
+        self.assertTrue(ld.equals(pd.Series([7, 6, 5, 3, 2, 1, 0, 0])))
+        ld2 = og.analysis.load_duration(ts, trim_zeros=True)
+        self.assertTrue(ld2.equals(pd.Series([7, 6, 5, 3, 2, 1])))
+
 
 if __name__ == '__main__':
     unittest.main()