首页 > 解决方案 > 索引超出数组范围

问题描述

我正在尝试以多种方式解决此问题,但仍然存在标题中的错误。

我想提一下,此代码是为生成光子包并计算下一步、偏转和在此期间的损失(以能量计)而创建的。

发射错误(第 58 行)绘图(X(i,1:r(i)),Y(i,1:r(i)),'-');

这是代码:

close all;
clc;

%% parameters
Nph=100;
mua=99.9;
mus=0.01;
g=0.01;

R=10000; %resolution
t=0:1/R:1;
n=1:1:Nph;

%% distributions
dS=zeros(size(t));
theta=zeros(size(t));
for i=1:1:R+1
    dS(i)=step(t(i),mua,mus);
    theta(i)=deflection(t(i),g);
end

%% throw Nph photons with no energy losses
Xn=zeros([Nph R+2]);
Yn=zeros([Nph R+2]);
for i=1:1:Nph
    phi=rand*2*pi;
    [Xn(i,:),Yn(i,:)]=throwPhotonNoLoss(phi,R,mua,mus,g);
end

% throw Nph photons with energy losses
X=zeros([Nph R+2]);
Y=zeros([Nph R+2]);
r=zeros([1 Nph]);
W=zeros([1 Nph]);
for i=1:1:Nph
    phi=rand*2*pi;
    [X,Y,r,W]=throwPhoton(phi,R,mua,mus,g);
end

%% plot
%figure(1), plot(t,dS,'-'), title('step');
%figure(2), plot(t,theta,'-'), title('deflection');

str='photons='+string(Nph)+', R='+string(R)+', mua='+string(mua)+', mus='+string(mus)+', g='+string(g);

figure(3), hold on, title(['Paths (No Loss): ', str]);
for i=1:1:Nph
    plot(Xn(i,:),Yn(i,:),'-');
end

figure(4), hold on, title(['Final positions (No Loss): ', str]);
for i=1:1:Nph
    plot(Xn(i,R+2),Yn(i,R+2),'.');
end

figure(5), hold on, title(['Paths: ', str]);
for i=1:1:Nph
    plot(X(i,1:r(i)),Y(i,1:r(i)),'-');
end

figure(6), hold on, title(['Final positions: ', str]);
for i=1:1:Nph
    plot(X(i,r(i)),Y(i,r(i)),'.');
end

figure(7), plot(n,r,'+'), title('r');
figure(8), plot(n,W,'+'), title('W');

功能代码:deflection.m:

function th=deflection(e,g)

num=1-g*g;
den=1-g+2*g*e;
f=1/(2*g);
th=f*(1+g*g-num/den);

end

ThrowPhoton.m

function [x,y,r,W]=throwPhoton(phi,R,mua,mus,g)

W=1.0;
a=mua/(mua+mus);
% dS=zeros([1 R+1]);
% theta=zeros([1 R+1]);
% x=zeros([1 R+2]);
% y=zeros([1 R+2]);

x(1)=0;
y(1)=0;
i=1;
while W>=0.1 %&& i~=R+1
    e1=rand;
    e2=rand;
    dS(i)=step(e1,mua,mus);
    theta(i)=deflection(e2,g);
    x(i+1)=x(i)+dS(i)*cos(theta(i)+phi);
    y(i+1)=y(i)+dS(i)*sin(theta(i)+phi);
    W=W*a;
    i=i+1;
end
r=i;
W;

end

ThrowPhotonNoLoss.m

function [x,y]=throwPhotonNoLoss(phi,R,mua,mus,g)

dS=zeros([1 R+1]);
theta=zeros([1 R+1]);
x=zeros([1 R+2]);
y=zeros([1 R+2]);

x(1)=0;
y(1)=0;
for i=1:1:R+1
    e1=rand;
    e2=rand;
    dS(i)=step(e1,mua,mus);
    theta(i)=deflection(e2,g);
    x(i+1)=x(i)+dS(i)*cos(theta(i)+phi);
    y(i+1)=y(i)+dS(i)*sin(theta(i)+phi);
end

end

步骤.m

function dS=step(e,mua,mus)

if e==0
    e=0.01;
end
dS=-log(e)/(mua+mus);
end

预先感谢您的回复

标签: matlab

解决方案


推荐阅读