循环混淆

Question

我正在尝试创建一个 while 循环，该循环使用 Bethe-Bloch 方程计算重带电粒子的停止能力和射程。我在设置循环以在循环迭代时保存每个值时遇到问题。

我的想法是，我有一个方程停止功率 (-dEdx)，当给定初始能量值时，它会输出一个值，表示在 dx 上损失了多少能量。

在这里的示例代码中，初始能量 = 250 (MeV)，我正在尝试计算每个 dx 的新能量。这应该看起来像 E = 250、249.99、249.96 等 我不是 100% 确定循环是否有任何其他问题，但我主要无法弄清楚如何保持 E 的每个单独值，以便我可以绘制能量在总距离 x 上的变化。

该公式似乎只为 E 吐出 1 个值，并且没有填充所有其他列。

任何帮助将不胜感激！

clear
clc

c = 2.998e10;
pm = 938.272; % proton mass
em = 0.510991; % electron mass
re = 2.818e-13; % classical electron radius
z = 1; % charge of proton (+)
na = 6.022e23; % avagadros number
rho = 1;
dx = 0.01;
x = 0:dx:350;    % 350mm = 35cm
n = na.*10./18.*rho;
E = x*0;
E(1) = 250;

while E > 0
    gam = E./pm + 1;
    beta = sqrt(gam.^2-1)./gam;
    dEdx = (4.*pi.*n.*z^2.*em.*re.^2)./beta.*(log(2.*em.*c.^2.*beta.^2./(75.*(1-beta.^2)))-beta.^2);
    dE = -dEdx.*dx;
    E(x(E)) = E(x(E)) + dE;
end

在每次迭代中保留 gam,beta,dEdx,dE,E 的所有值可能很有用。例如 E = [250, 249.99, 249.96, ..., 0]

提前致谢！

Answer 1

你的问题在E(x(E)) = E(x(E)) + dE;这一行

当您运行循环时，E=250 并且 E 中只有 1 个元素，因此 x(E) 是某个数字（x 的第 250 个元素），然后 E(x(E)) 是 E 的元素那个索引。通常，这会导致错误，因为 x(E) 在您的示例中不能保证是整数，而且它只是没有按照您希望的方式执行。为了实现你的目标，你需要一个索引变量：

index=1;
E(index)=250;

while E(index) > 0
    gam = E(index)/pm + 1;
    beta = sqrt(gam^2-1)/gam;
    dEdx = (4.*pi.*n.*z^2.*em.*re.^2)./beta.*(log(2.*em.*c.^2.*beta.^2./(75.*(1-beta.^2)))-beta.^2);
    dE = -dEdx*dx;
    E(index+1) = E(index) + dE;
    index=index+1;
end

您可能还想解决另一个不一致的问题。您提前设置 x 具有固定的最大值，但在计算 E 时您没有强制执行该最大值。如果您想使用 while 循环并确保获得 E 的完全衰减，您不应该提前定义 x时间，而是在你知道它需要多长时间之后在最后创建它。或者，如果您只想达到 x 的最大值，您可以使用 for 循环而不是 while 循环，并且您将内置索引变量。

Answer 2

对于任何偶然发现这个问题并查看类似问题的人，这里是更新后的代码，感谢回答问题的 John，看起来运行良好。

clear
clc

c = 2.998e10;
pm = 938.272; % proton mass
em = 0.510991; % electron mass
re = 2.818e-13; % classical electron radius
z = 1; % charge of proton (+)
na = 6.022e23; % avagadros number
rho = 1;
dx = 0.01;
x = 0:dx:350;    % Expected Maximum range: 350mm = 35cm
n = na.*10./18.*rho;

dEdx = x*0;
dE = x*0;
E = x*0;

i = 1;
E(i) = 250;

while E(i) > 0
    gam = E./pm + 1;
    beta = sqrt(gam.^2-1)./gam;
    dEdx = ((4.*pi.*n.*z^2.*em.*re.^2)./beta).*(log(2.*em.*c.^2.*beta.^2./...
    (75.*(1-beta.^2)))-beta.^2);

    dE(i) = -dEdx(i).*dx;
    E(i+1) = E(i) + dE(i);
    i = i + 1;
end