#!/usr/bin/perl

#
# Perl script to analyse the bright star catalogue and 
# compare the angular distribution of nearest neighbours to 
# randomly placed stars.
#

#
# Made by Robert Izzard for the Binary Stars lecture course (astro8501/6944)
# at the University of Bonn 2012.
# While this Perl script is copyrighted, nothing in it is worth copying!
#

# required modules
use strict;
use Math::Trig;
use Math::Trig ':radial';
use Math::Trig ':pi';
use Math::Trig 'great_circle_distance';
use POSIX qw(log10);

# lists of positions of the stars from the catalogue and randomly
# placed stars
my @stars;
my @random;

# output files (for plotting)
open(STARS,">stars_xyz.dat");
open(RAND,">random_xyz.dat");

# input file (RA,Dec from the bright star catalogue)
open(ALL,"<all.dat")||die;
while(<ALL>)
{
    #if(rand()>0.9) # if testing use this to reduce dataset
    {
	{
	    # get bright star positions in RA, DEC (J2000) 
	    # and convert to radians
	    s/^\s+//o;
	    my ($ra,$dec)=split(/\s+/,$_,3);
	    
	    # convert RA/Dec angles > radians
	    my $theta=$ra/24.0*pi2; # azimuth (pi2=2.0*pi)
	    my $phi=($dec+90.0)/180.0*pi; # polar (start from zero)
	    
	    # logging
	    printf STARS "%g %g %g\n",spherical_to_cartesian(1.0, $theta, $phi);
 
	    # save data in @stars array
	    push(@stars,[$theta,$phi]);
	}


	{
	    # generate random points on a sphere
	    my $z=(rand()-0.5)*2.0; # z [-1:1]
	    my $phi=pi2*rand(); # phi [0:2pi]

	    # we have z,phi : need theta so general cartesian x,y
	    my $sq=sqrt(1.0-$z*$z);
	    my $x=$sq*cos($phi);
	    my $y=$sq*sin($phi);
	    
	    # hence spherical co-ordinate
	    my $rho; # always 1
	    my $theta;
	    ($rho, $theta, $phi) = 
		cartesian_to_spherical($x,$y,$z);

	    # logging
	    printf RAND "%g %g %g\n",spherical_to_cartesian($rho, $theta,$phi);

	    # save data
	    push(@random,[$theta,$phi]);
	}
    }
}

close STARS;
close RAND;
close ALL;

print "data set up complete\n";

# bin in histograms and output:

# bright stars
my $h=angdist(\@stars,'bright stars');
output('angdists.dat',$h,$#stars-1);

# random 'stars'
$h=angdist(\@random,'random stars');
output('random.dat',$h,$#stars-1);

exit;

############################################################
# subroutines...
############################################################

sub output
{

    # output data for lecture plots
    my $file=shift;
    my $h=shift;
    my $n=shift;
    print "Output $file\n";

    my $cum=0.0;
    open(DATA,'>'.$file)||die;
    my $norm=1.0/$n;
    foreach my $bin (sort {$a<=>$b} keys %$h)
    {
	$cum+=$$h{$bin};
	printf DATA "%g %g\n",$bin,$cum*$norm;
    }
    close DATA;
}

sub angdist
{
    # calculate angular distribution

    # input
    my $array=shift; # pointer to the array of stars
    my $label=shift; # label (for output to the screen only)
    my $n=$#{$array}; # number of stars (+1)

    print "Stats for $label\n";

    # calculate angular separations
    my $thresh=1e-8;
    my $binwidth=0.1;
    my %histogram;

    # loop over stars to find the closest star to each
    for(my $i=0;$i<$n;$i++)
    {
	my @star=@{$$array[$i]};
	my $minangsep=1e100; # start huge
	my $nearest;
	
	# loop over all other stars (the slow bit)
	for(my $j=0;$j<$n;$j++)
	{
	    next if($i==$j); # *other* stars!

	    # find angular separation
	    # see e.g. http://www.astro.uu.nl/~strous/AA/en/reken/grootcirkel.html (eq 4)
	    my $angsep=great_circle_distance(@star,@{$$array[$j]});
	    
	    # if > minimum threshold and < previous smallest : save
	    if(($angsep>$thresh)&&($angsep<$minangsep))
	    {
		$minangsep=$angsep;
		$nearest=$j;
	    }
	}
	
	if(defined($nearest))
	{
	    # convert to degrees
	    $minangsep *= 360.0/pi2;

	    #print "star $i/$n : nearest is $nearest at $minangsep degrees\n";

	    # log10
	    $minangsep=log10($minangsep);

	    # bin
	    $minangsep=int($minangsep/$binwidth+0.5)*$binwidth;
	    
	    # save in the histogram
	    $histogram{$minangsep}++;
	}
	else
	{
	    print "Error : failed to find nearest star ($i)\n";
	}
    }
    
    # return histogram data
    return \%histogram;
}