diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..363095d
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,11 @@
+test/test_windfield.asv
+test/.mat
+test/test_windshear.asv
+Simulations/2021_9T_Data/setup.mlx
+Simulations/2021_9T_Data/setup.mlx
+after_init_simulation_T.mat
+flowfield_xyz_1_steps.mat
+flowfield_xyz_2_steps.mat
+input_T_1_steps.mat
+input_T_2_steps.mat
+input_setUpTmpWFAndRun_2_steps.mat
diff --git a/FLORIDynCL/FLORIDynCL.m b/FLORIDynCL/FLORIDynCL.m
index d229726..d07a5af 100644
--- a/FLORIDynCL/FLORIDynCL.m
+++ b/FLORIDynCL/FLORIDynCL.m
@@ -53,6 +53,9 @@
 SimTime = Sim.StartTime;
 for it = 1:Sim.nSimSteps
     Sim.SimStep = it;
+    turbulence = T.States_T(201, 3);
+    display("Info: 1: "+turbulence);
+
     % ========== PREDICTION ==========
     % Iterate OPs and states
     T = iterateOPs(T,Sim,paramFLORIS,paramFLORIDyn);
@@ -68,11 +71,23 @@
     T.intOPs = interpolateOPs(T);
     
     %   ======= Set up temp Wind farms & run FLORIS
+    if it == 2
+        filename = "input_setUpTmpWFAndRun_"+Sim.nSimSteps+"_steps.mat";
+        save(filename, 'T', 'paramFLORIS', 'Wind');
+        display("Saved: "+filename);
+    end
     [tmpM,T] = setUpTmpWFAndRun(T,paramFLORIS,Wind);
+    turbulence = T.States_T(201, 3);
+    display("Info: 4: "+turbulence);
+
     M((it-1)*T.nT+1:it*T.nT,2:4)    = tmpM;
     M((it-1)*T.nT+1:it*T.nT,1)      = SimTime;
     T.States_T(T.StartI,3)          = tmpM(:,2);
     M_int{it}                       = T.red_arr;
+    turbulence = T.States_T(201, 3);
+    format long e
+    display(tmpM(:,2));
+    display("Info: 5: "+turbulence);
     
     % ========== Get Wind field variables & correct values
     [T,Wind] = correctVel(T,Wind,SimTime,paramFLORIS,tmpM);
diff --git a/Simulations/2021_9T_Data/setup.mlx b/Simulations/2021_9T_Data/setup.mlx
index 434fe7f..03b068c 100644
Binary files a/Simulations/2021_9T_Data/setup.mlx and b/Simulations/2021_9T_Data/setup.mlx differ
diff --git a/Visualization/PlotFlowField.m b/Visualization/PlotFlowField.m
index 54e4708..8525595 100644
--- a/Visualization/PlotFlowField.m
+++ b/Visualization/PlotFlowField.m
@@ -22,6 +22,11 @@
 function [Vis] = PlotFlowField(T,~,Wind,Sim,Vis,paramFLORIDyn,paramFLORIS,SimTime)
 %PLOTFLOWFIELD Generate full flow field plot
 % exportgraphics(gcf,'9T_flow.pdf','ContentType','vector')
+
+filename = "input_T_"+Sim.nSimSteps+"_steps.mat";
+save(filename, 'T');
+display("Saved: "+filename);
+
 %% Preallocate field
 nM = countMeasurements(Vis.FlowField.Data);
 nM = 3;
@@ -42,6 +47,14 @@
     Z = getMeasurements(X,Y,nM,T.posNac(1,3),T,paramFLORIS,Wind,Vis);
 end
 
+out = "v_min: " + min(min(Z(:,:,3)));
+out = sprintf('%s\n', out);
+out = out + "v_max: " + max(max(Z(:,:,3)));
+display(out);
+filename = "flowfield_xyz_"+Sim.nSimSteps+"_steps.mat";
+save(filename, 'X', 'Y', 'Z');
+display("Saved: "+filename);
+
 %% Generate & save plots
 % figure
 % contourf(X,Y,Z(:,:,1),40,'LineColor','none')
diff --git a/bin/count_loc b/bin/count_loc
new file mode 100755
index 0000000..087f80c
--- /dev/null
+++ b/bin/count_loc
@@ -0,0 +1,4 @@
+# Copyright (c) 2025 Uwe Fechner
+# SPDX-License-Identifier: BSD-3-Clause
+
+scc -x csv -x toml ..
\ No newline at end of file
diff --git a/main.m b/main.m
index 70f9d14..b936438 100644
--- a/main.m
+++ b/main.m
@@ -20,8 +20,9 @@
 % ======================================================================= %
 
 %% MainFLORIDyn Center-Line model
-% Improved FLORIDyn approach over the gaussian FLORIDyn model
+% % Improved FLORIDyn approach over the gaussian FLORIDyn model
 % Link folder with simulation data
+tic()
 pathToSimulation = '2021_9T_Data';
 
 %% Load data from the simulation
@@ -30,6 +31,7 @@
 
 % Get the settings for the wind field, visualization, controller and Sim.
 [Wind, Vis, Sim, Con] = setup();
+% Sim.Dyn.OPiteration='IterateOPs_average';
 
 % Add according functions to the search path
 addFLORISPaths;
@@ -48,47 +50,55 @@
 %% ====== Init simulation
 % Run initial conditions until no more change happens
 T = initSimulation(T,Wind,Sim,Con,Vis,paramFLORIDyn,paramFLORIS);
-
+Vis.FlowField.Plot.Online = false;
 %% ============ RUN SIMULATION ============
 % T := Simulation state (OP states, Turbine states, wind field states(OPs))
 % M := Measurements from the simulation (T_red, T_addedTI, T_Ueff, T_pow)
+toc()
 tic
+% Sim.EndTime = 20000;
+Sim.nSimSteps = 2;
+
+filename = "after_init_simulation_T.mat";
+save(filename, 'T');
+display("Saved: "+filename);
+
 [T,M,Vis,Mint] = FLORIDynCL(T,Wind,Sim,Con,Vis,paramFLORIDyn,paramFLORIS);
 t = toc;
 disp(['Sec. per sim. step: ' num2str(t/Sim.nSimSteps) ' with '...
     num2str(T.nT) ' turbine(s), total sim. time: ' num2str(t) ' s.'])
 %% Plotting & visualization
 
-% Measurements
-if Vis.Msmnts.Output
-    PlotMeasurements(T,M,Wind,Sim,Vis,paramFLORIDyn,paramFLORIS);
-    pause(0.1)
-end
+% % Measurements
+% if Vis.Msmnts.Output
+%     PlotMeasurements(T,M,Wind,Sim,Vis,paramFLORIDyn,paramFLORIS);
+%     pause(0.1)
+% end
 
 % Flow Field
 if Vis.FlowField.Plot.Post
-    PlotFlowField(T,M,Wind,Sim,Vis,paramFLORIDyn,paramFLORIS)
-    pause(0.1)
+    PlotFlowField(T,M,Wind,Sim,Vis,paramFLORIDyn,paramFLORIS);
+    pause(0.1);
 end
 
-% 3D Flow Field (.vtk for ParaView)
-if Vis.FlowField3D.Generate
-    gen3DFlowField(T,M,Wind,Sim,Vis,paramFLORIDyn,paramFLORIS);
-end
-
-% Final state to skip initialization
-if Sim.SaveFinalState
-    save([Sim.PathToSim 'T_final.mat'],'T')
-end
+% % 3D Flow Field (.vtk for ParaView)
+% if Vis.FlowField3D.Generate
+%     gen3DFlowField(T,M,Wind,Sim,Vis,paramFLORIDyn,paramFLORIS);
+% end
 
-% Save movie if online plotting was activated
-if Vis.FlowField.Plot.Online
-    Vis.Film.InProgress = false;
-    % Write saved frames as a movie to the simulation folder
-    v = VideoWriter(Vis.Film.MovFileEffU);
-    v.FrameRate = 20;
-    v.Quality   = 100;
-    open(v)
-    writeVideo(v,Vis.Film.FrmFileEffU)
-    close(v)
-end
+% % Final state to skip initialization
+% if Sim.SaveFinalState
+%     save([Sim.PathToSim 'T_final.mat'],'T')
+% end
+% 
+% % Save movie if online plotting was activated
+% if Vis.FlowField.Plot.Online
+%     Vis.Film.InProgress = false;
+%     % Write saved frames as a movie to the simulation folder
+%     v = VideoWriter(Vis.Film.MovFileEffU);
+%     v.FrameRate = 20;
+%     v.Quality   = 100;
+%     open(v)
+%     writeVideo(v,Vis.Film.FrmFileEffU)
+%     close(v)
+% end
diff --git a/test/create_random_numbers.m b/test/create_random_numbers.m
new file mode 100644
index 0000000..ea3dc76
--- /dev/null
+++ b/test/create_random_numbers.m
@@ -0,0 +1,8 @@
+% a vector.
+clear all
+rng(1234);
+N=1000;
+
+vec=randn(1,N);
+filename = 'test/randn.mat';
+save(filename, 'vec')
\ No newline at end of file
diff --git a/test/test_single_turbine.m b/test/test_single_turbine.m
new file mode 100644
index 0000000..5a2f5ee
--- /dev/null
+++ b/test/test_single_turbine.m
@@ -0,0 +1,49 @@
+clear all
+addpath('./WindField/Velocity_Interpolation_wErrorCov/');
+
+% Simple linear data: speed = 5 at t=0, speed = 15 at t=10
+times = [0.0, 10.0];
+speeds = [5.0, 15.0];
+data = [times(:), speeds(:)];  % Column-wise concatenation
+
+% Simple 2x2 identity Cholesky: noise still generated but uncorrelated
+chol = chol(eye(2));  % chol returns upper triangular by default in MATLAB
+
+% Define a WindVelMatrix function equivalent
+wind = WindVelMatrix(data, chol);
+
+% Test 5: Single turbine (scalar iT input)
+iT_scalar = 1;
+vel_scalar = getWindSpeedT(wind, iT_scalar, 5.0);
+
+% assert(length(vel_scalar) == 1);
+% assert(abs(vel_scalar(1) - 10.0) <= 3.0);
+
+
+%% --- Supporting functions ---
+
+function wind = WindVelMatrix(data, chol)
+    % Store data and chol in a struct (MATLAB equivalent of a Julia type)
+    wind.Data = data;
+    wind.CholSig = chol;
+end
+
+% function vel = getWindSpeedT(wind, iT, tQuery)
+%     % Simple linear interpolation to get speed at time tQuery for turbine iT
+%     % Assuming wind.data(:,1) = times, wind.data(:,2) = speeds
+% 
+%     times = wind.data(:,1);
+%     speeds = wind.data(:,2);
+% 
+%     % Interpolate speed linearly at time tQuery
+%     speed_t = interp1(times, speeds, tQuery, 'linear');
+% 
+%     % Add noise using chol matrix (noise generation simplified)
+%     noise = wind.chol * randn(2,1);
+% 
+%     % For this example just take the first element noise as perturbation
+%     speed_noisy = speed_t + noise(1);
+% 
+%     % Return speed as scalar (for single turbine)
+%     vel = speed_noisy;
+% end
diff --git a/test/test_windfield.m b/test/test_windfield.m
new file mode 100644
index 0000000..5ccc2f4
--- /dev/null
+++ b/test/test_windfield.m
@@ -0,0 +1,136 @@
+%     dir_mode = Direction_Constant()
+%     WindDir = 270
+%     iT = [1, 2, 3]
+%     phi = getWindDirT(dir_mode, WindDir, iT, nothing)
+%     for ph in phi
+%         @test ph ≈ 270.0
+%     end
+
+clear all
+addpath('./WindField/Direction_Constant');
+WindDir = 270;
+iT = [1, 2, 3];
+phi = getWindDirT(WindDir, iT);
+
+tol = 1e-8; % set a tolerance for floating-point comparison
+for i = 1:length(phi)
+    ph = phi(i);
+    assert(abs(ph - 270.0) < tol);
+end
+rmpath('./WindField/Direction_Constant')
+
+%     dir_mode = Direction_Constant_wErrCov()
+%     struct WindDirType
+%         Data::Float64
+%         CholSig::Matrix{Float64}
+%     end
+%     WindDir = WindDirType(270.0, cholesky(Matrix{Float64}(I, 3, 3)).L)
+%     iT = [1, 2, 3]
+%     phi = getWindDirT(dir_mode, WindDir, iT)
+%     result = [269.6402710931765, 271.0872084924286, 269.5804103830612]
+%     for (i, ph) in pairs(phi)
+%         @test ph ≈ result[i]
+%     end
+
+clear all
+addpath('./WindField/Direction_Constant_wErrorCov');
+rng(1234);
+
+% Cholesky factor of 3x3 identity
+L = chol(eye(3), 'lower');
+
+% WindDir structure
+WindDir.Data = 270.0;
+WindDir.CholSig = L;
+
+iT = [1, 2, 3]';
+
+% Expected result
+result = [269.0527533560426, 270.54014974707036, 269.78339785902375];
+
+phi = getWindDirT(WindDir, iT);
+for i = 1:length(phi)
+    assert(abs(phi(i) - result(i)) < 1e-8, 'Test failed at index %d', i);
+end
+rmpath('./WindField/Direction_Constant_wErrorCov');
+
+%   dir_mode = Direction_EnKF_InterpTurbine()
+%   # Suppose WindDir is a matrix where each row is [time, phi_T0, phi_T1, ...]
+%     WindDir = [
+%         0.0  10.0  20.0
+%         1.0  12.0  22.0
+%         2.0  14.0  24.0
+%     ]
+%     phi = getWindDirT_EnKF(dir_mode, WindDir, 1, 0.5)
+%     @test phi ≈ 11.0
+
+clear all
+addpath('./WindField/Direction_EnKF_InterpTurbine');
+WindDir = [
+    0.0  10.0  20.0
+    1.0  12.0  22.0
+    2.0  14.0  24.0
+];
+phi = getWindDirT_EnKF(WindDir, 1, 0.5);
+assert(abs(phi - 11.0) < 1e-8, 'Test failed');
+rmpath('./WindField/Direction_EnKF_InterpTurbine');
+
+% dir_mode = Direction_Interpolation_wErrorCov()
+% 
+% # Example wind direction data (time, phi)
+% wind_data = [
+%     0.0  10.0;
+%     5.0  20.0;
+%     10.0 30.0
+% ]
+% 
+% # Example Cholesky factor (for 2 turbines)
+% chol_sig = cholesky([1.0 0.5; 0.5 1.0]).L
+% 
+% # Create WindDir instance
+% WindDir = WindDirMatrix(wind_data, chol_sig)
+% 
+% # Example turbine indices (for 2 turbines)
+% iT = [1, 2]
+% 
+% # Example time
+% t = 7.0
+% 
+% # Call the function
+% phi = getWindDirT(dir_mode, WindDir, iT, t)
+% @test size(phi) == (2,1)
+% @test phi[1] ≈ 24.92903352636283
+% @test phi[2] ≈ 24.363944731838128
+
+clear all
+addpath('./WindField/Direction_Interpolation_wErrorCov');
+
+% Wind direction data: [time, direction]
+wind_data = [
+    0.0, 10.0;
+    5.0, 20.0;
+    10.0, 30.0
+];
+
+% Cholesky factor (lower triangular)
+chol_sig = chol([1.0 0.5; 0.5 1.0], 'lower');
+
+% Structure to hold data
+WindDir.Data = wind_data;
+WindDir.CholSig = chol_sig;
+
+% Turbine indices and time
+iT = [1, 2]';
+t = 7.0;
+
+phi = getWindDirT(WindDir, iT, t);
+
+% Test
+assert(isequal(size(phi), [2,1]));
+assert(abs(phi(1) - 25.84752589085377) < 1e-6);
+assert(abs(phi(2) - 25.140544918823198128) < 1e-6);
+
+rmpath('./WindField/Direction_Interpolation_wErrorCov');
+
+disp('All tests passed.');
+
diff --git a/test/test_windshear.m b/test/test_windshear.m
new file mode 100644
index 0000000..47c0f8a
--- /dev/null
+++ b/test/test_windshear.m
@@ -0,0 +1,29 @@
+% shear_mode = Shear_Interpolation()
+% WindShear = [10 0.8; 20 0.9; 30 1.0]  # Example data
+% z = [5, 15, 25, 35]                   # Heights to interpolate at
+% shear = getWindShearT(shear_mode, WindShear, z)
+% @test shear[1] ≈ 0.8
+% @test shear[2] ≈ 0.85
+% @test shear[3] ≈ 0.95
+% @test shear[4] ≈ 1.0
+
+clear all
+addpath('./WindField/Shear_Interpolation');
+% Example data
+WindShear = [10 0.8; 20 0.9; 30 1.0];
+
+% Heights to interpolate at
+z = [5, 15, 25, 35];
+
+% Interpolation (linear by default)
+shear = getWindShearT(WindShear, z);
+
+% Tests (using assert with a tolerance)
+assert(abs(shear(1) - 0.8) < 1e-6)
+assert(abs(shear(2) - 0.85) < 1e-6)
+assert(abs(shear(3) - 0.95) < 1e-6)
+assert(abs(shear(4) - 1.0) < 1e-6)
+rmpath('./WindField/Shear_Interpolation')
+
+disp('All tests passed.');
+