Planet Orbit- Dormand Prince Algorithm

5 views (last 30 days)
Danielle Curtin on 1 Dec 2019
Answered: Vimal Rathod on 5 Dec 2019
function PlanetOrbit
[S,M,AU,C1,C2,R]=StateInitialization;
dt=3600;
t=0;
DT=datenum(2019,12,1,0,0,0);
SF=1.157407407407407e-05; %keeps track of date
plotState(S,AU,C1,C2,R,DT);
eps=0.000001; %error allowance in one step calculation.
t0=0; %initiallization of variables.
y0=0;
h0=0.1;
while(t0<=tf) %tf is the ending time of the calculation.
k1=h*f(t0,y0);
k2=h*f(t0 + (1/5)*h, y0 + (1/5)*k1);
k3=h*f(t0 + (3/10)*h, y0+(3/40)*k1+(9/40)*k4);
k4=h*f(t0+(4/5)*h, y0 + (44/45)*k1-(56/12)*k2+(32/9)*k3);
k5=h*f(t0+(8/9)*k1, y0+(19372/6561)*k1-(25360/2187)*k2+(64448/6561)*k3-(212/729)*k4);
k6=h*f(t0+h, y0+(9017/3168)*k1-(355/33)*k2-(46732/5247)*k3+(49/176)*k4-(5103/18656)*k5);
k7=h*f(t0+h,y0+(35/384)*k1+(500/1113)*k3+(125/192)*k4-(2187/6784)*k5+(11/84)*k6);
y1 = y0 + 35/384*k1 + 500/113*k3 + 125/192*k4 - 2187/6784*k5 +11/84*k6;
z1 = y0 + 5179/57600*k1 + (7571/16695)*k3+(393/640)*k4-(92097/339200)*k5+(187/2100)*k6+(1/40)*k7; %to estimate the error.
err = abs(z1-y1); %error estimation
s=pow(eps*h0/(2*err),1/5);
h1=s*h0; %optimal time interval. used in the next step.
if(h1<hmin)
h1=hmin; %confine time interval between hmin and hmax.
else if(h1>hmax)h1=hmax;
t0=t0+h0;
y0=y1;
h0=h1;
end
end
end
for j=1:50000
S=S+dt*physics(S,M);
t=t+dt;
DT=DT+SF*dt;
if mod(j,10)==0
plotState(S,AU,C1,C2,R,DT);
end
end
end
function plotState(S,AU,C1,C2,R,DT)
NB=round(length(S)/6);
clf;
hold on;
for j=1:NB
bi=6*(j-1)+1;
plot3(S(bi)/AU,S(bi+1)/AU,S(bi+2)/AU,C1(j,:),'MarkerSize',R(j),'MarkerFaceColor',C2(j));
end
axis([-2,2,-2,2,-2,2])
title(datestr(round(DT,4)))
drawnow;
end
I get the error message:
>> PlanetOrbit
Index exceeds matrix dimensions.
Error in PlanetOrbit>plotState (line 59)
plot3(S(bi)/AU,S(bi+1)/AU,S(bi+2)/AU,C1(j,:),'MarkerSize',R(j),'MarkerFaceColor',C2(j));
Error in PlanetOrbit (line 12)
plotState(S,AU,C1,C2,R,DT);
The other functions used are:
function [S,M,AU,C1,C2,R]=StateInitialization
% S state of system (x,y,z,vx,vy,vz) meters and m/s long vector
% M list of masses (sun,earth, etc.)
% AU astronomical unit(meters) earth to sun
% C1 keeps track of outer circle, fill interior blue
% C2 outer color
% R how big dots appear
AU = 149597870700;
C1(1,:)='yo'; C2(1)='y';
C1(2,:)='bo'; C2(2)='b';
C1(3,:)='mo'; C2(3)='m';
C1(4,:)='go'; C2(4)='g';
C1(5,:)='ro'; C2(5)='r';
C1(6,:)='ko'; C2(6)='k';
C1(7,:)='co'; C2(7)='c';
C1(8,:)='wo'; C2(8)='w';
R=[12;8;8;5];
S=zeros(6,1);
%Sun
X=0; Y=0; Z=0;
VX=0; VY=0; VZ=0;
j=1; bi=6*(j-1)+1;
S(bi)=1000*X; S(bi+1)=1000*Y; S(bi+2)=1000*Z;
S(bi+3)=1000*VX; S(bi+4)=1000*VY; S(bi+5)=1000*VZ;
M(1)=1.989*1e30;
%Earth/Moon
X = 5.476756046031971E+07; Y = 1.369854390984182E+08; Z =-6.418494023710489E+03;
VX=-2.814628809603848E+01; VY= 1.094601049527264E+01; VZ=-4.343632416126120E-04;
j=j+1; bi=6*(j-1)+1;
S(bi)=1000*X; S(bi+1)=1000*Y; S(bi+2)=1000*Z;
S(bi+3)=1000*VX; S(bi+4)=1000*VY; S(bi+5)=1000*VZ;
M(2)=5.972*1e24+7.34767309*1e22;
%Venus
X = 7.690918093481298E+07; Y =-7.695356068929358E+07; Z =-5.494117392182905E+06;
VX= 2.454395327263311E+01; VY= 2.462723606903867E+01; VZ=-1.078441904000677E+00;
j=j+2; bi=6*(j-1)+1;
S(bi)=1000*X; S(bi+1)=1000*Y; S(bi+2)=1000*Z;
S(bi+3)=1000*VX; S(bi+4)=1000*VY; S(bi+5)=1000*VZ;
M(3)=4.867*1e24;
%Apophis
X = 1.006188535413713E+08; Y = 7.025826519631593E+07; Z =-1.349533159624398E+06;
VX=-1.511473674904992E+01; VY= 3.112188699288047E+01; VZ=-2.016175343620677E+00;
j=j+3; bi=6*(j-1)+1;
S(bi)=1000*X; S(bi+1)=1000*Y; S(bi+2)=1000*Z;
S(bi+3)=1000*VX; S(bi+4)=1000*VY; S(bi+5)=1000*VZ;
M(4)=26.99*1e9;
%Saturn
X = 5.454260468832656E+08; Y = -1.399004902083351E+09; Z = 2.611933568101764E+06;
VX= 8.479463517851789E+00; VY= 3.483787999150237E+00; VZ=-3.985176203231436E-01;
j=j+4; bi=6*(j-1)+1;
S(bi)=1000*X; S(bi+1)=1000*Y; S(bi+2)=1000*Z;
S(bi+3)=1000*VX; S(bi+4)=1000*VY; S(bi+5)=1000*VZ;
M(5) = 5.6834*1e26;
%Jupiter
X = 4.421300143790728E+07; Y =-7.824664674728345E+08; Z = 2.260798916549563E+06;
VX= 1.290169373776935E+01; VY= 1.355135395534917E+00; VZ=-2.942549827353074E-01;
j=j+5; bi=6*(j-1)+1;
S(bi)=1000*X; S(bi+1)=1000*Y; S(bi+2)=1000*Z;
S(bi+3)=1000*VX; S(bi+4)=1000*VY; S(bi+5)=1000*VZ;
M(6) = 1.89813*1e27;
%Neptune
X = 4.371602845443687E+09; Y =-9.668331843241746E+08; Z =-8.085109646501470E+07;
VX= 1.151818638737934E+00; VY= 5.342596132371595E+00; VZ=-1.370365801570574E-01;
j=j+6; bi=6*(j-1)+1;
S(bi)=1000*X; S(bi+1)=1000*Y; S(bi+2)=1000*Z;
S(bi+3)=1000*VX; S(bi+4)=1000*VY; S(bi+5)=1000*VZ;
M(7) = 1.02413*1e26;
%Uranus
X = 2.437829683225929E+09; Y = 1.688160852693148E+09; Z =-2.530390129291379E+07;
VX=-3.914948553525651E+00; VY= 5.282609057004567E+00; VZ= 7.012996969005503E-02;
j=j+7; bi=6*(j-1)+1;
S(bi)=1000*X; S(bi+1)=1000*Y; S(bi+2)=1000*Z;
S(bi+3)=1000*VX; S(bi+4)=1000*VY; S(bi+5)=1000*VZ;
M(8) = 8.6813*1e25;
end
and
function dzdt = physics(z,t)
dzdt=0*z;
dzdt(3) = z(1) * cos(z(2));
dzdt(4) = z(1) * sin(z(2));
dzdt(2) = -9.81/z(1) * cos(z(2));
D = ((0.5) * (1.29)* (0.72))/2 * (z(1)^2);
dzdt(1) = -D/(80) - 9.81*sin(z(2));
end
function[T] = GravAcc(M1, P1, P2)
T = P2-P1;
T = -T*(6.67408*1e-11)*M1/(norm(T)^3);
end

1 Comment

Image Analyst on 1 Dec 2019
What does this show in the command window:
whos S
whos AU
whos C1
whos C2
whos R
whos DT
bi
Chances are your index is more than the length of them, like S does not have bi+2 elements in it.

Vimal Rathod on 5 Dec 2019
You could try debugging your code by putting a breakpoint on the line of error and check if the concerned variables are defined and have expected values.