Beefy Boxes and Bandwidth Generously Provided by pair Networks
Pathologically Eclectic Rubbish Lister
 
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 radius: $earth_radius meter\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 +ngle_rad) )); } 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 "Adding earth radius...\n"; $per += $earth_radius; $ap += $earth_radius; 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... Adding earth radius... Orbital parameters for Source orbit: Periapsis: 6871000 m Apoapsis: 7871000 m Inclination: 0 deg Converting to meter... Adding earth radius... 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=");'

Log In?
Username:
Password:

What's my password?
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 wandering the Monastery: (7)
As of 2019-02-22 21:27 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?
    I use postfix dereferencing ...









    Results (119 votes). Check out past polls.

    Notices?