#!/usr/bin/perl
use PDL;
use PDL::NiceSlice;

$pi = 4*(atan 1.0);

sub create_sample {
	my $n = shift;
	my ($x, $y, $z, $w, $t);

	$z = 2.0 * random($n) - 1.0;
	$t = 2.0 * $pi * random($n);
	$w = sqrt(1 - $z*$z);
	$x = $w * cos($t);
	$y = $w * sin($t);

	return ($x, $y, $z);
}

$npoints = 1000;
$n = $ARGV[0];
$torusangle = $ARGV[1];

(defined $n) or $n=40;
print "Calculating statistics for N=$n\n";

(defined $torusangle) or $torusangle = 40;
print "Opening angle: $torusangle degrees.\n";
$torusangle *= $pi/180.0;

($x, $y, $z) = create_sample($npoints);
$psz = sqrt($x*$x+$y*$y);
$theta = $pi/2.0 - atan(abs $z/$psz);
$i = which ($theta >= $torusangle);
$n_nlrg = $i->nelem;

while ($n_nlrg % $n) {
	($tmpx, $tmpy, $tmpz) = create_sample(10);
# PROBLEM SEEMS TO OCCUR IN THIS AREA
########################################
	$x = pdl($x,$tmpx)->flat;
	$y = pdl($y,$tmpy)->flat;
	$z = pdl($z,$tmpz)->flat;
########################################
	$psz = sqrt($x*$x+$y*$y);
	$theta = $pi/2.0 - atan(abs $z/$psz);
	$i = which ($theta >= $torusangle);
	$n_nlrg = $i->nelem;
}

$npoints = $x->nelem;
print $npoints," objects created.\n";