Asked by Alexandru Lapusneanu
on 20 Dec 2017

%%%%%%%%%%%%Valori pentru Rcsc

%%%%Pozitiile si vitezele pe cele 3 axe

y0(1,1)= 743322.3616 ;

y0(2,1)= -6346021.219 ;

y0(3,1)= -3394131.349 ;

y0(4,1)= 5142.38067;

y0(5,1)= 4487.44895 ;

y0(6,1)= -7264.00872;

%%%%Timpul

tspan=[0 :8640];

%%%%Masa(kg) si aria suprafetei satelitului (m^2)

m = 217 ; %320;

A = 1.2; %8;

%%%%Metoda Runge-Kutta de ordin 4

h=1;

y = zeros(6, tspan(end)/h);

y(:,1) = y0;

for i=1:(tspan(end)/h)

H=sqrt(y(1,i)^2+y(2,i)^2+y(3,i)^2);

k_1 = proiectia(tspan(i), y(:,i), H, m, A, y(4:6, i));

k1=double(k_1);

k_2 = proiectia(tspan(i)+0.5*h, y(:,i)+0.5*h*k_1, H, m, A, y(4:6, i));

k2=double(k_2);

k_3 = proiectia((tspan(i)+0.5*h), (y(:,i)+0.5*h*k_2), H, m, A, y(4:6, i));

k3=double(k_3);

k_4 = proiectia((tspan(i)+h),(y(:,i)+k_3*h), H, m, A, y(4:6, i));

k4=double(k_4);

y(:,i+1) = double(y(:,i) + (1/6)*(k1+2*k2+2*k3+k4)*h);

end

%%%Distanta satelitului

Rcsc = ((y(1,:).^2 + y(2,:).^2 + y(3,:).^2).^0.5);

n=50;

%plot(tspan,Rcsc)

%%Textured 3D Earth example

%

% Ryan Gray

% 8 Sep 2004

% Revised 9 March 2006, 31 Jan 2006, 16 Oct 2013

%%Options

space_color = 'k';

npanels = 180; % Number of globe panels around the equator deg/panel = 360/npanels

alpha = 1; % globe transparency level, 1 = opaque, through 0 = invisible

GMST0 = []; % Don't set up rotatable globe (ECEF)

%GMST0 = 4.89496121282306; % Set up a rotatable globe at J2000.0

% Earth texture image

% Anything imread() will handle, but needs to be a 2:1 unprojected globe

% image.

image_file = 'https://upload.wikimedia.org/wikipedia/commons/thumb/c/cd/Land_ocean_ice_2048.jpg/1024px-Land_ocean_ice_2048.jpg';

% Mean spherical earth

erad = 6371008.7714; % equatorial radius (meters)

prad = 6371008.7714; % polar radius (meters)

erot = 7.2921158553e-5; % earth rotation rate (radians/sec)

%%Create figure

figure('Color', space_color);

orbit=animatedline(y(1,:),y(2,:),y(3,:),'Color','y');

hold on;

% Turn off the normal axes

set(gca, 'NextPlot','add', 'Visible','off');

axis equal;

axis auto;

% Set initial view

view(0,30);

axis vis3d;

%%Create wireframe globe

% Create a 3D meshgrid of the sphere points using the ellipsoid function

[x, y, z] = ellipsoid(0, 0, 0, erad, erad, prad, npanels);

globe = surf(x, y, -z, 'FaceColor', 'none', 'EdgeColor', 0.5*[1 1 1]);

%%Texturemap the globe

% Load Earth image for texture map

cdata = imread(image_file);

% Set image as color data (cdata) property, and set face color to indicate

% a texturemap, which Matlab expects to be in cdata. Turn off the mesh edges.

set(globe, 'FaceColor', 'texturemap', 'CData', cdata, 'FaceAlpha', alpha, 'EdgeColor', 'none');

I'm sorry for the appearance of the code!

Answer by Meysam Mahooti
on 9 Mar 2019

Edited by Meysam Mahooti
on 9 Mar 2019

Sign in to comment.

Answer by Walter Roberson
on 20 Dec 2017

Alexandru Lapusneanu
on 21 Dec 2017

%%%%%%%%%%%%Valori pentru Rcsc

tspan=[0 :86400];

figure;

orbit=animatedline('MaximumNumPoints',length(tspan));

%%%%Pozitiile si vitezele pe cele 3 axe

y0(1,1)= 743322.3616 ;

y0(2,1)= -6346021.219 ;

y0(3,1)= -3394131.349 ;

y0(4,1)= 5142.38067;

y0(5,1)= 4487.44895 ;

y0(6,1)= -7264.00872;

%%%%Timpul

%%%%Masa(kg) si aria suprafetei satelitului (m^2)

m = 217 ; %320;

A = 1.2; %8;

%%%%Metoda Runge-Kutta de ordin 4

h=1;

y = zeros(6, tspan(end)/h);

y(:,1) = y0;

for i=1:(tspan(end)/h)

H=sqrt(y(1,i)^2+y(2,i)^2+y(3,i)^2);

k_1 = proiectia(tspan(i), y(:,i), H, m, A, y(4:6, i));

k1=double(k_1);

k_2 = proiectia(tspan(i)+0.5*h, y(:,i)+0.5*h*k_1, H, m, A, y(4:6, i));

k2=double(k_2);

k_3 = proiectia((tspan(i)+0.5*h), (y(:,i)+0.5*h*k_2), H, m, A, y(4:6, i));

k3=double(k_3);

k_4 = proiectia((tspan(i)+h),(y(:,i)+k_3*h), H, m, A, y(4:6, i));

k4=double(k_4);

y(:,i+1) = double(y(:,i) + (1/6)*(k1+2*k2+2*k3+k4)*h);

addpoints(orbit,y(1,:),y(2,:),y(3,:));

drawnow

end

Walter Roberson
on 21 Dec 2017

Alexandru Lapusneanu
on 21 Dec 2017

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.