Issue with the cross product of my r and v vectors.

So I wrote code where you take an r and v vector and it will give you back 6 orbital elements. I was given a bunch of different r and v vectors dependent on time and I'm supposed to create 6 graphs of the orbital elements based on these r and v vectors. I found equation estimates of the components of the r and v vectors and put those into my r and v vectors in the code. When I run it I get an error for the cross product saying A and B must be of length 3 in the dimension in which the cross product is taken. I also know the values are wrong for the magnitudes. I think the code is adding up all of these different vectors magnitudes but I need them to just graph the individual values over the time interval. And I need them to take the individual cross products of the different vectors and then graph those individual values. Any suggestions?

t=linspace(0,600,100);
r=[-0.0037.*t.^2-7.5378.*t+1269.5,-0.0048.*t.^2+0.6894.*t.^1+6293.9,-0.0001.*t.^2+3.3758.*t.^1+1668.6];% in km
rmag=(sum(r.^2))^0.5;
v=[6*10^(-6).*t.^2-0.0108.*t.^1-7.1346,-8*10^(-8).*t.^2-0.0093.*t.^1+0.6499,-2*10^(-6).*t.^2+0.0012.*t.^1+3.1943];% in km/s
vmag=(sum(v.^2))^0.5;
h=cross(r,v);% in km^2/s
hmag=(sum(h.^2))^0.5;
k=[0,0,1];
n=cross(k,h);
nmag=(sum(n.^2))^0.5;
u=3.986*10^5;% in km^3/s^2
e=(1/u)*((vmag^2-u/rmag)*r-((dot(r,v))*v));
emag=(sum(e.^2))^0.5;%eccentricity(e)
a=hmag^2/(u*(1-emag^2));%semimajor axis in km
h_k=h(1,3);
i=acos(h_k/hmag);%inclination in radians
n_i=n(1,1);
omega=acos(n_i/nmag);%longitude of ascending node in radians
w=acos(dot(n,e)/(nmag*emag));%argument of periapsis in radians
v_0=acos(dot(e,r)/(emag*rmag));%true anomaly at epoch in radians

Answer 1

[James Tursa]

You need to make your r and v matrices of 3-element column vectors, not one long string of number all in a row. That means putting ; instead of , for delimiters. And then downstream in your code you need to replace all of your matrix operators * and / and ^ with the element-wise operators .* and ./ and .^ for things to work. Finally, you should double check those last three lines. The acos( ) function only gives results in a 180 deg range, but you probably need to use atan2( ) for these since they need to be 360 deg angles. E.g.,

t=linspace(0,600,100);
r=[-0.0037.*t.^2-7.5378.*t+1269.5;-0.0048.*t.^2+0.6894.*t.^1+6293.9;-0.0001.*t.^2+3.3758.*t.^1+1668.6];% in km
rmag=(sum(r.^2)).^0.5;
v=[6*10^(-6).*t.^2-0.0108.*t.^1-7.1346;-8*10^(-8).*t.^2-0.0093.*t.^1+0.6499;-2*10^(-6).*t.^2+0.0012.*t.^1+3.1943];% in km/s
vmag=(sum(v.^2)).^0.5;
h=cross(r,v);% in km^2/s
hmag=(sum(h.^2)).^0.5;
k=repmat([0;0;1],1,size(r,2));
n=cross(k,h);
nmag=(sum(n.^2)).^0.5;
u=3.986*10^5;% in km^3/s^2
e=(1/u)*((vmag.^2-u./rmag).*r-((dot(r,v)).*v));
emag=(sum(e.^2)).^0.5;%eccentricity(e)
a=hmag.^2./(u*(1-emag.^2));%semimajor axis in km
h_k=h(3,:);
i=acos(h_k./hmag);%inclination in radians
n_i=n(1,:);
% DOUBLE CHECK THESE NEXT THREE LINES ... THEY LOOK WRONG TO ME
omega=acos(n_i./nmag);%longitude of ascending node in radians
w=acos(dot(n,e)./(nmag.*emag));%argument of periapsis in radians
v_0=acos(dot(e,r)./(emag.*rmag));%true anomaly at epoch in radians

[Son of Coul]

Thank you. Such a simple mistake I should've known