How draw an ellipse by acquiring points in 3D

Hi. I have solved an Orbital Mechanics question about orbiting and acquired three points about the orbits, for example R1=(-9.5i 6.04j 3.1k)*10^3 and R2=(-9.6i 6.03j 3.4k)*10^3 as you see, the orbits are ellipse. i want to show it in cartesian axis. how can i do this..? should i use plot3 formation ? would you please help me please?

[José-Luis]

I don't follow. What are i, j, and k?

[KSSV]

i, j and k might be unit vectors along x,y,z respectively.

[José-Luis]

So you have the coordinates of two points. Two points are not enough to define a unique ellipse.

Yes, you could use plot3() to visualize your points.

[Adam]

Comment by Taha Smith moved from flag:

yeah, But in plot3() I should write z as x and y i mean: t=0:1:2*pi x=cos(xt) y=sin(xt) z=t while z parameter has a measure in this problem and it is 3.1*10^3 in first orbit and 3.4*10^3 in second one. if there is another formation for plot3(), would you please give me ? thanks

[Adam]

@Taha Smith. Please reply with comments rather than flagging someone's comment. Flags are used for reporting spam or other inappropriate content.

[Morteza.Moslehi]

Oh sorry, I didn't know.. you're right

Answer 1

[Morteza.Moslehi]

answer is: step 1:

function [r,v] = coes_RV(h,inc,RA,e,omega,theta)   
%From the orital elements, the position and velocity vectors are is found. 
%The Inputs are inputed the following matter: 
%     h is angular momentum in km^2/s 
%     inc is inclination in radians  
%     RA is right ascension of ascending node in radians 
%     e is eccelntricity (unitless) 
%     omega is argument of perigee in radians 
%     theta is true anomaly in radians 
%Equations (EQN) used come from Curtis       
mu = 398600; %Gravitational parameter of Earth (km^3/s^2)   
%Now following Algorithm 4.2 from Curtis, I use EQN. 4.37 to find r in 
%perifocal frame, and EQN 4.38 to find v in perifocal plane   
rperi = h^2/mu/(1+e*cos(theta))*[cos(theta); sin(theta); 0]; % (km)  
vperi = mu/h.*[-sin(theta); e + cos(theta); 0]; % (km/s)   
%Next use EQN 4.44 to find transformation matrix from perifocl to 
%geocentric equatorial coordinates   

Qperi=[cos(RA)*cos(inc)*sin(omega)+cos(RA)*cos(omega) -sin(RA)*cos(inc)*cos(omega)-cos(RA)*sin(omega) sin(RA)*sin(inc);
   cos(RA)*cos(inc)*sin(omega)+sin(RA)*cos(omega) cos(RA)*cos(inc)*cos(omega)-sin(RA)*sin(omega) -cos(RA)*sin(inc);
   sin(inc)*sin(omega) sin(inc)*cos(omega) cos(inc)];

%Finally use, EQN 4.46 to find the r and V in the geocentric equatorial 
%plane  
r=Qperi*rperi; %in km   
v=Qperi*vperi; %in km/s   
%Convert r and v into row vectors   
r = r'; %km  
v = v'; %km/s

step 2:

function [accel] = orbit(t,y)   
%y is the initial conditions, and time is in seconds   
%Set constant values 
mu = 398600; %Gravitational parameter of Earth (km^3/s^2)   
%Pull out the initial conditions to components 
rx=y(1); %km 
ry=y(2); %km 
rz=y(3); %km 
vx=y(4); %km/s
vy=y(5); %km/s 
vz=y(6); %km/s
%Normalize the position vector for futre use
R=norm([rx, ry, rz]);   %Find acceleration from the position vector
ax=-mu*rx/R^3; %km/s^2 
ay=-mu*ry/R^3; %km/s^2 
az=-mu*rz/R^3; %km/s^2   
%Set up new conditions after t seconds 
accel = [vx; vy; vz; ax; ay; az]; 

step 3:

%Senior Project   
  close all 
  clear all 
  clc   
  %First state initial orbital elements in order to find initial r and v in 
  %geocentric equatorial coordinates 

    %Set constants 
    Re = 6378; %Radius of Earth (km) 
    mu = 398600; %Gravitational Parameter of Earth (km^3/s^2)
    %Angular Momentum 
    h0 = 133929.0857; %km^2/s   
    %Inclination 
    inc0 = 45*pi/180; %radians   
    %Right Ascension of Ascending Node 
    RA0 = 0; %radians   
    %Eccentricity 
    e0 = 0;   
    %Argument of Perigee 
    omega0 = 0*pi/180; %radians   
    %True Anomaly 
    theta0 = 30*pi/180; %radians   
    %Now find r in km and v in km/s 
    [r0,v0] = coes_RV(h0,inc0,RA0,e0,omega0,theta0);   
    %Find Range 
    range0 = norm(r0) - Re; %in km   
    %Find Period of orbit 
    T0 = 2*pi/sqrt(mu)*(norm(r0))^(1.5); %in seconds  
    %Now these are the initial conditions and time span 
    y0 = [r0 v0]; t0 = [0 172800]; %in seconds 
    % t0 = [0 T0]; %in seconds   
    %Use ode45 to get orbit 
    [t,y] = ode45('orbit', t0, y0);   
    % %Plot Orbit 
    % plot3(y(:,1), y(:,2), y(:,3))   
    % %Simulation 
    % for i = 1:length(t) 
    %     comet3(y(1:i,1), y(1:i,2), y(1:i,3)) 
    %     drawnow 
    % end   
    %Simulation 
    for i = 1:length(t)    
        plot3(y(1:i,1), y(1:i,2), y(1:i,3))    
        drawnow
    end

give thanks to Nancy Teresa Cabrera too.

[Matthew Worker]

Hello please could you tell me what y(1),y(2),y(3),y(4) stand for?

Answer 2

[José-Luis]

Adapt to your needs:

z = [0.5,1];
fx = @(x) sin(x);
fy = @(x) cos(x);

t = linspace(0,2*pi,101);

plot3(fx(t),fy(t),repmat(z(1),1,numel(t)));
hold on
plot3(fx(t),fy(t),repmat(z(2),1,numel(t)));

[Morteza.Moslehi]

@José-Luis at first i should say thank you so much for your kind. But look at here. my problem is that, how I can draw a graph with three points: i have written this code to acquire these three pints:

if true
  % A=[cos(pi/2.769231) sin(pi/2.769231) 0;-sin(pi/2.769231) cos(pi/2.769231) 0;0 0 1];
B=[1 0 0;0 cos(pi/7.2) sin(pi/7.2);0 -sin(pi/7.2) cos(pi/7.2)];
C=[cos(pi/36) sin(pi/36) 0;-sin(pi/36) cos(pi/36) 0;0 0 1];
F=A*B*C
G=inv(F)
h=69788.206382
mi=398600
e=0.15
f=(h^2)/(mi*(1+(e*cos(pi/2.4))))
g=[cos(pi/2.4);sin(pi/2.4);0]
R=f*g
J=G*R
end
i want to plot "J" to give me a ellipse.
is it possible?

[José-Luis]

I don't understand what you mean.

[Morteza.Moslehi]

@José-Luis thanks for your favor there. i found it

[Morteza.Moslehi]

in following this solution, i want to change parameters of radius in Cartesian coordinate to spherical coordinate. i know i should use this code:

[teta,phi,r] = cart2sph (x,y,z)

right? how can I plot " teta,phi,r" in spherical coordinate? would you please help me?

i want to demonstrate a spherical and show this orbit on that.

regards