#!/usr/bin/perl	

use lib $ENV{HOME}.'/lib/perl';
use Param;

my $param = Param::Parse
(
	{'info','einstein modell'},
	'quanta',1000,'q','Energie des Systems in Portionen',
	'impulse',2000,'g','Photonenimpulszustände (hier fließt Volumen und Frequenzbreite ein)',
	'atome',1000,'n','Anzahl Atome gesamt',
	'grund',undef,'G','Anzahl Atome im Grundzustand (wenn undefiniert, dann wir mit wahrscheinlichstem Wert angefangen)',
	'zeit',1000,'t','Anzahl Zeitschritte'
);

my $nq = $param->{quanta};
my $g = $param->{impulse};
my $n = $param->{atome};
my $ng = defined $param->{grund} ? $param->{grund} : ProbVal();

print STDERR "Parameter:\n nq = $nq; g = $g; n = $n; ng = $ng\n";
print STDERR "\$\\bar n_g = ", &ProbVal(), "\$\n";


open(DAT, ">tmp.dat") or die;

print DAT "#\$n_g = $ng\$\n";
print DAT "#Parameter:\n# nq = $nq; g = $g; n = $n; ng = $ng\n";
print DAT "#\$\\bar n_g = ", &ProbVal(), "\$\n";
print DAT "#Zeit\t\$n_g\$\n";

for(my $t = 0; $t < $param->{zeit}; ++$t)
{
	my $r = rand();
	if($r > 0.5)
	{
		$r = 2*$r - 1; # ($r-0.5)/0.5
		my $pa = $ng * ($nq - $n + $ng);
		my $pe = ($n-$ng) * ($g + $nq - $n + $ng);
		if( $pa+$pe > 0 and $r < $pa/($pa+$pe))
		{
			$ng -= 1 unless $pa == 0;
		}
		else{ $ng += 1 unless $pe == 0; }
	}
	print DAT "$t\t$ng\n";
}

close(DAT);
open(DAT, "|gnuplot -persist -");
print DAT "unset key\nset xlabel 'Zeit'\nset yrange[0:",$n+1,"]\nset ylabel 'Atome im Grundzustand'\nset grid\nplot 'tmp.dat' w l\nq\n";
close(DAT);

sub ProbVal
{
	my $p = (3*$n-2*$nq-$g);
	my $n_1 = 0.25*( $p + sqrt($p**2 + 8*$n*($g+$nq-$n)) );
	my $n_2 = 0.25*( $p - sqrt($p**2 + 8*$n*($g+$nq-$n)) );
	$nm = ($n_1 >= 0 and  $nq - $n + $n_1 >= 0 ? $n_1 : $n_2);
	return($nm - int($nm) >= 0.5 ? int($nm) + 1 : int($nm)  );
}