#!/afs/athena/contrib/perl5/perl

unshift(@INC, "/mit/mkgray/lib");

use Plotter;

require Gaussian;
import Gaussian;

$N = 100;
$alpha = $ARGV[0];
$alpha = .1 unless $alpha;
srand(17);

$pl = new Plotter(500,200);
$pl->addPlot("p", 0);

for $i (0..($N-1)) {
    $j[$i] = &gauss_rand;
}

for $i (0..($N-1)) {
    if(rand() > .5){
	$spin[$i] = 1;
    }
    else{
	$spin[$i] = -1;
    }
}

$emin = &emin;
print("Emin: $emin\n");
$oldenergy = &energy;
while(1){
    $t++;
    $flip = int(100*rand());
    $spin[$flip] = ($spin[$flip]==1) ? -1:1;
    $energy = &energy;
    if($energy > $oldenergy) {
	$p = exp(-1*$alpha*$t*($energy-$oldenergy));
	$c = rand();
	if($p < $c){
	    $spin[$flip] = ($spin[$flip]==1) ? -1:1;
	}
    }
    $oldenergy = &energy;
    $pl->plot("p", -1*(($oldenergy-$emin)/$emin));
    $pl->update();
#    print("$oldenergy ($emin)\n");
}

sub emin {
    my $emin = 0;
    for $i (0..($N-1)){
	$emin+= abs($j[$i]);
    }
    return -1*$emin;
}

sub energy {
    my $energy=0;
    for $i (0..($N-2)){
	$energy += $j[$i]*$spin[$i]*$spin[$i+1];
    }
    $energy += $j[$N-1]*$spin[0]*$spin[$N-1];
    return -1*$energy;
}