Creating a step function for multiple array indexes.

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;% N: Max number of site, L: Length of  box
for n = 1:N % generate nth site
    
    x(n) = rand*L; % generate a random x coordinate
    
    y(n) = rand*L; % generate a random y coordinate
    
    E = randraw( 'pareto',[1*10^19 3], N ); % energy
    
    % Count it if it is in the circle % Count it if it is in the circle
    d(n) = sqrt(x(n).^2+y(n).^2); % calculate the distance from the centre of the circle
    theta(n) = atan2(y(n),x(n)); % compute the angle made
    

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;
                end
            end
        end
        end
    end

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.

Answer 1

[Mrunmayee Gaikwad]

As ‘d’ is a 1xN array, the floor function in Line 3 of the if loop (steps = floor(d)) will operate on each element of ‘d’ array, creating 1xN ‘steps’ array. The next if statement compares the whole 1xN ‘steps’ array to ‘L/2’, which returns a 1xN logical array (containing ‘0’s and ‘1’s). The if statement will be executed only if all the elements in the logical array are 1s.

The same is true for the line (EN = E/2) in the if loop.

Assuming the if loop is used inside the main for loop, use steps = floor(d(n)); to get step value for the distance meeting the condition. Similarly, use EN = E(n)/2. Here, ‘steps’ and ‘EN’ will be scalar variables, which would make it easier for comparison.

[Nate Taylor]

Thank you!

I think this method will work but I ended up creating a seperate function:

function E = stepenergy(E0,d)
    global L
if d > L/2
   E = 0;
else  
    steps = floor(d);
        for i = 1:steps
            if rand < 0.1
                E0 = E0/2;               
            end
        end
E = E0;
end

Which I then call in an if loop:

 if d(i) <= L/2 % only apply the step function to specified distances
         sites = sites +1;% increase counter for number of potential sitesd 
         E(i) = stepenergy(E0(i),d(i)); % call step function for energy reduction
      else E(i) = 0; % ensure all other energy values are reduced to zero so not to be included
      end