diff --git a/python_scripts/example/from_vicondar_to_cs_files.py b/python_scripts/example/from_vicondar_to_cs_files.py
new file mode 100644
index 0000000000000000000000000000000000000000..c7d3cc39d589d8e5d4857cd1b1720db48df2f1f0
--- /dev/null
+++ b/python_scripts/example/from_vicondar_to_cs_files.py
@@ -0,0 +1,215 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sat Mar 14 09:49:17 2020
+
+@author: davcon
+
+This script reads output files from the VICONDAR framework and
+create inputs to the constrained simulation framework.
+
+CS: stands for Constrained simulation framework
+
+"""
+
+import numpy as np
+import scipy.io
+import os
+from scipy.interpolate import interp1d
+import numpy as np
+import matplotlib.pyplot as plt
+import os
+import scipy.io as sio
+import scipy
+import scipy.optimize as optimization
+from scipy.interpolate import interp1d,interp2d
+from scipy.interpolate import griddata
+import pickle
+from scipy.optimize import minimize
+import fnmatch
+from datetime import datetime, timedelta
+import re
+from mann_parameters import *
+from mann_parameters import var2ae
+import pandas as pd
+
+
+# User-defined inputs
+PATH_VICONDAR_OUTPUT = './'
+
+# Turbulence box inputs
+dy, dz, Ny, Nz, Nx = 6.5, 6.5, 32, 32, 8192
+ly, lz = dy*Ny, dz*Nz
+Yvalues = dy/2 + np.arange(0,(Ny),1)*dy
+Zvalues = dz/2 + np.arange(0,(Nz),1)*dz
+Tsim = 600 # Simulation time in seconds
+
+# Arbritary lidar scanning patterns
+LidarSys = 'circ7p' # Name of the lidar
+Rrotor = 89.5 # Turbine parameters (rotor)
+hh = 119 # Turbine parameters (hub-height)
+frac = Rrotor* 0.85
+phi = np.deg2rad(60)
+lidar_7pcirc = np.array([-frac,-frac*np.cos(phi),0,frac*np.cos(phi),frac,frac*np.cos(phi),-frac*np.cos(phi)]),\
+ np.array([0,-frac*np.sin(phi),0,frac*np.sin(phi),0,-frac*np.sin(phi),frac*np.sin(phi)])
+lidar_7pcirc = lidar_7pcirc + Yvalues[int(len(Yvalues)/2)]
+idx_y = [np.argmin(np.abs(lidar_7pcirc[0][i]- Yvalues)) for i in range(len(lidar_7pcirc[0]))]
+idx_z = [np.argmin(np.abs(lidar_7pcirc[1][i]- Yvalues)) for i in range(len(lidar_7pcirc[0]))]
+cc = Yvalues[int(len(Yvalues)/2)] # center of the box
+y_lidar,z_lidar = Yvalues[idx_y]-cc,Zvalues[idx_z]-cc
+
+# Mann parameters
+L_mann = 29.4
+Gamma_Mann = 3.9
+
+
+for filename in os.listdir(PATH_VICONDAR_OUTPUT):
+ if fnmatch.fnmatch(filename, '*.mat'):
+ print(filename)
+
+ # Get wind field statitistcs and turbulence seed number etc. from the filename
+ temp=re.findall(r'\d+', filename)
+
+ # Shear
+ alpha = temp[0]
+
+ # Mean wind speed at hub-height
+ U_amb = int(temp[2])
+
+ # Ambient turbulence
+ TI_amb = temp[3]
+
+ # Turbulence seed
+ seed = temp[4]
+ tu_seed = int(seed)
+
+ # Probe volume size of the lidar
+ vsize = int(temp[6])
+
+ # Vertical wind profile
+ Uz = int(U_amb)*(Zvalues/hh)**(int(alpha)/100)
+
+ # Load VICONDAR outputs
+ test_mat=scipy.io.loadmat(PATH_VICONDAR_OUTPUT + '/%s' %filename)
+ u_lidar = test_mat['uval'][0].reshape(-1)
+ t_lidar = test_mat['time'][0].reshape(-1)
+
+
+ # Center the lidar pattern to the turbulence box
+ y_lidar = Yvalues[int(len(Yvalues)/2)] + y_lidar[~np.isnan(y_lidar)]
+ z_lidar = Zvalues[int(len(Zvalues)/2)] + z_lidar[~np.isnan(z_lidar)]
+
+ # Number of points per scan
+ ppt_scans = np.shape(u_lidar)[0]
+
+ # Number of scan within the time period (e.g., 10-min period)
+ nr_scans = np.shape(u_lidar[0].reshape(-1))[0]
+
+ # Initialize variables
+ u_lidar_tot = np.zeros(shape=(ppt_scans,nr_scans))
+ t_lidar_tot = np.zeros_like(u_lidar_tot); y_lidar_tot = np.zeros_like(u_lidar_tot)
+ z_lidar_tot = np.zeros_like(u_lidar_tot)
+
+ # Reshape vectors according to VICONDAR notation
+ for i in range(ppt_scans):
+ u_lidar_tot[i,:]= u_lidar[i].reshape(-1)
+ t_lidar_tot[i,:]= t_lidar[i].reshape(-1)
+ y_lidar_tot[i,:]= z_lidar[i].reshape(-1)
+ z_lidar_tot[i,:]= y_lidar[i].reshape(-1)
+
+ # Interpolate the vertical ambient wind speed at the location of the constraints
+ u_amb_z = interp1d(Zvalues,Uz)(z_lidar)
+
+ # Compute wind speed residuals to be input to the constrained simulation framework
+ u_lidar_res = np.asarray([(u_lidar_tot[i_coord,:] -u_amb_z[i_coord]) for i_coord in range(ppt_scans)])
+
+ # Reshape vectors according to the CS convention
+ t_cnst_array = (t_lidar_tot.T).reshape(-1)
+ u_res_cnst_array = (u_lidar_res.T).reshape(-1)
+ z_cnst_array = (z_lidar_tot.T).reshape(-1)
+ y_cnst_array = (y_lidar_tot.T).reshape(-1)
+
+ # Create a single array with all info to be input to the CS framework
+ constraint_array =np.zeros(shape=(len(t_cnst_array),4))
+ for iC in range(0,len(t_cnst_array)):
+ constraint_array[iC,:] = np.array([t_cnst_array[iC], y_cnst_array[iC],z_cnst_array[iC],u_res_cnst_array[iC]])
+
+ # Fit ae (of the Mann model) from TI statistics
+ var_u = ((int(TI_amb)/100)*int(U_amb))**2
+ ae = var2ae(var_u, L=29.4, G=3.9)
+
+ idSim_idx = 0
+
+ # Name of the 'wake/simulation' case
+ iWake = 'av'
+
+ # ==============================================
+ # Create the CS file for simulations
+
+ constraint_cells=[]
+ constraint_cells=['[%.4f, %.4f, %.4f, %.4f]' %( t_cnst_array[iC], y_cnst_array[iC],z_cnst_array[iC],u_res_cnst_array[iC]) for iC in range(0,len(t_cnst_array))]
+
+ title = "./cs_seed%s_%s_%s_vol%d" %(seed,iWake,LidarSys,vsize)
+
+ # Add .py to the end of the script.
+ title = title + '.py'
+
+ # Convert all letters to lower case.
+ title = title.lower()
+
+ # Remove spaces from the title.
+ title = title.replace(' ', '_')
+
+ # Create a file that can be written to.
+ filename = open(title, 'w')
+
+ # Write the data to the file.
+ filename.write('\n#Script for imposing constraints in a turbulence box, seed no.%s' %tu_seed)
+ filename.write('\n#INPUT TO CONSTRAINING PROCEDURE\n')
+ filename.write('\nimport pickle')
+ filename.write('\nimport pandas as pd')
+ filename.write('\nimport numpy as np')
+ filename.write('\nimport matplotlib.pyplot as plt')
+ filename.write('\nimport os')
+ filename.write('\nimport scipy.io')
+ filename.write('\nfrom scipy import stats')
+ filename.write('\nimport scipy.io as sio')
+ filename.write('\nfrom MannTensor_lib import MannTensor,TrapezoidalSum2D')
+ filename.write('\nimport cmath')
+ filename.write('\nfrom scipy.interpolate import interp1d')
+ filename.write('\nfrom ConstrainTurbulenceBox_U_S import ConstrainTurbulenceBox_U_S\n')
+ filename.write('\nfrom ConstrainTurbulenceBox_U_S import ConstrainTurbulenceBox_U_S\n')
+
+ filename.write('\n#Initialization of variables:')
+ filename.write('\nInput={} \n')
+ filename.write("\nInput['Udatafile']=" + "'./turb%s_ti%s_av_setbu.bin'" %(tu_seed,int(TI_amb)))
+ filename.write("\nInput['Vdatafile']=" + "'turb%s_av_setbv.bin'" %(tu_seed))
+ filename.write("\nInput['Wdatafile']=" + "'turb%s_av_setbw.bin'" %(tu_seed))
+ filename.write("\nInput['LidarSys']=" + "'%s'" %LidarSys)
+ filename.write("\nInput['CrossSpectrumFilenames']=" + "'./MannCrossSpectrum_Gamma'")
+ filename.write("\nInput['CalculateSpectrum']=" + '%d' %2)
+ filename.write("\nInput['SpectrumSteps']=" + '%d' %128)
+ filename.write("\nInput['UseBoxCovariance']=" + '%d' %0)
+ filename.write("\nInput['Gamma']=" + '%.2f' %Gamma_Mann) # Mann spectral velocity tensor parameters
+ filename.write("\nInput['L']=" + '%.4f' %L_mann) # Mann spectral velocity tensor parameters
+ filename.write("\nInput['AlphaEpsilon']=" + '%.4f' %ae) # Mann spectral velocity tensor parameters
+ filename.write("\nInput['Nx']=" + '%d' %Nx) # Mann turbulence box dimensions
+ filename.write("\nInput['Ny']=" + '%d' %Ny) # Mann turbulence box dimensions
+ filename.write("\nInput['Nz']=" + '%d' %Nz) # Mann turbulence box dimensions
+ filename.write("\nInput['T']=" + '%d' %Tsim) # Simulation time in seconds
+ filename.write("\nInput['Umean']=" + '%d' %U_amb) # Ambient wind speed at hub-height
+ filename.write("\nInput['BoxWidth']=" + '%d' %(Ny*dy)) # in meters
+ filename.write("\nInput['BoxHeight']=" + '%d' %(Nz*dz)) # in meters
+ filename.write("\nInput['SaveToFile']=" + '%d' %1)
+ filename.write("\nInput['PrintDetailedOutput']=" + '%d' %0)
+ filename.write("\nInput['CrossSpectrumComponent1']=" + '%d' %1)
+ filename.write("\nInput['CrossSpectrumComponent2']=" + '%d' %1)
+ filename.write("\nInput['OutputFileName']=" + "'./turb%s_ti%s_av_setb_%s_c_vol%du.bin'" %(tu_seed,int(TI_amb),LidarSys,vsize))
+ filename.write("\nInput['Constraints']= [")
+ filename.write(", \n".join(constraint_cells))
+ filename.write("]\n")
+ filename.write("\nprint('Start applying constraints:')")
+ filename.write("\nConstrainTurbulenceBox_U_S(Input)")
+
+ # Close the file after writing to it.
+ filename.close()
+
\ No newline at end of file