matlab - 如何在 Matlab 中循环显示和绘制数据?
问题描述
我试图模拟两颗卫星的位置并确定它们是否会发生碰撞。为此,我使用了 Rk4 方法。如您所见,在循环中,我有一个drawnow
命令,它在每次迭代之后在 3d 图上添加一个点。我感兴趣的是找到一种方法来表示 3d 绘图,直到碰撞条件满足(我有消息的地方),并在.txt
文件中获取值,而不是在循环中。当我使用 return 命令并尝试绘制或添加文本时,我有一个错误“向量必须是相同的长度”并且不允许我获取数据。发生这种情况是因为让我们说 sat1=1X200 而 tspan=1X201。是否有任何其他命令可以用来代替 return 或更好的解决方案?
这是我的代码:
clear; close all; clc;
figure(1);
orbitsat1 = animatedline('Color','r');
msat1 = 124;
msat2 = 234;
Asat1 = 2.2;
Asat2 = 3.2;
sat10(1) = 453322.3616;
sat10(2) = -2346021.219;
sat10(3) = -7894131.349;
sat10(4) = 2142.38067;
sat10(5) = 7487.44895;
sat10(6) = -9864.00872;
sat10 = sat10';
sat20(1) = 543322.3616;
sat20(2) = -3246021.219;
sat20(3) = -8794131.349;
sat20(4) = 1242.38067;
sat20(5) = 4787.44895;
sat20(6) = -8964.00872;
sat20 = sat20';
tspan = 0:200;
secpday = 60*60*24;
a1 = 2018;
la1 = 1;
z1 = 2;
o1 = 12;
m1 = 0;
s1 = 0;
n1 = datenum(a1,la1,z1,o1,m1,s1);
n1 = n1 + tspan/secpday;
[an,luna,zi,ora,min,sec] = datevec(n1);
h = 1;
sat1 = zeros(6, tspan(end)/h);
sat1(:,1) = sat10;
sat2 = zeros(6, tspan(end)/h);
sat2(:,1) = sat20;
for i = 1:tspan(end)/h
k_1 = simsat1(tspan(i), sat1(:,i), msat1, Asat1, sat1(4:6,i));
k_2 = simsat1(tspan(i)+0.5*h, sat1(:,i)+0.5*h*k_1,msat1, Asat1, sat1(4:6,i));
k_3 = simsat1(tspan(i)+0.5*h, sat1(:,i)+0.5*h*k_2, msat1, Asat1, sat1(4:6,i));
k_4 = simsat1(tspan(i)+h, sat1(:,i)+k_3*h, msat1, Asat1, sat1(4:6,i));
k_12 = simsat2(tspan(i), sat2(:,i), msat2, Asat2, sat2(4:6,i));
k_22 = simsat2(tspan(i)+0.5*h, sat2(:,i)+0.5*h*k_12, msat2, Asat2, sat2(4:6,i));
k_32 = simsat2(tspan(i)+0.5*h, sat2(:,i)+0.5*h*k_22, msat2, Asat2, sat2(4:6,i));
k_42 = simsat2(tspan(i)+h, sat2(:,i)+k_32*h, msat2, Asat2, sat2(4:6,i));
sat2(:,i+1) = sat2(:,i) + (1/6)*(k_12+2*k_22+2*k_32+k_42)*h;
sat1(:,i+1) = sat1(:,i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;
addpoints(orbitsat1, sat1(1,i), sat1(2,i), sat1(3,i));
drawnow update;
hold off;
Rorb = sqrt(sum(sat1(1:3,i).^2));
Rcsc = sqrt(sum(sat2(1:3,i).^2));
dif = Rorb-Rcsc;
if (dif >0.5e10 && dif < 2e5) && i>350
Message = sprintf('Data:\nan:%d\nluna:%d\nzi:%d\nora:%d\nminut:%d\nsecunda:%d',an(i),luna(i),zi(i),ora(i),min(i),sec(i));
msgbox(Message)
msgbox('Coliziune!','Coliziune','warn',Message);
return
end
end
Pozx = sat1(1,:);
Pozy = sat1(2,:);
Pozz = sat1(3,:);
Tbody = [an luna zi ora min sec Pozx Pozy Pozz]';
twobody = fopen('nobody.txt','w');
fprintf(twobody,'Rezultate Twobody (metri) \n\n');
fprintf(twobody,' An luna zi ora min sec pozitia pe x pozitia pe y pozitia pe z viteza pe x viteza pe y viteza pe z\n\n');
fprintf(twobody,'---------------------------------------------------------------------------------------------------------------------------------------------\n%6d-%3d-%3d\t%3d:%3d:%3d\t\t\t\t%12.6f\t%12.6f\t\t%12.6f\t\t%12.6f\t%\n',Tbody);
fclose(twobody);
simsat1 在哪里
function sat1prim = simsat1(t,sat1,msat1,Asat1,vit)
miu = 398600.4418e9;
magn = sum(sat1(1:3).^2)^(3/2);
sat1prim = zeros(6,1);
sat1prim(1:3) = sat1(4:6);
sat1prim(4:6) = -miu.*sat1(1:3)./magn;
end
和 simsat2
function sat2prim = simsat2(t,sat2,msat2,Asat2,vit)
miu = 398600.4418e9;
magn = sum(sat2(1:3).^2)^(3/2);
sat2prim = zeros(6,1);
sat2prim(1:3) = sat2(4:6);
sat2prim(4:6) = -miu.*sat2(1:3)./magn;
end
解决方案
好的; 如果我理解正确,这应该是你想要的。
我已将Tbody
变量移到“for”循环中,并且在每个时间步发生时将其写入文件。此外,您搜索的关键字是break
,不是return
。break
只是让您脱离当前的 for 或 while 循环(仅上升一级),同时return
将控制权返回给调用函数。
我稍微改变了初始化时间向量的方式,使其更清晰。您对 "sat1=1X200 while tspan=1X201" 的问题来自这样一个事实,即您的tspan = 0:200
which 是201点,而tpsan(end)
value 仅为200。这很好并且很有意义,因为您将初始位置包含在向量中。基本上确保他们都有正确的点数。我通过定义一个 timestep dt=1
hour 然后 number of timesteps来做到这一点nt=200
。这样我的时间向量就变成tspan = 0:nt*dt
了,它会有点nt+1
。for i = 1:(nt-1)
然后我可以像你一样迭代。
我在整个过程中添加了注释,我鼓励您在代码中也这样做。反正; 这里是:
clear; close all; clc;
warning('on'); % Turn on warnings (Because they may be off)
figure(1);
orbitsat1 = animatedline('Color','r');
%% Initial conditions
msat1 = 124; msat2 = 234; Asat1 = 2.2; Asat2 = 3.2;
% Initialise satellites
sat10 = [453322.3616 -2346021.219 -7894131.349 2142.38067 7487.44895 -9864.00872];
sat20 = [543322.3616 -3246021.219 -8794131.349 1242.38067 4787.44895 -8964.00872];
dt = 1; % On hour timestep
nt = 200; % Number of timesteps
tspan = (0:nt)*dt; % Time from 0 to 200 timesteps in hours
% Make the time vectors
n1 = datenum(2018,1,2,12,0,0); % Starting point (serial date number in days)
n1 = n1 + tspan/24; % Timespan from hours -> days
[an,luna,zi,ora,min,sec] = datevec(n1);
%% Initialise containers
% Satelites
% nt+1 points because (:,1) is the initial condition
sat1 = zeros(6,nt+1); sat1(:,1) = sat10;
sat2 = zeros(6,nt+1); sat2(:,1) = sat20;
% Initialise Tbody here
Tbody = zeros(9,nt+1); Tbody(1:6,:) = [an; luna; zi; ora; min; sec];
%% Open the file
f_twobody = fopen('nobody.txt','w');
fprintf2(f_twobody,'Rezultate Twobody (metri)\n');
fprintf2(f_twobody,'Numar\tAn\tluna\tzi\tora\tmin\tsec\tpozitia_pe_x\tpozitia_pe_y\tpozitia_pe_z\n');
fprintf2(f_twobody,'%06d\t%6d\t%3d\t%3d\t%3d\t%3d\t%3d\t%12.6f\t%12.6f\t%12.6f\n',0,Tbody(:,1));
%% Go though all the timesteps
for i = 1:(nt-1)
% Calculate timestep
k_1 = simsat1(tspan(i), sat1(:,i), msat1, Asat1, sat1(4:6,i));
k_2 = simsat1(tspan(i)+0.5*dt, sat1(:,i)+0.5*dt*k_1,msat1, Asat1, sat1(4:6,i));
k_3 = simsat1(tspan(i)+0.5*dt, sat1(:,i)+0.5*dt*k_2, msat1, Asat1, sat1(4:6,i));
k_4 = simsat1(tspan(i)+dt, sat1(:,i)+k_3*dt, msat1, Asat1, sat1(4:6,i));
k_12 = simsat2(tspan(i), sat2(:,i), msat2, Asat2, sat2(4:6,i));
k_22 = simsat2(tspan(i)+0.5*dt, sat2(:,i)+0.5*dt*k_12, msat2, Asat2, sat2(4:6,i));
k_32 = simsat2(tspan(i)+0.5*dt, sat2(:,i)+0.5*dt*k_22, msat2, Asat2, sat2(4:6,i));
k_42 = simsat2(tspan(i)+dt, sat2(:,i)+k_32*dt, msat2, Asat2, sat2(4:6,i));
sat2(:,i+1) = sat2(:,i) + (1/6)*(k_12+2*k_22+2*k_32+k_42)*dt;
sat1(:,i+1) = sat1(:,i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*dt;
% Put the calculated positions into Tbody
Tbody(7:9,i+1) = [sat1(1,i+1) sat1(2,i+1) sat1(3,i+1)]';
% Write next line to the file
fprintf2(f_twobody,'%06d\t%6d\t%3d\t%3d\t%3d\t%3d\t%3d\t%12.6f\t%12.6f\t%12.6f\n',i,Tbody(:,i+1));
% Update plot
addpoints(orbitsat1, sat1(1,i), sat1(2,i), sat1(3,i));
drawnow update; hold off;
% Check for collision
Rorb = sqrt(sum(sat1(1:3,i+1).^2));
Rcsc = sqrt(sum(sat2(1:3,i+1).^2));
dif = abs(Rorb-Rcsc);
%if (dif >0.5e10 && dif < 2e5) && i>350
if sat1(1,i)>sat2(1,i) % Some kindof collision condition
warning('Coliziune! %d/%d/%d %d:%d:%d',Tbody(1:6,i+1));
fprintf(f_twobody,'Coliziune!\n');
break % Break out of the for loop (use this instead of 'return')
end
end
fclose(f_twobody); % Close the file
function sat1prim = simsat1(t,sat1,msat1,Asat1,vit)
miu = 398600.4418e9;
magn = sum(sat1(1:3).^2)^(3/2);
sat1prim = zeros(6,1);
sat1prim(1:3) = sat1(4:6);
sat1prim(4:6) = -miu.*sat1(1:3)./magn;
end
function sat2prim = simsat2(t,sat2,msat2,Asat2,vit)
miu = 398600.4418e9;
magn = sum(sat2(1:3).^2)^(3/2);
sat2prim = zeros(6,1);
sat2prim(1:3) = sat2(4:6);
sat2prim(4:6) = -miu.*sat2(1:3)./magn;
end
function fprintf2(f,varargin)
% Write both to stdout (console) and file
fprintf(varargin{:}); fprintf(f,varargin{:});
end
推荐阅读
- windows - 为什么在 cygwin 运行时,shell 中的字符串比较不起作用
- java - 从特定模式中拆分字符串
- firebase - 只有一个帐户的 Firebase 身份验证
- xml - 如何有条件地在 xpath 中包含元素?
- dataframe - 如何在 PySpark 中进行范围查找和搜索
- x86 - Does ASLR affect the maps file?
- python - GeoDjango Distance() Annotation on Related Model
- cuda - 我应该让 cublas 处理全局并在不同的(主机)函数中重用它们吗?
- javascript - need something like event.target.focus() in React
- android - Force Token Refresh in Xamarin Forms Android, Cloud Messaging