#!/usr/bin/perl -w 
use PDL;
use PDL::Graphics::PGPLOT::Window;
use PDL::NiceSlice;
use PDL::Bad;
use PDL::Math;
use strict;
 
my $pi=2*acos(0);

my $p=0.7;
my $q=1-$p;
my $nmax=100;

my $N=pdl [$nmax];


my $n=sequence($nmax+1);

my $P=$p**$n*$q**($N-$n)*(&fact($N)/(&fact($n)*&fact($N-$n)));

my $ncont=sequence(100*$nmax+1)/100;
my $Pgauss=(1/sqrt(2*$pi*$nmax*$p*$q))*exp(-($ncont-$p*$N)**2/(2*$N*$p*$q));

#my $Pgauss=(1/sqrt(2*$pi*$nmax*$p*$q))*exp(-($n-$p*$N)**2/(2*$N*$p*$q));


my $win=PDL::Graphics::PGPLOT::Window->new(Device => '/ps', Aspect => 0.6);

$win->line($ncont,$Pgauss);
$win->hold;
#$win->points($n,$P,{SYMBOL => 'STAR'});
$win->points($n,$P);

#while (!<>) {}
$win -> close;

sub fact {
    my $x=rint(shift(@_)) ;
    my $res=$x ->copy;
    foreach my $i (0..($x->nelem-1)) { 
	if ($x($i)==0) {
	    $res($i) .= 1;
	    next;
	}
	while ($x($i)>1) {
	    $x($i)--;
	    $res($i) .= $res($i)*$x($i);
	}
    }
    return $res;
}