good chemistry is complicated,and a little bit messy -LW PerlMonks

Rough Delta-V estimate for Hohmann transfer orbits

by cavac (Curate)
 on Nov 20, 2018 at 11:30 UTC ( #1226059=CUFP: print w/replies, xml ) Need Help??

I've been interested in space stuff as far back as i can remember. Quite a while back, Scott Manley had a series on the Youtubes explaining how to calculate orbital changes and i implemented this in Javascript on my Blog.

But i always wanted a command line tool in Perl, so lets get to implementing it. It's not the nicest piece of code and it isn't optimizing the transfer burn (check the different options when to do the burns to minimize Delta-V), nor does it try for things like three or (more) burns/multiple intermediate orbits. But it gives a rough idea of "how bad did they miss the target orbit".

Here is the code (long, so in a "readmore" tag):

```#!/usr/bin/env perl

use strict;
use warnings;

# magic numbers
my \$km = 1000; # 1 km = 1000 meter
my \$pi = 3.14159265359;
my \$earth_radius = 6371 * \$km; # earth radius at equator in meters
my \$M = 5.97219 * (10 ** 24); # earth mass in kg
my \$G = 6.67384 * (10 ** -11); # gravitational constant
my \$mu = \$M * \$G;

print "Constants:\n";
print "   Earth mass: \$M kg\n";
print "   Gravity: \$G\n";
print "   mu: \$mu\n";
print "\n";

#my @av = ('1000/e0.5/0', '1000/1000/10');

my %source = convertOrbitDefinition('Source', shift @ARGV);
my %target = convertOrbitDefinition('Target', shift @ARGV);

my \$totaldelta = 0;

print "Calculating semi major axis...\n";
\$source{semi} = (\$source{per} + \$source{ap}) / 2;
\$target{semi} = (\$target{per} + \$target{ap}) / 2;

print "Defining Hohmann orbit...\n";
my %hohmann = (
per => \$source{per},
ap => \$target{per},
incl => 0,
);
\$hohmann{semi} = (\$hohmann{per} + \$hohmann{ap}) / 2;

print "Orbital parameters for Hohmann orbit:\n";
print "   Periapsis: \$hohmann{per} m\n";
print "   Apoapsis: \$hohmann{ap} m\n";
print "   Inclination: \$hohmann{incl} deg\n";
print "\n";

print "Calculating first burn...\n";

# Speed at source orbit periapsis
my \$source_per_speed = sqrt(abs(\$mu * (2/\$source{per} - 1/\$source{semi
+})));

# Required speed at transfer orbit periapsis (at the same point in orb
+it as the source orbit periapsis)
my \$hohmann_per_speed = sqrt(abs(\$mu * (2/\$hohmann{per} - 1/\$hohmann{s
+emi})));

# Speed difference
my \$first_burn = abs(\$hohmann_per_speed - \$source_per_speed);

\$totaldelta += \$first_burn;

print "Data for burn at source orbit periapsis:\n";
print "   Speed before burn: \$source_per_speed m/s\n";
print "   Speed after burn: \$hohmann_per_speed m/s\n";
print "   Required Delta-V: \$first_burn m/s\n";
print "\n";

print "Calculating second burn...\n";
#
# Speed at source orbit periapsis
my \$target_per_speed = sqrt(abs(\$mu * (2/\$target{per} - 1/\$target{semi
+})));

# Required speed at transfer orbit apoapsis (at the same point in orbi
+t as the target orbit periapsis)
my \$hohmann_ap_speed = sqrt(abs(\$mu * (2/\$hohmann{ap} - 1/\$hohmann{sem
+i})));

my \$second_burn;
my \$change_angle_deg = abs(\$source{incl} - \$target{incl});
if(\$change_angle_deg == 0) {
# No inclination change, just calculate the difference
\$second_burn = abs(\$target_per_speed - \$hohmann_ap_speed);
} else {
print "Need to include plane change in second burn...\n";

# Convert to radians
my \$change_angle_rad = \$change_angle_deg * \$pi / 180;
\$second_burn = sqrt(abs((\$hohmann_ap_speed ** 2) + (\$target_per_sp
+eed ** 2) - 2 * \$hohmann_ap_speed * \$target_per_speed * cos(\$change_a
}

print "Data for burn at Hohmann orbit apoapsis:\n";
print "   Speed before burn: \$hohmann_ap_speed m/s\n";
print "   Speed after burn: \$target_per_speed m/s\n";
if(\$change_angle_deg != 0) {
print "   Inclination before burn: \$source{incl} deg\n";
print "   Inclination after burn: \$target{incl} deg\n";
print "   Inclination change: \$change_angle_deg deg\n";
}
print "   Required Delta-V: \$second_burn m/s\n";
print "\n";
\$totaldelta += \$second_burn;

print "Total Delta-V requirement: \$totaldelta m/s\n";

exit(0);

sub convertOrbitDefinition {
my (\$name, \$rawstring) = @_;

my (\$per, \$ap, \$incl) = split/\//, \$rawstring; # Periapsis, Apoaps
+is, Inclination

if(\$ap =~ /^e/) {
print \$name, " orbit is using eccentricity notation, convertin
+g...\n";
\$ap =~ s/^e//g;
my \$semi = \$per; # Semi major axis
my \$ecc = \$ap; # Eccentricity
\$ap = \$semi * (1 + \$ecc);
\$per = \$semi * (1 - \$ecc);
}
if(\$per > \$ap) {
print "Swapping periapsis and apoapsis for \$name orbit...\n";
my \$tmp = \$per;
\$per = \$ap;
\$ap = \$tmp;
}

print "Converting to meter...\n";
\$per *= \$km;
\$ap *= \$km;

print "Orbital parameters for \$name orbit:\n";
print "   Periapsis: \$per m\n";
print "   Apoapsis: \$ap m\n";
print "   Inclination: \$incl deg\n";
print "\n";

return (
ap => \$ap,
per => \$per,
incl => \$incl,
);
}

The basic commandline is perl lithobreak.pl 1000/e0.5/0 1000/1000/10

First argumnent is the source orbit (initial orbit), the second is the target orbit. You can define each orbit in two ways, either using periapsis/apoapsis/inclination or semimajor-axis/eccentricity/inclination. All measurements in kilometers above earth surface. This example uses both:

1. Source orbit of 1000/e0.5/0: Semi major axis of 1000km, "e" defined eccentricity mode with an eccentricity of 0.5, and an inclination of 0 degrees. Internally, this gets converted to a periapsis of 500km and an apoapsis of 1500km. Could have also been written as 500/1500/0
2. Target orbit of 1000/1000/10: Circular orbit of 1000km with a 10 degree inclination. Could have also been written as 1000/e0/10.

Note: when using ap/per mode, don't worry about using the wrong order, it gets sorted out internally.

And here is the result of our calculation:

```Constants:
Earth radius: 6371000 meter
Earth mass: 5.97219e+24 kg
Gravity: 6.67384e-11
mu: 398574405096000

Source orbit is using eccentricity notation, converting...
Converting to meter...
Orbital parameters for Source orbit:
Periapsis: 6871000 m
Apoapsis: 7871000 m
Inclination: 0 deg

Converting to meter...
Orbital parameters for Target orbit:
Periapsis: 7371000 m
Apoapsis: 7371000 m
Inclination: 10 deg

Calculating semi major axis...
Defining Hohman orbit...
Orbital parameters for Hohman orbit:
Periapsis: 6871000 m
Apoapsis: 7371000 m
Inclination: 0 deg

Calculating first burn...
Data for burn at source orbit periapsis:
Speed before burn: 7870.39409917514 m/s
Speed after burn: 7748.85334888779 m/s
Required Delta-V: 121.540750287353 m/s

Calculating second burn...
Need to include plane change in second burn...
Data for burn at Hohman orbit apoapsis:
Speed before burn: 7223.22227109049 m/s
Speed after burn: 7353.45599234394 m/s
Inclination before burn: 0 deg
Inclination after burn: 10 deg
Inclination change: 10 deg
Required Delta-V: 1277.04850388302 m/s

Total Delta-V requirement: 1398.58925417037 m/s

Yeah, 1.3km/s required Delta-V. Not good. Better call our Kerbals to build us a new sat, this one isn't going to make it...

perl -e 'use MIME::Base64; print decode_base64("4pmsIE5ldmVyIGdvbm5hIGdpdmUgeW91IHVwCiAgTmV2ZXIgZ29ubmEgbGV0IHlvdSBkb3duLi4uIOKZqwo=");'

Create A New User
Node Status?
node history
Node Type: CUFP [id://1226059]
Front-paged by haukex
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others making s'mores by the fire in the courtyard of the Monastery: (6)
As of 2019-10-17 09:05 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?
In 2019 the site I miss most is:

Results (42 votes). Check out past polls.

Notices?