How can I create a circle of points where radius of a given point depends on a Gaussian distribution and point's value, & points are uniformly distributed around the circle?

Hello,

I am trying to create a circle of points where the radius of each point is based on a Gaussian distribution and the point's initial value (representing a photon's energy). The initial value is based on a probability (commented as "6 MV photon spectrum energy fluence fractions" below).

Instructions state: "Once you pick a certain photon energy from the given 6 MV spectrum, set its initial position at a radius 0 ≤ r ≤ 5cm according to a Gaussian distribution of photon energy as a function of radius, with the maximum on the beam axis. Use a Gaussian distribution with a height of 6 MeV, mean of 0, and a sigma of 1 cm." (Note: Beam axis is along z-axis, where x = 0 and y = 0.)

I tried to create a length x_photon such that it would be a random fraction of the radius r_circ (so, x_photon would fall randomly on the range [0, r_circ]), and then calculate y using the Pythagorean theorem. However, this has resulted in no photons having a position with y = 0, as visible on the scatter plot.

(This is a part of a larger script, but I think that everything below is what is relevant to this issue.)

Thanks!

clear

%%%%%%%%%%%%%%%%%% 6 MV photon spectrum energy fluence fractions
% 0.25 MeV
nmin_025 = 1./100000;
nmax_025 = 2480./100000;

% 0.5 MeV
nmin_050 = 2481./100000;
nmax_050 = 15000./100000;

% 0.75 MeV
nmin_075 = 15001./100000;
nmax_075 = 27290./100000;

% 1 MeV
nmin_100 = 27291./100000;
nmax_100 = 37590./100000;

% 1.25 MeV
nmin_125 = 37591./100000;
nmax_125 = 46310./100000;

% 1.5 MeV
nmin_150 = 46311./100000;
nmax_150 = 53760./100000;

% 1.75 MeV
nmin_175 = 53761./100000;
nmax_175 = 60140./100000;

% 2 MeV
nmin_200 = 60141./100000;
nmax_200 = 65680./100000;

% 2.25 MeV
nmin_225 = 65681./100000;
nmax_225 = 70460./100000;

% 2.5 MeV
nmin_250 = 70461./100000;
nmax_250 = 74630./100000;

% 2.75 MeV
nmin_275 = 74631./100000;
nmax_275 = 78290./100000;

% 3 Mev
nmin_300 = 78291./100000;
nmax_300 = 81510./100000;

% 3.25 MeV
nmin_325 = 81511./100000;
nmax_325 = 84330./100000;

% 3.5 MeV
nmin_350 = 84331./100000;
nmax_350 = 86860./100000;

% 3.75 MeV
nmin_375 = 86861./100000;
nmax_375 = 89090./100000;

% 4 MeV
nmin_400 = 89091./100000;
nmax_400 = 91060./100000;

% 4.25 MeV
nmin_425 = 91061./100000;
nmax_425 = 92790./100000;

% 4.5 MeV
nmin_450 = 92791./100000;
nmax_450 = 94330./100000;

% 4.75 MeV
nmin_475 = 94331./100000;
nmax_475 = 95670./100000;

% 5 MeV
nmin_500 = 95671./100000;
nmax_500 = 96840./100000;

% 5.25 MeV
nmin_525 = 96841./100000;
nmax_525 = 97850./100000;

% 5.5 MeV
nmin_550 = 97851./100000;
nmax_550 = 98710./100000;

% 5.75 MeV
nmin_575 = 98711./100000;
nmax_575 = 99420./100000;

% 6 MeV
nmin_600 = 99421./100000;
nmax_600 = 100000./100000;


%%%%%%%%%%%%%%%%%% pre-allocations
nhistories = 100000; % number of histories to run
x_photon = zeros(1,nhistories); % pre-allocate x position
y_photon = zeros(1,nhistories); % pre-allocate y position
z_photon = zeros(1,nhistories); % pre-allocate z position
E_photon = zeros(1,nhistories); % pre-allocate energies of photons for each history

%%%%%%%%%%%%%%%%%% Gaussian for choosing radius
a = 6; % height of Gaussian (MeV)
mu = 0; % mean
sig = 1; % std dev
r_circ = 0:0.001:5; % radius limits of circle on surface
gauss_r = a.*exp(-((r_circ - mu).^2)./(2.*(sig.^2))); % gaussian
% Note: I'm not certain whether r_circ should go from -5 to 5 or 0 to 5; 
% it doesn't seem to affect the result though.


%%%%%%%%%%%%%%%%%% For j loop, over histories
tic
for j = 1:nhistories % run j loop through histories
    
    % first get initial value (energy) based on probability
    n = rand().*nhistories;
    if n < nmin_050.*nhistories
            E_photon(j) = 0.25; % energy in MeV
        elseif (n > nmax_025.*nhistories) && (n < nmin_075.*nhistories)
            E_photon(j) = 0.50; % energy in MeV
        elseif (n > nmax_050.*nhistories) && (n < nmin_100.*nhistories)
            E_photon(j) = 0.75; % energy in MeV
        elseif (n > nmax_075.*nhistories) && (n < nmin_125.*nhistories)
            E_photon(j) = 1.00; % energy in MeV
        elseif (n > nmax_100.*nhistories) && (n < nmin_150.*nhistories)
            E_photon(j) = 1.25; % energy in MeV
        elseif (n > nmax_125.*nhistories) && (n < nmin_175.*nhistories)
            E_photon(j) = 1.50; % energy in MeV
        elseif (n > nmax_150.*nhistories) && (n < nmin_200.*nhistories)
            E_photon(j) = 1.75; % energy in MeV
        elseif (n > nmax_175.*nhistories) && (n < nmin_225.*nhistories)
            E_photon(j) = 2.00; % energy in MeV
        elseif (n > nmax_200.*nhistories) && (n < nmin_250.*nhistories)
            E_photon(j) = 2.25; % energy in MeV
        elseif (n > nmax_225.*nhistories) && (n < nmin_275.*nhistories)
            E_photon(j) = 2.50; % energy in MeV
        elseif (n > nmax_250.*nhistories) && (n < nmin_300.*nhistories)
            E_photon(j) = 2.75; % energy in MeV
        elseif (n > nmax_275.*nhistories) && (n < nmin_325.*nhistories)
            E_photon(j) = 3.00; % energy in MeV
        elseif (n > nmax_300.*nhistories) && (n < nmin_350.*nhistories)
            E_photon(j) = 3.25; % energy in MeV
        elseif (n > nmax_325.*nhistories) && (n < nmin_375.*nhistories)
            E_photon(j) = 3.50; % energy in MeV
        elseif (n > nmax_350.*nhistories) && (n < nmin_400.*nhistories)
            E_photon(j) = 3.75; % energy in MeV
        elseif (n > nmax_375.*nhistories) && (n < nmin_425.*nhistories)
            E_photon(j) = 4.00; % energy in MeV
        elseif (n > nmax_400.*nhistories) && (n < nmin_450.*nhistories)
            E_photon(j) = 4.25; % energy in MeV
        elseif (n > nmax_425.*nhistories) && (n < nmin_475.*nhistories)
            E_photon(j) = 4.50; % energy in MeV
        elseif (n > nmax_450.*nhistories) && (n < nmin_500.*nhistories)
            E_photon(j) = 4.75; % energy in MeV
        elseif (n > nmax_475.*nhistories) && (n < nmin_525.*nhistories)
            E_photon(j) = 5.00; % energy in MeV
        elseif (n > nmax_500.*nhistories) && (n < nmin_550.*nhistories)
            E_photon(j) = 5.25; % energy in MeV
        elseif (n > nmax_525.*nhistories) && (n < nmin_575.*nhistories)
            E_photon(j) = 5.50; % energy in MeV
        elseif (n > nmax_550.*nhistories) && (n < nmin_600.*nhistories)
            E_photon(j) = 5.75; % energy in MeV
        else 
            E_photon(j) = 6.00; % energy in MeV
    end

%%%%%%%%%%%%%%%%%% Get photon inital positions on phantom surface

    % get index of radius based on photon energy
    [~,idx]=min(abs(gauss_r - E_photon(j)));
    % (the reason gauss_r hass so many more values than E_photon is that
    % E_photon can change in a later part of the script that I have
    % excluded here)
    
    % randomly sample whether x and y values are positive or negative
    Rpmx = rand();
    if Rpmx >= 0.5
        pmx = 1;
    else
        pmx = -1;
    end

    Rpmy = rand();
    if Rpmy >= 0.5
        pmy = 1;
    else
        pmy = -1;
    end

    % randomly sample x from between 0 and r
    x_frac = randi([0, 4294967296])/4294967296.0; % what fraction of r_circ is the length x_photon
    x_photon(j) = pmx.*x_frac.*abs(r_circ(idx)); % (cm)

    % get y from pythagorean theorem
    y_photon(j) = pmy.*sqrt((r_circ(idx)).^2 - (x_photon(j)).^2); % (cm)

    % z = 0 cm at surface
    z_photon(j) = 0;

end % end for j over histories
position_time = toc

%%%%%%%%%%%%%%%%%% plot resulting positions of photons on phantom surface
scatter3(x_photon,y_photon,z_photon)
xlabel('x')
ylabel('y')

Answer 1

[Joan B]

I changed from randomly sampling a fraction for x to randomly sampling an angle. So I now have the lines:

theta_quarter = 0:0.01:90; % degrees

and

% get index of radius based on energy
    [~,idx]=min(abs(gauss_r - E_photon(j)));
    
    % randomly sample whether x and y values are positive or negative
    Rpmx = rand();
    if Rpmx >= 0.5
        pmx = 1;
    else
        pmx = -1;
    end

    Rpmy = rand();
    if Rpmy >= 0.5
        pmy = 1;
    else
        pmy = -1;
    end
    
    % choose a random angle between 0 and 90 degrees
    R_theta = ceil(rand().*length(theta_quarter));

    % get x and y using cos and sin of selected theta_quarter
    x_photon(j) = pmx.*r_circ(idx).*cosd(theta_quarter(R_theta));
    y_photon(j) = pmy.*r_circ(idx).*sind(theta_quarter(R_theta));