Hi I really hope someone can help me!
I am creating a 2D model for a university project to approximate the density of intergalatic cosmic ray sources within 100 Mpc from Earth. I plan to repeat the model and average the results. The cosmic ray sources I want to investigate have arrival energies above 10^19 eV. The model shall assume that starting energies are randomly distributed with an intensity following a power-law disribution with α = 3.1 and that as the CR propogates toward Earth (other directions are not considered), for every 10Mpc travelled there is a 1/10 chance that the energy of that CR will be reduced by 50%. The goal of the model is to determine the average number of sources of CR from the population, thier angles, energy density, and average density.
To begin, I generate random position vectors within a square of length L = 200, generate an array of random energies using randraw power-distribution function, and compute the vector distance to the centre of the square, as well as angle (note: N will be much larger once I have the model working correclty):
N= 2; L= 200; counts= 0; source = 0; sites=0;
for n = 1:N
x(n) = rand*L;
y(n) = rand*L;
E = randraw( 'pareto',[1*10^19 3], N );
d(n) = sqrt(x(n).^2+y(n).^2);
theta(n) = atan2(y(n),x(n));
Next, I wish to only consider values with a distance d(n) <= 100 which I also wish to count, so I create an if loop and create a counter of values called 'sites'. Inside my if loop I want to take each value of d(n) that meets my distance requirement, and step the energy reduction probability for every 1Mpc. As the energy reduction is based on distance, each entry of d(n) that meets the distance condition will have a unique integar value for the number of steps needed to be taken for the probability energy reduction determined by taking the floor of that distance value. Now, I need to apply the unique number of steps to each distance value whereby at each step a random number p is generated such that if p < 0.1 the energy related to that corresponding distance value is reduced by 50%. Each step repeats the same process of assigning a random number and potentially reducing the energy value. The resulting energy at the completion of the loop is the "arrival energy" mentioned above. Note: I have not yet attempted to code for tasks beyond this point. I then need to identify how many sites arrive with energies above 10^19 eV and count them (I plan to do this using a simple if loop and new counter called "source"), as well as sum their corresonding angles, compute and plot the resultant amplitude and phase for all iterations.
My attempt at the if loop:
if d(n) <= L/2
sites = sites+1;
steps = floor(d);
if steps <= L/2
for i = 1:steps
p(i,n) = rand;
if p(i,n) < 0.1
EN = E/2;
if EN >= 10^19
source = source+1;
Problem 1: I would like to call the unique energy values for only those distance values that meet my distance condition. I have tested this using a low N value and guaranteed values for d by generating x and y values randomly but multiplied by L/2 rather than L. In this way every energy value calulated will be reduced if the probability function meets the condition of p <0.1 so I need a way to exclude the energy values associated with those distances that do not meet the condition, or conversely, only inlcude those energy values associated to distance values that do meet the condition.
Problem 2: I need to use the unique step value for each unique distance value. My current code creates an array called "steps" which does successfully attribute and integar value to each distance, however, it calculates the steps for all distance values not only those that meet the distance requirement. Furthermore, I can't tell if each number of steps is being applied to the correct index (if at all).
I will be very grateful for any help.