#!/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 orbit as the source orbit periapsis)
my $hohmann_per_speed = sqrt(abs($mu * (2/$hohmann{per} - 1/$hohmann{semi})));

# 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 orbit as the target orbit periapsis)
my $hohmann_ap_speed = sqrt(abs($mu * (2/$hohmann{ap} - 1/$hohmann{semi})));

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_speed ** 2) - 2 * $hohmann_ap_speed * $target_per_speed * cos($change_angle_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, Apoapsis, Inclination

    if($ap =~ /^e/) {
        print $name, " orbit is using eccentricity notation, converting...\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,
    );
}



##</code><code>##

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