Using a function with input parameters in ode45

Hi,

Im having a lot of trouble getting some code for an orbital calculation to work and I don't understand why I am getting a particular error. I think its telling me that I am trying to look at a place in my matrix that does not exist, but I don't see how that would be.

??? Attempted to access cartesian(:,1); index out of bounds because
numel(cartesian)=1.

Error in ==> accel_func at 7
AccTB(4:6,1)=-mu/(r^3)*cartesian(1:3);

Error in ==> @(cartesian,velocity)accel_func(cartesian,velocity)

Help much appreciated. The m file:

clear all;
a=6978.14;
e=0;
n=sqrt(398600/(a^3));
M=0;
E=M;
Theta=(cos(E)-e)/(1-e*cos(E));
r=a*(1-e*cos(E));
p=a(1-e^2);
vx=(sqrt(398600/p))*(-sin(Theta));
vy=(sqrt(398600/p))*(e+cos(Theta));
nodem=[cos(0) -sin(0) 0; sin(0) cos(0) 0; 0 0 1];
incm=[1 0 0;0 cos(98*pi/180) -sin(98*pi/180);0 sin(98*pi/180) cos(98*pi/180)];
argm=[cos(0) -sin(0) 0;sin(0) cos(0) 0;0 0 1];
cartesian=nodem*incm*argm*[r*cos(Theta);r*sin(Theta);0];
velocity=nodem*incm*argm*[vx;vy;0];
final=zeros(6,1);

orb_opt=odeset('RelTol',1*10^-11,'AbsTol',[1*10^-11 1*10^-11 1*10^-11 1*10^-11 1*10^-11 1*10^-11]);
tof=5*24*3600;
[t xdot]=ode45(@(cartesian,velocity) accel_func(cartesian,velocity),[0 tof], cartesian, orb_opt)

and the function definition file:

function final=accel_func(cartesian,velocity)
Re=6378.1363;
J2=0.0010826269;
mu=398600;
final=zeros(6,1);
r=norm(cartesian);
AccTB(4:6,1)=-mu/(r^3)*cartesian(1:3);
temp1=(cartesian(3)/r)^2;
temp=5*temp1;
coeff2=-3*J2*mu*Re^2/(2*r^5);
accel_J2X=coeff2*cartesian(1)*(1-temp);
accel_J2Y=coeff2*cartesian(2)*(1-temp);
accel_J2Z=coeff2*cartesian(3)*(3-temp);
Earthrate=7.9292115*10^-5;
rho=5*10^-4;
Cd=2.2;
A=15*10^-6;
mass=400;
vrel=[velocity(1)+Earthrate*cartesian(2) velocity(2)-Earthrate*cartesian(1) velocity(3)];
Accel_dragX=-.5*rho*Cd*A/mass*norm(vrel)*vrel(1);
Accel_dragY=-.5*rho*Cd*A/mass*norm(vrel)*vrel(2);
Accel_dragZ=-.5*rho*Cd*A/mass*norm(vrel)*vrel(3);

final(1)=velocity(1);
final(2)=velocity(2);
final(3)=velocity(3);
final(4)=AccTB(4)+accel_J2X+Accel_dragX;
final(5)=AccTB(5)+accel_J2Y+Accel_dragY;
final(6)=AccTB(6)+accel_J2Z+Accel_dragZ;

Answer 1

[Matt Tearle]

The problem here seems to be that you're mixing up your representation of your dependent variables. The rate equation function accel_func should take two inputs: t and z. In your case z will be a 6-element vector of dependent variables. The return should be the corresponding derivatives. Your equations are autonomous, so the first input t will go unused. You can replace it with ~ if you prefer.

Given the equations, I think z should be [cartesian;velocity]. Hence:

[t, z] = ode45(@accel_func,[0 tof],[cartesian;velocity],orb_opt);

And

function final = accel_func(~,z)
cartesian = z(1:3);
velocity = z(4:6);
  etc