Zeros array outputting the last data set from loop instead of all iterations

%% Problem Inputs

M0 = 0.558505;

a = 26562;

e = 0.74105;

i = 63.4;

w = 270 ;

Om = 260;

tol = 10^-8;

mu = 398600;

P = (2*pi/sqrt(mu))*a^(3/2);

h= sqrt(mu*a*(1-e^2));

%% Calculations

% Time vector

t0 = 0;

t = linspace(t0, P, 1000);

% Initializating variables for storage

M = zeros(1, length(t));

E = zeros(1, length(t));

theta = zeros(1, length(t));

rECI = zeros(3, length(t));

vECI = zeros(3, length(t));

X = zeros(6, length(t));

R3w = [cosd(w) -sind(w) 0; sind(w) cosd(w) 0; 0 0 1]; %rotation matrix

R1 = [1 0 0; 0 cosd(i) -sind(i); 0 sind(i) cosd(i)]; %roation matrix

R3Om = [cosd(Om) -sind(Om) 0; sind(Om) cosd(Om) 0; 0 0 1]; %rotation matrix

for idx = 1:1:length(t)

time = t(idx);

M(idx) = M0 + (2*pi/P)*(time - t0);

E(idx) = kepler(M(idx), e, tol);

beta = e/(1 + sqrt(1-e*e));

theta(idx) = E(idx) + 2*atan((beta*sin(E(idx)))/(1-beta*cos(E(idx))));

% theta(idx) = 2*atan( sqrt((1+e)/(1-e)) * tan(E(idx)/2) );

%%conversion of Orbital elements to ECI

% Formula for perifocal r vector s

r_p = (h^2/mu)*1/1+e*cosd(theta(idx));

rx = r_p*cosd(theta(idx)); % perifocal x coordinate

ry = r_p*sind(theta(idx)); % perifocal y coordinate

r = [rx; ry; 0]; % Perifocal r coordinates

v_p = mu/h; % velocity vector formula perifocal

vx= v_p*(-1*sind(theta(idx))); % Velocity x coordinate

vy = v_p*(e+cosd(theta(idx))); % velocity y coordinate

v = [vx; vy; 0]; %velocity vector perifocal

rECI = R3Om*R1*R3w*r; % multiplying r perifocal vector by rotation matrices

vECI = R3Om*R1*R3w*v; % multiplying v perifocal vector by rotation matrices

X = [rECI; vECI];

end

figure(10);clf;

plot3(rECI(1,:),rECI(2,:),rECI(3,:), '.');

hold on;

quiver3(rECI(1,1),rECI(2,1),rECI(3,1),vECI(1,1),vECI(2,1),vECI(3,1),1, 'r');

keyboard;

%% Function Definition

function [E] = kepler(M,e,tol)

E = M; %Initializes the first guess as the value for the mean anomoly

f=@(E) M-E+e*sin(E); %keplers equation rearranged to equal 0

fdiff=@(E) -1+e*cos(E); %First differential of keplers equation

E = E - f(E)/fdiff(E); %newtons methods

y = f(E);

while abs(y) > tol

E = E - y/fdiff(E);

y = f(E);

end

end

So the problem here is with the rECI and vECI values, it is not saving the values for each iteration and only saves the last one. This is causing my plot to be a dot with a vector arrow. Here is a picture of where the problem is occuring

Answer 1

[VBBV]

rECI(:,idx) = R3Om*R1*R3w*r; % multiplying r perifocal vector by rotation matrices
vECI(:,idx) = R3Om*R1*R3w*v; % multiplying v perifocal vector by rotation matrices        

[VBBV]

Use the for loop index

[Rowan McLoughlin]

thanks, saved my desk from further abuse