#!/usr/bin/perl
use strict;

use Data::Dumper;
use Math::Trig;

my $report = 0;

my @sign_list = qw (Aries Taurus Gemini Leo Cancer Virgo Libra Scorpio Sagittarius Capricorn Aquarius Pisces);

#vectors by name
my %pos_of;

# read in all points
foreach my $point_file ("axes.txt", "points.txt") {
    open(POINTS, $point_file) || die "can't open $point_file $!";
    while (my $point_line = <POINTS>) {
	chomp;
	my ($name, $x, $y, $z) = $point_line =~ /(.+)\s+x=(\S+)\s+y=(\S+)\s+z=(\S+)\s_*$/;
	$pos_of{$name} = [$x, $y, $z];
    }
    close(POINTS);
}

# optional: delete Mercury
delete $pos_of{Mercury};

# positions of axes in 
my @axis = qw(ecliptic_x normal vernal);
my @univ_to_ecliptic; 
foreach my $i (0 .. 2) {
    foreach my $j (0 .. 2) {
	$univ_to_ecliptic[$i][$j] = $pos_of{$axis[$j]}[$i];
    }
}
# multiply vector[i] by matrix[i] to get ecliptic-frame version

foreach my $name (keys %pos_of) {
    $pos_of{$name} = transformvector($pos_of{$name}, \@univ_to_ecliptic);
}

# now we have everything in true mars ecliptic coordinates

# Ascendant is where the horizon intersects the ecliptic
$pos_of{Ascendant} = normalize(cross($pos_of{normal}, $pos_of{horizon_normal}));
$pos_of{Descendant} = vsmul($pos_of{Ascendant}, -1);

# Medium Coeli meridian goes through pole and horizon_normal
$pos_of{mc_meridian} = normalize(cross($pos_of{pole}, $pos_of{horizon_normal}));
# Medium Coeli is intersection of mc_meridian and ecliptic
$pos_of{MediumCoeli} = normalize(cross($pos_of{mc_meridian}, $pos_of{normal}));
$pos_of{ImumCoeli} = vsmul($pos_of{MediumCoeli}, -1);

# get houses
my @mc_as_house = house_set('MediumCoeli');
$pos_of{HouseCusp10} = @mc_as_house[0];
$pos_of{HouseCusp11} = @mc_as_house[1];
$pos_of{HouseCusp12} = @mc_as_house[2];
$pos_of{HouseCusp1} = @mc_as_house[3];

my @ic_as_house = house_set('ImumCoeli');
$pos_of{HouseCusp4} = @ic_as_house[0];
$pos_of{HouseCusp3} = @ic_as_house[1];
$pos_of{HouseCusp2} = @ic_as_house[2];
$pos_of{HouseCusp1} = @ic_as_house[3];

foreach my $i (5 .. 9) {
    my $reverse = $i + 6;
    $reverse -= 12 if $reverse > 12;
    $pos_of{"HouseCusp$i"} = vsmul($pos_of{"HouseCusp$reverse"}, -1);
}

my %angle_of;
foreach my $point (keys %pos_of) {
    my $pos = $pos_of{$point};
    $angle_of{$point} = rad2deg(atan2($pos->[0], $pos->[2]));
}

# All the points we care about (planets, ascendant, mc) have ucfirst names
my @planets = grep {/^[A-Z][a-z]/} keys %angle_of;

# find aspects
my @point_list = grep {$_ !~ /HouseCusp|ImumCoeli|Descendant/} @planets;
# brighter planets conjunct easier: 
# Sol and Moon from Earth, Sol and both Moons from Mars.
my %brighter = (Sol => 1, Phobos => 1, Deimos => 1);
my @aspect_list;
my @aspect_type_list = (
    ['conjunction', 0, 8, 10],
    ['opposition', 180, 8],
    ['trine', 120, 8],
    ['square', 90, 8],
    ['sextile', 60, 3],
    ['quincunx', 150, 3],
    ['semisquare', 45, 2],
    ['sesquiquadrate', 135, 2],
    ['semisextile', 30, 2],
    ['quintile', 72, 2],
    ['biquintile', 144, 2],
    # ['none', 90, 180],
);
for (my $i = 0; $i < @point_list; $i++) {
    for (my $j = $i + 1; $j < @point_list; $j++) {
	my $name_1 = $point_list[$i];
	my $name_2 = $point_list[$j];
	my $angle_1 = $angle_of{$name_1};
	my $angle_2 = $angle_of{$name_2};
	my $separation = angle_separation($angle_1, $angle_2);
	$separation = abs $separation;
	# $separation <= 180
	ASPECT: foreach my $aspect_type (@aspect_type_list) {
	    my ($aspect_name, $aspect_angle, $aspect_error, $bright_error) = 
		@$aspect_type;
	    if ($bright_error && ($brighter{$name_1} || $brighter{$name_2})) {
		$aspect_error = $bright_error;
	    }
	    if (my $orb = near($separation, $aspect_angle, $aspect_error)) {
		push @aspect_list, [$name_1, $name_2, $aspect_name, $orb];
		last ASPECT;
	    }
	}
    }
}

# check for retrograde
#  Retro-able points are aspect-able points, minus the sun, moons,
# and location-based points
my @no_retro_list = qw(Sol Deimos Phobos MediumCoeli Ascendant);
my %no_retro_check_for = ( map {($_ => 1)} @no_retro_list);
my %in_retrograde;
foreach my $point (@point_list) {
    if ($no_retro_check_for{$point}) {
	$in_retrograde{$point} = 0;
    }
    else {
	my $pre_motion = angle_separation(
	    $angle_of{"yesterday_$point"}, $angle_of{$point});
	my $post_motion = angle_separation(
	    $angle_of{$point}, $angle_of{"tomorrow_$point"});
	if ($post_motion > 0 && $pre_motion > 0) {
	    $in_retrograde{$point} = 0;
	} 
	elsif ($post_motion < 0 && $pre_motion < 0) {
	    $in_retrograde{$point} = 1;
	}
	else {
	    print STDERR "retrograde check for $point inconclusive\n";
	}
    }
}


if ($report) {
    foreach my $planet ( sort {$angle_of{$a} <=> $angle_of{$b}} @planets) {
	my $angle = $angle_of{$planet};
	if ($angle < 0) {
	    $angle += 360;
	} 
	my $sign = $sign_list[$angle / 30];
	my $degrees = int ($angle % 30);
	my $minutes = ($angle * 60) % 60;
	my $retrograde = $in_retrograde{$planet} ? 'Retrograde': undef;
	print "$planet: $degrees&deg; $sign $minutes' $retrograde<br/>\n";
    }
    foreach my $aspect_set (sort {$a->[2] cmp $b->[2]} @aspect_list) {
	my ($first, $second, $aspect, $orb) = @$aspect_set;
	$orb = abs($orb);
	my $degrees = int ($orb % 30);
	my $minutes = ($orb * 60) % 60;
	print "$first - $aspect - $second (orb $degrees&deg; $minutes')<br/>\n";
    }

}
else {
    # print angles
    print "/AngleOf <<\n";
    foreach my $planet ( sort {$angle_of{$a} <=> $angle_of{$b}} @planets) {
	print "    /$planet $angle_of{$planet}\n";
    }
    print ">> def \n";

    print "/IsRetrograde <<\n";
    foreach my $planet (keys %in_retrograde) {
	my $bool = $in_retrograde{$planet} ? 'true' : 'false';
	print "    /$planet $bool\n";
    }
    print ">> def\n";

    print "/AspectList [\n";
    foreach my $aspect_set (sort {$a->[2] cmp $b->[2]} @aspect_list) {
	my $orb = pop @$aspect_set;
	print "    [", join (" ", (map "/$_", @$aspect_set)), " $orb ]\n";
    }
    print "] def\n";
}

sub angle_separation {
    my ($angle_1, $angle_2) = @_;
    my $separation = $angle_2 - $angle_1;
    while ($separation < 0) {
	$separation += 360;
    }
    if (abs $separation > 180) {
	$separation -= 360;
    }
    return $separation;
}

sub near {
    my $number_a = shift;
    my $number_b = shift;
    my $error = shift;
    my $diff = $number_a - $number_b;
    if (abs($diff) <= $error) {
	return $diff || '0 but true';
    }
    else {
	return;
    } 
}

sub transformvector{
    my $vec = shift;
    my $transform_matrix = shift;
    my $new_vec = [0, 0, 0];
    foreach my $i (0 .. 2) {
	foreach my $j (0 .. 2) {
	    $new_vec->[$j] += $vec->[$i] * $transform_matrix->[$i][$j];
	}
    }
    return $new_vec;
}


sub vlength {
    my $vec = shift;
    my $sum = 0;
    foreach my $i (@$vec) {
	$sum += $i * $i;
    }
    return 0 if $sum <= 0;
    my $len = sqrt $sum;
    return $len;
}
sub vsmul {
    my $vec = shift;
    my $sca = shift;
    return [ map {$_ * $sca} @$vec];
}
sub normalize {
    my $vec = shift;
    my $len = vlength($vec);
    return $vec if $len <= 0;
    my $norm = vsmul($vec, 1 / $len);
    return $norm;
}
sub cross {
    my $a_vec = shift;
    my $b_vec = shift;
    my $x_vec = [];
    foreach my $i (0..2) {
	my $j = ($i + 1) % 3;
	my $k = ($i + 2) % 3;
	$x_vec->[$i] = 
	    $a_vec->[$j] * $b_vec->[$k] - $a_vec->[$k] * $b_vec->[$j];
    }
    return $x_vec;
}

sub vprint {
    my $vec = shift;
    my $out = sprintf("%0.6g %0.6g %0.6g (%0.6g)\n", @$vec, vlength($vec));
    return $out;
}

sub house_set {
    my $set_name = shift;
    # figure out houses
    # 0. transfer positions we need into local storage
    my %ec_pos_of;
    foreach my $name (qw(pole normal horizon_normal), $set_name) {
	$ec_pos_of{$name} = $pos_of{$name};
    }

    # 1. get equatorial coordinates: y is pole, z is mc x pole, x is y x z.
    $ec_pos_of{equatorial_z} = 
	normalize(cross($ec_pos_of{$set_name}, $ec_pos_of{pole}));
    $ec_pos_of{equatorial_x} = 
	cross($ec_pos_of{pole}, $ec_pos_of{equatorial_z});
    my @equatorial_axis = qw(equatorial_x pole equatorial_z);
    my @ecliptic_to_equatorial; 
    foreach my $i (0 .. 2) {
	foreach my $j (0 .. 2) {
	    $ecliptic_to_equatorial[$i][$j] = 
		$ec_pos_of{$equatorial_axis[$j]}[$i];
	}
    }

    # 2. transform reference point and horizon
    my %eq_pos_of;
    foreach my $name (qw(horizon_normal normal), $set_name) {
	$eq_pos_of{$name} = 
	    transformvector($ec_pos_of{$name}, \@ecliptic_to_equatorial);
    }

    # 3. get horizon_normal angle for ref point
    # this is the most complicated bit - we want the normal for the horizon
    # that goes through the rp.  Since we're in equatorial coordinates, staying
    # at a constant latitude means a constant y.  The horizon_normal is 90 deg
    # from the rp, so it's on the great circle normal to the rp.  We chose the
    # reference frame such that we put the rp on the xy plane to simplify 
    # the math.  I did all the arithmetic on a whiteboard, so you'll have to do 
    # them yourself or trust me.
    { # scope block
	my $Hy = $eq_pos_of{horizon_normal}[1];
	my @MC = @{$eq_pos_of{$set_name}};
	my $Rx = -1 * ($MC[1] * $Hy) / $MC[0];
	# two solutions - which do we want?
	my $Rz = 1 * sqrt(1 - $Hy**2 - $Rx**2);
	$eq_pos_of{first_house_normal} = [$Rx, $Hy, $Rz];
    }

    # sanity check
    # my $pseudo_mc = 
	# normalize(cross($eq_pos_of{normal}, $eq_pos_of{first_house_normal}));
    # print "pseudo_mc = ", vprint($pseudo_mc);
    # print "real_mc = ", vprint($eq_pos_of{$set_name});

    # 4. get house_normal angles
    my @house_normal_angle;
    my $eq_1h_pos = $eq_pos_of{first_house_normal};
    $house_normal_angle[0] = rad2deg(atan2($eq_1h_pos->[0], $eq_1h_pos->[2]));
    my $eq_hn_pos = $eq_pos_of{horizon_normal};
    $house_normal_angle[3] = rad2deg(atan2($eq_hn_pos->[0], $eq_hn_pos->[2]));
    foreach my $base (0,  3) {
	while ($house_normal_angle[$base] < 0) {
	    $house_normal_angle[$base] += 360;
	}
    }
    my $start = $house_normal_angle[0];
    my $end = $house_normal_angle[3];
    while ($end < $start) {
	$end += 360;
    }
    if (($end - $start) > 180) {
	$end -= 360;
    }
    my $step = ($end - $start) / 3;
    $house_normal_angle[1] = $start + $step;
    $house_normal_angle[2] = $end - $step;
    foreach my $angle (@house_normal_angle) {
	$angle -= 360 if $angle > 180;
    }

    # 5. get normal vectors
    my $xz_len = sqrt(1 - abs($eq_hn_pos->[1]) ** 2);
    my @house_normal_vector;
    foreach my $i ( 0..3) {
	my $normal_angle = deg2rad($house_normal_angle[$i]);
	$house_normal_vector[$i] = [ 
	    sin($normal_angle) * $xz_len, 
	    $eq_hn_pos->[1], 
	    cos($normal_angle) * $xz_len 
	];
    } 

    # 6. transform back to ecliptic coordinates
    my @equatorial_to_ecliptic;
    foreach my $i (0 .. 2) {
	foreach my $j (0 .. 2) {
	    $equatorial_to_ecliptic[$i][$j] = $ecliptic_to_equatorial[$j][$i];
	}
    }
    my @ecliptic_house_normal_vector = map {
	transformvector($_, \@equatorial_to_ecliptic)
    } @house_normal_vector;
    # print map {vprint($_)} @ecliptic_house_normal_vector;
    # print vprint($pos_of{horizon_normal});

    # 7. get ecliptic intersections
    my @ecliptic_house_vector;
    foreach my $i (0..3) {
	my $cusp_num = ($i + 10) % 12;
	$ecliptic_house_vector[$i] =  
	    normalize(
		cross($pos_of{normal}, $ecliptic_house_normal_vector[$i]));
    }
    return @ecliptic_house_vector;
}