How to respresent and plot data while on loop in Matlab?

Question

What I'm trying to to is to simulate the position of two satellites and to determine if they will collide. For that I'm using a Rk4 method. As you can see, in the loop I have the drawnow command that adds a point on a 3d graph, after avery iteration. What I'm interested in is to find a way to represent the 3d plot until the condition of collision meets (where I have the message) and also to get the values on the .txt file, while in the loop not outside. When I use the return command, and try to plot or to add text, I have an error "vectors must be the same length" and doesn't let me obtain the data. That happens because lets say sat1=1X200 while tspan=1X201. Is there any other command that I could use instead of return or a better solution?

Here is my code:

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);

where 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

and 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

Answer 1

Allright; If I understand correctly this should be what you want.

I've moved the Tbody variable inside of the ‘for’ loop and I'm writing each timestep to the file as it happens. Also, the keyword you were searching for is break, not return. break just breaks you out of the current for or while loop (up only one level) while return returns the control to the invoking function.

I've changed slightly the way to initialise the time vector to make it more clear. Your problem with "sat1=1X200 while tspan=1X201" came from the fact that your tspan = 0:200 which is 201 points, while tpsan(end) value is only 200. This is fine and makes sense as you're including your initial position in your vectors. Basically make sure all of them have the right number of points. I've done this by defining a timestep dt=1 hour and then number of timesteps nt=200. This way my time vector becomes tspan = 0:nt*dt and it will have nt+1 points. I can then iterate for i = 1:(nt-1) just like you did.

I've added comments throughout which I'd encourage you to do too in your code. Anyway; here it is:

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