How can I plot a Hohmann Transfer Orbit in MATLAB using ode45?

[SpaceGuy17]

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!

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

 

[mmwave]

Thank you for your interest in asteroid mining. I am excited to see more people exploring the potential of this field. I have reviewed your function for asteroid mining and would like to provide some feedback.

Firstly, I appreciate the use of initial conditions and constants in your code. This allows for easy modification and adaptation to different scenarios.

Next, I noticed that you have calculated the delta-v for the transfer orbit and the asteroid orbit separately. While this approach is valid, I would suggest considering a more comprehensive approach that takes into account the gravitational assist from Earth's gravity during the transfer orbit. This can significantly reduce the amount of delta-v required for the transfer.

Additionally, I would suggest incorporating the effects of solar radiation pressure and other perturbations in your calculations. These factors can have a significant impact on the trajectory of the spacecraft and should not be overlooked.

Finally, I would recommend including a function for optimizing the trajectory and minimizing delta-v to maximize the efficiency of the mission. This can be achieved through techniques such as Lambert's problem or genetic algorithms.

Overall, your function for asteroid mining is a good starting point, but I would encourage you to further explore and refine your code to account for all relevant factors. Thank you for your contribution to the field of asteroid mining.