R=[-2436.45; -2436.45; 6891.0379];
V=[5.088611; -5.088611; 0];
x0=[R;V];
x0=x0*1000;
dt=1;
tVector=linspace(0,86400,2);
xfinal=rk4_2(x0,dt,tVector);
format longG
disp(xfinal(end,:))
%random noise
sz=size(xfinal(:,1:3));
error= 1/sqrt(3);
random=randn(sz)*error;
stateserror(:,1:3)=random+xfinal(:,1:3);
stateserror(:,4:6)=xfinal(:,4:6);
whos stateserror
stateserror is 2 x 6
%% Gibbs and Herrick loop
i=10;
Notice that i starts at 10
k=1;
tmax=20*60;
while i<=tmax/2
r1=stateserror(1,1:3);
That is okay because you are accessing fixed row 1 (out of 2)
r2=stateserror(i+1,1:3);
but i is 10 so that is asking for row 11 (out of 2)
r3=stateserror(i+2+1,1:3);
and that would ask for row 13 (out of 2)
r1magnitude=(r1);
r2magnitude=(r2);
r3magnitude=(r3);
z12=cross(r1,r2);
z31=cross(r3,r1);
z23=cross(r2,r3);
%Angels
angle12=acosd(dot(r1,r2)/(norm(r1*norm(r2))));
angle23=acosd(dot(r2,r3)/(norm(r2*norm(r3))));
angle13=acosd(dot(r1,r3)/(norm(r1*norm(r3))));
%Gibbs
N=r1magnitude.*z23+r2magnitude.*z21+r3magnitude.z12;
Nmag=norm(B);
D=z12+z23+z31;
Dmag=norm(D);
S=(r2magnitude-r2magnitude).*r1+(r3magnitude-r1magnitude).*r2+(r1magnitude-r2magnitude).*r3;
B=cross(D,r2);
Lg=sqrt(mu/(Nmag*Dmag));
v2G=Lg.*B/r2m+Lg.*S;
v2Gm=norm(v2gG);
statesm=norm(states(i,4:6));
Gibb_error=abs((V2Gm-statesm)/statesm)*100;
%Herrick
t31=2*i;
t32=i;
t21=i;
v2HG=-t32*((1/(t21*t31))+(mu/(12*r1magnitude^3)).*r1+(t32-t21)*(1/(t21*t32))+(mu/(12*r2magnitude^3)).*r2+t21*(1/(t32*t31))+(mu/(12*r3magnitude^3)).*r3);
v2HGm=norm(v2HG);
HGibb_error=abs((v2HGm-statesm)/statesm)*100;
angle(k,1)=angle13;
k=k+1;
i=i+10;
end
and nothing in your code grows stateserror, so it has to be pre-created to be 3 more than the largest possible i value you can get (so, 20*60/2 + 3 == 603)
Perhaps you need a longer time vector?
