Путаница

Question

Я пытаюсь создать цикл while, который вычисляет тормозную способность и диапазон тяжелой заряженной частицы, используя уравнение Бете-Блоха. У меня проблемы с настройкой цикла для сохранения каждого значения во время итерации цикла.

Идея состоит в том, что у меня есть уравнение «Остановочная сила» (-dEdx), которое выплевывает значение, сколько энергии теряется в течение dx, когда дано начальное значение энергии.

В приведенном здесь примере кода начальная энергия = 250 (МэВ), и я пытаюсь рассчитать новую энергию для каждого дх. Это должно выглядеть примерно так E = 250, 249,99, 249,96 и т. Д. Я не уверен на 100%, есть ли у цикла какие-либо другие проблемы, но я в основном не могу понять, как сохранить каждое отдельное значение E, чтобы я мог построить график того, что происходит с энергией на всем расстоянии, х.

Формула, по-видимому, выделяет только 1 значение для E и не заполняет все остальные столбцы.

Любая помощь будет принята с благодарностью!

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

Для тех, кто сталкивается с этим вопросом и рассматривает аналогичную проблему, вот обновленный код, который, кажется, отлично работает, благодаря Джону, который ответил на вопрос.

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

Answer 2

Ваша проблема в линии E(x(E)) = E(x(E)) + dE;

Когда вы запускаете свой цикл, E = 250 и в E есть только 1 элемент, поэтому x (E) - это некоторое число (250-й элемент x), а затем 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, и у вас будет встроенная индексная переменная.