diff --git a/Matlab App v5/BOS Pattern Description.png b/Matlab App v5/BOS Pattern Description.png
new file mode 100644
index 0000000..c36e1c5
Binary files /dev/null and b/Matlab App v5/BOS Pattern Description.png differ
diff --git a/Matlab App v5/Error_lightfield_dialog.mlapp b/Matlab App v5/Error_lightfield_dialog.mlapp
new file mode 100644
index 0000000..34be2f8
Binary files /dev/null and b/Matlab App v5/Error_lightfield_dialog.mlapp differ
diff --git a/Matlab App v5/Figures.pptx b/Matlab App v5/Figures.pptx
new file mode 100644
index 0000000..a2adb0a
Binary files /dev/null and b/Matlab App v5/Figures.pptx differ
diff --git a/Matlab App v5/Location_dialog.mlapp b/Matlab App v5/Location_dialog.mlapp
new file mode 100644
index 0000000..95e3c08
Binary files /dev/null and b/Matlab App v5/Location_dialog.mlapp differ
diff --git a/Matlab App v5/MIRAGE.mlapp b/Matlab App v5/MIRAGE.mlapp
new file mode 100644
index 0000000..2dc16b9
Binary files /dev/null and b/Matlab App v5/MIRAGE.mlapp differ
diff --git a/Matlab App v5/Warning_Dialog.mlapp b/Matlab App v5/Warning_Dialog.mlapp
new file mode 100644
index 0000000..5022306
Binary files /dev/null and b/Matlab App v5/Warning_Dialog.mlapp differ
diff --git a/Matlab App v5/function_circular_aperture_cutoff_heisenberg_diffraction.m b/Matlab App v5/function_circular_aperture_cutoff_heisenberg_diffraction.m
new file mode 100644
index 0000000..4f3c018
--- /dev/null
+++ b/Matlab App v5/function_circular_aperture_cutoff_heisenberg_diffraction.m	
@@ -0,0 +1,89 @@
+function [lightfield,del_theta_x,del_theta_y] = function_circular_aperture_cutoff(lightfield, circular_aperture_diameter, wavelength, diffraction_type, low_lim,upp_lim)
+% Function to apply a circular aperture cutoff to a lightfield
+% INPUTS:
+%   - lightfield: Structure containing ray positions and directions
+%   - circular_aperture_diameter: Diameter of the circular aperture
+% OUTPUT:
+%   - new_lightfield: Updated lightfield structure with rays outside the aperture set to NaN
+
+% Start timer for performance measurement
+tic
+
+% Create a logical mask that identifies rays outside the aperture radius
+% Radius is half of the aperture diameter
+mask = sqrt(lightfield.position.x.^2 + lightfield.position.y.^2) > (circular_aperture_diameter / 2);
+
+% % Copy the original lightfield to preserve unchanged rays
+% lightfield = lightfield;
+
+% Apply the mask to all position and direction components
+% Rays outside the aperture are set to NaN in the new lightfield structure
+lightfield.position.x(mask) = nan;
+lightfield.position.y(mask) = nan;
+lightfield.position.z(mask) = nan;
+lightfield.direction.dx(mask) = nan;
+lightfield.direction.dy(mask) = nan;
+lightfield.direction.dz(mask) = nan;
+
+lightfield.position.x = lightfield.position.x(~isnan(lightfield.position.x));
+lightfield.position.y = lightfield.position.y(~isnan(lightfield.position.y));
+lightfield.position.z = lightfield.position.z(~isnan(lightfield.position.z));
+lightfield.direction.dx = lightfield.direction.dx(~isnan(lightfield.direction.dx));
+lightfield.direction.dy = lightfield.direction.dy(~isnan(lightfield.direction.dy));
+lightfield.direction.dz = lightfield.direction.dz(~isnan(lightfield.direction.dz));
+
+if diffraction_type ~= "none"
+    %Calculate distance to each aperture edge
+    aperture_radius = circular_aperture_diameter/2;
+    x_distance = sqrt(aperture_radius^2-lightfield.position.y.^2)-abs(lightfield.position.x); %dimensionalized distance to edge
+    y_distance = sqrt(aperture_radius^2-lightfield.position.x.^2)-abs(lightfield.position.y); %dimensionalized distance to edge
+    del_theta_x = zeros(size(lightfield.position.x)); %Initialize vector for change in theta_x due to diffraction
+    del_theta_y = zeros(size(lightfield.position.y)); %Initialize vector for change in theta_y due to diffraction
+    k = 2*pi/wavelength;
+
+    %Calculate norm distribution characteristics
+    sigma_x = atan(1./(2*x_distance*k));
+    sigma_y = atan(1./(2*y_distance*k));
+
+    %Diffraction model parameters
+    if string(diffraction_type) == 'dynamic limited'
+        %Calculate angle of defelction due to diffraction
+        lower_limit_x = -sigma_x;
+        upper_limit_x = -lower_limit_x;
+        lower_limit_y = -sigma_y;
+        upper_limit_y = -lower_limit_y;
+        %Modeling diffraction in the Monte-Carlo ray-trace environment,
+        %Coffey 1998
+    elseif string(diffraction_type) == "-pi/2 to pi/2 limited"
+        lower_limit_x = -pi/2*ones(size(sigma_x));
+        upper_limit_x = -lower_limit_x;
+        lower_limit_y = -pi/2*ones(size(sigma_y));
+        upper_limit_y = -lower_limit_y;
+    else
+        lower_limit_x = low_lim*ones(size(sigma_x));
+        upper_limit_x = upp_lim*ones(size(sigma_x));
+        lower_limit_y = lower_limit_x;
+        upper_limit_y = upper_limit_x;
+    end %End of type of diffraction
+
+    %Calculate angle due to diffraction
+    parfor counter = 1:numel(lightfield.position.x) % For all rays
+        %Generate values within the truncated range
+        del_theta_x(counter) = truncate_gaussian(0, sigma_x(counter), lower_limit_x(counter), upper_limit_x(counter));
+        del_theta_y(counter) = truncate_gaussian(0, sigma_y(counter), lower_limit_y(counter), upper_limit_y(counter));
+    end
+    %Edge diffraction in Monte Carlo ray tracing, Freniere et al 1999
+    % DOI:10.1117/12.363773
+    
+
+
+    %Update new directions in the x and y component by converting from
+    %spherical to cartesian
+
+    lightfield.direction.dx = lightfield.direction.dx + tan(del_theta_x);
+    lightfield.direction.dy = lightfield.direction.dy + tan(del_theta_y);
+    %Display completion message and stop timer
+    fprintf("Lightfield cutoff.");
+end %End of diffraction check
+toc
+end
\ No newline at end of file
diff --git a/Matlab App v5/function_convert_rho_field_to_n_field.m b/Matlab App v5/function_convert_rho_field_to_n_field.m
new file mode 100644
index 0000000..1725df2
--- /dev/null
+++ b/Matlab App v5/function_convert_rho_field_to_n_field.m	
@@ -0,0 +1,11 @@
+function n_field = function_convert_rho_field_to_n_field(rho,gladstone_dale)
+%Convert rho to n field
+
+%% Load data
+
+load("sample-density.mat");
+rho = rho(:,:,1);
+
+n = gladstone_dale*rho+1;
+end
+
diff --git a/Matlab App v5/function_generate_image_2.m b/Matlab App v5/function_generate_image_2.m
new file mode 100644
index 0000000..1be745f
--- /dev/null
+++ b/Matlab App v5/function_generate_image_2.m	
@@ -0,0 +1,72 @@
+function [rendered_image] = function_generate_image_2(lightfield, x_pixel_number, y_pixel_number, pixel_pitch)
+
+% Remove invalid rays
+invalid = isnan(lightfield.position.x) | isnan(lightfield.position.y);
+lightfield.direction.dx(invalid) = [];
+lightfield.direction.dy(invalid) = [];
+lightfield.position.x(invalid) = [];
+lightfield.position.y(invalid) = [];
+
+tic
+rendered_image = zeros(y_pixel_number, x_pixel_number); %rows by columns
+
+% Compute subpixel locations relative to center of sensor
+x_center = x_pixel_number / 2;
+y_center = y_pixel_number / 2;
+
+% Convert positions to pixel space
+x_pixel_f = (lightfield.position.x / pixel_pitch) + x_center + 1;
+y_pixel_f = (lightfield.position.y / pixel_pitch) + y_center + 1;
+
+% Compute
+alpha = atan(sqrt((lightfield.direction.dx./lightfield.direction.dz).^2+(lightfield.direction.dy./lightfield.direction.dz).^2));
+intensity = cos(alpha).^4;
+
+% Get surrounding pixel indices
+x0 = floor(x_pixel_f);
+y0 = floor(y_pixel_f);
+x1 = x0 + 1;
+y1 = y0 + 1;
+
+% Calculate weights for bilinear interpolation
+wx = x_pixel_f - x0;
+wy = y_pixel_f - y0;
+w00 = (1 - wx) .* (1 - wy); % top-left
+w10 = wx .* (1 - wy);       % top-right
+w01 = (1 - wx) .* wy;       % bottom-left
+w11 = wx .* wy;             % bottom-right
+
+% Filter valid rays that hit the sensor
+valid = x0 >= 1 & x1 <= x_pixel_number & y0 >= 1 & y1 <= y_pixel_number;
+x0 = x0(valid); x1 = x1(valid);
+y0 = y0(valid); y1 = y1(valid);
+w00 = w00(valid); w10 = w10(valid); w01 = w01(valid); w11 = w11(valid);
+intensity = intensity(valid);
+
+%Reshape vector
+pixel_indices = [
+    sub2ind([y_pixel_number, x_pixel_number], y0, x0),  % top-left
+    sub2ind([y_pixel_number, x_pixel_number], y0, x1),  % top-right
+    sub2ind([y_pixel_number, x_pixel_number], y1, x0),  % bottom-left
+    sub2ind([y_pixel_number, x_pixel_number], y1, x1)   % bottom-right
+];
+weights = [
+    intensity .* w00;
+    intensity .* w10;
+    intensity .* w01;
+    intensity .* w11
+];
+
+% Accumulate all contributions using accumarray
+accumulated = accumarray(pixel_indices(:), weights(:), [y_pixel_number * x_pixel_number, 1]);
+rendered_image = reshape(accumulated, [y_pixel_number, x_pixel_number]);
+
+% Flip vertically to match image coordinate convention
+rendered_image = rendered_image(end:-1:1, end:-1:1);
+
+
+fprintf("Image rendered with bilinear interpolation. ");
+toc
+
+figure; imshow(rendered_image(2:end-1,2:end-1), []);
+end
diff --git a/Matlab App v5/function_generate_image_with_splat.m b/Matlab App v5/function_generate_image_with_splat.m
new file mode 100644
index 0000000..4e4e765
--- /dev/null
+++ b/Matlab App v5/function_generate_image_with_splat.m	
@@ -0,0 +1,59 @@
+function [] = function_generate_image_with_splat(lightfield, x_pixel_number, y_pixel_number, pixel_pitch)
+
+% Remove invalid rays
+valid_mask = ~isnan(lightfield.position.x) & ~isnan(lightfield.position.y);
+x = lightfield.position.x(valid_mask);
+y = lightfield.position.y(valid_mask);
+dx = lightfield.direction.dx(valid_mask);
+dy = lightfield.direction.dy(valid_mask);
+
+% Floating-point pixel positions
+x_fp = (x / pixel_pitch) + (x_pixel_number / 2);
+y_fp = (y / pixel_pitch) + (y_pixel_number / 2);
+
+% Compute intensity from direction magnitude
+intensity = max(0, min(1, 1 ./ (1 + sqrt(dx.^2 + dy.^2))));
+
+% Bilinear splatting
+x0 = floor(x_fp);
+y0 = floor(y_fp);
+dx_frac = x_fp - x0;
+dy_frac = y_fp - y0;
+
+% Clip to valid ranges (x0+1 and y0+1 must also be within bounds)
+valid = x0 >= 1 & x0 < x_pixel_number & y0 >= 1 & y0 < y_pixel_number;
+
+x0 = x0(valid);
+y0 = y0(valid);
+dx_frac = dx_frac(valid);
+dy_frac = dy_frac(valid);
+intensity = intensity(valid);
+
+rendered_image = zeros(y_pixel_number, x_pixel_number);
+
+% Accumulate contributions into 4 neighboring pixels
+for i = 1:length(x0)
+    x1 = x0(i) + 1;
+    y1 = y0(i) + 1;
+    
+    % Contributions to four neighbors
+    w00 = (1 - dx_frac(i)) * (1 - dy_frac(i));
+    w01 = (1 - dx_frac(i)) * dy_frac(i);
+    w10 = dx_frac(i) * (1 - dy_frac(i));
+    w11 = dx_frac(i) * dy_frac(i);
+
+    rendered_image(y0(i),  x0(i))  = rendered_image(y0(i),  x0(i))  + intensity(i) * w00;
+    rendered_image(y1,     x0(i))  = rendered_image(y1,     x0(i))  + intensity(i) * w01;
+    rendered_image(y0(i),  x1)     = rendered_image(y0(i),  x1)     + intensity(i) * w10;
+    rendered_image(y1,     x1)     = rendered_image(y1,     x1)     + intensity(i) * w11;
+end
+
+% Flip image vertically and horizontally
+rendered_image = flipud(fliplr(rendered_image));
+
+fprintf("Image rendered. ");
+toc
+
+figure; imshow(rendered_image, []);
+
+end
\ No newline at end of file
diff --git a/Matlab App v5/function_generate_lightfield.asv b/Matlab App v5/function_generate_lightfield.asv
new file mode 100644
index 0000000..b9fc20e
--- /dev/null
+++ b/Matlab App v5/function_generate_lightfield.asv	
@@ -0,0 +1,100 @@
+function lightfield = function_generate_lightfield(number_of_light_rays_per_location, scale, dot_width, dot_height, spacing_x, spacing_y, num_column, num_row, z_location, cone_half_angle,use_gpu)
+
+num_location_per_dot_x = round(dot_width/scale);
+lightsource.location.x = zeros(num_column*num_location_per_dot_x);
+x_values = (scale * (1:num_location_per_dot_x)') + (dot_width + spacing_x) * (0:num_column-1);
+lightsource.location.x = reshape(x_values, [], 1);
+lightsource.location.x = lightsource.location.x-(lightsource.location.x(end)-lightsource.location.x(1))/2; %center around (0,y)
+
+num_location_per_dot_y = round(dot_height/scale);
+lightsource.location.y = zeros(num_row*num_location_per_dot_y);
+y_values = (scale * (1:num_location_per_dot_y)') + (dot_height + spacing_y) * (0:num_row-1);
+lightsource.location.y = reshape(y_values, [], 1);
+lightsource.location.y = lightsource.location.y-(lightsource.location.y(end)-lightsource.location.y(1))/2; %center around (x,0)
+
+%Used to address issue of scale being approximately equal to pixel pitch.
+%Perturb lightsource locations so gridded image doesn't appear.
+
+% lightsource.location.x = lightsource.location.x + 1*scale * (rand(size(lightsource.location.x)) - 0.5);
+% lightsource.location.y = lightsource.location.y + 1*scale * (rand(size(lightsource.location.y)) - 0.5);
+% lightsource.location.x = unifrnd(-1,1,size(lightsource.location.x))*scale+lightsource.location.x;
+% lightsource.location.y = unifrnd(-1,1,size(lightsource.location.y))*scale+lightsource.location.y;
+
+lightsource.location.z = z_location; %Initialize position
+%% Generate uniform lightfield
+
+% Collimated light should have very small divergence/diffusion angle
+cone_half_angle_radians = deg2rad(cone_half_angle);
+tic
+[temp_x_mesh,temp_y_mesh] = meshgrid(lightsource.location.x,lightsource.location.y);
+temp_x_flat = temp_x_mesh(:);
+mu = number_of_light_rays_per_location;
+for counter = 1:size(temp_x_flat,1)
+    num_lightray_array(counter) = min(max(round(normrnd(mu,(2*mu-1)/6)),1),2*mu);
+end
+total_lightray = sum(num_lightray_array);
+
+for i = 1:M
+    for j = 1:N
+        % Calculate the number of rays for the current position
+        num_rays = min(max(round(normrnd(mu, (2*mu - 1) / 6)), 1), 2 * mu);
+
+        % Assign the x, y positions to the arrays
+        x_positions = temp_x_mesh(i,j) * ones(1, num_rays);
+        y_positions = temp_y_mesh(i,j) * ones(1, num_rays);
+        
+        % Store the x and y positions in the preallocated arrays
+        lightfield_position_x(current_index : current_index + num_rays - 1) = x_positions;
+        lightfield_position_y(current_index : current_index + num_rays - 1) = y_positions;
+
+        % Update the current index
+        current_index = current_index + num_rays;
+    end
+end
+
+
+% for counter = 1:number_of_light_rays_per_location
+%     lightfield.position.x(:,:,counter) = temp_x_mesh;
+%     lightfield.position.y(:,:,counter) = temp_y_mesh;
+% end
+
+
+
+
+
+
+
+clearvars temp_x_mesh temp_y_mesh;
+lightfield.position.x = lightfield.position.x(:);
+lightfield.position.y = lightfield.position.y(:);
+lightfield.position.z = lightsource.location.z*ones(size(lightfield.position.y));
+
+lightfield.direction.random_cone_angle = (rand(size(lightfield.position.x))-1/2)*cone_half_angle_radians*2;
+lightfield.direction.random_circular_angle = (rand(size(lightfield.position.x))-1/2)*pi*2;
+lightfield.direction.dx = sin(lightfield.direction.random_cone_angle).*cos(lightfield.direction.random_circular_angle);
+lightfield.direction.dy = sin(lightfield.direction.random_cone_angle).*sin(lightfield.direction.random_circular_angle);
+lightfield.direction.dz = (-1)*ones(size(lightfield.direction.dy));
+
+clear lightfield.direction.random_cone_angle lightfield.direction.random_circular_angle
+
+if use_gpu
+    lightfield.position.x = gpuArray(lightfield.position.x);
+    lightfield.position.y = gpuArray(lightfield.position.y);
+    lightfield.position.z = gpuArray(lightfield.position.z);
+    lightfield.direction.dx = gpuArray(lightfield.direction.dx);
+    lightfield.direction.dy = gpuArray(lightfield.direction.dy);
+    lightfield.direction.dz = gpuArray(lightfield.direction.dz);
+end
+
+
+
+% if show_ray_tracing_results
+    % figure;
+    % plot3(lightfield.position.x(1:37:end),lightfield.position.y(1:37:end),lightfield.position.z(1:37:end),'.');
+    % % axis equal;
+    % grid on;
+    % hold on;
+% end
+
+fprintf("Lightfield generated. Number of light rays = %d. ", size(lightfield.position.x,1));
+toc
\ No newline at end of file
diff --git a/Matlab App v5/function_generate_lightfield.m b/Matlab App v5/function_generate_lightfield.m
new file mode 100644
index 0000000..d17f227
--- /dev/null
+++ b/Matlab App v5/function_generate_lightfield.m	
@@ -0,0 +1,101 @@
+function lightfield = function_generate_lightfield(number_of_light_rays_per_location, scale, dot_width, dot_height, spacing_x, spacing_y, num_column, num_row, z_location, cone_half_angle,use_gpu)
+
+num_location_per_dot_x = round(dot_width/scale);
+lightsource.location.x = zeros(num_column*num_location_per_dot_x);
+x_values = (scale * (1:num_location_per_dot_x)') + (dot_width + spacing_x) * (0:num_column-1);
+lightsource.location.x = reshape(x_values, [], 1);
+lightsource.location.x = lightsource.location.x-(lightsource.location.x(end)-lightsource.location.x(1))/2; %center around (0,y)
+
+num_location_per_dot_y = round(dot_height/scale);
+lightsource.location.y = zeros(num_row*num_location_per_dot_y);
+y_values = (scale * (1:num_location_per_dot_y)') + (dot_height + spacing_y) * (0:num_row-1);
+lightsource.location.y = reshape(y_values, [], 1);
+lightsource.location.y = lightsource.location.y-(lightsource.location.y(end)-lightsource.location.y(1))/2; %center around (x,0)
+
+lightsource.location.z = z_location; %Initialize position
+%% Generate uniform lightfield
+
+% Collimated light should have very small divergence/diffusion angle
+cone_half_angle_radians = deg2rad(cone_half_angle);
+tic
+[temp_x_mesh,temp_y_mesh] = meshgrid(lightsource.location.x,lightsource.location.y);
+
+% % Flatten the mesh grids
+% temp_x_flat = temp_x_mesh(:);
+% temp_y_flat = temp_y_mesh(:);
+% 
+% % Calculate the number of rays per position using normal distribution
+% mu = number_of_light_rays_per_location;
+% num_lightray_array = zeros(size(temp_x_flat));  % Preallocate
+
+% % Generate rays per position using normal distribution
+% for counter = 1:size(temp_x_flat, 1)
+%     num_lightray_array(counter) = min(max(round(normrnd(mu, (2 * mu - 1) / 6)), 1), 2 * mu);
+% end
+
+% % Calculate the total number of light rays
+% total_lightray = sum(num_lightray_array);
+% 
+% % Preallocate arrays to store light ray positions
+% lightfield.position.x = zeros(1, total_lightray);
+% lightfield.position.y = zeros(1, total_lightray);
+% 
+% % Initialize current index
+% current_index = 1;
+% 
+% % Loop through each (x, y) position and assign the light ray positions
+% for counter = 1:size(temp_x_flat, 1)
+%     num_rays = num_lightray_array(counter);
+% 
+%     % Assign x and y positions for the current number of rays
+%     x_positions = temp_x_flat(counter) * ones(1, num_rays);
+%     y_positions = temp_y_flat(counter) * ones(1, num_rays);
+% 
+%     % Store the positions in the preallocated arrays
+%     lightfield.position.x(current_index : current_index + num_rays - 1) = x_positions;
+%     lightfield.position.y(current_index : current_index + num_rays - 1) = y_positions;
+% 
+%     % Update the current index
+%     current_index = current_index + num_rays;
+% end
+
+
+for counter = 1:number_of_light_rays_per_location
+    lightfield.position.x(:,:,counter) = temp_x_mesh;
+    lightfield.position.y(:,:,counter) = temp_y_mesh;
+end
+
+clearvars temp_x_mesh temp_y_mesh;
+lightfield.position.x = lightfield.position.x(:);
+lightfield.position.y = lightfield.position.y(:);
+lightfield.position.z = lightsource.location.z*ones(size(lightfield.position.y));
+
+lightfield.direction.random_cone_angle = (rand(size(lightfield.position.x))-1/2)*cone_half_angle_radians*2;
+lightfield.direction.random_circular_angle = (rand(size(lightfield.position.x))-1/2)*pi*2;
+lightfield.direction.dx = sin(lightfield.direction.random_cone_angle).*cos(lightfield.direction.random_circular_angle);
+lightfield.direction.dy = sin(lightfield.direction.random_cone_angle).*sin(lightfield.direction.random_circular_angle);
+lightfield.direction.dz = (-1)*ones(size(lightfield.direction.dy));
+
+clear lightfield.direction.random_cone_angle lightfield.direction.random_circular_angle
+
+if use_gpu
+    lightfield.position.x = gpuArray(lightfield.position.x);
+    lightfield.position.y = gpuArray(lightfield.position.y);
+    lightfield.position.z = gpuArray(lightfield.position.z);
+    lightfield.direction.dx = gpuArray(lightfield.direction.dx);
+    lightfield.direction.dy = gpuArray(lightfield.direction.dy);
+    lightfield.direction.dz = gpuArray(lightfield.direction.dz);
+end
+
+
+
+% if show_ray_tracing_results
+    % figure;
+    % plot3(app.lightfield.position.x(1:37:end),app.lightfield.position.y(1:37:end),app.lightfield.position.z(1:37:end),'.');
+    % % axis equal;
+    % grid on;
+    % hold on;
+% end
+
+fprintf("Lightfield generated. Number of light rays = %d. ", size(lightfield.position.x,1));
+toc
\ No newline at end of file
diff --git a/Matlab App v5/function_generate_lightfield_sobol.m b/Matlab App v5/function_generate_lightfield_sobol.m
new file mode 100644
index 0000000..370c2e0
--- /dev/null
+++ b/Matlab App v5/function_generate_lightfield_sobol.m	
@@ -0,0 +1,65 @@
+function lightfield = function_generate_lightfield_sobol(num_rays_per_dot, dot_width, dot_height, spacing_x, spacing_y, num_column, num_row, z_location, cone_half_angle, use_gpu)
+% FUNCTION_GENERATE_LIGHTFIELD_V2 generates a lightfield using Sobol sequence sampling
+% with evenly spaced dot layout and a cone angle for angular divergence.
+
+tic
+
+% -------------------------------------------------------------------------
+% Dot Grid Definition
+dot_centers_x = (0:num_column-1) * (dot_width + spacing_x);
+dot_centers_y = (0:num_row-1) * (dot_height + spacing_y);
+[dot_grid_x, dot_grid_y] = meshgrid(dot_centers_x, dot_centers_y);
+
+dot_centers = [dot_grid_x(:), dot_grid_y(:)];
+num_dots = size(dot_centers, 1);
+total_rays = num_dots * num_rays_per_dot;
+
+% Center the dot grid around (0,0)
+dot_centers(:,1) = dot_centers(:,1) - mean(dot_centers(:,1));
+dot_centers(:,2) = dot_centers(:,2) - mean(dot_centers(:,2));
+
+% -------------------------------------------------------------------------
+% Generate Sobol sequence (vectorized)
+sobol = sobolset(2, 'Skip', 1e3, 'Leap', 1e2);
+sobol = scramble(sobol, 'MatousekAffineOwen');
+samples = net(sobol, total_rays);  % Generate all at once
+
+% Map to local dot rectangles
+x_local = samples(:,1) * dot_width  - dot_width/2;
+y_local = samples(:,2) * dot_height - dot_height/2;
+
+% Assign dot centers (repeat each one num_rays_per_dot times)
+dot_indices = repelem(1:num_dots, num_rays_per_dot)';
+dot_x = dot_centers(dot_indices, 1);
+dot_y = dot_centers(dot_indices, 2);
+
+% Final global positions
+lightfield.position.x = dot_x + x_local;
+lightfield.position.y = dot_y + y_local;
+lightfield.position.z = z_location * ones(total_rays, 1);
+
+% -------------------------------------------------------------------------
+% Generate Diverging Directions (Cone Distribution)
+cone_half_angle_radians = deg2rad(cone_half_angle);
+cos_theta = (1 - cos(cone_half_angle_radians)) * rand(total_rays,1) + cos(cone_half_angle_radians);
+theta = acos(cos_theta);
+phi = 2 * pi * rand(total_rays,1);
+
+lightfield.direction.dx = sin(theta) .* cos(phi);
+lightfield.direction.dy = sin(theta) .* sin(phi);
+lightfield.direction.dz = -cos(theta); 
+
+% -------------------------------------------------------------------------
+% Move to GPU (if requested)
+if use_gpu
+    lightfield.position.x = gpuArray(lightfield.position.x);
+    lightfield.position.y = gpuArray(lightfield.position.y);
+    lightfield.position.z = gpuArray(lightfield.position.z);
+    lightfield.direction.dx = gpuArray(lightfield.direction.dx);
+    lightfield.direction.dy = gpuArray(lightfield.direction.dy);
+    lightfield.direction.dz = gpuArray(lightfield.direction.dz);
+end
+
+fprintf("Lightfield generated using Sobol. Number of light rays = %d. ", total_rays);
+toc
+end
diff --git a/Matlab App v5/function_generate_lightfield_v2.asv b/Matlab App v5/function_generate_lightfield_v2.asv
new file mode 100644
index 0000000..8f3cd64
--- /dev/null
+++ b/Matlab App v5/function_generate_lightfield_v2.asv	
@@ -0,0 +1,68 @@
+function lightfield = function_generate_lightfield_v2(num_rays_per_dot, dot_width, dot_height, spacing_x, spacing_y, num_column, num_row, z_location, cone_half_angle, use_gpu)
+% FUNCTION_GENERATE_LIGHTFIELD_V2 generates a lightfield using Latin Hypercube Sampling
+% with evenly spaced dot layout and a cone angle for angular divergence.
+
+% -------------------------------------------------------------------------
+% Dot Grid Definition
+tic
+dot_centers_x = (0:num_column-1) * (dot_width + spacing_x);
+dot_centers_y = (0:num_row-1) * (dot_height + spacing_y);
+[dot_grid_x, dot_grid_y] = meshgrid(dot_centers_x, dot_centers_y);
+
+dot_centers = [dot_grid_x(:), dot_grid_y(:)];
+num_dots = size(dot_centers, 1);
+
+% Center the dot grid around (0,0)
+dot_centers(:,1) = dot_centers(:,1) - mean(dot_centers(:,1));
+dot_centers(:,2) = dot_centers(:,2) - mean(dot_centers(:,2));
+
+% -------------------------------------------------------------------------
+% Preallocate and Sample Positions Using LHS
+total_rays = num_dots * num_rays_per_dot;
+lightfield.position.x = zeros(total_rays, 1);
+lightfield.position.y = zeros(total_rays, 1);
+lightfield.position.z = z_location * ones(total_rays, 1);
+
+x_cell = cell(num_dots, 1);
+y_cell = cell(num_dots, 1);
+
+parfor k = 1:num_dots
+    lhs_samples = lhsdesign(num_rays_per_dot, 2);  % LHS in [0,1]^2
+    x_local = lhs_samples(:,1) * dot_width  - dot_width/2;
+    y_local = lhs_samples(:,2) * dot_height - dot_height/2;
+
+    x_cell{k} = dot_centers(k,1) + x_local;
+    y_cell{k} = dot_centers(k,2) + y_local;
+    fprintf('%')
+end
+
+% Concatenate results after parallel loop
+lightfield.position.x = vertcat(x_cell{:});
+lightfield.position.y = vertcat(y_cell{:});
+lightfield.position.z = z_location * ones(num_dots * num_rays_per_dot, 1);
+
+% -------------------------------------------------------------------------
+% Generate Diverging Directions (Cone Distribution)
+cone_half_angle_radians = deg2rad(cone_half_angle);
+
+random_cone_angle = (rand(total_rays,1) - 0.5) * 2 * cone_half_angle_radians;
+random_circular_angle = (rand(total_rays,1) - 0.5) * 2 * pi;
+
+lightfield.direction.dx = sin(random_cone_angle) .* cos(random_circular_angle);
+lightfield.direction.dy = sin(random_cone_angle) .* sin(random_circular_angle);
+lightfield.direction.dz = -cos(random_cone_angle);  % Or simply -1 for narrow cones
+
+% -------------------------------------------------------------------------
+% Move to GPU (if requested)
+if use_gpu
+    lightfield.position.x = gpuArray(lightfield.position.x);
+    lightfield.position.y = gpuArray(lightfield.position.y);
+    lightfield.position.z = gpuArray(lightfield.position.z);
+    lightfield.direction.dx = gpuArray(lightfield.direction.dx);
+    lightfield.direction.dy = gpuArray(lightfield.direction.dy);
+    lightfield.direction.dz = gpuArray(lightfield.direction.dz);
+end
+
+fprintf("Lightfield generated. Number of light rays = %d. ", total_rays);
+toc
+end
diff --git a/Matlab App v5/function_generate_lightfield_v2.m b/Matlab App v5/function_generate_lightfield_v2.m
new file mode 100644
index 0000000..93060e5
--- /dev/null
+++ b/Matlab App v5/function_generate_lightfield_v2.m	
@@ -0,0 +1,53 @@
+function lightfield = function_generate_lightfield_v2(num_rays_per_dot, dot_width, dot_height, spacing_x, spacing_y, num_column, num_row, z_location, cone_half_angle)
+
+% -------------------------------------------------------------------------
+% Dot Grid Definition
+tic
+dot_centers_x = (0:num_column-1) * (dot_width + spacing_x);
+dot_centers_y = (0:num_row-1) * (dot_height + spacing_y);
+[dot_grid_x, dot_grid_y] = meshgrid(dot_centers_x, dot_centers_y);
+
+dot_centers = [dot_grid_x(:), dot_grid_y(:)];
+num_dots = size(dot_centers, 1);
+total_rays = num_dots * num_rays_per_dot;
+
+% Center the dot grid around (0,0)
+dot_centers(:,1) = dot_centers(:,1) - mean(dot_centers(:,1));
+dot_centers(:,2) = dot_centers(:,2) - mean(dot_centers(:,2));
+
+% Generate purely random samples for all rays
+total_rays = num_dots * num_rays_per_dot;
+
+% Generate random samples in [0, 1]
+samples = rand(total_rays, 2);  % [x, y]
+
+% Map to local dot areas
+x_local = samples(:,1) * dot_width  - dot_width/2;
+y_local = samples(:,2) * dot_height - dot_height/2;
+
+% Repeat dot centers to match total_rays
+dot_indices = repelem(1:num_dots, num_rays_per_dot)';
+dot_x = dot_centers(dot_indices, 1);
+dot_y = dot_centers(dot_indices, 2);
+
+
+
+% Final global positions
+lightfield.position.x = dot_x + x_local;
+lightfield.position.y = dot_y + y_local;
+lightfield.position.z = z_location * ones(total_rays, 1);
+
+% Generate Diverging Directions (Cone Distribution)
+cone_half_angle_radians = deg2rad(cone_half_angle);
+
+cos_theta = (1 - cos(cone_half_angle_radians)) * rand(total_rays,1) + cos(cone_half_angle_radians);
+theta = acos(cos_theta);  % in [0, cone_half_angle]
+phi = 2 * pi * rand(total_rays,1);  % in [0, 2*pi]
+
+lightfield.direction.dx = sin(theta) .* cos(phi);
+lightfield.direction.dy = sin(theta) .* sin(phi);
+lightfield.direction.dz = -cos(theta); 
+
+fprintf("Lightfield generated. Number of light rays = %d. ", total_rays);
+toc
+end
diff --git a/Matlab App v5/function_knife_edge_cutoff.m b/Matlab App v5/function_knife_edge_cutoff.m
new file mode 100644
index 0000000..444f4e0
--- /dev/null
+++ b/Matlab App v5/function_knife_edge_cutoff.m	
@@ -0,0 +1,112 @@
+function [lightfield] = function_knife_edge_cutoff(lightfield, cutoff_location, wavelength, diffraction_type, knife_type, low_lim, upp_lim)
+
+if knife_type == "Top Knife Edge" %Remove lightrays from y = (-Inf, cutoff_location]
+    mask = lightfield.position.y <= cutoff_location;
+    lightfield = remove_lightrays(lightfield,mask);
+    y_distance = abs(lightfield.position.y - cutoff_location);
+elseif knife_type == "Bottom Knife Edge" %Remove lightrays from y = [cutoff_location, Inf)
+    mask = lightfield.position.y >= cutoff_location;
+    lightfield = remove_lightrays(lightfield,mask);
+    y_distance = abs(lightfield.position.y - cutoff_location);
+elseif knife_type == "Left Knife Edge" %Remove lightrays from x = (-Inf, cutoff_location]
+    mask = lightfield.position.x <= cutoff_location;
+    lightfield = remove_lightrays(lightfield,mask);
+    x_distance = abs(lightfield.position.x - cutoff_location);
+else %Remove lightrays from x = [cutoff_location,Inf)
+    mask = lightfield.position.x >= cutoff_location;
+    lightfield = remove_lightrays(lightfield,mask);
+    x_distance = abs(lightfield.position.x - cutoff_location);
+end
+
+%Check which direction for diffraction implementation
+x_exist = exist("x_distance",'var');
+y_exist = exist("y_distance",'var');
+
+if diffraction_type ~= "none"
+    k = 2*pi/wavelength;
+    if x_exist
+        del_theta_x = zeros(size(lightfield.position.x)); %Initialize vector for change in theta_x due to diffraction
+        sigma_x = atan(1./(2*x_distance*k)); %Calculate norm distribution characteristics
+    end
+    if y_exist
+        del_theta_y = zeros(size(lightfield.position.y)); %Initialize vector for change in theta_y due to diffraction
+        sigma_y = atan(1./(2*y_distance*k)); %Calculate norm distribution characteristics
+    end
+
+    %Diffraction model parameters
+    if string(diffraction_type) == 'dynamic limited'
+        %Calculate angle of defelction due to diffraction
+        if x_exist
+            lower_limit_x = -sigma_x;
+            upper_limit_x = -lower_limit_x;
+        end
+        if y_exist
+            lower_limit_y = -sigma_y;
+            upper_limit_y = -lower_limit_y;
+        end
+        %Modeling diffraction in the Monte-Carlo ray-trace environment,
+        %Coffey 1998
+    elseif string(diffraction_type) == "-pi/2 to pi/2 limited"
+        if x_exist
+            lower_limit_x = -pi/2*ones(size(sigma_x));
+            upper_limit_x = -lower_limit_x;
+        end
+        if y_exist
+            lower_limit_y = -pi/2*ones(size(sigma_y));
+            upper_limit_y = -lower_limit_y;
+        end
+    else
+        if x_exist
+            lower_limit_x = low_lim*ones(size(sigma_x));
+            upper_limit_x = upp_lim*ones(size(sigma_x));
+        end
+        if y_exist
+            lower_limit_y = low_lim*ones(size(sigma_y));
+            upper_limit_y = upp_lim*ones(size(sigma_y));
+        end
+    end %End of type of diffraction
+
+    %Calculate angle due to diffraction
+    parfor counter = 1:numel(lightfield.position.x) % For all rays
+        %Generate values within the truncated range
+        if x_exist
+            del_theta_x(counter) = truncate_gaussian(0, sigma_x(counter), lower_limit_x(counter), upper_limit_x(counter));
+        end
+        if y_exist
+            del_theta_y(counter) = truncate_gaussian(0, sigma_y(counter), lower_limit_y(counter), upper_limit_y(counter));
+        end
+    end
+
+    %Edge diffraction in Monte Carlo ray tracing, Freniere et al 1999
+    % DOI:10.1117/12.363773
+
+    %Update new directions in the x and y component by converting from
+    %spherical to cartesian
+    if x_exist
+        lightfield.direction.dx = lightfield.direction.dx + tan(del_theta_x);
+    end
+    if y_exist
+        lightfield.direction.dy = lightfield.direction.dy + tan(del_theta_y);
+    end
+    %Display completion message and stop timer
+    fprintf("Lightfield cutoff.");
+end %End of diffraction check
+
+end
+
+function [lightfield] = remove_lightrays(lightfield,mask)
+%Remove lightrays that are blocked by a mask from memory
+lightfield.position.x(mask) = nan;
+lightfield.position.y(mask) = nan;
+lightfield.position.z(mask) = nan;
+lightfield.direction.dx(mask) = nan;
+lightfield.direction.dy(mask) = nan;
+lightfield.direction.dz(mask) = nan;
+
+lightfield.position.x = lightfield.position.x(~isnan(lightfield.position.x));
+lightfield.position.y = lightfield.position.y(~isnan(lightfield.position.y));
+lightfield.position.z = lightfield.position.z(~isnan(lightfield.position.z));
+lightfield.direction.dx = lightfield.direction.dx(~isnan(lightfield.direction.dx));
+lightfield.direction.dy = lightfield.direction.dy(~isnan(lightfield.direction.dy));
+lightfield.direction.dz = lightfield.direction.dz(~isnan(lightfield.direction.dz));
+end
\ No newline at end of file
diff --git a/Matlab App v5/function_light_propagation_to_camera.m b/Matlab App v5/function_light_propagation_to_camera.m
new file mode 100644
index 0000000..85cebe8
--- /dev/null
+++ b/Matlab App v5/function_light_propagation_to_camera.m	
@@ -0,0 +1,81 @@
+function lightfield = function_light_propagation_to_camera(lightfield, optical_components, wavelength,density_object)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%Function to propagate through all components defined in the system by the
+%user
+
+%optical_components is app.optical_layout.components
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+components_list = optical_components.list; %initialize components list
+num_step = 1;
+%Loop through all optical components
+for counter = 2:max(size(components_list)-1) %Loop through all components that arent the camera or the lightsource
+    str = char(components_list(counter)); %Initialize character array to test which component to propagate through
+    field_name = components_list(counter); %Initialize field name to call component properties
+    
+    %Conduct linear propagation to current optical component
+    lightfield = function_linear_propagation(lightfield,optical_components.z_vector(counter)*1000);
+
+    %Conduct optical propagation through a thin lens
+    if str(1:9) == "Thin_Lens"
+
+        %Set object parameters
+        object_param.f = optical_components.(field_name).focal_length*1000;
+        object_param.d = optical_components.(field_name).diameter*1000;
+
+        %Conduct propagation
+        lightfield = function_optical_object_propagation_v2(lightfield,'thin lens',object_param);
+
+    %Conduct optical propagation through an aperture
+    elseif str(1:8) == "Aperture"
+        low_limit = [];
+        upp_limit = [];
+        if isfield(optical_components.(field_name),'diffraction_low_limit')
+            low_limit = optical_components.(field_name).diffraction_low_limit;
+            upp_limit = optical_components.(field_name).diffraction_upp_limit;
+        end
+        a = max(size(lightfield.direction.dx));
+        type = string(optical_components.(field_name).type);
+        if type == "Circular"
+            lightfield = function_circular_aperture_cutoff_heisenberg_diffraction(lightfield, optical_components.(field_name).cutoff_value, wavelength, optical_components.(field_name).diffraction_model, low_limit, upp_limit);
+        else
+            lightfield = function_knife_edge_cutoff(lightfield,optical_components.(field_name).cutoff_value, wavelength, optical_components.(field_name).diffraction_model, type, low_limit,upp_limit);
+        end
+        b = max(size(lightfield.direction.dx));
+        fprintf('\n Percent Cutoff from aperture = %0.2f',(1-b/a)*100);
+    %Conduct optical propagation through the density gradient
+    elseif str == "Density_Gradient"
+        if density_object.GPU
+            lightfield = function_trace_through_density_grad_rk4_sharma_gpu_v3(lightfield, density_object.loc, density_object);
+        elseif density_object.grid == 'Gridded Data'
+            lightfield = function_trace_through_density_grad_rk4_sharma_gpu_while_loop(lightfield, density_object.loc, density_object);
+        else
+            lightfield = function_trace_through_density_grad_rk4_sharma_nongrid(lightfield, density_object.loc, density_object);
+        end
+
+        if density_object.save %If save lightfield after tracing through density gradient    
+            save(sprintf("%s\\%s.mat",density_object.save_folder, density_object.save_file),"lightfield",'-v7.3');
+        end
+    %Conduct optical propagation through a thick lens
+    else 
+
+        %Set object parameters
+        object_param.d = optical_components.(field_name).diameter*1000;
+        object_param.R1 = optical_components.(field_name).r1*1000;
+        object_param.R2 = optical_components.(field_name).r2*1000;
+        object_param.t = optical_components.(field_name).thickness*1000;
+        object_param.n = optical_components.(field_name).refractive_index;
+
+        %Conduct propagation
+        lightfield = function_optical_object_propagation_v2(lightfield,'thick lens',object_param); 
+
+    end %End of component type determination if statement
+    %figure; quiver(lightfield.position.x(1:num_step:end),lightfield.position.y(1:num_step:end),lightfield.direction.dx(1:num_step:end),lightfield.direction.dy(1:num_step:end))
+    if density_object.progress
+        figure(100)
+        plot3(lightfield.position.x(1:num_step:end),lightfield.position.y(1:num_step:end),lightfield.position.z(1:num_step:end),'.');
+        hold on
+    end
+end %End of optical component loop
+
+end
\ No newline at end of file
diff --git a/Matlab App v5/function_linear_propagation.m b/Matlab App v5/function_linear_propagation.m
new file mode 100644
index 0000000..5bb410e
--- /dev/null
+++ b/Matlab App v5/function_linear_propagation.m	
@@ -0,0 +1,32 @@
+function [lightfield] = function_linear_propagation(lightfield, destination_z)
+% Function for linear propagation of rays to a new z-position
+% INPUTS:
+%   - lightfield: Structure containing ray positions and directions
+%   - destination_z: Target z-coordinate for propagation
+% OUTPUTS:
+%   - new_lightfield: Updated lightfield after propagation
+%   - travel: Structure containing the dx, dy, and dz of travel
+
+% Start timer for performance measurement
+tic
+
+% Calculate the distance to propagate along z-axis
+travel_dz = destination_z - lightfield.position.z;
+
+% Pre-compute the scaling factor to apply to x and y direction components
+scale_factor = travel_dz ./ (lightfield.direction.dz);
+
+% Calculate dx and dy travel using the scale factor
+travel_dx = lightfield.direction.dx .* scale_factor;
+travel_dy = lightfield.direction.dy .* scale_factor;
+
+% Update the new lightfield structure
+% lightfield = lightfield;
+lightfield.position.x = lightfield.position.x + travel_dx;
+lightfield.position.y = lightfield.position.y + travel_dy;
+lightfield.position.z = lightfield.position.z + travel_dz;
+
+% Display completion message with destination z and stop timer
+fprintf("Lightfield linearly propagated (to z = %f). ", destination_z);
+toc
+end
\ No newline at end of file
diff --git a/Matlab App v5/function_optical_object_propagation_v2.m b/Matlab App v5/function_optical_object_propagation_v2.m
new file mode 100644
index 0000000..6718c29
--- /dev/null
+++ b/Matlab App v5/function_optical_object_propagation_v2.m	
@@ -0,0 +1,77 @@
+function [lightfield] = function_optical_object_propagation(lightfield,object_type,object_param,n0)
+% function_optical_object_propagation Propagate a lightfield through an
+% optical component
+%
+% A = function_optical_object_propagation(lightfield,type,object)
+% outputs the lightfield after tracing rays through an optical component
+% 
+% lightfield is the incoming lightfield
+%
+% object type => 'thin lens', 'thick lens' and 'mirror' are
+% supported in this current version
+%
+% object_param provides the object parameters such as diameter, focal length, etc
+%
+% n0 
+
+
+%% Remove lightrays that don't hit the surface of optical component
+
+mask = sqrt(lightfield.position.x.^2 + lightfield.position.y.^2) > (object_param.d / 2);
+
+% Apply the mask to all position and direction components
+% Rays outside the aperture are set to NaN in the new lightfield structure
+lightfield.position.x(mask) = nan;
+lightfield.position.y(mask) = nan;
+lightfield.position.z(mask) = nan;
+lightfield.direction.dx(mask) = nan;
+lightfield.direction.dy(mask) = nan;
+lightfield.direction.dz(mask) = nan;
+
+lightfield.position.x = lightfield.position.x(~isnan(lightfield.position.x));
+lightfield.position.y = lightfield.position.y(~isnan(lightfield.position.y));
+lightfield.position.z = lightfield.position.z(~isnan(lightfield.position.z));
+lightfield.direction.dx = lightfield.direction.dx(~isnan(lightfield.direction.dx));
+lightfield.direction.dy = lightfield.direction.dy(~isnan(lightfield.direction.dy));
+lightfield.direction.dz = lightfield.direction.dz(~isnan(lightfield.direction.dz));
+
+%% Conduct Propagation through component
+
+%Propogation through optical component
+if strcmpi(object_type,"thin lens") ||  strcmpi(object_type,"mirror") %Thin lens and mirror propagation
+
+    fprintf("Conducting propagation through a %s.\n", lower(object_type))
+    multiplier = -1/object_param.f;
+    lightfield.direction.dx = lightfield.direction.dx + multiplier * lightfield.position.x;
+    lightfield.direction.dy = lightfield.direction.dy + multiplier * lightfield.position.y;
+
+elseif strcmpi(object_type,"thick lens") %Thick lens propagation
+    %fprintf("Conducting propagation through a %s %s.\n", upper(object_param.material), lower(object_type))
+    if nargin < 4
+        n0 = 1; % Assume air if not specified
+    end
+
+    %Construct matrix for thick lens
+    D1 = (object_param.n-n0)/object_param.R1;
+    D2 = (n0-object_param.n)/(-1*object_param.R2);
+    A11 = 1-D2*object_param.t/object_param.n;
+    A12 = -D1-D2-D1*D2*object_param.t/object_param.n;
+    A21 = object_param.t/object_param.n;
+    A22 = 1-D1*object_param.t/object_param.n;
+    %fprintf("Effective focus of lens is %f mm.\n",(-1/A12)/1000);
+    
+    %Calculate angle with respect to z plane
+    alpha_0x = lightfield.direction.dx;
+    alpha_0y = lightfield.direction.dy;
+
+    alpha_exit_x = A11*alpha_0x+A12*lightfield.position.x;
+    lightfield.position.x = A21*alpha_0x+A22*lightfield.position.x;
+    lightfield.direction.dx = alpha_exit_x;
+    alpha_exit_y = A11*alpha_0y+A12*lightfield.position.y;
+    lightfield.direction.dy = alpha_exit_y;
+    lightfield.position.y = A21*alpha_0y+A22*lightfield.position.y;
+    lightfield.position.z = lightfield.position.z-object_param.t;
+    
+end
+
+end
\ No newline at end of file
diff --git a/Matlab App v5/function_trace_sharma.m b/Matlab App v5/function_trace_sharma.m
new file mode 100644
index 0000000..e833ac0
--- /dev/null
+++ b/Matlab App v5/function_trace_sharma.m	
@@ -0,0 +1,111 @@
+function lightfield = function_trace_sharma(lightfield, start_z, density_object)
+% Implements Sharma et al. (1982) ray tracing through graded-index media
+% using 4th-order Runge-Kutta integration
+
+% Extract and scale grid
+x = density_object.rho.x * 1000;
+y = density_object.rho.y * 1000;
+z = (density_object.rho.z + start_z) * 1000;
+z_step = abs(density_object.step_size); % Arclength step size (mm)
+
+% Refractive index field
+rho = density_object.rho.rho;
+n = (density_object.fluid_glasstone * rho + 1);
+
+% Compute gradient of n
+[dn_dx, dn_dy, dn_dz] = gradient(n, ...
+    mean(diff(x(:,1,1))), ...
+    mean(diff(y(1,:,1))), ...
+    mean(diff(z(1,1,:))));
+
+% Create interpolants
+F_n = griddedInterpolant(x, y, z, n, lower(density_object.interp_method));
+F_gradx = griddedInterpolant(x, y, z, dn_dx, lower(density_object.interp_method));
+F_grady = griddedInterpolant(x, y, z, dn_dy, lower(density_object.interp_method));
+F_gradz = griddedInterpolant(x, y, z, dn_dz, lower(density_object.interp_method));
+
+% Initialize ray position and direction
+pos = cat(4, lightfield.position.x, lightfield.position.y, lightfield.position.z);
+dir = cat(4, lightfield.direction.dx, lightfield.direction.dy, lightfield.direction.dz);
+
+% Normalize direction
+dir = dir ./ vecnorm(dir,2,4);
+
+% Compute number of steps
+z_extent = abs(z(end) - z(1));
+num_steps = round(z_extent / z_step);
+
+% RK4 Integration
+for counter = 1:num_steps
+    % k1
+    n1 = F_n(pos(:,:,:,1), pos(:,:,:,2), pos(:,:,:,3));
+    gradn1 = cat(4, F_gradx(pos(:,:,:,1), pos(:,:,:,2), pos(:,:,:,3)), ...
+                  F_grady(pos(:,:,:,1), pos(:,:,:,2), pos(:,:,:,3)), ...
+                  F_gradz(pos(:,:,:,1), pos(:,:,:,2), pos(:,:,:,3)));
+    t1 = dir;
+    tdotgrad1 = sum(t1 .* gradn1, 4);
+    dt1 = (gradn1 - tdotgrad1 .* t1) ./ n1;
+    k1_pos = t1;
+    k1_dir = dt1;
+
+    % k2
+    pos2 = pos + 0.5 * z_step * k1_pos;
+    dir2 = dir + 0.5 * z_step * k1_dir;
+    dir2 = dir2 ./ vecnorm(dir2,2,4);
+    n2 = F_n(pos2(:,:,:,1), pos2(:,:,:,2), pos2(:,:,:,3));
+    gradn2 = cat(4, F_gradx(pos2(:,:,:,1), pos2(:,:,:,2), pos2(:,:,:,3)), ...
+                  F_grady(pos2(:,:,:,1), pos2(:,:,:,2), pos2(:,:,:,3)), ...
+                  F_gradz(pos2(:,:,:,1), pos2(:,:,:,2), pos2(:,:,:,3)));
+    t2 = dir2;
+    tdotgrad2 = sum(t2 .* gradn2, 4);
+    dt2 = (gradn2 - tdotgrad2 .* t2) ./ n2;
+    k2_pos = t2;
+    k2_dir = dt2;
+
+    % k3
+    pos3 = pos + 0.5 * z_step * k2_pos;
+    dir3 = dir + 0.5 * z_step * k2_dir;
+    dir3 = dir3 ./ vecnorm(dir3,2,4);
+    n3 = F_n(pos3(:,:,:,1), pos3(:,:,:,2), pos3(:,:,:,3));
+    gradn3 = cat(4, F_gradx(pos3(:,:,:,1), pos3(:,:,:,2), pos3(:,:,:,3)), ...
+                  F_grady(pos3(:,:,:,1), pos3(:,:,:,2), pos3(:,:,:,3)), ...
+                  F_gradz(pos3(:,:,:,1), pos3(:,:,:,2), pos3(:,:,:,3)));
+    t3 = dir3;
+    tdotgrad3 = sum(t3 .* gradn3, 4);
+    dt3 = (gradn3 - tdotgrad3 .* t3) ./ n3;
+    k3_pos = t3;
+    k3_dir = dt3;
+
+    % k4
+    pos4 = pos + z_step * k3_pos;
+    dir4 = dir + z_step * k3_dir;
+    dir4 = dir4 ./ vecnorm(dir4,2,4);
+    n4 = F_n(pos4(:,:,:,1), pos4(:,:,:,2), pos4(:,:,:,3));
+    gradn4 = cat(4, F_gradx(pos4(:,:,:,1), pos4(:,:,:,2), pos4(:,:,:,3)), ...
+                  F_grady(pos4(:,:,:,1), pos4(:,:,:,2), pos4(:,:,:,3)), ...
+                  F_gradz(pos4(:,:,:,1), pos4(:,:,:,2), pos4(:,:,:,3)));
+    t4 = dir4;
+    tdotgrad4 = sum(t4 .* gradn4, 4);
+    dt4 = (gradn4 - tdotgrad4 .* t4) ./ n4;
+    k4_pos = t4;
+    k4_dir = dt4;
+
+    % Update
+    pos = pos + (z_step/6) * (k1_pos + 2*k2_pos + 2*k3_pos + k4_pos);
+    dir = dir + (z_step/6) * (k1_dir + 2*k2_dir + 2*k3_dir + k4_dir);
+    dir = dir ./ vecnorm(dir,2,4); % normalize
+
+    if mod(counter,10) == 0
+        fprintf("RK4 Sharma step %d / %d\n", counter, num_steps);
+    end
+end
+
+% Store back into 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);
+
+end
diff --git a/Matlab App v5/function_trace_through_density_grad.asv b/Matlab App v5/function_trace_through_density_grad.asv
new file mode 100644
index 0000000..e88e61b
--- /dev/null
+++ b/Matlab App v5/function_trace_through_density_grad.asv	
@@ -0,0 +1,128 @@
+function [lightfield] = function_trace_through_density_grad(lightfield, start_z, density_object)
+
+% x = linspace(-20*1024/0.6562/2,20*1024/0.6562/2,size(rho,2));
+% y = linspace(-20*1024/0.6562/2,20*1024/0.6562/2,size(rho,1));
+% z = linspace(end_z,start_z,size(rho,3));
+%z = linspace(990,1000,size(rho,3));
+%z_step = (z(1)-z(2));
+
+%Initialize variables
+x = density_object.rho.x*1000;
+y = density_object.rho.y*1000;
+z = (density_object.rho.z+start_z)*1000;
+z_step = -density_object.step_size;
+num_step = (z(end)-z(1))/z_step;
+rho = density_object.rho.val;
+
+
+rho = 10*rho;
+%gladstone_dale_air = 0.225e-3; % kg/m^3
+
+%Calculate n field
+n = (app.density_object.fluid_glasstone*rho+1);
+
+%Interpolate based off grid or no grid
+if density_object.grid == "Gridded Data"
+    if min(size(x))== 1
+        [x,y,z] = ndgrid(x,y,z);
+    end
+    F = griddedInterpolant(X,Y,Z,n,lower(density_object.interp_method));
+    
+else
+    F = scatteredInterpolant(X(:),Y(:),Z(:),n(:),lower(density_object.interp_method));
+end
+
+%for counter = 1:size(z, 2)
+for counter = 1:num_step;
+    tic
+    if ~exist("lightfield.n.original", 'var')
+        lightfield.n.original = 1.0003 * ones(size(lightfield.position.x)); %Refractive index of air is 1.0003
+    else
+        lightfield.n.original = lightfield.n.inquiry;
+    end
+
+    % RK4 Step 1
+    k1x = lightfield.direction.dx ./ lightfield.direction.dz;
+    k1y = lightfield.direction.dy ./ lightfield.direction.dz;
+    k1z = 1; % dz_step assumed constant
+
+    % Inquiry position for k2
+    inquiry.position.x = lightfield.position.x + 0.5 * z_step * k1x;
+    inquiry.position.y = lightfield.position.y + 0.5 * z_step * k1y;
+    inquiry.position.z = lightfield.position.z + 0.5 * z_step * k1z;
+
+    % lightfield.n.inquiry = interp3(y, x, z, n, ...
+    %     inquiry.position.y, ...
+    %     inquiry.position.x, ...
+    %     inquiry.position.z, "linear", nan);
+    % lightfield.n.inquiry(isnan(lightfield.n.inquiry)) = lightfield.n.original(isnan(lightfield.n.inquiry));
+    lightfield.n.inquiry = F(inquiry.position.y,inquiry.position.x,inquiry.position.z);
+    lightfield.n.inquiry(isnan(lightfield.n.inquiry)) = lightfield.n.original(isnan(lightfield.n.inquiry));
+
+    k2x = (lightfield.direction.dx .* lightfield.n.original ./ lightfield.n.inquiry) ./ lightfield.direction.dz;
+    k2y = (lightfield.direction.dy .* lightfield.n.original ./ lightfield.n.inquiry) ./ lightfield.direction.dz;
+    k2z = 1;
+
+    % Inquiry position for k3
+    inquiry.position.x = lightfield.position.x + 0.5 * z_step * k2x;
+    inquiry.position.y = lightfield.position.y + 0.5 * z_step * k2y;
+    inquiry.position.z = lightfield.position.z + 0.5 * z_step * k2z;
+
+    % lightfield.n.inquiry = interp3(y, x, z, n, ...
+    %     inquiry.position.y, ...
+    %     inquiry.position.x, ...
+    %     inquiry.position.z, "linear", nan);
+    % lightfield.n.inquiry(isnan(lightfield.n.inquiry)) = lightfield.n.original(isnan(lightfield.n.inquiry));
+
+    lightfield.n.inquiry = F(inquiry.position.y,inquiry.position.x,inquiry.position.z);
+    lightfield.n.inquiry(isnan(lightfield.n.inquiry)) = lightfield.n.original(isnan(lightfield.n.inquiry));
+
+    k3x = (lightfield.direction.dx .* lightfield.n.original ./ lightfield.n.inquiry) ./ lightfield.direction.dz;
+    k3y = (lightfield.direction.dy .* lightfield.n.original ./ lightfield.n.inquiry) ./ lightfield.direction.dz;
+    k3z = 1;
+
+    % Inquiry position for k4
+    inquiry.position.x = lightfield.position.x + z_step * k3x;
+    inquiry.position.y = lightfield.position.y + z_step * k3y;
+    inquiry.position.z = lightfield.position.z + z_step * k3z;
+
+    % lightfield.n.inquiry = interp3(y, x, z, n, ...
+    %     inquiry.position.y, ...
+    %     inquiry.position.x, ...
+    %     inquiry.position.z, "linear", nan);
+    % lightfield.n.inquiry(isnan(lightfield.n.inquiry)) = lightfield.n.original(isnan(lightfield.n.inquiry));\
+    lightfield.n.inquiry = F(inquiry.position.y,inquiry.position.x,inquiry.position.z);
+    lightfield.n.inquiry(isnan(lightfield.n.inquiry)) = lightfield.n.original(isnan(lightfield.n.inquiry));
+
+    k4x = (lightfield.direction.dx .* lightfield.n.original ./ lightfield.n.inquiry) ./ lightfield.direction.dz;
+    k4y = (lightfield.direction.dy .* lightfield.n.original ./ lightfield.n.inquiry) ./ lightfield.direction.dz;
+    k4z = 1;
+
+    % Combine RK4 steps for positions
+    lightfield.position.x = lightfield.position.x + (z_step / 6) * (k1x + 2*k2x + 2*k3x + k4x);
+    lightfield.position.y = lightfield.position.y + (z_step / 6) * (k1y + 2*k2y + 2*k3y + k4y);
+    lightfield.position.z = lightfield.position.z + (z_step / 6) * (k1z + 2*k2z + 2*k3z + k4z);
+
+    % Update direction using Snell's law
+    lightfield.direction.dx = lightfield.direction.dx .* lightfield.n.original ./ lightfield.n.inquiry;
+    lightfield.direction.dy = lightfield.direction.dy .* lightfield.n.original ./ lightfield.n.inquiry;
+    lightfield.direction.dz = lightfield.direction.dz; % Adjust if dz_step changes with refractive index
+
+    fprintf('For Loop Counter = %d ', counter);
+    toc
+end %End of rk4 step preocess
+
+clearvars inquiry;
+
+% Exiting the density field
+lightfield.n.original = lightfield.n.inquiry;
+lightfield.n.inquiry = 1.003*ones(size(lightfield.position.x));
+% Use Snell's law to update the direction
+lightfield.direction.dx = lightfield.direction.dx .* lightfield.n.original ./ lightfield.n.inquiry;
+lightfield.direction.dy = lightfield.direction.dy .* lightfield.n.original ./ lightfield.n.inquiry;
+% Propagate rays by a step
+lightfield.position.x = lightfield.position.x + lightfield.direction.dx ./ lightfield.direction.dz * z_step;
+lightfield.position.y = lightfield.position.y + lightfield.direction.dy ./ lightfield.direction.dz * z_step;
+lightfield.position.z = lightfield.position.z + z_step;
+
+end
\ No newline at end of file
diff --git a/Matlab App v5/function_trace_through_density_grad.m b/Matlab App v5/function_trace_through_density_grad.m
new file mode 100644
index 0000000..89c9af6
--- /dev/null
+++ b/Matlab App v5/function_trace_through_density_grad.m	
@@ -0,0 +1,136 @@
+function [lightfield] = function_trace_through_density_grad(lightfield, start_z, density_object)
+tic
+% x = linspace(-20*1024/0.6562/2,20*1024/0.6562/2,size(rho,2));
+% y = linspace(-20*1024/0.6562/2,20*1024/0.6562/2,size(rho,1));
+% z = linspace(end_z,start_z,size(rho,3));
+%z = linspace(990,1000,size(rho,3));
+%z_step = (z(1)-z(2));
+
+%Initialize variables
+x = density_object.rho.x*1000;
+y = density_object.rho.y*1000;
+z = (density_object.rho.z+start_z)*1000;
+z_step = -density_object.step_size;
+num_step = abs((z(end)-z(1))/z_step);
+rho = density_object.rho.rho;
+
+
+%rho = 10*rho;
+%gladstone_dale_air = 0.225e-3; % kg/m^3
+
+%Calculate n field
+n = (density_object.fluid_glasstone*rho+1);
+
+%Interpolate based off grid or no grid
+if density_object.grid == "Gridded Data"
+    if min(size(x))== 1
+        [x,y,z] = ndgrid(x,y,z);
+    end
+    F = griddedInterpolant(x,y,z,n,lower(density_object.interp_method));
+else
+    F = scatteredInterpolant(x(:),y(:),z(:),n(:),lower(density_object.interp_method));
+end
+
+lightfield.n.original = 1.0003 * ones(size(lightfield.position.x));
+lightfield.n.inquiry = lightfield.n.original;
+
+%for counter = 1:size(z, 2)
+for counter = 1:num_step;
+    tic
+    % if ~exist("lightfield.n.original", 'var')
+    %     lightfield.n.original = 1.0003 * ones(size(lightfield.position.x)); %Refractive index of air is 1.0003
+    % else
+    %     lightfield.n.original = lightfield.n.inquiry;
+    % end
+
+    % RK4 Step 1
+    k1x = lightfield.direction.dx ./ lightfield.direction.dz;
+    k1y = lightfield.direction.dy ./ lightfield.direction.dz;
+    k1z = 1; % dz_step assumed constant
+
+    % Inquiry position for k2
+    inquiry.position.x = lightfield.position.x + 0.5 * z_step * k1x;
+    inquiry.position.y = lightfield.position.y + 0.5 * z_step * k1y;
+    inquiry.position.z = lightfield.position.z + 0.5 * z_step * k1z;
+
+    % lightfield.n.inquiry = interp3(y, x, z, n, ...
+    %     inquiry.position.y, ...
+    %     inquiry.position.x, ...
+    %     inquiry.position.z, "linear", nan);
+    % lightfield.n.inquiry(isnan(lightfield.n.inquiry)) = lightfield.n.original(isnan(lightfield.n.inquiry));
+    lightfield.n.inquiry = F(inquiry.position.x,inquiry.position.y,inquiry.position.z);
+    lightfield.n.inquiry(isnan(lightfield.n.inquiry)) = lightfield.n.original(isnan(lightfield.n.inquiry));
+
+    k2x = (lightfield.direction.dx .* lightfield.n.original ./ lightfield.n.inquiry) ./ lightfield.direction.dz;
+    k2y = (lightfield.direction.dy .* lightfield.n.original ./ lightfield.n.inquiry) ./ lightfield.direction.dz;
+    k2z = 1;
+
+    % Inquiry position for k3
+    inquiry.position.x = lightfield.position.x + 0.5 * z_step * k2x;
+    inquiry.position.y = lightfield.position.y + 0.5 * z_step * k2y;
+    inquiry.position.z = lightfield.position.z + 0.5 * z_step * k2z;
+
+    % lightfield.n.inquiry = interp3(y, x, z, n, ...
+    %     inquiry.position.y, ...
+    %     inquiry.position.x, ...
+    %     inquiry.position.z, "linear", nan);
+    % lightfield.n.inquiry(isnan(lightfield.n.inquiry)) = lightfield.n.original(isnan(lightfield.n.inquiry));
+
+    lightfield.n.inquiry = F(inquiry.position.x,inquiry.position.y,inquiry.position.z);
+    lightfield.n.inquiry(isnan(lightfield.n.inquiry)) = lightfield.n.original(isnan(lightfield.n.inquiry));
+
+    k3x = (lightfield.direction.dx .* lightfield.n.original ./ lightfield.n.inquiry) ./ lightfield.direction.dz;
+    k3y = (lightfield.direction.dy .* lightfield.n.original ./ lightfield.n.inquiry) ./ lightfield.direction.dz;
+    k3z = 1;
+
+    % Inquiry position for k4
+    inquiry.position.x = lightfield.position.x + z_step * k3x;
+    inquiry.position.y = lightfield.position.y + z_step * k3y;
+    inquiry.position.z = lightfield.position.z + z_step * k3z;
+
+    % lightfield.n.inquiry = interp3(y, x, z, n, ...
+    %     inquiry.position.y, ...
+    %     inquiry.position.x, ...
+    %     inquiry.position.z, "linear", nan);
+    % lightfield.n.inquiry(isnan(lightfield.n.inquiry)) = lightfield.n.original(isnan(lightfield.n.inquiry));\
+    lightfield.n.inquiry = F(inquiry.position.x,inquiry.position.y,inquiry.position.z);
+    lightfield.n.inquiry(isnan(lightfield.n.inquiry)) = lightfield.n.original(isnan(lightfield.n.inquiry));
+
+    k4x = (lightfield.direction.dx .* lightfield.n.original ./ lightfield.n.inquiry) ./ lightfield.direction.dz;
+    k4y = (lightfield.direction.dy .* lightfield.n.original ./ lightfield.n.inquiry) ./ lightfield.direction.dz;
+    k4z = 1;
+
+    % Combine RK4 steps for positions
+    lightfield.position.x = lightfield.position.x + (z_step / 6) * (k1x + 2*k2x + 2*k3x + k4x);
+    lightfield.position.y = lightfield.position.y + (z_step / 6) * (k1y + 2*k2y + 2*k3y + k4y);
+    lightfield.position.z = lightfield.position.z + (z_step / 6) * (k1z + 2*k2z + 2*k3z + k4z);
+
+    % Update direction using Snell's law
+    lightfield.direction.dx = lightfield.direction.dx .* lightfield.n.original ./ lightfield.n.inquiry;
+    lightfield.direction.dy = lightfield.direction.dy .* lightfield.n.original ./ lightfield.n.inquiry;
+    lightfield.direction.dz = lightfield.direction.dz .* lightfield.n.original ./ lightfield.n.inquiry;
+
+    %Normalize Components
+    % ray_norm = sqrt(lightfield.direction.dx.^2+lightfield.direction.dy.^2+lightfield.direction.dz.^2);
+    % lightfield.direction.dx = lightfield.direction.dx./ray_norm;
+    % lightfield.direction.dy = lightfield.direction.dy./ray_norm;
+    % lightfield.direction.dz = lightfield.direction.dz./ray_norm;
+
+    fprintf('For Loop Counter = %d ', counter);
+    toc
+end %End of rk4 step preocess
+
+clearvars inquiry;
+
+% Exiting the density field
+lightfield.n.original = lightfield.n.inquiry;
+lightfield.n.inquiry = 1.003*ones(size(lightfield.position.x));
+% Use Snell's law to update the direction
+lightfield.direction.dx = lightfield.direction.dx .* lightfield.n.original ./ lightfield.n.inquiry;
+lightfield.direction.dy = lightfield.direction.dy .* lightfield.n.original ./ lightfield.n.inquiry;
+% Propagate rays by a step
+lightfield.position.x = lightfield.position.x + lightfield.direction.dx ./ lightfield.direction.dz * z_step;
+lightfield.position.y = lightfield.position.y + lightfield.direction.dy ./ lightfield.direction.dz * z_step;
+lightfield.position.z = lightfield.position.z + z_step;
+toc
+end
\ No newline at end of file
diff --git a/Matlab App v5/function_trace_through_density_grad_rk4_sharma_gpu_v3.m b/Matlab App v5/function_trace_through_density_grad_rk4_sharma_gpu_v3.m
new file mode 100644
index 0000000..c06433f
--- /dev/null
+++ b/Matlab App v5/function_trace_through_density_grad_rk4_sharma_gpu_v3.m	
@@ -0,0 +1,134 @@
+function [lightfield] = function_trace_through_density_grad_rk4_sharma_gpu_v3(lightfield, start_z, density_object)
+% GPU RK4 ray tracer using finite-difference gradients and interp3
+tic
+interp_method = density_object.interp_method;
+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+start_z) * 1000;
+rho = permute(density_object.rho.rho, [2,1,3]);
+n = density_object.fluid_glasstone * rho + 1;
+
+% Move to gpu
+[x, y, z] = deal(gpuArray(x), gpuArray(y), gpuArray(z));
+n_grid = gpuArray(n);
+
+
+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 = density_object.partition;
+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('\nNumber 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')
+        break;
+    end
+end
+fprintf('\nLightrays passed through density gradient in %d total steps\n',counter)
+
+% 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/Matlab App v5/function_trace_through_density_grad_rk4_sharma_nongrid.m b/Matlab App v5/function_trace_through_density_grad_rk4_sharma_nongrid.m
new file mode 100644
index 0000000..2184186
--- /dev/null
+++ b/Matlab App v5/function_trace_through_density_grad_rk4_sharma_nongrid.m	
@@ -0,0 +1,107 @@
+function [lightfield] = function_trace_through_density_grad_rk4_sharma_nongrid(lightfield, start_z, density_object)
+
+% 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+start_z) * 1000;
+rho = permute(density_object.rho.rho, [2,1,3])
+n = density_object.fluid_glasstone * rho + 1;
+x_min = min(x);
+x_max = max(x);
+y_min = min(y);
+y_max = max(y);
+
+Fn = scatteredInterpolant({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;
+while sum(is_active) > 0
+    idx = find(is_active);
+    if isempty(idx); break; end
+    
+    n0 = Fn(pos(idx,1), pos(idx,2), pos(idx,3));
+    gradn0 = finite_gradient(Fn, pos(idx,:), delta);
+    delta_t = step_size;
+    Tn = dir(idx,:) .* n0;
+    %Calculate A
+    A = gradn0 .* n0 .* step_size;
+
+    %Update Position and Calculate B
+    pos1 = pos(idx,:) + (0.5 * Tn + 1/8 * A) .* delta_t;
+    n1 = Fn(pos1(:,1), pos1(:,2), pos1(:,3));
+    gradn1 = finite_gradient(Fn, pos1, delta);
+    B = gradn1 .* n1 .* delta_t;
+
+    %Update Position and Calculate C
+    pos2 = pos(idx,:) + (Tn + 0.5 * B) .* delta_t;
+    n2 = Fn(pos2(:,1), pos2(:,2), pos2(:,3));
+    gradn2 = finite_gradient(Fn, pos2, delta);
+    C = gradn2 .* n2 .* delta_t;
+
+    %Clalculate 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(:);
+    dir(idx,:) = normalize_vector(Tn_new ./ n0);
+
+    % 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 | temp_pos(idx,3) < z_min;
+    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;
+    dir(invalid,:) = NaN;
+    pos(invalid,:) = NaN;
+    pos(idx,:) = temp_pos(idx,:);
+    curr_active_count = sum(is_active);
+    if curr_active_count <= prev_active_count-100
+        fprintf('\n Number of lightrays remaining = %d', curr_active_count);
+        prev_active_count = curr_active_count;
+    end
+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/Matlab App v5/function_trace_through_density_grad_rk4_sharma_while_loop.m b/Matlab App v5/function_trace_through_density_grad_rk4_sharma_while_loop.m
new file mode 100644
index 0000000..7100d56
--- /dev/null
+++ b/Matlab App v5/function_trace_through_density_grad_rk4_sharma_while_loop.m	
@@ -0,0 +1,107 @@
+function [lightfield] = function_trace_through_density_grad_rk4_sharma_while_loop(lightfield, start_z, density_object)
+
+% 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+start_z) * 1000;
+rho = permute(density_object.rho.rho, [2,1,3])
+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;
+while sum(is_active) > 0
+    idx = find(is_active);
+    if isempty(idx); break; end
+    
+    n0 = Fn(pos(idx,1), pos(idx,2), pos(idx,3));
+    gradn0 = finite_gradient(Fn, pos(idx,:), delta);
+    delta_t = step_size;
+    Tn = dir(idx,:) .* n0;
+    %Calculate A
+    A = gradn0 .* n0 .* step_size;
+
+    %Update Position and Calculate B
+    pos1 = pos(idx,:) + (0.5 * Tn + 1/8 * A) .* delta_t;
+    n1 = Fn(pos1(:,1), pos1(:,2), pos1(:,3));
+    gradn1 = finite_gradient(Fn, pos1, delta);
+    B = gradn1 .* n1 .* delta_t;
+
+    %Update Position and Calculate C
+    pos2 = pos(idx,:) + (Tn + 0.5 * B) .* delta_t;
+    n2 = Fn(pos2(:,1), pos2(:,2), pos2(:,3));
+    gradn2 = finite_gradient(Fn, pos2, delta);
+    C = gradn2 .* n2 .* delta_t;
+
+    %Clalculate 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(:);
+    dir(idx,:) = normalize_vector(Tn_new ./ n0);
+
+    % 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 | temp_pos(idx,3) < z_min;
+    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;
+    dir(invalid,:) = NaN;
+    pos(invalid,:) = NaN;
+    pos(idx,:) = temp_pos(idx,:);
+    curr_active_count = sum(is_active);
+    if curr_active_count <= prev_active_count-100
+        fprintf('\n Number of lightrays remaining = %d', curr_active_count);
+        prev_active_count = curr_active_count;
+    end
+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/Matlab App v5/truncate_gaussian.m b/Matlab App v5/truncate_gaussian.m
new file mode 100644
index 0000000..a1b4fdb
--- /dev/null
+++ b/Matlab App v5/truncate_gaussian.m	
@@ -0,0 +1,9 @@
+% Truncated Gaussian function
+function val = truncate_gaussian(mu, sigma, lower, upper)
+    while true
+        val = normrnd(mu, sigma); % Generate a normal random value
+        if val >= lower && val <= upper
+            break; % Accept the value if within range
+        end
+    end
+end
\ No newline at end of file