gnuplot histogram bins divided by volume

Question

I am simulating points in a sphere volume with radius 1. I generated 1.000.000 monte-carlo based points in this volume. To make a gnuplot histogram i calculated the length of each vector (every vector length is between 0 and 1). With 100 bins the histogram looks like: gnuplot data histogram.
If someone is wondering why there no points greater than 0.91 are generated, i also dont know, but this is not the question here.

This is my gnuplot Code:

n=100 #number of intervals
max=1.0 #max value
min=0.0 #min value
width=(max-min)/n #interval width

#function used to map a value to the intervals
hist(x,width)=width*floor(x/width)+width/2.0

#settings
set xlabel "Radius"
set ylabel "Primarys/Intervall"
set xrange [-0.1:1.1] 
set yrange [0:32000]
set boxwidth width*0.8
set style fill solid 0.5 #fillstyle
set tics out nomirror

#plot
plot "primaryPosition(1).csv" u (hist($1,width)):(1.0) smooth freq w boxes lc rgb"green"

In general: A Volume grows by r^3 to Radius r.
In my histrogram every spherical shell is one bin and the bin number is 100. So, as the bin number increases, the volume of each sperical shell grows cubically (with r^3). From this point of view, the histogram looks good.
But what i want to do is to plot the density of points per volume: points/shellvolume.
This should be a linear distribution from the center of the sphere to its border.
How can i tell gnuplot to divide each bin by its corresponding volume, which depends on the outer and the inner radius of each spherical shell?
The formula is: (4/3)pi(R^3-r^3) with R outer and r inner radius a shell.

theozh theozh · Accepted Answer · 2020-07-07T06:16:12

The following example creates some random test data (should be 20'000 equally distributed random points). One possibility would be that you first you create your histogram data via binning into a table and then you divide it by the volume of the shell.

By the way, the volume of a sphere shell is (4./3)*pi*(R**3-r**3), not the formula you've given. And why are you setting max < min? Maybe you want to fine tune the binning to your exact needs.

Code:

### histogram normalized by sphere shell volume
reset session
set view equal xyz

# create some test data
set print $Data
    do for [i=1:20000] {
        x = rand(0)*2-1
        y = rand(0)*2-1
        z = rand(0)*2-1
        r = sqrt(x**2 + y**2 + z**2)
        if (r <= 1) { print sprintf("%g %g %g %g",x,y,z,r) }
    }
set print

n = 100                 # number of intervals
min = 0.0               # max value
max = 1.0               # min value
myWidth=(max-min)/n     # interval width
bin(x)=myWidth*floor(x/myWidth)

ShellVolume(r) = (4./3)*pi*((r+myWidth)**3-r**3)
set boxwidth myWidth absolute

set table $Histo
    plot $Data u (bin($4)):(1) smooth freq
unset table

set multiplot layout 2,1
    plot $Histo u 1:2 w boxes ti "Occurrences"

    plot $Histo u 1:($2/ShellVolume($1)) w boxes ti "Density"
unset multiplot
### end of code

Result:

gnuplot histogram bins divided by volume

1 Answers