r=rMaxin steps of
drand display the total number per bin over the total radius.
Conceptually, this is fairly simple. Fortunately, this is also a fairly easy program too. If we have a radius of 15 units (for the galaxy I was modeling, this would be 15 kiloparsecs), the easiest thing to do would make an integer array with 15 cells so that the bin size is 1 unit. Then you could just take the floor of the radial position and get the distribution quickly:
histare appropriately defined as integers.
Alternatively, since Fortran returns integer values for integer math, dividing by 1 can accomplish the same thing:
Upon testing these methods on a the same set of 10 million real values,
floortakes about 3.9 seconds compared to 3.8 for dividing by 1. Similar results are found for using
ceilingand dividing by 1 and adding 1.
If we want smaller bins, say
dr=0.5, we might think we would run into the problem of determining which half-unit this point belongs to in order to add it to the correct bin--this would basically involve something like finding the remainder of
real(floor(radius(i)))-radius(i)and seeing if it is greater or lesser than
0.5and then putting it into the appropriate bin.
However, there is far more simple method of doing this that requires us to first do two things:
(1) Set the histogram lower index to 0 and the upper index to number of bins - 1 (i.e.,
(2) Multiply the radius by 2:
allocate(hist(0:31)) do i=1,iMax
radius(i)=0.75, multiplying this by 2 gives 1, which is the correct bin. If
radius(i)=4.25, then multiplying by 2 gives 8, which is also the correct bin. Note that we could have left the indices as they were and just add 1 to the multiplication of 2, but this method adds a needless (and potentially costly) extra mathematical operation.
So, the next time you have to create a histogram of your data, recall this easy-binning method, which can easily be extended to larger or smaller bin sizes.
Comments always welcome!