#!/usr/bin/perl

my $occommand="oc";

my %cmdline; #this will contain all command line options, hashed by the option name
for (my $i=0; $i<@ARGV; $i++) {
    my @opt = split /=/, $ARGV[$i];
    if ( @opt == 1 ) {$opt[1]=1}
    $cmdline{$opt[0]}=$opt[1];
}

if ( exists $cmdline{"-h"} ) {
    print "Usage: phasenorm -tdb1=tdbfile1 -tdb2=tdbfile2 -n=nb_of_samples -T0=min_temperature -T1=max_temperature -e=element,element,...\n";
    print "";
    exit;
}

for my $tag (("-n","-T0","-T1","-tdb1","-tdb2","-e")) {
    if ( ! exists $cmdline{$tag} ) {
	print "Missing " . $tag  . ". Use -h for help\n";
	exit;
    }
}

sub capitalize_first {
    return uc(substr($_,0,1)) . lc(substr($_,1));
}

sub is_same_phase {
    my ($pha,$phb)=@_;
    $pha =~ s/_AUTO.*//g ;
    $phb =~ s/_AUTO.*//g ;
    if ($pha eq $phb) {return 1;} else {return 0;}
}

#main

my $maxn=$cmdline{"-n"};
my $tempmin=$cmdline{"-T0"};
my $tempmax=$cmdline{"-T1"};
my @tdbfiles=($cmdline{"-tdb1"},$cmdline{"-tdb2"});
my @elem=split /[^a-zA-Z]/ , $cmdline{"-e"};
@elem=map(capitalize_first,@elem);
my @varelem=@elem;
my $filelem=shift @varelem;
my $myseed=srand();
my $nbproc=1;
if (exists $cmdline{"-np"}) {
  $nbproc=$cmdline{"-np"}
}
for (my $fi=0; $fi<2; $fi++) {
    srand($myseed);
    for (my $i=0; $i<$nbproc; $i++) {
	open(my $fh,">","tmp_${i}_${fi}.OCM");
	print $fh "read tdb $tdbfiles[$fi]\n";
	print $fh join(" ",@elem) . "\n\n\n";
	#    print $fh "set num\n50\n1e-6\n";
	#    print $fh "set adv dense\n";
	print $fh "set bit global 2\n";
	print $fh "set bit global 8\n";
	my $tiny=2e-4;
	for (my $n=0; $n<$maxn/$nbproc; $n++) {
	    my %x;
	    do {
		$x{$filelem}=1;
		for my $e (@varelem) {
		    $x{$e}=$tiny+(1-2*$tiny)*rand(1);
		    $x{$filelem}-=$x{$e};
		}
	    } while ($x{$filelem} < 0);
	    my $temp=$tempmin+($tempmax-$tempmin)*rand(1);
	    print $fh "set cond t=" . $temp . " p=1e5 n=1 " . join("",(map("x(" . $_ . ")=" . sprintf("%1.4f",$x{$_}) . " " , @varelem))) . "\n";
	    print $fh "calc equil\n\n";
	    print $fh "list\n\n\n";
	}
	print $fh "fin\n";
	close($fh);
    }
}
for (my $fi=0; $fi<2; $fi++) {
    my $ocpara="";
    for (my $i=0; $i<$nbproc; $i++) {
	$ocpara = $ocpara . "$occommand < tmp_${i}_${fi}.OCM > tmp_${i}_${fi}.oco & ";
    }
    system($ocpara . " wait");
    my $catcmd="cat ";
    for (my $i=0; $i<$nbproc; $i++) {
	$catcmd=$catcmd . "tmp_${i}_${fi}.oco ";
    }
    system($catcmd . " > alltmp_${fi}.oco");
    for (my $i=0; $i<$nbproc; $i++) {
	unlink("tmp_${i}_${fi}.oco");
    }
}

my %phaseindex;
my @revphaseindex;
my $curphaseindex=0;
my @bothphfrac;
for (my $fi=0; $fi<2; $fi++) {
    open(my $oh,"<","alltmp_${fi}.oco");
    my $bad=0;
    my @phfrac;
    my $nbpt=0;
    while (<$oh>) {
	if ( index($_,"not a valid equilibrium") != -1 || index($_,"Error") != -1 ) {$bad=1;}
	my @c=split /\s+/ , $_;
	if ( $c[0] eq "Name" ) {
	    $nbpt++;
	    my @allfrac;
	    my $nbphase=0;
	    while (1) {
		my $line=<$oh>;
		my @c=split /\s+/, $line;
		if ( $c[0] eq "" || $c[0] =~ /OC.*/ ) {last;}
		my $phase=$c[0];
		$phase =~ s/\.*$//g;
#		$phase =~ s/_AUTO.*//g;
		if ( ! exists $phaseindex{$phase} ) {
		    $curphaseindex++;
		    $phaseindex{$phase}=$curphaseindex;
		    $revphaseindex[$curphaseindex]=$phase;
		}
		$allfrac[$phaseindex{$phase}]=$c[2];
		print $phase," ",$c[2]," ";
		$line=<$oh>;
		$line=<$oh>;
	    }
	    if ( $bad == 0 ) {
		for (my $i=1; $i<=$curphaseindex ; $i++) {
		    $allfrac[$i]=$allfrac[$i]+0;
		}
		print "\n";
		for (my $s=1; $s<=$curphaseindex ; $s++) {
		    for (my $i=1; $i<=$curphaseindex ; $i++) {
			for (my $j=$i+1; $j<=$curphaseindex ; $j++) {
			    if (is_same_phase($revphaseindex[$i],$revphaseindex[$j])) {
				if ($allfrac[$i]<$allfrac[$j]) {
				    my $tmp=$allfrac[$j];
				    $allfrac[$j]=$allfrac[$i];
				    $allfrac[$i]=$tmp;
				}
			    }
			}
		    }
		}
#		for (my $i=1; $i<=$curphaseindex ; $i++) {
#		    print $allfrac[$i]," ";
#		}
#		print "\n";
		$phfrac[$nbpt]=\@allfrac;
	    }
	    else {
		@allfrac=();
		$phfrac[$nbpt]=\@allfrac;
	    }
	    $bad=0;
	}
    }
    $bothphfrac[$fi]=\@phfrac;
}

my $num=0;
my $den=0;
my $nbpt=@{@bothphfrac[0]};
for (my $i=1; $i<$nbpt; $i++) {
    my @v1=@{@{$bothphfrac[0]}[$i]};
    my @v2=@{@{$bothphfrac[1]}[$i]};
    if ( @v1 != 0 && @v2 != 0) {
	for (my $p=1; $p<=$curphaseindex; $p++) {
	    $num+=abs($v1[$p]-$v2[$p]);
	    print $v1[$p]+0,"-",$v2[$p]+0," ";
	    $den++;
	}
	print "\n";
    }
}

print $num/$den , " " , $den ,"\n";
