diff --git a/wafo/doc/rainflow_example.py b/wafo/doc/rainflow_example.py
new file mode 100644
index 0000000..105136b
--- /dev/null
+++ b/wafo/doc/rainflow_example.py
@@ -0,0 +1,178 @@
+#from tutor_init import *
+import numpy as np
+import logging
+import matplotlib.pyplot as plt
+import itertools
+# import sys
+import itertools
+log = logging.getLogger(__name__)
+log.setLevel(logging.DEBUG)
+
+MARKERS = ('o', 'x', '+', '.', '<', '>', '^', 'v')
+
+
+def set_windows_title(title='', log=None):
+    pass
+
+def plot_varying_symbols(x, y, color='red', size=5):
+    """
+    Create a plot with varying symbols
+    Parameters
+    ----------
+    x : numpy array with x data of the points
+    y : numpy array with y data of the points
+    color : color of the symbols
+
+    Returns
+    -------
+
+    """
+    markers = itertools.cycle(MARKERS)
+    for q, p in zip(x, y):
+        plt.plot(q, p, marker=markers.next(), linestyle='', color=color, markersize=size)
+
+def damage_vs_S(S, beta, K):
+    """
+    calculate the damage 1/N for a given stress S
+    Parameters
+    ----------
+    S       : Stress [Pa]
+    beta    : coefficient, typically 3
+    K       : constant
+
+    Returns
+    -------
+
+    """
+    return K * np.power(S, beta)
+
+# Section 4.3.1 Crossing intensity
+#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+import wafo.data as wd
+import wafo.objects as wo
+import wafo.misc as wm
+
+xx_sea = wd.sea()
+
+Tlength = xx_sea[-1, 0] - xx_sea[0, 0]
+beta = 3
+K1 = 6.5e-31
+Np = 200
+Tp = Tlength/Np
+A = 100e6
+log.info("setting sin wave with Tp={} and T={}".format(Tp, Tlength))
+Nc = 1.0 / damage_vs_S(A, beta, K1)
+damage = float(Np)/float(Nc)
+log.info("budget at S={} N={}: damage = {} ".format(A, Nc, damage))
+#xx_sea[:, 1] = A * np.cos(2 * np.pi * xx_sea[:, 0]/Tp)
+xx_sea[:, 1] *= 500e6
+
+log.info("loaded sea time series {}".format(xx_sea.shape))
+ts = wo.mat2timeseries(xx_sea)
+
+tp = ts.turning_points()
+mM = tp.cycle_pairs(kind='min2max')
+Mm = tp.cycle_pairs(kind='max2min')
+lc = mM.level_crossings(intensity=True)
+T_sea = ts.args[-1]-ts.args[0]
+
+#for i in dir(mM):
+#    print(i)
+
+
+ts1 = wo.mat2timeseries(xx_sea[:,:])
+tp1 = ts1.turning_points()
+sig_tp = ts.turning_points(h=0, wavetype='astm')
+try:
+    sig_cp = sig_tp.cycle_astm()
+    log.info("Successfully used cycle_astm")
+except AttributeError:
+    log.warning("Could use cycle_astm")
+    sig_cp = None
+tp1 = ts1.turning_points()
+tp2 = ts1.turning_points(wavetype='Mw')
+mM1 = tp1.cycle_pairs(kind='min2max')
+Mm1 = tp1.cycle_pairs(kind='max2min')
+
+tp_rfc = tp1.rainflow_filter(h=100e6)
+mM_rfc = tp_rfc.cycle_pairs()
+try:
+    mM_rfc_a = tp1.cycle_astm()
+except AttributeError:
+    mM_rfc_a = None
+tc1 = ts1.trough_crest()
+min_to_max = True
+rfc_plot = True
+if min_to_max:
+    m1, M1 = mM1.get_minima_and_maxima()
+    i_min_start = 0
+else:
+    m1, M1 = Mm1.get_minima_and_maxima()
+    i_min_start = 2
+
+m_rfc, M_rfc = mM_rfc.get_minima_and_maxima()
+#m_rfc_a, M_rfc_a = mM_rfc_a.get_minima_and_maxima()
+ts1.plot('b-')
+if rfc_plot:
+    plot_varying_symbols(tp_rfc.args[0::2], m_rfc, color='red', size=10)
+    plot_varying_symbols(tp_rfc.args[1::2], M_rfc, color='green', size=10)
+else:
+    plot_varying_symbols(tp.args[i_min_start::2], m1, color='red', size=10)
+    plot_varying_symbols(tp.args[1::2], M1, color='green', size=10)
+
+set_windows_title("Sea time series", log)
+
+plt.figure()
+plt.subplot(122),
+mM.plot()
+plt.title('min-max cycle pairs')
+plt.subplot(121),
+mM_rfc.plot()
+
+title = 'Rainflow filtered cycles'
+plt.title(title)
+set_windows_title(title)
+
+
+# Min-max and rainflow cycle distributions
+#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# import wafo.misc as wm
+ampmM_sea = mM.amplitudes()
+ampRFC_sea = mM_rfc.amplitudes()
+plt.figure()
+title = "s_n_curve"
+set_windows_title(title)
+S = np.linspace(1e6,1000e6)
+plt.loglog(S, damage_vs_S(S, beta, K1))
+plt.figure()
+plt.subplot(121)
+stress_range = (1,1e9)
+n_bins = 100
+wm.plot_histgrm(ampmM_sea, bins=n_bins, range=stress_range)
+plt.xlim(stress_range)
+ylim = plt.gca().get_ylim()
+plt.title('min-max amplitude distribution')
+plt.subplot(122)
+if sig_cp is not None:
+    wm.plot_histgrm(sig_cp[:, 0],bins=n_bins, range=stress_range)
+    plt.gca().set_ylim(ylim)
+    title = 'Rainflow amplitude distribution'
+    plt.title(title)
+    plt.semilogy
+    set_windows_title(title)
+
+    hist, bin_edges = np.histogram(sig_cp[:, 0], bins=n_bins, range=stress_range)
+
+    plt.figure()
+    title = "my_bins"
+    plt.title(title)
+    plt.title(title)
+    set_windows_title(title)
+    plt.semilogy
+    plt.bar(bin_edges[:-1], hist, width=stress_range[1]/n_bins)
+
+    print("damage min/max : {}".format(mM_rfc.damage([beta], K1)))
+
+    damage_rfc = K1 * np.sum(sig_cp[:, 0] ** beta)
+    print("damage rfc : {}".format(damage_rfc))
+plt.show('hold')