4D spherical shell heat map

I have written a code (attached) to plot particles striking a sphere of radius r – I can produce the image in example_1. I divide the shell into segments and calculate the number of particles falling within each.
I would like to represent the data as a smooth 3d heatmap over the spherical shell. I can produce the image shown in example_2 using scatter3 and a color bar, but I would like a smooth surface. I have tried to follow the example of my previous question using an isosurface –
however ux != uy != uz and I have the "Number of elements must not change" error.
I have also tried meshgrid(), but the size of photons results in an array that excedes maximum size preference.

Is there an easier way to smooth this data?

photons = readmatrix(‘input.txt’);

r = 70; % in unit mm

prop = zeros(length(photons),3);

for i = 1:length(photons)

c(i) = ((photons(i,1)*photons(i,1)) + (photons(i,2)*photons(i,2)) + (photons(i,3)*photons(i,3)) – (r*r));
b(i) = 2*((photons(i,1)*photons(i,7)) + (photons(i,2)*photons(i,8)) + (photons(i,3)*photons(i,9)));
a(i) = ((photons(i,7)*photons(i,7)) + (photons(i,8)*photons(i,8)) + (photons(i,9)*photons(i,9)));

t(i) = (-b(i) + sqrt(b(i).^2 – 4*a(i)*c(i)))/(2*a(i));

prop(i,1) = photons(i,1) + t(i)*photons(i,7);
prop(i,2) = photons(i,2) + t(i)*photons(i,8);
prop(i,3) = photons(i,3) + t(i)*photons(i,9);

prop_sph(i,3) = sqrt(prop(i,1).^2 + prop(i,2).^2 + prop(i,3).^2); % calculate spherical radius (should be r)

prop_sph(i,1) = acos((prop(i,3))/(prop_sph(i,3))); % calculate spherical theta (elevation)

if prop(i,1) > 0 && prop(i,2) >= 0
prop_sph(i,2) = atan(prop(i,2)/prop(i,1)); % calculate spherical phi (azimuth)
end

if prop(i,1) > 0 && prop(i,2) < 0
prop_sph(i,2) = atan(prop(i,2)/prop(i,1)) + 2*pi;
end

if prop(i,1) < 0 && prop(i,2) >=0
prop_sph(i,2) = atan(prop(i,2)/prop(i,1)) + pi;
end

if prop(i,1) < 0 && prop(i,2) < 0
prop_sph(i,2) = (atan(prop(i,2)/prop(i,1)) – pi) + 2*pi ;
end

if prop(i,1) == 0 && prop(i,2) > 0
prop_sph(i,2) = pi / 2;
end

if prop(i,1) == 0 && prop(i,2) < 0
prop_sph(i,2) = (-pi / 2) + 2*pi;
end

if prop(i,1) == 0 && prop(i,2) == 0
prop_sph(i,2) = NaN;
end

end

% plot the photon distribution in Cartesian coordinates.
% scatter3(prop(:,1),prop(:,2),prop(:,3), 10,’k’,’filled’,’MarkerFaceAlpha’,.9);

% sort photons into bins of size d_th*d_phi in spherical coordinates

max_th = pi;
max_phi = 2*pi;

d_th = max_th / 100;
d_phi = max_phi / 100;

l_th = round((max_th / d_th) * (max_phi / d_phi));

theta = (0:d_th:max_th)’;
phi = (0:d_phi:max_phi)’;

mat = [zeros(l_th,1)’; zeros(l_th,1)’; zeros(l_th,1)’;]’;

n = 1;

for j = 1:length(phi)
for i = 1:length(theta)

mat(n,1) = theta(i);
mat(n,2) = phi(j);

n = n+1;
end
end

for w = 1:length(mat)
for k = 1:length(prop_sph)
if (((prop_sph(k,1) >= mat(w,1)) && (prop_sph(k,1) < mat(w,1) + d_th)) && ((prop_sph(k,2) >= mat(w,2)) && (prop_sph(k,2) < mat(w,2) + d_phi)))
mat(w,3) = mat(w,3) + 1;
end
end
end

% ——- convert back to cartesian coordinates and include color data

for i = 1:length(mat)

prop_map_fin(i,1) = prop_sph(i,3)*sin(mat(i,1) + d_th/2)*cos(mat(i,2) + d_phi/2);
prop_map_fin(i,2) = prop_sph(i,3)*sin(mat(i,1) + d_th/2)*sin(mat(i,2) + d_phi/2);
prop_map_fin(i,3) = prop_sph(i,3)*cos(mat(i,1) + d_th/2);
prop_map_fin(i,4) = mat(i,3);

% X = prop_sph(:,3).*sin(prop_sph(:,1)).*cos(prop_sph(:,2));
% Y = prop_sph(:,3).*sin(prop_sph(:,1)).*sin(prop_sph(:,2));
% Z = prop_sph(:,3).*cos(prop_sph(:,1));

end

% — visualize ——————————————————–

%surf(prop_map_fin(:,1),prop_map_fin(:,2),prop_map_fin(:,3),prop_map_fin(:,4))

scatter3(prop_map_fin(:,1),prop_map_fin(:,2),prop_map_fin(:,3), 10,prop_map_fin(:,4), ‘filled’);
cb = colorbar; % create and label the colorbar
cb.Label.String = ”;

% xslice = [5 9.9]; % define the cross sections to view
% yslice = 3;
% zslice = ([-3 0]);

%slice(prop_map_fin(:,1), prop_map_fin(:,2), prop_map_fin(:,3), prop_map_fin(:,4), xslice, yslice, zslice)

%
T = array2table(prop_map_fin);

x = T{:, 1}; y = T{:, 2}; z = T{:, 3}; c2 = T{:, 4};
ux = unique(x); uy = unique(y); uz = unique(z);

data_sorted = sortrows(T, 1:4);
v = reshape(data_sorted{:, 4}, length(ux),length(uy),length(uz));

isosurface(uz, uy, ux, v, 0);

%[X2,Y2,Z2] = meshgrid(prop_map_fin(:,1), prop_map_fin(:,2), prop_map_fin(:,3));

function [x y] = GetCircle(r, h, k, a, b)
t = linspace(a, b, 40);
x = r*cos(t) + h;
y = r*sin(t) + k;
endI have written a code (attached) to plot particles striking a sphere of radius r – I can produce the image in example_1. I divide the shell into segments and calculate the number of particles falling within each.
I would like to represent the data as a smooth 3d heatmap over the spherical shell. I can produce the image shown in example_2 using scatter3 and a color bar, but I would like a smooth surface. I have tried to follow the example of my previous question using an isosurface –
however ux != uy != uz and I have the "Number of elements must not change" error.
I have also tried meshgrid(), but the size of photons results in an array that excedes maximum size preference.