#!/usr/bin/perl -w
#-- FILE:   randtree.pl
#-- CODER:  KUBO Takuya (kubo@ees.hokudai.ac.jp)
#-- 
#-- (C)opyright 2001 KUBO Takuya (kubo@ees.hokudai.ac.jp).
#--    All rights are reserved.
#--    Do not modify without amending copyright,
#--    distribute freely but retain this paragraph.
#-- 
#-- CHANGE: Sun Aug 26 16:10:59 2001
#   MODIFY: Sun Aug 26 01:36:41 2001

use strict;
use vars qw( $PI %parameter );
$PI = 3.141592;
%parameter = (
	"time_step_max" => 99,
	"start_branching" => 9,
	"length_modifier" => 0.7,
	"branch_angle_range" => 0.5,
	"theta_modifier_range" => 5.0e-3,
	"branching_rate_modifier" => 0.4
);
# initialization for the primary terminal
sub rand_zero {
	my ( $range ) = @_;
	return $range * ( rand( 1.0 ) - 0.5 );
}
my $primary_term = {
	"x" => 1.0,
	"y" => 1.0,
	"theta" => 0.5 * $PI,
	"length" => 1.0,
	"trajectory" => [ [ 1.0, 1.0 ] ],
	"branching_rate" => 0.1,
	"theta_modifier" => 1.0 + rand_zero( 1.0e-3 )
};
# branching from parent
sub child_of {
	my ( $parent_term ) = ( @_ );
	my $child_term = {
		"x" => $$parent_term{ "x" },
		"y" => $$parent_term{ "y" },
		"length" => $$parent_term{ "length" }
			* $parameter{ "length_modifier" },
		"trajectory" =>
			[ [ $$parent_term{ "x" } , $$parent_term{ "y" } ] ],
		"branching_rate" => $$parent_term{ "branching_rate" }
			* $parameter{ "branching_rate_modifier" }
	};
	$$child_term{ "theta" } = $$parent_term{ "theta" }
		* ( 1 + rand_zero( $parameter{ "branch_angle_range" } ) );
	$$child_term{ "theta_modifier" }
		= 1 + rand( $parameter{ "theta_modifier_range" } )
			* ( ( $$child_term{ "theta" } > 0.5 * $PI ) ? 1 : -1 );
	return $child_term;
}

# main loop
warn "# Starting randtree ...\n";
my %copy_primary = %$primary_term;
my @terminals;
push @terminals, \%copy_primary;
foreach my $time_step ( 0 .. $parameter{ "time_step_max" } ) {
	my @new_terminals;
	foreach my $term ( @terminals ) {
		if (
			rand( 1.0 ) < $$term{ "branching_rate" }
			and $time_step >= $parameter{ "start_branching" }
		) {
			push @new_terminals, child_of( $term );
		}
	}
	foreach my $new_term ( @new_terminals ) {
		push @terminals, $new_term;
	}
	foreach my $term ( @terminals ) {
		$$term{ "x" } += $$term{ "length" } * cos( $$term{ "theta" } );
		$$term{ "y" } += $$term{ "length" } * sin( $$term{ "theta" } );
		push @{ $$term{ "trajectory" } },
			[ $$term{ "x" }, $$term{ "y" } ];
		$$term{ "theta" } *= $$term{ "theta_modifier" };
	}
}

# output
my ( $output_terminal, $output_trajectory )
	= ( "randtree_terminal.txt", "randtree_trajectory.txt" );
warn "# terminal locations\t-> \"$output_terminal\"\n";
open TERMINAL, "> $output_terminal";
foreach my $term ( @terminals ) {
	printf TERMINAL "%.2e\t%.2e\t%.2e\n",
		$$term{ "x" }, $$term{ "y" }, $$term{ "theta" };
}
close TERMINAL;
warn "# term trajecctories\t-> \"$output_trajectory\"\n";
open TRAJECTORY, "> $output_trajectory";
foreach my $term ( @terminals ) {
	my $xy;
	while ( $xy = shift @{ $$term{ "trajectory" } } ) {
		printf TRAJECTORY "%.2e\t%.2e\n", $$xy[ 0 ], $$xy[ 1 ];
	}
	print TRAJECTORY "\n";
}
close TRAJECTORY;

# generating a gnuplot script
my ( $x_min, $x_max, $y_min, $y_max ) = (
	$$primary_term{ "x" } - 0.5
		* $$primary_term{ "length" } * $parameter{ "time_step_max" },
	$$primary_term{ "x" } + 0.5
		* $$primary_term{ "length" } * $parameter{ "time_step_max" },
	$$primary_term{ "y" },
	$$primary_term{ "y" }
		+ $$primary_term{ "length" }
			* ( $parameter{ "time_step_max" } + 5 )
);
my $gnuplot_file = "randtree.gp";
warn "# gunuplot driver\t-> \"$gnuplot_file\"\n";
open GNUPLOT, "> $gnuplot_file";
my $png = "randtree.png";
print GNUPLOT <<EOF;
set term png color
set size square 0.5
set output '$png'
set nokey
set noxtics
set noytics
set noborder
set xrange [$x_min:$x_max]
set yrange [$y_min:$y_max]

plot '$output_trajectory' u 1:2 w l 3, \\
	'$output_terminal' u 1:2 w p 1
print "# $png is generated";
EOF
close GNUPLOT;

exit 0;
