$max_runs = 100; # Number of trials to average over per configuration

( $min_edge, $max_edge ) = ( 5, 65 ); # Min and max lattice size

for $c ( 1..15 ) {
  $cross_section = 0.01 * $c;

  for $edge ( $min_edge..$max_edge ) {

    ( $qq, $tt, $gg, $t2, $g2 ) = ( 0., 0, 0, 0, 0 );
    for $q ( 0..$max_runs ) {
      $qq ++;
      ( $t, $g ) = run_simulation( $edge, $cross_section );
      $tt += $t;
      $gg += $g;
      $t2 += $t*$t;
      $g2 += $g*$g;
    }

    $t  = $tt/$qq;		# Avg final timestep
    $g  = $gg/$qq;		# Avg final remaining atoms
    $gd = $g2/$qq - $g*$g;	# Std dev of remaining atoms

    $gamma = $edge * $cross_section; # Scaled "coverage" parameter

    $out  = "$edge\t$cross_section\t$gamma\t$t\t";
    $out .= $g/$edge**2 . "\t" . sqrt($gd)/$edge**2 . "\n";

    print $out;
    print STDERR $out;
  }

  print "\n\n";
  print STDERR '='x50, "\n";
}

# Exits when there are no more excited atoms, 
# returns timestep and number of atoms in groundstate
sub run_simulation {
  my ( $edge, $cross_section ) = @_;

  my $crystal = make_crystal( $edge );
  my $excited = [ [ int($edge/2), int($edge/2) ] ]; # Excited seed atom

  my ( $t, $r ) = ( 0, 0 );
  while( 1 ) {
    $t++;

    $excited = explode( $crystal, $excited, $cross_section );

    $r = analyze_crystal( $crystal );

#     print $t, "\t", scalar( @$excited ), "\t|\t";
#     for my $k ( 'groundstate', 'excited', 'decayed' ) {
#       print $r->{$k} || 0, "\t";
#     }
#     print "\n";

    if( scalar( @$excited ) == 0 ) { last; }
  }

  return ( $t, $r->{groundstate} || 0 );
}

sub explode {
  my ( $crystal, $excited, $cross_section ) = @_;

  my @newly_excited = ();

  foreach my $e ( @$excited ) {   
    my ( $ex, $ey ) = @$e;

    $crystal->[$ex][$ey] = 'decayed';
  
#   my $f1 = propagate_neutron( $crystal, $ex, $ey, $cross_section );
#   my $f2 = propagate_neutron( $crystal, $ex, $ey, $cross_section );

    # Generate exactly two neutrons per decaying atom...
    my $f1 = propagate_fast( $crystal, $ex, $ey, $cross_section );
    my $f2 = propagate_fast( $crystal, $ex, $ey, $cross_section );

    if( defined $f1 ) { push @newly_excited, $f1; }
    if( defined $f2 ) { push @newly_excited, $f2; }
  }

  return \@newly_excited;
}

# For the excited atom given, shoot out a neutron in a random direction.
# If the neutron passes within the "cross-section" of another one, this
# atom will be "excited" and decay in the next time step. If there are 
# multiple candidates, select the one closest to the atom decaying in the
# current time step.
# Note: This procedure iterates over the whole crystal, and counts TWO
# directions for the neutron: forward and backward (this is arguably a
# bug).
sub propagate_neutron {
  my ( $crystal, $ex, $ey, $cross_section ) = @_;

  my ( $s, $c ) = random_direction();

  my @ee = ();
  my ( $n, $rr ) = ( scalar( @$crystal )-1, 0 );

  my $r2 = 2.0*$n*$n; # length of crystal diagonal (squared)
  for my $y ( 0..$n ) {
    for my $x ( 0..$n ) {
      unless( $crystal->[$x][$y] eq 'groundstate' ) { next; }
      
      my $beta = - $s * ( $x - $ex ) + $c * ( $y - $ey );
      if( abs( $beta ) < $cross_section ) {
	$rr = ( $x-$ex )**2 + ( $y-$ey )**2; 

	if( $rr < $r2 ) { 
	  @ee = ( $x, $y ); # keep the ($x,$y) with the smallest $rr
	  $r2 = $rr;
	}
      }
    }
  }

  if( scalar( @ee ) > 0 ) {
    $crystal->[$ee[0]][$ee[1]] = 'excited';
    return \@ee; 
  }
  return;
}

# The same as above, but using a much faster algorithm, since it does
# not iterate over the entire crystal, but only over those atoms which
# are reasonable close to the generated neutron. Also, takes into account
# only a single direction for the neutron. (No forward/backward.)
sub propagate_fast {
  my ( $crystal, $ex, $ey, $cross_section ) = @_;

  my $n = scalar( @$crystal )-1;
  my ( $s, $c ) = random_direction();

  if( abs($s) < abs($c) ) { # step in x-direction 
    $direction = 'X';

    $dx = ( $c > 0 ? 1 : -1 );
    $dy = $s/abs($c);

  } else { # step in y-direction
    $direction = 'Y';

    $dx = $c/abs($s);
    $dy = ( $s > 0 ? 1 : -1 );
  }

  my @ee = ();
  my ( $rr, $r2 ) = ( 0, 2.0*$n*$n ); # length of crystal diagonal (squared)   
  
  for( $x=$ex, $y=$ey; ; $x+=$dx, $y+=$dy ) {
    if( $x < 0 || $y < 0 || $x > $n || $y > $n ) { last; }

    if( $direction eq 'X' ) {
      ( $x1, $x2 ) = ( $x, $x );                 # Stepping in x-direction
      ( $y1, $y2 ) = ( int( $y ), int( $y+1 ) ); # The atom "above" and "below"
    } else {
      ( $x1, $x2 ) = ( int( $x ), int( $x+1 ) ); # same for y-direction
      ( $y1, $y2 ) = ( $y, $y );
    }

    # The two nearest candidates...
    unless( $crystal->[$x1][$y1] eq 'groundstate' ) { next; }
    unless( $crystal->[$x2][$y2] eq 'groundstate' ) { next; }
    
    # Calculate beta for the first candidate...
    my $beta = - $s * ( $x1 - $ex ) + $c * ( $y1 - $ey );
    if( abs( $beta ) < $cross_section ) {
      $rr = ( $x1-$ex )**2 + ( $y1-$ey )**2; 

      if( $rr < $r2 ) { 
	@ee = ( $x1, $y1 );
	$r2 = $rr;
      }
    }
    
    # ... and for the second.
    $beta = - $s * ( $x2 - $ex ) + $c * ( $y2 - $ey );
    if( abs( $beta ) < $cross_section ) {
      $rr = ( $x2-$ex )**2 + ( $y2-$ey )**2; 

      if( $rr < $r2 ) { 
	@ee = ( $x2, $y2 );
	$r2 = $rr;
      }
    }
  }

  if( scalar( @ee ) > 0 ) {
    $crystal->[$ee[0]][$ee[1]] = 'excited';
    return \@ee; 
  }
  return;
}

# Returns a hash-ref with counts of the atoms in all states: 
# groundstae, excited, decayed
sub analyze_crystal {
  my ( $crystal ) = @_;

  my %results;

  foreach my $row ( @$crystal ) {
    foreach my $atom ( @$row ) {
      $results{ $atom } ++;
    }
  }

  return \%results;
}

sub make_crystal {
  my ( $n ) = @_;

  my @crystal = ();
  for my $x ( 0..$n-1 ) {
    push @crystal, [];
    for my $y ( 0..$n-1 ) {
      push @{$crystal[$x]}, 'groundstate';
    }
  }

  return \@crystal;
}

sub random_direction {
  my $pi  = 3.14159265358979323846;
  my $phi = rand( 2.0*$pi );
  return ( sin( $phi ), cos( $phi ) );
}