use strict;
use Math::Complex qw(:trig);

sub sqr($) {
  my $x = shift;
  $x * $x;
}

my @spheres = ();
if (0) {
  my $N = 3;
  for my $i (0..$N-1) {
    push @spheres, { 'x' => 0.5 * cos($i*2*pi()/$N),
                     'y' => 0,
		     'z' => 0.5 * sin($i*2*pi()/$N),
		     'radius' => 0.7 };
  }
} else {
  my $r0 = 0.3;
  my $D = 1.3;
#  srand 1;
  for my $i (1..6000) {
      my $rnd = rand(1);
      my $radius = $rnd > 0 ? $r0 * exp(1-1/$rnd) : 0;
      push @spheres, { 'x' => 30 * (rand(1) - 0.5),
		       'y' => $radius * (rand(2)-1),
		       'z' => 30 * (rand(1) - 0.5),
		       'radius' => $radius };
  }
}

sub absv ($) {
  my $p = shift;
  sqrt(sqr($p->{'x'}) +
       sqr($p->{'y'}) +
       sqr($p->{'z'}));
}

sub addv ($$) {
  my $p1 = shift;
  my $p2 = shift;
  { 'x' => $p2->{'x'} + $p1->{'x'},
    'y' => $p2->{'y'} + $p1->{'y'},
    'z' => $p2->{'z'} + $p1->{'z'} };
}

sub subv ($$) {
  my $p1 = shift;
  my $p2 = shift;
  { 'x' => $p2->{'x'} - $p1->{'x'},
    'y' => $p2->{'y'} - $p1->{'y'},
    'z' => $p2->{'z'} - $p1->{'z'} };
}

sub mulv ($$) {
  my $p1 = shift;
  my $p2 = shift;
  $p2->{'x'} * $p1->{'x'} +
  $p2->{'y'} * $p1->{'y'} +
  $p2->{'z'} * $p1->{'z'};
}

sub mulvs ($$) {
  my $p = shift;
  my $m = shift;
  { 'x' => $p->{'x'} * $m,
    'y' => $p->{'y'} * $m,
    'z' => $p->{'z'} * $m };
}

sub distv ($$) {
  my $p1 = shift;
  my $p2 = shift;
  absv(subv($p1, $p2));
}


open(BUBBLES, ">bubbles.inc") or die "bubbles.inc: $!";
print BUBBLES "#declare Bubbles = object { union {\n";
for my $i (0..$#spheres) {
    my $s = $spheres[$i];
    my @planes = ();
    for my $j (0..$#spheres) {
        next if $i == $j;
        my $r1 = $spheres[$i]->{'radius'};
        my $r2 = $spheres[$j]->{'radius'};
        my $d = distv($spheres[$i], $spheres[$j]);
#        print "$i,$j: r1=$r1, r2=$r2, d=$d\n";
        if ($d > abs($r1 - $r2) &&
            $d < $r1 + $r2) {
#            print "spheres $i and $j intersect\n";
            my $d1 = - (sqr($r2) - sqr($r1) - sqr($d))/(2 * $d);
            my $v = subv($spheres[$i], $spheres[$j]);
            my $h = addv($spheres[$i], mulvs($v, $d1/$d));
            my $n = mulvs($v, -1/$d);
	    $n->{'d'} = mulv($n, $h);
            push @planes, $n;
        }
    }
    if(@planes) {
        print "sphere $i intersects with ", scalar @planes, " sphere(s).\n";
	print BUBBLES "    difference {\n        sphere { <$s->{x}, $s->{y}, $s->{z}>, $s->{radius} }\n";
        foreach my $n (@planes) {
	    print BUBBLES "        plane { <$n->{x}, $n->{y}, $n->{z}>, $n->{d} }\n";
     	}
	print BUBBLES "    }\n";
    } else {
        print BUBBLES "    sphere {\n        <$s->{x}, $s->{y}, $s->{z}>, $s->{radius}\n    }\n";
    }
}
print BUBBLES "    }\n    hollow\n    no_shadow\n}\n";
close BUBBLES;