diff --git a/python_scripts/example/mann_functions.py b/python_scripts/example/mann_functions.py
new file mode 100644
index 0000000000000000000000000000000000000000..7776dec4466e8fa612ba22082ad499670c203386
--- /dev/null
+++ b/python_scripts/example/mann_functions.py
@@ -0,0 +1,318 @@
+'''
+Created on 03/06/2014
+
+@author: MMPE
+@ adapt from wetb: DAVCON 05/02/2019
+'''
+
+
+import os
+from scipy.interpolate import RectBivariateSpline
+import numpy as np
+from spectra import spectra, logbin_spectra, plot_spectra, detrend_wsp
+
+# Look-Up Tables
+sp1, sp2, sp3, sp4 = np.load("C:/Users/davcon/Desktop/0_WakeIndicators/Mann_model/mann_spectra_data.npy")
+yp = np.arange(-3, 3.1, 0.1)
+xp = np.arange(0, 5.1, 0.1)
+RBS1 = RectBivariateSpline(xp, yp, sp1)
+RBS2 = RectBivariateSpline(xp, yp, sp2)
+RBS3 = RectBivariateSpline(xp, yp, sp3)
+RBS4 = RectBivariateSpline(xp, yp, sp4)
+
+
+
+def get_mann_model_spectra(ae, L, G, k1):
+ """Mann model spectra
+
+ Parameters
+ ----------
+ ae : int or float
+ Alpha epsilon^(2/3) of Mann model
+ L : int or float
+ Length scale of Mann model
+ G : int or float
+ Gamma of Mann model
+ k1 : array_like
+ Desired wave numbers
+
+ Returns
+ -------
+ uu : array_like
+ The u-autospectrum of the wave numbers, k1
+ vv : array_like
+ The v-autospectrum of the wave numbers, k1
+ ww : array_like
+ The w-autospectrum of the wave numbers, k1
+ uw : array_like
+ The u,w cross spectrum of the wave numbers, k1
+ """
+ xq = np.log10(L * k1)
+ yq = (np.zeros_like(xq) + G)
+ f = L ** (5 / 3) * ae
+ uu = f * RBS1.ev(yq, xq)
+ vv = f * RBS2.ev(yq, xq)
+ ww = f * RBS3.ev(yq, xq)
+ uw = f * RBS4.ev(yq, xq)
+ return uu, vv, ww, uw
+
+
+
+
+def _local_error(x, k1, uu, vv, ww=None, uw=None):
+
+ ae, L, G = x
+ val = 10 ** 99
+ if ae >= 0 and G >= 0 and G <= 5 and L > 0 and np.log10(k1[0] * L) >= -3 and np.log10(k1[0] * L) <= 3:
+ tmpuu, tmpvv, tmpww, tmpuw = get_mann_model_spectra(ae, L, G, k1)
+ val = np.sum((k1 * uu - k1 * tmpuu) ** 2)
+ if vv is not None:
+ val += np.sum((k1 * vv - k1 * tmpvv) ** 2)
+ if ww is not None:
+ val += np.sum((k1 * ww - k1 * tmpww) ** 2) + np.sum((k1 * uw - k1 * tmpuw) ** 2)
+ return val
+
+def fit_mann_model_spectra(k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2, min_bin_count=2, start_vals_for_optimisation=(0.01, 50, 3.3), plt=False):
+ """Fit a mann model to the spectra
+
+ Bins the spectra, into logarithmic sized bins and find the mann model parameters,
+ that minimize the error between the binned spectra and the Mann model spectra
+ using an optimization function
+
+ Parameters
+ ----------
+ k1 : array_like
+ Wave numbers
+ uu : array_like
+ The u-autospectrum of the wave numbers, k1
+ vv : array_like, optional
+ The v-autospectrum of the wave numbers, k1
+ ww : array_like, optional
+ The w-autospectrum of the wave numbers, k1
+ uw : array_like, optional
+ The u,w cross spectrum of the wave numbers, k1
+ log10_bin_size : int or float, optional
+ Bin size (log 10, based)
+ start_vals_for_optimization : (ae, L, G), optional
+ - ae: Alpha epsilon^(2/3) of Mann model\n
+ - L: Length scale of Mann model\n
+ - G: Gamma of Mann model
+
+ Returns
+ -------
+ ae : int or float
+ Alpha epsilon^(2/3) of Mann model
+ L : int or float
+ Length scale of Mann model
+ G : int or float
+ Gamma of Mann model
+
+ Examples
+ --------
+ >>> sf = sample_frq / u_ref
+ >>> u,v,w = # u,v,w wind components
+ >>> ae, L, G = fit_mann_model_spectra(*spectra(sf, u, v, w))
+ >>> u1,v1 = # u,v wind components
+ >>> ae, L, G = fit_mann_model_spectra(*spectra(sf, u, v))
+ """
+ from scipy.optimize import fmin
+ x = fmin(_local_error, start_vals_for_optimisation, logbin_spectra(k1, uu, vv, ww, uw, log10_bin_size, min_bin_count), disp=False)
+
+ if plt:
+ if not hasattr(plt, 'plot'):
+ import matplotlib.pyplot as plt
+# plot_spectra(k1, uu, vv, ww, uw, plt=plt)
+# plot_mann_spectra(*x, plt=plt)
+ ae, L, G = x
+ plot_fit(ae, L, G, k1, uu, vv, ww, uw, log10_bin_size=log10_bin_size, plt=plt)
+ plt.title('ae:%.3f, L:%.1f, G:%.2f' % tuple(x))
+ plt.xlabel('Wavenumber $k_{1}$ [$m^{-1}$]')
+ plt.ylabel(r'Spectral density $k_{1} F(k_{1})/U^{2} [m^2/s^2]$')
+ plt.legend()
+ plt.show()
+ return x
+
+def residual(ae, L, G, k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2):
+ """Fit a mann model to the spectra
+
+ Bins the spectra, into logarithmic sized bins and find the mann model parameters,
+ that minimize the error between the binned spectra and the Mann model spectra
+ using an optimization function
+
+ Parameters
+ ----------
+ ae : int or float
+ Alpha epsilon^(2/3) of Mann model
+ L : int or float
+ Length scale of Mann model
+ G : int or float
+ Gamma of Mann model
+ k1 : array_like
+ Wave numbers
+ uu : array_like
+ The u-autospectrum of the wave numbers, k1
+ vv : array_like, optional
+ The v-autospectrum of the wave numbers, k1
+ ww : array_like, optional
+ The w-autospectrum of the wave numbers, k1
+ uw : array_like, optional
+ The u,w cross spectrum of the wave numbers, k1
+ log10_bin_size : int or float, optional
+ Bin size (log 10, based)
+ start_vals_for_optimization : (ae, L, G), optional
+ - ae: Alpha epsilon^(2/3) of Mann model\n
+ - L: Length scale of Mann model\n
+ - G: Gamma of Mann model
+
+ Returns
+ -------
+ residual : array_like
+ rms of each spectrum
+ """
+ k1_sp = np.array([sp for sp in logbin_spectra(k1, uu, vv, ww, uw, log10_bin_size) if sp is not None])
+ bk1, sp_meas = k1_sp[0], k1_sp[1:]
+ sp_fit = np.array(get_mann_model_spectra(ae, L, G, bk1))[:sp_meas.shape[0]]
+ return np.sqrt(((bk1 * (sp_meas - sp_fit)) ** 2).mean(1))
+
+
+def var2ae(variance, spatial_resolution, N, L, G):
+ """Fit alpha-epsilon to match variance of time series
+
+ Parameters
+ ----------
+ variance : array-like
+ variance of u vind component
+ spatial_resolution : int, float or array_like
+ Distance between samples in meters
+ - For turbulence boxes: 1/dx = Nx/Lx where dx is distance between points,
+ Nx is number of points and Lx is box length in meters
+ - For time series: Sample frequency / U
+ N : int
+ L : int or float
+ Length scale of Mann model
+ G : int or float
+ Gamma of Mann model
+
+ Returns
+ -------
+ ae : float
+ Alpha epsilon^(2/3) of Mann model that makes the energy of the model equal to the varians of u
+ """
+
+ k1 = np.logspace(1,10,1000)/100000000
+ def get_var(uu):
+ return np.trapz(2 * uu[:], k1[:])
+
+ v1 = get_var(get_mann_model_spectra(0.1, L, G, k1)[0])
+ v2 = get_var(get_mann_model_spectra(0.2, L, G, k1)[0])
+ ae = (variance - v1) / (v2 - v1) * .1 + .1
+ return ae
+
+
+
+def fit_ae(spatial_resolution, u, L, G, plt=False):
+ """Fit alpha-epsilon to match variance of time series
+
+ Parameters
+ ----------
+ spatial_resolution : int, float or array_like
+ Distance between samples in meters
+ - For turbulence boxes: 1/dx = Nx/Lx where dx is distance between points,
+ Nx is number of points and Lx is box length in meters
+ - For time series: Sample frequency / U
+ u : array-like
+ u vind component
+ L : int or float
+ Length scale of Mann model
+ G : int or float
+ Gamma of Mann model
+
+ Returns
+ -------
+ ae : float
+ Alpha epsilon^(2/3) of Mann model that makes the energy of the model equal to the varians of u
+ """
+ #if len(u.shape) == 1:
+ # u = u.reshape(len(u), 1)
+# if min_bin_count is None:
+# min_bin_count = max(2, 6 - u.shape[0] / 2)
+# min_bin_count = 1
+ def get_var(k1, uu):
+ l = 0 #128 // np.sqrt(u.shape[1])
+ return np.trapz(2 * uu[l:], k1[l:])
+
+ k1, uu = spectra(spatial_resolution, u)[:2]
+ v = get_var(k1,uu)
+ v1 = get_var(k1, get_mann_model_spectra(0.1, L, G, k1)[0])
+ v2 = get_var(k1, get_mann_model_spectra(0.2, L, G, k1)[0])
+ ae = (v - v1) / (v2 - v1) * .1 + .1
+# print (ae)
+#
+# k1 = spectra(sf, u)[0]
+# v1 = get_var(*logbin_spectra(k1, get_mann_model_spectra(0.1, L, G, k1)[0], min_bin_count=min_bin_count)[:2])
+# v2 = get_var(*logbin_spectra(k1, get_mann_model_spectra(0.2, L, G, k1)[0], min_bin_count=min_bin_count)[:2])
+# k1, uu = logbin_spectra(*spectra(sf, u), min_bin_count=2)[:2]
+# #variance = np.mean([detrend_wsp(u_)[0].var() for u_ in u.T])
+# v = get_var(k1, uu)
+# ae = (v - v1) / (v2 - v1) * .1 + .1
+# print (ae)
+ if plt is not False:
+ if not hasattr(plt, 'plot'):
+ import matplotlib.pyplot as plt
+ plt.semilogx(k1, k1 * uu, 'b-', label='uu')
+ k1_lb, uu_lb = logbin_spectra(*spectra(sf, u), min_bin_count=1)[:2]
+
+ plt.semilogx(k1_lb, k1_lb * uu_lb, 'r--', label='uu_logbin')
+ muu = get_mann_model_spectra(ae, L, G, k1)[0]
+ plt.semilogx(k1, k1 * muu, 'g', label='ae:%.3f, L:%.1f, G:%.2f' % (ae, L, G))
+ plt.legend()
+ plt.xlabel('Wavenumber $k_{1}$ [$m^{-1}$]')
+ plt.ylabel(r'Spectral density $k_{1} F(k_{1}) [m^2/s^2]$')
+ plt.show()
+ return ae
+
+
+def plot_fit(ae, L, G, k1, uu, vv=None, ww=None, uw=None, mean_u=1, log10_bin_size=.2, plt=None):
+# if plt is None:
+# import matplotlib.pyplot as plt
+ plot_spectra(k1, uu, vv, ww, uw, mean_u, log10_bin_size, plt)
+ plot_mann_spectra(ae, L, G, "-", mean_u, plt)
+
+
+
+def plot_mann_spectra(ae, L, G, style='-', u_ref=1, plt=None, spectra=['uu', 'vv', 'ww', 'uw']):
+ if plt is None:
+ import matplotlib.pyplot as plt
+ mf = 10 ** (np.linspace(-4, 1, 1000))
+ muu, mvv, mww, muw = get_mann_model_spectra(ae, L, G, mf)
+ plt.title("ae: %.3f, L: %.2f, G:%.2f"%(ae,L,G))
+ if 'uu' in spectra: plt.semilogx(mf, mf * muu * 10 ** 0 / u_ref ** 2, 'r' + style)
+ if 'vv' in spectra: plt.semilogx(mf, mf * mvv * 10 ** 0 / u_ref ** 2, 'g' + style)
+ if 'ww' in spectra: plt.semilogx(mf, mf * mww * 10 ** 0 / u_ref ** 2, 'b' + style)
+ if 'uw' in spectra: plt.semilogx(mf, mf * muw * 10 ** 0 / u_ref ** 2, 'm' + style)
+
+
+#
+#if __name__ == "__main__":
+# import gtsdf
+# from geometry import wsp_dir2uv
+# #from wetb import wind
+# import matplotlib.pyplot as plt
+# # """Example of fitting Mann parameters to a "series" of a turbulence box"""
+# l = 1800
+# nx = 8192; ny, nz = 32,32; lx=0.2197265625;ly,lz=1,1;
+# sf = (nx / l)
+# s=1;
+# fn = r'./turb100%d%%s.bin'%s
+# u, v, w = [np.fromfile(fn % uvw, np.dtype('<f'), -1).reshape(nx , ny*nz) for uvw in ['u', 'v', 'w'] ]
+# ae, L, G = fit_mann_model_spectra(*spectra(sf, u, v, w), plt= True)
+#
+# u, v, w = [np.fromfile(fn % uvw, np.dtype('<f'), -1).reshape(nx , ny,nz) for uvw in ['u', 'v', 'w'] ]
+# # """ Y and Z vectors """
+# Y = np.linspace(-ly/2., ly/2., ny)
+# Z = np.linspace(-lz/2., lz/2., nz)
+# plt.close('all')
+# plt.ion()
+#
+
+
\ No newline at end of file