问题描述
在我下面的代码中，我有向量 r 和 v，r = [r1,0] v = [0,sqrt(mu/r1),0] 定义了航天器的初始位置和速度。我想在 r(double dot)=-(mu)*r/rmag^3 上使用 ode45，它定义了一个二体问题。我的目标是绘制飞船从地球到小行星的霍曼转移图，任何帮助将不胜感激！
function Asteroid_Mining
clc
%Initial conditions
g0 = 9.81; %gravity (m/s^2)
p = 1.225; %atmospheric density at sea level (kg/m3)
Re = 6378; %radius of Earth (km)
Ra = 7.431e7; %distance of Bennu from Earth in (km) [August 2023]
G = 6.674e-11/1e9; % Gravitational constant (km3/kg.s2)
mu = 3.986e5; %gravitational parameter (km3/s2)
Me = 5.972e24; %mass of Earth (kg)
Ma = 7.8e10; %Mass of asteroid (kg)
Ms = 2110; %launch mass of spacecraft (kg)
Isp = 230; %specific impulse (s)
tspan = [0 10000]; %time at 0s
%% Relative to a non-rotating Earth-centred Cartesian coordinate system [0 0 0]
ar = [Re+Ra,0];
av = [0,10000,10000];
r1 = Re+500;
r2 = r1+Ra;
r = [r1,0]; %initial position of the spacecraft
v = [0,0]; %initial veLocity of the spacecraft
rv = [r;v];
h = cross(r,v);
rdoth = cross(v,h)*-1;
rmag = norm(r);
e = (rdoth/mu) - (r/rmag);
emag = norm(e); %eccentricity
%Calculating delta-v for transfer orbit
ve1 = sqrt(mu/r1);
h2 = sqrt(2*mu)*sqrt((r1*r2)/(r1+r2));
ve2 = h2/r1;
%Calculating delta-v for Asteroid Orbit
va2 = h2/r2;
va3 = sqrt(mu/r2);
deltaV = norm(ve2-ve1) + norm(va3-va2);
deltaVc = deltaV*1e3;
%mass of propellant consumed
deltaM_M0 = 1 - emag^-(deltaVc/(IsP*g0));
deltaM = deltaM_M0*Ms;
[tout,rvout] = ode45(@Orbit,tspan,rv);
xout = rvout(:,1);
yout = rvout(:,2);
zout = rvout(:,3);
vxout = rvout(:,4);
vyout = rvout(:,4);
vzout = rvout(:,4);
function drvt = Orbit(t,state)
r = [r1,0]; %initial position of the spacecraft
v = [0,0]; %initial veLocity of the spacecraft
rv = [r,v];
x = rv(1);
y = rv(2);
z = rv(3);
vx = rv(4);
vy = rv(5);
vz = rv(6);
r = sqrt(x^2+y^2+z^2);
dxdt = x;
dydt = y;
dzdt = z;
dvxdt = -mu*x*(r^-3);
dvydt = -mu*y*(r^-3);
dvzdt = -mu*z*(r^-3);
drvt = zeros(size(state));
drvt(1) = dxdt;
drvt(2) = dydt;
drvt(3) = dzdt;
drvt(4) = dvxdt;
drvt(5) = dvydt;
drvt(6) = dvzdt;
end
end
解决方法
这对你有用吗？
function Numeric_Kepler
G = 6.674e-11/1e9; %gravitational constant
M_Earth = 5.972e24; %mass of Earth
Re = 6378; %radius of Earth (km)
GM_Earth = G*M_Earth;
dt = 70; %time step in sec
mu = 3.986e5; %gravitational parameter (km3/s2)
r1 = Re+500; %distance between initial position of spacecraft and Earth
%Ra = 7.431e7; %distance of Bennu from Earth in (km) [August 2023]
r_0 = [r1,0]; %initial position of the spacecraft
v_0 = [0,sqrt(mu/r1),0]; %initial uncorrected velocity of the spacecraft
dv0 = [0,1.3,0]; %velocity correction
r0 = r_0; %initial position of spacecraft
v0 = v_0 + dv0; %corrected initial velocity
%initialization of the dynamical variables
r_ = r_0;
v_ = v_0;
r = r0;
v = v0;
%keepoint the plots in the same window
hold on
grid on
% setting up the window's shape by definig the coordinate axes' ranges
r2 = 3*r1;
axis([ -r2 r2 -r2 r2])
% generating the dynamics and potting it at the same time,% simple animation style
plot(0,'ro')
for j = 1:1000
%propagation along spacraft's uncorrected orbit
[r_,v_] = Kepler_flow(r_,v_,dt);
plot(r_(1),r_(2),'bo')
%propagation along spacraft's corrected orbit
[r,v] = Kepler_flow(r,v,dt);
plot(r(1),r(2),'ro')
pause(0.1)
end
function a = Gravity(x)
a = - ( (GM_Earth) / norm(x)^3 ) * x;
end
function [r_1,v_1] = Kepler_step(r,dt)
r_1 = r + (dt/2) * v;
v_1 = v + dt * Gravity(r_1);
r_1 = r_1 + (dt/2) * v_1;
end
function [r_1,v_1] = Kepler_flow(r,dt)
[r_1,dt/2);
[r_1,v_1] = Kepler_step(r_1,v_1,dt);
[r_1,dt/2);
end
end
这是另一个版本，它具有积分器的四阶收敛，与显式 Runge-Kutta 4 相同，但该积分器是时间可逆的和 symplecitc（准守恒系统能量）。对于动态的长时间传播，这种积分器非常稳定、一致且相当准确。
function Numeric_Kepler
G = 6.674e-11/1e9; %gravitational constant
M_Earth = 5.972e24; %mass of Earth
Re = 6378; %radius of Earth (km)
GM_Earth = G*M_Earth;
mu = 3.986e5; %gravitational parameter (km3/s2)
% integration parameters
w0 = nthroot(2,3);
w1 = 1 / (2 - w0);
w0 = - w0*w1;
c = [w1,w0+w1,w1]/2;
d = [w1,w0,w1];
r1 = Re+500; %distance between initial position of spacecraft and Earth
%Ra = 7.431e7; %distance of Bennu from Earth in (km) [August 2023]
dt = 70; %time step in sec
r_0 = [r1; 0; 0]; %initial position of the spacecraft
v_0 = [0; sqrt(mu/r1); 0]; %initial uncorrected velocity of the spacecraft
dv0 = [0; 1.5; 0]; %velocity correction
r0 = r_0; %initial position of spacecraft
v0 = v_0 + dv0; %corrected initial velocity
%initialization of the dynamical variables
r_ = r_0;
v_ = v_0;
v = v0;
r = r0;
%keeping the plots in the same window
hold on
grid on
% setting up the window's shape by defining the coordinate axes' ranges
r2 = 3*r1;
axis([ -r2 r2 -r2 r2])
% generating the dynamics and potting it at the same time,'ro')
for j = 1:300
%propagation along spacecraft's uncorrected orbit
[r_,v_] = Kepler_flow_4(r_,'bo')
%propagation along spacecraft's corrected orbit
[r,v] = Kepler_flow_4(r,v_1] = Kepler_flow_4(r,dt)
r_1 = r;
v_1 = v;
for i = 1:3
r_1 = r_1 + c(i)*dt * v_1;
v_1 = v_1 + d(i)*dt * Gravity(r_1);
end
r_1 = r_1 + c(4)*dt * v_1;
end
end