diff --git a/Validation/Luneberg/Error Convergence Analysis.csv b/Validation/Luneberg/Error Convergence Analysis.csv
new file mode 100644
index 0000000..6511ed8
--- /dev/null
+++ b/Validation/Luneberg/Error Convergence Analysis.csv	
@@ -0,0 +1,21 @@
+Grid Resolution,Avg normalized error,,Number of lightrays,GPU clock,CPU clock
+25,0.0257,,1.00E+03,1.12,0.4
+50,0.0107,,5.00E+03,1.15,0.65
+75,0.0065,,1.00E+04,1.41,1.05
+100,0.0049,,5.00E+04,1.67,6.33
+125,0.0037,,1.00E+05,1.7,11.02
+150,0.0031,,5.00E+05,4.52,95.45
+175,0.0025,,1.00E+06,11.86,192.8
+200,0.0022,,5.00E+06,68.7,1026.16
+225,0.0019,,1.00E+07,153.95,1819.67
+250,0.0017,,,,
+275,1.50E-03,,,,
+300,0.0013,,,,
+325,0.0012,,,,
+350,0.0011,,,,
+375,0.001,,,,
+400,9.62E-04,,,,
+425,8.93E-04,,,,
+450,8.25E-04,,,,
+475,7.70E-04,,,,
+500,7.26E-04,,,,
diff --git a/Validation/Luneberg/Error Convergence Analysis.xlsx b/Validation/Luneberg/Error Convergence Analysis.xlsx
new file mode 100644
index 0000000..8c2adbf
Binary files /dev/null and b/Validation/Luneberg/Error Convergence Analysis.xlsx differ
diff --git a/Validation/Luneberg/Luneberg Error.png b/Validation/Luneberg/Luneberg Error.png
new file mode 100644
index 0000000..de3ae75
Binary files /dev/null and b/Validation/Luneberg/Luneberg Error.png differ
diff --git a/Validation/Luneberg/Luneberg time elapsed.png b/Validation/Luneberg/Luneberg time elapsed.png
new file mode 100644
index 0000000..85927c0
Binary files /dev/null and b/Validation/Luneberg/Luneberg time elapsed.png differ
diff --git a/Validation/Luneberg/function_trace_through_density_grad_rk4_sharma_gpu.m b/Validation/Luneberg/function_trace_through_density_grad_rk4_sharma_gpu.m
new file mode 100644
index 0000000..7081154
--- /dev/null
+++ b/Validation/Luneberg/function_trace_through_density_grad_rk4_sharma_gpu.m
@@ -0,0 +1,136 @@
+function [lightfield] = function_trace_through_density_grad_rk4_sharma_gpu(lightfield, start_z, density_object)
+% GPU-native RK4 ray tracer using finite-difference gradients and interp3
+% No use of griddedInterpolant or gather (fully GPU-resident)
+
+tic
+interp_method = 'linear';
+
+% === Setup Grid and Refractive Index ===
+step_size = density_object.step_size;  % [um]
+delta = step_size;                     % finite difference step
+
+x = density_object.rho.x * 1000;       % mm → um
+y = density_object.rho.y * 1000;
+z = density_object.rho.z * 1000;
+rho = density_object.rho.rho;
+n_cpu = density_object.fluid_glasstone * rho + 1;
+
+[x, y, z] = deal(gpuArray(x), gpuArray(y), gpuArray(z));
+n_grid = gpuArray(n_cpu);              % move n to GPU
+
+x_min = min(x); x_max = max(x);
+y_min = min(y); y_max = max(y);
+z_min = min(z); z_max = max(z);
+
+% === Ray Initialization ===
+pos = gpuArray([lightfield.position.x, lightfield.position.y, lightfield.position.z]); % [um]
+dir = gpuArray([lightfield.direction.dx, lightfield.direction.dy, lightfield.direction.dz]);
+
+num_rays = size(pos,1);
+is_active = true(num_rays,1);
+prev_active_count = num_rays + 1;
+
+counter = 1;
+max_iter = 10e7;
+chunk_size = 1.0e6;
+
+while true || counter < max_iter
+    idx = find(is_active);
+    if isempty(idx); break; end
+    
+
+    for c = 1:ceil(length(idx)/chunk_size)
+
+        % Get chunk indices
+        chunk_start = (c-1)*chunk_size + 1;
+        chunk_end = min(c*chunk_size, length(idx));
+        chunk_idx = idx(chunk_start:chunk_end);
+
+        % Stage A
+        [n0, gradn0] = interpolate_n_and_gradient_gpu(x, y, z, n_grid, pos(chunk_idx,:), delta, interp_method);
+        delta_t = step_size ./ n0;
+        Tn = dir(chunk_idx,:) .* n0;
+        A = gradn0 .* n0 .* delta_t;
+
+        % Stage B
+        %pos1 = pos;
+        pos1(chunk_idx,:) = pos(chunk_idx,:) + (0.5 * Tn + 1/8 * A) .* delta_t;
+        [n1, gradn1] = interpolate_n_and_gradient_gpu(x, y, z, n_grid, pos1(chunk_idx,:), delta, interp_method);
+        clearvars pos1
+        B = gradn1 .* n1 .* delta_t;
+
+        % Stage C
+        %pos2 = pos;
+        pos2(chunk_idx,:) = pos(chunk_idx,:) + (Tn + 0.5 * B) .* delta_t;
+        [n2, gradn2] = interpolate_n_and_gradient_gpu(x, y, z, n_grid, pos2(chunk_idx,:), delta, interp_method);
+        clearvars pos2
+        C = gradn2 .* n2 .* delta_t;
+
+        % Final direction update
+        Tn_new = Tn + (1/6) * (A + 4*B + C);
+        %temp_pos = pos;
+        temp_pos(chunk_idx,:) = pos(chunk_idx,:) + (Tn + 1/6 * (A + 2*B)) .* delta_t;
+
+        % Final direction normalization
+        n_new = interp3(x, y, z, n_grid, temp_pos(chunk_idx,1), temp_pos(chunk_idx,2), temp_pos(chunk_idx,3), interp_method, 1.0);
+        dir(chunk_idx,:) = normalize_vector_gpu(Tn_new ./ n_new);
+
+        % Deactivate invalid rays
+        invalid = temp_pos(chunk_idx,1) < x_min | temp_pos(chunk_idx,1) > x_max | ...
+            temp_pos(chunk_idx,2) < y_min | temp_pos(chunk_idx,2) > y_max | ...
+            temp_pos(chunk_idx,3) < z_min | temp_pos(chunk_idx,3) > z_max;
+
+        nan_mask = any(isnan(temp_pos(chunk_idx,:)), 2) | ...
+            any(isnan(dir(chunk_idx,:)), 2) | ...
+            isnan(n0) | n0 <= 0;
+
+        is_active(chunk_idx(invalid | nan_mask)) = false;
+
+        %Position update
+        pos(chunk_idx,:) = temp_pos(chunk_idx,:);
+        clearvars temp_pos
+        curr_active_count = sum(is_active);
+        if curr_active_count ~= prev_active_count
+            fprintf('\n Number of lightrays remaining = %d', curr_active_count);
+            prev_active_count = curr_active_count;
+        end
+    end
+    counter = counter + 1;
+    if counter == max_iter
+        fprintf('\ncounter threshold met\n')
+    end
+end
+
+
+% === Output to CPU ===
+pos = gather(pos);
+dir = gather(dir);
+
+lightfield.position.x = pos(:,1);
+lightfield.position.y = pos(:,2);
+lightfield.position.z = pos(:,3);
+lightfield.direction.dx = dir(:,1);
+lightfield.direction.dy = dir(:,2);
+lightfield.direction.dz = dir(:,3);
+
+toc
+end
+
+function [n_val, grad_n] = interpolate_n_and_gradient_gpu(x, y, z, n_grid, pos, delta, interp_method)
+X = pos(:,1); Y = pos(:,2); Z = pos(:,3);
+
+n_val = interp3(x, y, z, n_grid, X, Y, Z, interp_method, 1.0);
+
+grad_n = zeros(size(pos), 'like', pos);
+grad_n(:,1) = (interp3(x, y, z, n_grid, X+delta, Y, Z, interp_method, 1.0) - ...
+               interp3(x, y, z, n_grid, X-delta, Y, Z, interp_method, 1.0)) ./ (2*delta);
+grad_n(:,2) = (interp3(x, y, z, n_grid, X, Y+delta, Z, interp_method, 1.0) - ...
+               interp3(x, y, z, n_grid, X, Y-delta, Z, interp_method, 1.0)) ./ (2*delta);
+grad_n(:,3) = (interp3(x, y, z, n_grid, X, Y, Z+delta, interp_method, 1.0) - ...
+               interp3(x, y, z, n_grid, X, Y, Z-delta, interp_method, 1.0)) ./ (2*delta);
+end
+
+function norm_vec = normalize_vector_gpu(v)
+mag = sqrt(sum(v.^2, 2));
+norm_vec = v ./ mag;
+end
\ No newline at end of file
diff --git a/Validation/Luneberg/function_trace_through_density_grad_rk4_sharma_while_loop.m b/Validation/Luneberg/function_trace_through_density_grad_rk4_sharma_while_loop.m
new file mode 100644
index 0000000..9cdc9d2
--- /dev/null
+++ b/Validation/Luneberg/function_trace_through_density_grad_rk4_sharma_while_loop.m
@@ -0,0 +1,131 @@
+function [lightfield] = function_trace_through_density_grad_rk4_sharma_while_loop(lightfield, start_z, density_object)
+% Finite-difference gradient RK4 ray tracer with dynamic stopping
+% Based on Sharma et al., but uses finite-difference gradients (no precomputed dn/dx)
+
+% INPUTS:
+%   lightfield: structure containing initial ray positions and directions
+%   start_z: starting z-plane [mm]
+%   density_object: structure with refractive index interpolant (Fn) and grid info
+
+% OUTPUT:
+%   lightfield: structure with updated ray positions and directions
+
+tic
+
+% Setup
+step_size = density_object.step_size; % [um]
+delta = step_size; % [um] finite difference step for gradient estimation
+
+x = density_object.rho.x * 1000; % Convert mm to um
+y = density_object.rho.y * 1000;
+z = density_object.rho.z * 1000;
+rho = density_object.rho.rho;
+n = density_object.fluid_glasstone * rho + 1;
+x_min = min(x);
+x_max = max(x);
+y_min = min(y);
+y_max = max(y);
+
+Fn = griddedInterpolant({x, y, z}, n, lower(density_object.interp_method));
+%Fn.ExtrapolationMethod = 'nearest'; 
+z_min = min(z);
+z_max = max(z);
+
+% Ray position and direction initialization
+pos = [lightfield.position.x, lightfield.position.y, lightfield.position.z]; % [um]
+dir = [lightfield.direction.dx, lightfield.direction.dy, lightfield.direction.dz]; % unit vectors
+
+num_rays = size(pos,1);
+is_active = true(num_rays,1);
+prev_active_count = num_rays + 1;
+counter = 1;
+while sum(is_active) > 0%0.01*num_rays
+    idx = find(is_active);
+
+    n0 = Fn(pos(idx,1), pos(idx,2), pos(idx,3));
+    gradn0 = finite_gradient(Fn, pos(idx,:), delta);
+    delta_t = step_size ./ n0;
+    Tn = dir(idx,:) .* n0;
+    %Calculate A
+    A = gradn0 .* n0 .* delta_t;
+
+    %Update Position and Calculate B
+    % pos1 = pos(idx,:) + (0.5 * Tn + 1/8 * A) .* delta_t;
+
+    pos1(idx,1) = pos(idx,1) + (0.5 * Tn(:,1) + 1/8 * A(:,1)) .* delta_t(:);
+    pos1(idx,2) = pos(idx,2) + (0.5 * Tn(:,2) + 1/8 * A(:,2)) .* delta_t(:);
+    pos1(idx,3) = pos(idx,3) + (0.5 * Tn(:,3) + 1/8 * A(:,3)) .* delta_t(:);
+    n1 = Fn(pos1(idx,1), pos1(idx,2), pos1(idx,3));
+    gradn1 = finite_gradient(Fn, pos1(idx,:), delta);
+    % disp(size(n1))
+    % disp(size(gradn1))
+    B = gradn1 .* n1 .* delta_t(:);
+
+    %Update Position and Calculate C
+    %pos2 = pos(idx,:) + (Tn + 0.5 * B) .* delta_t;
+    pos2 = pos;
+    pos2(idx,1) = pos(idx,1) + (Tn(:,1) + 0.5 * B(:,1)) .* delta_t(:);
+    pos2(idx,2) = pos(idx,2) + (Tn(:,2) + 0.5 * B(:,2)) .* delta_t(:);
+    pos2(idx,3) = pos(idx,3) + (Tn(:,3) + 0.5 * B(:,3)) .* delta_t(:);
+    n2 = Fn(pos2(idx,1), pos2(idx,2), pos2(idx,3));
+    gradn2 = finite_gradient(Fn, pos2(idx,:), delta);
+    C = gradn2 .* n2 .* delta_t(:);
+
+
+    %Calculate new positions and new directions
+    Tn_new = Tn + (1/6) * (A + 4*B + C);
+    temp_pos(idx,1) = pos(idx,1) + (Tn(:,1) + 1/6 * (A(:,1) + 2*B(:,1))) .* delta_t(:);
+    temp_pos(idx,2) = pos(idx,2) + (Tn(:,2) + 1/6 * (A(:,2) + 2*B(:,2))) .* delta_t(:);
+    temp_pos(idx,3) = pos(idx,3) + (Tn(:,3) + 1/6 * (A(:,3) + 2*B(:,3))) .* delta_t(:);
+    n_new = Fn(temp_pos(idx,1), temp_pos(idx,2), temp_pos(idx,3));
+    dir(idx,:) = normalize_vector(Tn_new ./ n_new);
+
+    % Deactivate rays that leave volume or encounter NaN
+    invalid = temp_pos(idx,1) < x_min | temp_pos(idx,1) > x_max | temp_pos(idx,2) < y_min | temp_pos(idx,2) > y_max;
+    exited = temp_pos(idx,3) > z_max;
+    nan_mask = any(isnan(temp_pos(idx,:)), 2) | any(isnan(dir(idx,:)), 2) | isnan(n0) | n0 <= 0;
+
+    % Mark these rays as inactive
+    is_active(idx(invalid | nan_mask | exited)) = false;
+    % is_active(idx(nan_mask | exited)) = false;
+    dir(invalid,:) = NaN;
+    pos(invalid,:) = NaN;
+    pos(idx,:) = temp_pos(idx,:);
+    
+    curr_active_count = sum(is_active);
+    if curr_active_count ~= prev_active_count
+        fprintf('\n Number of lightrays remaining = %d', curr_active_count);
+        prev_active_count = curr_active_count;
+    end
+    % if mod(counter,1000) == 0
+    % figure(1);
+    % % plot3(lightfield.position.x(1),lightfield.position.y(1),lightfield.position.z(1),'.')
+    % hold on
+    % plot3(pos(:,1), pos(:,2), pos(:,3),'.')
+    % end
+    counter = counter + 1;
+end
+
+% Update final lightfield
+lightfield.position.x = pos(:,1);
+lightfield.position.y = pos(:,2);
+lightfield.position.z = pos(:,3);
+lightfield.direction.dx = dir(:,1);
+lightfield.direction.dy = dir(:,2);
+lightfield.direction.dz = dir(:,3);
+
+toc
+end
+
+function grad_n = finite_gradient(Fn, pos, delta)
+N = size(pos,1);
+grad_n = zeros(N,3);
+grad_n(:,1) = (Fn(pos(:,1)+delta, pos(:,2), pos(:,3)) - Fn(pos(:,1)-delta, pos(:,2), pos(:,3))) ./ (2*delta);
+grad_n(:,2) = (Fn(pos(:,1), pos(:,2)+delta, pos(:,3)) - Fn(pos(:,1), pos(:,2)-delta, pos(:,3))) ./ (2*delta);
+grad_n(:,3) = (Fn(pos(:,1), pos(:,2), pos(:,3)+delta) - Fn(pos(:,1), pos(:,2), pos(:,3)-delta)) ./ (2*delta);
+end
+
+function norm_vec = normalize_vector(v)
+mag = sqrt(sum(v.^2,2));
+norm_vec = v ./ mag;
+end
diff --git a/Validation/Luneberg/luneberg_test.m b/Validation/Luneberg/luneberg_test.m
new file mode 100644
index 0000000..319e6c0
--- /dev/null
+++ b/Validation/Luneberg/luneberg_test.m
@@ -0,0 +1,83 @@
+clear
+clc
+close all
+
+tic
+
+%% Create Luneburg Lens
+R = 1;               % Radius of Luneburg lens [mm]
+N = 250;             % Number of points in x and y
+fluid_glasstone = 0.226e-3; % Gladstone-Dale constant of air
+
+x = linspace(-R, R, N);
+y = linspace(-R, R, N);
+z = linspace(-2*R-1e-6, 1e-6, N);
+Z_center = -R;
+rho = zeros(N, N, N); % preallocate output
+
+[Xs, Ys] = ndgrid(x, y);
+
+for k = 1:N
+    z_shifted = z(k) - Z_center;
+    r_slice = sqrt(Xs.^2 + Ys.^2 + z_shifted^2);
+
+    n_slice = ones(N, N); % background
+    inside = r_slice <= R;
+    n_slice(inside) = sqrt(2 - (r_slice(inside)/R).^2);
+
+    % Convert to density
+    rho(:,:,k) = (n_slice - 1) / fluid_glasstone;
+end
+
+% Package into density_object
+density_object.rho.x = x; %mm
+density_object.rho.y = y; %mm
+density_object.rho.z = z; %mm
+density_object.rho.rho = rho;
+density_object.fluid_glasstone = fluid_glasstone;
+density_object.step_size = 0.01*R*1000; % [um] step size for RK4
+density_object.interp_method = 'linear';
+start_z = -2*R;
+
+%% Create Lightfield 
+
+num_rays_total = 1e3; 
+
+% Generate random points inside a circular aperture of radius 0.99*R
+theta = 2*pi*rand(num_rays_total,1); % random angles [0, 2pi]
+r_rand = 0.99*R * sqrt(rand(num_rays_total,1));
+X0 = r_rand .* cos(theta); 
+Y0 = r_rand .* sin(theta);
+Z0 = -2*R*ones(size(X0));
+
+DX0 = zeros(size(X0)); 
+DY0 = zeros(size(X0));     
+DZ0 = ones(size(X0)); %all rays initially moving along +z
+
+% Create lightfield structure
+lightfield.position.x = X0(:) * 1000; 
+lightfield.position.y = Y0(:) * 1000;
+lightfield.position.z = Z0(:) * 1000;
+lightfield.direction.dx = DX0(:);
+lightfield.direction.dy = DY0(:);
+lightfield.direction.dz = DZ0(:);
+
+%% Step 3: Raytrace Through the Luneburg Lens
+lightfield = function_trace_through_density_grad_rk4_sharma_while_loop(lightfield, start_z, density_object);
+%lightfield = function_trace_through_density_grad_rk4_sharma_gpu(lightfield, start_z, density_object); 
+
+% figure;
+% plot3(lightfield.position.x/1000,lightfield.position.y/1000,lightfield.position.z,'.')
+% xlabel('x/R')
+% ylabel('y/R')
+
+%Calculate error
+distance = (sqrt(lightfield.position.x.^2+lightfield.position.y.^2)/(R*1000)); %Normalized by r
+alpha = atan2(lightfield.position.y,lightfield.position.x);
+alpha = mod(alpha,2*pi);
+distance = distance(~isnan(distance));
+spot_size = max(distance);
+avg_normalized_error = mean(distance);
+
+toc
+
diff --git a/Validation/Luneberg/plot_error_convergence.m b/Validation/Luneberg/plot_error_convergence.m
new file mode 100644
index 0000000..c0e747f
--- /dev/null
+++ b/Validation/Luneberg/plot_error_convergence.m
@@ -0,0 +1,30 @@
+%%Clean up
+clear
+clc
+close all
+
+A = csvread("Error Convergence Analysis.csv",1,0);
+B = A(:,4:6);
+A = A(:,1:2);
+
+figure
+loglog(A(:,1),A(:,2),'*-','LineWidth',2)
+xlabel('Number of grid points to define luneberg lens')
+ylabel('Average error (\Deltar/R)')
+grid on
+xlim([10 1000])
+ylim([1e-4 1e-1])
+%xticks(logspace(0, log10(500), 2)) 
+set(gca,'FontSize',16,'TickLabelInterpreter','Latex')
+set(gcf,'Position',[820   232   920   641])
+
+figure
+loglog(B(:,1),B(:,2),'*-','LineWidth',2)
+hold on
+loglog(B(:,1),B(:,3),'*-','LineWidth',2)
+xlabel('Number of rays for luneberg raytracing')
+ylabel('Elapsed time (s)')
+grid on
+legend('GPU','CPU')
+set(gca,'FontSize',16,'TickLabelInterpreter','Latex')
+set(gcf,'Position',[820   232   920   641])
\ No newline at end of file
diff --git a/Validation/Mixing Layer/CFD - BOS comparison qualitative.jpg b/Validation/Mixing Layer/CFD - BOS comparison qualitative.jpg
new file mode 100644
index 0000000..26743ef
Binary files /dev/null and b/Validation/Mixing Layer/CFD - BOS comparison qualitative.jpg differ
diff --git a/Validation/Mixing Layer/CFD - BOS comparison quantitative.jpg b/Validation/Mixing Layer/CFD - BOS comparison quantitative.jpg
new file mode 100644
index 0000000..1c80e5e
Binary files /dev/null and b/Validation/Mixing Layer/CFD - BOS comparison quantitative.jpg differ
diff --git a/Validation/Mixing Layer/CFD density grad.jpg b/Validation/Mixing Layer/CFD density grad.jpg
new file mode 100644
index 0000000..68240b2
Binary files /dev/null and b/Validation/Mixing Layer/CFD density grad.jpg differ
diff --git a/Validation/Mixing Layer/CFD-BOS Qualitative.png b/Validation/Mixing Layer/CFD-BOS Qualitative.png
new file mode 100644
index 0000000..7d103fc
Binary files /dev/null and b/Validation/Mixing Layer/CFD-BOS Qualitative.png differ
diff --git a/Validation/Mixing Layer/bos_image_1000_per.mat b/Validation/Mixing Layer/bos_image_1000_per.mat
new file mode 100644
index 0000000..7f2f7a1
Binary files /dev/null and b/Validation/Mixing Layer/bos_image_1000_per.mat differ
diff --git a/Validation/Mixing Layer/dir_estimation.m b/Validation/Mixing Layer/dir_estimation.m
new file mode 100644
index 0000000..969178c
--- /dev/null
+++ b/Validation/Mixing Layer/dir_estimation.m	
@@ -0,0 +1,129 @@
+%% Clean up
+%clear
+close all
+clc
+
+%% Load data
+ref = load("mixing_layer_ref_1000_per.mat");
+grad = load("bos_image_1000_per.mat");
+CFD_data = load("rho_mixing_001_centered_z_adjusted_dimensioned.mat");
+%% Estimate displacement using DIR
+tic
+%[DIR_results,~] = imregdeform(ref.image,grad.bos_image,'NumPyramidLevels',1,'GridSpacing',[8 8],'GridRegularization',0,'DisplayProgress',0); %Adjust parameters of num pyramid and spacing as you see fit
+%DIR_U = DIR_results(:,:,1)*20e-6; %20e-6 is pixel pitch of camera in m
+%DIR_V = DIR_results(:,:,2)*20e-6; %20e-6 is pixel pitch of camera in m
+fprintf("DIR estimation complete.")
+toc
+
+%Normalize gradient based off physical spacing of dimensions
+dx = abs(CFD_data.x(1)-CFD_data.x(2))/1000; %m
+dy = abs(CFD_data.y(1)-CFD_data.y(2))/1000; %m
+
+
+%Calculate estimated density gradient in y
+n0 = 1.0003; %Refractive index of air
+del_z = abs(max(CFD_data.z)-min(CFD_data.z))*1e-3; %Size of density object in m
+Z_D = (50e-3+del_z/2); %Standoff distance (distance from bos pattern to density grad) + 1/2 width of density object
+mag = 1000/1524; %Magnification of the optical system, f1/f2 for telecentric
+K = 0.000226; %Gladstone-Dale constant in kg/m^3
+drho_dy_BOS = n0*DIR_V/(Z_D*mag*del_z*K); %estimated density gradient (kg/m^4)
+
+%% Convert to physical coordinates 
+%Clean up extra space from raytracing from mismatch of sensor size to density
+%gradient size
+density_mag_x_size = ceil(max(CFD_data.x)*1000*mag/20); %Number of pixels in half sensor plane needed to show gradient
+density_mag_y_size = ceil(max(CFD_data.y)*1000*mag/20);
+excess_pixels_x = 8192/2-density_mag_x_size;
+excess_pixels_y = 1024/2-density_mag_y_size; %Negative value indicates small cutoff on density grad in y direction
+
+%Therefore, cutoff excess pixels in ray traced image in x dir, and cutoff
+%excess data in y dir for CFD data for 1 to 1 comparison.
+
+drho_dy_BOS_post = drho_dy_BOS(:,excess_pixels_x+1:end-excess_pixels_x);
+x_y_pixel_ratio = size(drho_dy_BOS_post,1)/size(drho_dy_BOS_post,2);
+excess_CFD_data_y = floor(size(CFD_data.rho,2)-size(CFD_data.rho,1)*x_y_pixel_ratio);
+%% Prepare CFD and Raytraced data for comparison
+[drdy_CFD_raw,drdx_CFD_raw,~] = gradient(CFD_data.rho); %dy,dx flipped since rho is originally 1200x200 and should be 200x1200.
+drdx_CFD_raw = mean(drdx_CFD_raw,3)/dx;
+drdy_CFD_raw = mean(drdy_CFD_raw,3)/dy;
+
+%Interpolate CFD to higher resolution grid
+x_original = linspace(0, 1, 1200);
+y_original = linspace(0, 1, 200);
+x_new = linspace(0, 1, size(drho_dy_BOS_post,2));
+y_new = linspace(0, 1, 1024+abs(excess_pixels_y*2));
+drdy_CFD_interp =  interp2(x_original, y_original, drdy_CFD_raw', x_new', y_new, 'linear');
+drdy_CFD_interp_post = drdy_CFD_interp(abs(excess_pixels_y)+1:end-abs(excess_pixels_y),:);
+
+%% Crop out initial region of mixing layer due to bad DIR estimation 
+drdy_CFD_cropped = drdy_CFD_interp_post(:,1:7200);
+drho_dy_BOS_cropped = drho_dy_BOS_post(:,1:7200);
+
+figure;
+subplot(2,1,1)
+%imshow(drdy_CFD_post',[])
+imshow(drdy_CFD_cropped,[])
+clim(max(abs(drho_dy_BOS_cropped),[],'all')*[-1,1]);
+title("CFD ground truth")
+ylabel("d\rho/dy")
+colorbar
+subplot(2,1,2)
+imshow(drho_dy_BOS_cropped,[])
+clim(max(abs(drho_dy_BOS_cropped),[],'all')*[-1,1]);
+colorbar;
+title("Estimated density gradient from raytracing")
+ylabel("d\rho/dy")
+
+%% Compare qualitatively and quantitatively
+
+%Qualitative reconstructions of density gradient
+figure;
+set(gcf,'Position',[1000         436        1066         442])
+subplot(2,1,1)
+imshow(drdy_CFD_interp_post,[])
+clim(max(abs(drho_dy_BOS_post),[],'all')*[-1,1]);
+title("CFD",'FontWeight','normal')
+ylabel("\partial\rho/\partialy (kg/m^{4})")
+colorbar
+set(gca,'FontSize',16,'TickLabelInterpreter','Latex')
+%set(gcf,"Position",[69         302        1752         308])
+axis equal
+
+%figure;
+subplot(2,1,2)
+imshow(drho_dy_BOS_post,[])
+clim(max(abs(drho_dy_BOS_post),[],'all')*[-1,1]);
+colorbar;
+title("BOS",'FontWeight','normal')
+ylabel("\partial\rho/\partialy (kg/m^{4})")
+set(gca,'FontSize',16,'TickLabelInterpreter','Latex')
+%set(gcf,"Position",[69         302        1752         308])
+axis equal
+
+%Calculate statistics between two datasets
+
+a = drdy_CFD_cropped(:)\drho_dy_BOS_cropped(:);
+y_fit = drdy_CFD_cropped(:)*a;
+
+% Compute R^2
+SS_res = sum((drho_dy_BOS_cropped(:) - y_fit).^2);      
+SS_tot = sum((drho_dy_BOS_cropped(:) - mean(drho_dy_BOS_cropped(:))).^2);
+accuracy = 1-abs(1-a)
+consistency = 1 - SS_res / SS_tot
+
+%Quantitative analysis of density gradient*
+figure;
+plot(drdy_CFD_cropped(:),drho_dy_BOS_cropped(:),'*');
+hold on
+XL = get(gca,'XLim');
+plot(XL,a*XL,'-','LineWidth',2)
+plot(XL,XL,'--','LineWidth',2)
+legend('CFD-BOS comparison','Best fit line','Perfect match line','Location','northwest')
+axis equal
+xlim([-300 100])
+ylim([-300 100])
+xlabel("\partial\rho/\partialy (kg/m^{4})")
+ylabel("\partial\rho/\partialy (kg/m^{4})")
+grid on
+set(gca,'FontSize',16,'TickLabelInterpreter','Latex')
+set(gcf,'Position',[820   232   920   641])