lowthrust_mars6.m

小推力航天器转移轨道程序

close all
clear all
clc
tic
global z

z = 1;

F = 200e-3; %N
Isp = 3000;
go = 9.81;

R_earth = 6378.1363; %km
mu_earth = 398600.4415; %km^3/s^2

a_earth = 149598023; %km
mu_sun = 1.32712428e11; %km^3/s^2
Mearth = 5.9742e24; %kg
au = 149597870; %km
Msun = 1.988435e30; %kg

a_mars = 1.52367934*au; %km
mu_mars = 4.305e4; %km^3/s^2
R_mars = 3397; %km

SOI_E = (Mearth/Msun)^(2/5) * au; %km
SOI_M = 5.78e5; %km

M_Init = 1500; %kg

%This assumes a GTO orbit
R_E_init = [au;0];
V_E_init = [0;(mu_sun/norm(R_E_init))^.5];

EarthAng = pi;
MarsAng = -1.63363191667473;
temp = R3(-MarsAng)*[a_mars;0;0];
R_M_init = temp(1:2);
V_vec = cross(temp,[0,0,-1]);
V_vec = V_vec/norm(V_vec);
V_M_init = (mu_sun/norm(R_M_init))^.5*V_vec(1:2)';

R_sc_temp = R3(EarthAng)*[R_earth + 300; 0; 0];
R_sc_init = R_E_init + R_sc_temp(1:2);
V_sc_temp = R3(EarthAng)*[0; (2*mu_earth/(R_earth + 300) - mu_earth/24582)^.5; 0];
V_sc_init = V_E_init + V_sc_temp(1:2);

xo = [M_Init;R_sc_init; V_sc_init; R_E_init; V_E_init; R_M_init; V_M_init];
  
Tspan = [0,4*60*60*24*365];
opts = odeset('MaxStep',5000,'RelTol',10^-4,'AbsTol',10^-4*ones(1,13));
[t,x] = ode23(@OrbCalc7,Tspan,xo,opts,F,Isp,go,mu_earth,mu_sun,mu_mars,SOI_E,SOI_M,a_mars,R_mars);

M = x(:,1);
R_sc = x(:,2:3);
V_sc = x(:,4:5);
R_E = x(:,6:7);
V_E = x(:,8:9);
R_M = x(:,10:11);
V_M = x(:,12:13);

R_sc_M = R_sc - R_M;
V_sc_M = V_sc - V_M;

figure(1)
plot(R_sc(:,1),R_sc(:,2),'b-',R_E(:,1),R_E(:,2),'g-',R_M(:,1),R_M(:,2),'r-')
title('Position with Respect to Sun')   
xlabel('km')
ylabel('km')
legend('S/C','Earth','Mars')
axis equal

R_sc_E = R_sc - R_E;
V_sc_E = V_sc - V_E;

figure(2)
plot(R_sc_E(:,1),R_sc_E(:,2))
hold on
circle([0,0],R_earth,200,'g');
title('Position with Respect to Earth')
xlabel('km')
ylabel('km')
axis equal

figure(3)
plot(R_sc_M(:,1),R_sc_M(:,2))
hold on
circle([0,0],R_mars,200,'g');
title('Position with Respect to Mars')
xlabel('km')
ylabel('km')
axis equal

[a,ecc,inc,O,w,nutemp,w_true] = elorb([R_sc_M(length(t),:) 0],[V_sc_M(length(t),:) 0],mu_mars);

nu2 = linspace(0,2*pi,1000);
R = randv(ones(1,length(nu2))*a,ones(1,length(nu2))*ecc,...
          zeros(1,length(nu2))*inc,zeros(1,length(nu2)),...
          ones(1,length(nu2))*w_true,nu2,mu_mars);

h = plot(R(1,:),R(2,:),'r')
set(h,'LineWidth',1.5)
axis equal
hold off

toc