diff --git a/python_scripts/from_vicondar_to_cs_files.py b/python_scripts/from_vicondar_to_cs_files.py
new file mode 100644
index 0000000000000000000000000000000000000000..eb6532fbdd56ff89b37b38df18e54a9dc6fd3610
--- /dev/null
+++ b/python_scripts/from_vicondar_to_cs_files.py
@@ -0,0 +1,207 @@
+# -*- 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
+
+
+# For each file in the folder
+PATH_VICONDAR_OUTPUT = './'
+
+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)*(z_ref/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)
+
+ # Gets coordinates of the turbulence box
+ dy, dz, Ny, Nz = 6.5, 6.5, 32, 32
+ ly, lz = dy*Ny, dz*Nz
+ Yvalues = dy/2 + np.arange(0,(ny),1)*dy
+ Zvalues = dz/2 + np.arange(0,(nz),1)*dz
+
+ # Gets coordinates of the lidar scanning patterns (arbitrary example)
+ LidarSys = 'circ7p' # Name of the lidar
+
+ Rrotor = 89.5
+ 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]))]
+ y_lidar,z_lidar = Yvalues[idx_y]-cc,Zvalues[idx_z]-cc
+
+ # 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(z_ref,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' %3.9) # Mann spectral velocity tensor parameters
+ filename.write("\nInput['L']=" + '%.4f' %29.4) # 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