#!/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,
    );
}
[download]