#!/usr/bin/perl -w

# Plots the bifurcation diagram of the logistic function.
# x -> r x (1-x)

use Getopt::Long; # defines GetOptions
use File::Basename; # defines basename

if ( !GetOptions("N1=i",\$N1, # number of iterations in advance
		 "N2=i",\$N2, # number of iterations to plot
		 "rs=f",\$rs, # start value of parameter r
		 "re=f",\$re, # end value of parameter r
		 "resolution=i",\$R, # number of parameters to plot
		 "help",\$help)
     || defined($help) || !defined($N1) || !defined($N2)
     || !defined($rs) || !defined($re) || !defined($R)  ) {
  $myname=basename($0);
  print STDERR "usage: $myname --N1=N1 --N2=N2 --rs=R_START re=R_END --resolution=R [--help]\n";
  exit 1;
}


# logistic function
sub logistic($) {
  my $x=shift @_;
  return $r*$x*(1-$x);
}


# Plot all limit sets for $r in [$rs,$re].
for ($r=$rs; $r<=$re; $r+= ($re-$rs)/$R ) {

  # We have to try several start values because there may be several
  # separate limit sets.
  for ($i=0; $i<5; ++$i) {
    $x=rand(); # random value in [0,1] as start value

    # Iterate the logistic function to get a stable limit set.
    for ($t=0; $t<$N1; ++$t) {
      $x = logistic($x);
    }

    # Iterate the logistic function and plot the limit set.
    for ($t=0; $t<$N2; ++$t) {
      print $r." ".$x."\n";
      $x = logistic($x);
    }
  }

}
