Add files via upload

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,,,,

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

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

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
+