gnuplot 直方图箱除以体积

Question

我正在模拟半径为 1 的球体体积中的点。我在该体积中生成了 1.000.000 个基于蒙特卡罗的点。为了制作 gnuplot 直方图，我计算了每个向量的长度（每个向量长度在 0 和 1 之间）。使用 100 个 bin，直方图如下所示： gnuplot data histogram.
如果有人想知道为什么没有生成大于 0.91 的点，我也不知道，但这不是这里的问题。

这是我的 gnuplot 代码：

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"

一般来说：体积增加 r^3 到半径 r。
在我的直方图中，每个球壳都是一个 bin，bin 编号是 100。因此，随着 bin 数量的增加，每个球壳的体积会呈立方增长（r^3）。从这个角度来看，直方图看起来不错。
但我想做的是绘制每体积点的密度：points/shellvolume。
这应该是从球体中心到其边界的线性分布。
我如何告诉 gnuplot 将每个 bin 除以其对应的体积，这取决于每个球壳的外半径和内半径？
公式为：(4/3)pi(R^3-r^3)，R 外径，r 内径为壳。

Answer 1

以下示例创建一些随机测试数据（应该是 20'000 个均匀分布的随机点）。 一种可能性是，您首先通过分箱创建直方图数据，然后将其除以外壳的体积。

顺便说一句，球壳的体积是(4./3)*pi*(R**3-r**3)，而不是你给出的公式。你为什么要设置max < min？也许您想根据您的确切需求微调分箱。

代码：

### 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

结果：