- #1

SpaceGuy17

- 2

- 0

- TL;DR Summary
- In my code below I have vectors r and v, r = [r1, 0, 0] v = [0, sqrt(mu/r1), 0] which define the initial position and velocity of the spacecraft. I want to use ode45 on r(double dot)=-(mu)*r/rmag^3 which defines a two body problem. My objective is to plot a Hohmann transfer of the spacecraft from Earth to an asteroid. In this problem I am neglecting the gravitational Sun on the transfer and treating Earth as the centre [0 0 0]. Any help will be greatly appreciated!

Matlab:

```
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, 0];
av = [0, 10000, 10000];
r1 = Re+500;
r2 = r1+Ra;
r = [r1, 0, 0]; %initial position of the spacecraft
v = [0, sqrt(mu/r1), 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, 0]; %initial position of the spacecraft
v = [0, sqrt(mu/r1), 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
```