While-Loop混淆

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中只有一个元素，所以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循环，并内置index变量。

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