1

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
Improve this question
7
  • If they're only 1 element different, why not just exclude the final tspan value for plotting? Commented Aug 13, 2018 at 16:18
  • That would work only in case of a collision , but if the satellites won t collide over a period of time that wouldn't be correct. Thats because sat1=201( in the case without collision) and tspan=201; you propouse tspan=tspan-1; and I would get the same error. Commented Aug 13, 2018 at 16:24
  • Your code example does not work. It's missing secundePeZi, msat1, Asat1 and then it throws an error when indexing simsat1. There may be error afterwards but I gave up trying by this point. Do the following; add close all; clear; clc; line at the very top of your code. Then make sure it runs and update your question. I'm guessing this code worked for you because you had all of these variables already in your workspace. To avoid this type of 'crosstalk' between the scripts I'd recommend always putting the close all; clear; clc; at the beginning of every piece of code. Commented Aug 13, 2018 at 17:36
  • You are wright. I've just uptade the code with A,m, and the functions. Commented Aug 13, 2018 at 17:46
  • Three things: 1. I don't understand your collision condition: dif >0.5*10^10 && dif < 2*10^5) && i>350. How can dif be simultaneously bigger than 5e9 and smaller than 2e5 ?? Are you sure you've got your <> signs correct. Also; are you sure that on the line before it shouldn't say dif = abs(Rorb-Rcsc);? 2. Your code as supplied does not actually produce the collision you speak of (nor the error on that matter). 3. Your functions take 5 input arguments but use only one - the sat1 and sat2's. The others are apparently superfluous but are used and varied in your code...(?) Commented Aug 13, 2018 at 23:37

1 Answer 1

Reset to default
2

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
Improve this answer
Sign up to request clarification or add additional context in comments.

2 Comments

That works very good. But what if I want to integrate for every second, assuming Tspan=200 seconds?Also, I see in the .txt file that sec(second) are from are sec+24s, and that it starts from 0 at the first timestep.Is there a way to put instead of 0 at first step the values for init cond and the seconds sec=sec+1?
I'm not sure whether your algorithm will work for that, but in principle you could just set dt = 1/3600 hours (i.e. 1s). As to the 'sec' entry in the file, it is what's generated by datevec. The first entry in the file it literally the initial conditions (there is an fprintf putting it there just before the loop). I've also noticed there were some errors in my code (I've assumed datevec has units of seconds but it in fact has units of days. Also brackets are hard...) so I've corrected that (see code above). Give it a go now. The entries in .txt should make more sense.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.