Hi iontrap,
I understand that you want to divide the sphere into a number of segments and calculate the number of "photons" falling in each of those segments. Then, you wish to color the segments using a color map according to the number of "photons" in them and get a smooth sphere as output.
Please follow the steps given below to achieve this:
1.) From your code, it seems you have calculated the x, y, and z coordinates of "photons" (stored in the "prop" matrix) from the "input.txt" file, and then you are converting them to spherical coordinates in the for loop. This step can be performed using the cart2sph function, as shown in the code snippet below.
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);
end
x = prop(:,1); y = prop(:,2); z = prop(:,3);
% Convert Cartesian coordinates to spherical coordinates
[azimuth, elevation, ~] = cart2sph(x, y, z);
2.) Divide the sphere into a number of surface segments with equal intervals. This can be done by dividing the elevation and azimuth as shown in the code snippet below:
% Define the number of segments in azimuth and elevation
nAzimuthSegments = 72; % 5-degree intervals
nElevationSegments = 36; % 5-degree intervals
% Find the range of azimuth and elevation
azimuthRange = linspace(-pi, pi, nAzimuthSegments+1);
elevationRange = linspace(-pi/2, pi/2, nElevationSegments+1);
3.) Calculate the number of "photons" in each of those segments.
% Initialize the count matrix
counts = zeros(nElevationSegments, nAzimuthSegments);
% Count photons in each segment
for i = 1:nElevationSegments
for j = 1: nAzimuthSegments
counts(i, j) = sum(azimuth >= azimuthRange(j) & azimuth < azimuthRange(j+1) & ...
elevation >= elevationRange(i) & elevation < elevationRange(i+1));
end
end
4.) Create a meshgrid to plot the surface for each segments. Convert this meshgrid to cartesian coordinates for plotting. Use sph2cart for this conversion.
% Generate a meshgrid for plotting that matches the dimensions of 'counts'
[theta, phi] = meshgrid(linspace(-pi, pi, nAzimuthSegments), linspace(-pi/2, pi/2, nElevationSegments));
% Convert to Cartesian coordinates for plotting
[xPlot, yPlot, zPlot] = sph2cart(theta, phi, r);
5.) Get the RGB triplets to color the segments according to the number of "photons" stored in "counts" variable for each segment.
% Normalize 'counts' to the range [0, 1] for colormap mapping
countsNormalized = counts / max(counts(:));
N = 256;
colormapMatrix = parula(N);
ind = round(countsNormalized * (N-1)) + 1;
% Initialize the RGB matrix
rgbMatrix = zeros(size(counts, 1), size(counts, 2), 3);
% Fill the RGB matrix using the colormap
for i = 1:size(counts, 1)
for j = 1:size(counts, 2)
rgbMatrix(i, j, :) = colormapMatrix(ind(i, j), :);
end
end
6.) Plot surface plot for each segment with the respective color. This can be done using the surf function as shown in the code snippet below.
figure;
surf(xPlot, yPlot, zPlot, rgbMatrix, 'EdgeColor', 'none');
colormap('parula'); % Use a heat map color scheme
colorbar("TickLabels", linspace(0,250,11 ));
axis equal;
title('Photon Distribution on Sphere');
Please refer to the following documentation links to learn more about cart2sph, sph2cart and surf functions.
- https://www.mathworks.com/help/matlab/ref/cart2sph.html
- https://www.mathworks.com/help/matlab/ref/sph2cart.html
- https://www.mathworks.com/help/matlab/ref/surf.html
Hope this helps!
