#!/usr/bin/perl

our $occommand="oc";
our $mkaxes="mkaxes";
our $simplexize="simplexize";
our $mpirun="mpirun";

my $configfile=$0;
$configfile =~ s/ocplotpd/.ocplotpd/g;
if ( ! -f $configfile) {
    $configfile=$ENV{"HOME"} . "/.ocplotpd";
    if ( ! -f $configfile ) {
	$configfile="";
    }
}
if ( length($configfile)>0 ) {
    require $configfile
}

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];
}

my $pi=3.1415927;

if ( exists $cmdline{"-h"} ) {
    print
	"A front-end to OpenCalphad for PLOTting Phase Diagrams.\n" .
	"Usage: ocplotpd [-ha] [-nt] [-tlf=fraction] [-1f] [-vtk] [-rmax=r_max] [-rmin=r_min] [-np=nb_of_proc] -tdb=tdbfile -n=nb_of_samples -T0=min_temperature -T1=max_temperature -e=element,element,...\n" .
	"\n" .
	"Mandatory parameters:\n" .
	" tdbfile: name of the Thermodynamic DataBase (TDB) file to use\n" .
	" element,element,... : a list of the elements (and associated phases) to extract from the TDB file\n" .
	" nb_of_samples: Number of phase equilibria sampled to generate the phase diagram\n" .
	" min_temperature: lower bound on the temperature axis\n" .
	" max_temperature: upper bound on the temperature axis (if min_temperature=max_temperature, a cross section is plotted)\n" .
	" \n" .
	"optional parameters:\n" .
	" -ha:   high accuracy calculations\n" .
	" -nt:   no tie lines\n" .
	" -tlf:  fraction of tie lines displayed (between 0.0 and 1.0)\n" .
	" -1f:   generate a single gnuplot file\n" .
	" -vtk:  generate vtk files (instead of gnuplot)\n" .
	" -rmax: max distance between points in triangulation\n" .
	" -rmin: discard points closer to each other than this cutoff\n" .
	" -np:   launch nb_of_proc copies of opencalphad\n" .
	" -noc:  do not run OpenCalphad - just create the graphs from the *.dat files\n" .
	" -nd:   general n-dimensional phase diagram (generate input files for ndviewer - experimental)\n" .
	" -debug:do not delete temporary files\n" .
	" -mpi:  use mpi for simplexize code (experimental)\n" .
	"";
    exit;
}

if (exists $cmdline{"-cut"} ) {
    my $cutfile=$cmdline{"-cut"};
    $cutfile =~ s/$/:/g;
    $cutfile =~ s/:+/:/g;
    $cutfile =~ s/:/\n/g;
    $cutfile =~ s/,/ /g;
    open(my $fh,">","cut.in");
    print $fh $cutfile;
    $cmdline{"-T0"}=0;
    $cmdline{"-T1"}=1;
}

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

if (system("echo fin > fin.OCM ; $occommand fin.OCM >& /dev/null")!=0) {
    print
	"Cannot find the OpenCalphad executable ${occommand}\n" .
	"Please install OpenCalphad from http://www.opencalphad.com/ and name the executable ${occommand} and make sure it is in your path\n";
    unlink("fin.OCM");
    exit;
}
unlink("fin.OCM");

my $pttoofar=0.2;
if ( exists $cmdline{"-rmax"} ){
    $pttoofar=$cmdline{"-rmax"};
}
if ( exists $cmdline{"-cc"} ){
    $pttoofar=$cmdline{"-cc"};
}

my $pttooclose=1e-2;
if ( exists $cmdline{"-rmin"} ){
    $pttooclose=$cmdline{"-rmin"};
}

my $tielinefrac=1.0;
if ( exists $cmdline{"-tlf"} ){
    $tielinefrac=$cmdline{"-tlf"};
}

sub pddist {
    my ($deltatemp,$pta,$ptb)=@_;
    if ($deltatemp == 0) {$deltatemp=1;}
    my @vpta=split /\s+/, $pta;
    my @vptb=split /\s+/, $ptb;
    if ($vpta[0] ne $vptb[0]) {
	return 100.;
    }
    my $dist=0;
    for (my $i=0; $i<@vpta-2; $i++) {
	my $delta=($vpta[3+$i]-$vptb[3+$i]);
	if ($i==0) {$delta=$delta/$deltatemp;}
	$dist+=$delta*$delta;
    }
    $dist=sqrt($dist);
    return $dist;
}

sub reorderpoints {
    my $deltatemp=shift;
    my $maxdist=shift;
    my @raw=split /\n/,$_[0];
    if (@raw==0) {return "\n";}
    my @ord;
    my @done;
    my $numdone=1;
    my $result;
    $ord[0]=$raw[0];
    $done[0]=1;
    while ($numdone < @raw ) {
	my $bestdist=1000.;
	my $bestpt=-1;
	my $bestend=-1;
	for (my $i=1; $i<@raw; $i++) {
	    if ( $done[$i] != 1 ) {
		my $curdistbeg=pddist($deltatemp,$ord[0],$raw[$i]);
		my $curdistend=pddist($deltatemp,$ord[@ord-1],$raw[$i]);
		if ($curdistbeg < $curdistend) {
		    if ($curdistbeg < $bestdist) {
			$bestdist=$curdistbeg;
			$bestpt=$i;
			$bestend=0;
		    }
		}
		else {
		    if ($curdistend < $bestdist) {
			$bestdist=$curdistend;
			$bestpt=$i;
			$bestend=@ord-1;
		    }
		}
	    }
	}
	$done[$bestpt]=1;
	$numdone++;
	if ($bestdist > $maxdist) {
	    $result = $result . "\n\n" . join("\n",@ord);
	    @ord=();
	    $ord[0]=$raw[$bestpt];
	}
	else {
	    if ($bestend==0) {
		unshift @ord, $raw[$bestpt];
	    }
	    else {
		push @ord, $raw[$bestpt];
	    }
	}
    }
    $result = $result . "\n\n" . join("\n",@ord) . "\n";
    return $result;
}

sub copy_file_to_handle {
    my $destfh=shift;
    my $orgfile=shift;
    open(my $fh,"<",$orgfile);
    while (<$fh>) {
	print $destfh $_;
    }
}

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

sub pick_phase {
    my $data=shift;
    my $ph1=shift;
    my $ph2=shift;
    my @lines=split /\n/,$data;
    my $outl=0;
    my @linesout;
    for (my $l=0; $l<@lines; $l++) {
	my $phases=$lines[$l];
	$phases =~ s/ +.*//g;
	$phases= ',' . $phases . ',';
	if ( index($phases,",$ph1,") > -1 && index($phases,",$ph2,") > -1 ) {
	    $linesout[$outl]=$lines[$l];
	    $outl++;
	}
    }
    if (@linesout==0) {
	return "";
    }
    else {
	return join("\n",@linesout)."\n";
    }
}

sub pick_phase1 {
    my $data=shift;
    my $ph2=shift;
    my @lines=split /\n/,$data;
    my $outl=0;
    my @linesout;
    for (my $l=0; $l<@lines; $l++) {
	my $phases=$lines[$l];
	$phases =~ s/ +.*//g;
	$phases= ',' . $phases . ',';
	if ( index($phases,",$ph2,") > -1 ) {
	    $linesout[$outl]=$lines[$l];
	    $outl++;
	}
    }
    if (@linesout==0) {
	return "";
    }
    else {
	return join("\n",@linesout)."\n";
    }
}

sub escape_bad {
  my $cln=shift;
  $cln =~ s/_/-/g;
  return $cln;
}

#main

my $maxn=$cmdline{"-n"};
my $tempmin=$cmdline{"-T0"};
my $tempmax=$cmdline{"-T1"};
my $tdbfile=$cmdline{"-tdb"};
my @elem=split /[^a-zA-Z]/ , $cmdline{"-e"};
@elem=map(capitalize_first,@elem);
my @varelem=@elem;
my $filelem=shift @varelem;
my $nbproc=1;
if (exists $cmdline{"-np"}) {
  $nbproc=$cmdline{"-np"}
}

my %phaseindex;
if (! exists $cmdline{"-noc"}) {
    for (my $i=0; $i<$nbproc; $i++) {
	open(my $fh,">","tmp${i}.OCM");

	#print $fh "new Y\n";
	print $fh "read tdb $tdbfile\n";
	print $fh join(" ",@elem) . "\n\n\n\n";
	if (exists $cmdline{"-ha"}) {
	    print $fh "set num\n400\n1e-7\n1e-4\n-0.005\n";
	    print $fh "set adv grid 2\n";
	}
	#print $fh "set bit global 2\n";
	print $fh "set bit global 8\n";
	print $fh "set bit global 18\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);
	    #my $temp=$tempmin+($tempmax-$tempmin)*$n/($maxn/$nbproc-1);
	    print $fh "set cond t=" . $temp . " p=1e5 n=1\n";
	    print $fh join("",(map("set cond x(" . $_ . ")=" . sprintf("%1.4f",$x{$_}) . "\n" , @varelem)));
	    print $fh "calc only_grid\n\n";
	    print $fh "calc no_global\n\n";
	    #print $fh "calc equil\n\n";
	    print $fh "list result 1\n";
	    #print $fh "delete equil 1\n";
	}
	print $fh "fin\n";
	close($fh);
    }
    my $ocpara="";
    for (my $i=0; $i<$nbproc; $i++) {
	# thanks to Algirdas Baskys for fixing the next line:
	$ocpara = $ocpara . "$occommand tmp${i}.OCM > tmp${i}.oco & ";
    }
    system($ocpara . " wait");
    my $catcmd="cat ";
    for (my $i=0; $i<$nbproc; $i++) {
	$catcmd=$catcmd . "tmp${i}.oco ";
    }
    system($catcmd . " > alltmp.oco");
    if ( ! exists $cmdline{"-debug"} ) {
	for (my $i=0; $i<$nbproc; $i++) {
	    unlink("tmp${i}.oco");
	    unlink("tmp${i}.OCM");
	}
    }
    {
	open(my $oh,"<","alltmp.oco");
	my $curphaseindex=0;
	my $temp=0;
	my $bad=0;
	my @ptbyphase;
	my $tielines;
	my $nbpt=0;
	while (<$oh>) {
	    if ( index($_,"not a valid equilibrium") != -1 || index($_,"Error") != -1 || index($_,"warning") != -1 ) {$bad=1; print "OpenCalphad message: " . $_}
	    my @c=split /\s+/ , $_;
	    if ( $c[0] eq "T=" ) {
		$temp=$c[1];
	    }
	    if ( $c[0] eq "Name" ) {
		$nbpt++;
		my @equil;
		my @allphase;
		my $nbphase=0;
		while (1) {
		    my $line=<$oh>;
		    my @c=split /\s+/, $line;
		    if ( $c[0] eq "" || $c[0] =~ /OC.*/ ) {
			if ($c[0] eq "") {
			    while (1) {
				my $line=<$oh>;
				if ( $line =~ /OC.*/  || eof($oh) ) {last;}
				if ( index($line,"would like to be stable") != -1 ) {$bad=1; print "OpenCalphad message: " . $line}
			    }
			}
			last;
		    }
		    my $phase=$c[0];
		    $phase =~ s/\.*$//g;
		    $phase =~ s/_AUTO.*//g;
		    $phase =~ s/\#[0-9]*//g;
		    if ( ! exists $phaseindex{$phase} ) {
			$curphaseindex++;
			$phaseindex{$phase}=$curphaseindex;
		    }
		    my %conc;
		    while (1) {
			$line=<$oh>;
			chomp $line;
			@c=split /\s+/, $line;
			if (@c == 0) {last;}
			for (my $i=1; $i<@c; $i+=2) {
			    $conc{$c[$i]}=$c[$i+1];
			}
		    }
		    $equil[$nbphase]=$phaseindex{$phase} . " " . $phase . " " . $temp . " " . join(" ",map($conc{uc $_},@elem));
		    $allphase[$nbphase]=$phaseindex{$phase};
		    $nbphase++;
		}
		if ( $bad == 0 ) {
		    if ($nbphase>1) {
			my $dotieline=(rand(1)<=$tielinefrac ? 1:0);
			my $multiphase=join(",",sort @allphase);
			for (my $i=0; $i<$nbphase; $i++) {
			    $ptbyphase[$allphase[$i]] = $ptbyphase[$allphase[$i]] . $multiphase . " " . $equil[$i] . "\n";
			    if ($dotieline) {
				$tielines = $tielines . $multiphase . " " . $equil[$i] . "\n";
			    }
			}
			if ($dotieline) {
			    $tielines = $tielines . "\n";
			}
		    }
		}
		$bad=0;
	    }
	}
	for my $ph (keys %phaseindex) {
	    if (length($ptbyphase[$phaseindex{$ph}]) == 0) {
		delete $phaseindex{$ph};
	    }
	}
	for my $ph (keys %phaseindex) {
	    open(my $phh,">","phase_" . $ph . ".dat");
	    if ( (@elem-1 + ($tempmin == $tempmax ? 0 : 1)) == 2 ) {
		my $orderedpt=reorderpoints($tempmax-$tempmin,$pttoofar,$ptbyphase[$phaseindex{$ph}]);
		print $phh $orderedpt;
	    }
	    else {
		#	    for (my $ph2=0; $ph2<=$curphaseindex; $ph2++) {
		#		for (my $ph1=0; $ph1<=$ph2; $ph1++ ) {
		#		    my $picked_phases=pick_phase($ptbyphase[$phaseindex{$ph}],$ph1,$ph2);
		#		    if (length($picked_phases)>0) {
		#			print $phh $picked_phases;
		#			print $phh "\n";
		#		    }
		#		}
		#	    }
		for (my $ph2=0; $ph2<=$curphaseindex; $ph2++) {
		    
		    my $picked_phases;
		    if ($ph2==$phaseindex{$ph}) {
			$picked_phases=pick_phase1($ptbyphase[$phaseindex{$ph}],"$ph2,$ph2");
		    }
		    else {
			$picked_phases=pick_phase1($ptbyphase[$phaseindex{$ph}],$ph2);
		    }
		    if (length($picked_phases)>0) {
			print $phh $picked_phases;
			print $phh "\n";
		    }
		}
	    }
	    close($phh);
	}
	{
	    open(my $tih,">","tielines.dat");
	    print $tih $tielines;
	    close($tih);
	}
    }
}
else {
    open(my $fh, 'ls -1 phase_*.dat |');
    while (my $line=<$fh>) {
	chomp $line;
	$line =~ s/^phase_//g;
	$line =~ s/.dat$//g;
	print $line,"\n";
	$phaseindex{$line}="present";
    }
}

{
    my $gnh;
    if ( ( ( ! exists $cmdline{"-vtk"} ) && (! exists $cmdline{"-nd"}) ) || @elem==2 || ( @elem==3 && $tempmin == $tempmax) ) {
	open($gnh,">","ocplotpd.gnu");
	print $gnh "";
	print $gnh "set key outside\n";
    }
    my $tielinecolor='lc rgbcolor "light-gray"';
    if (exists $cmdline{"-nt"}) {
	$tielinecolor='lc rgbcolor "white"';
    }
    my $tiefile="tielines.dat";
    if (exists $cmdline{"-1f"}) {
	$tiefile="-";
    }
    if (@elem==2) {
	print $gnh "set xlabel '$elem[1]'\n";
	print $gnh "set ylabel 'T (K)'\n";
	print $gnh "plot [0:1] [] ";
	print $gnh "'" . $tiefile ."' u 6:4 w l $tielinecolor notitle";
	for my $ph (keys %phaseindex) {
	    print $gnh ", \\\n";
	    my $filename="phase_" . $ph . ".dat";
	    if (exists $cmdline{"-1f"}) {
		$filename="-";
	    }
	    print $gnh "'" . $filename . "' u 6:4 w l t '" . escape_bad($ph) . "'";
	}
    }
    elsif (@elem==3 && (! exists $cmdline{"-nd"}) ) {
	if ($tempmin == $tempmax) {
	    print $gnh "unset border\nunset tics\n";
	    print $gnh "set size square 1,1\n";
	    print $gnh "set lmargin 4\n";
	    print $gnh "set rmargin 2\n";
	    print $gnh "set tmargin 4\n";
	    print $gnh "set bmargin 1\n";
	    print $gnh "set label 'T=$tempmin K' at 0.05,0.866\n";
	    print $gnh "set label '$elem[0]' at -0.05,0\n";
	    print $gnh "set label '$elem[1]' at 1.05,0\n";
	    print $gnh "set label '$elem[2]' at 0.5,0.866+0.05\n";
	    print $gnh "set arrow from 0,0 to 1,0 nohead\n";
	    print $gnh "set arrow from 1,0 to 0.5,0.866 nohead\n";
	    print $gnh "set arrow from 0.5,0.866 to 0,0 nohead\n";
	    print $gnh "plot [0:1] [0:0.866] ";
	    print $gnh "'" . $tiefile . "' u (\$6+0.5*\$7):(0.866*\$7) w l $tielinecolor notitle";
	    for my $ph (keys %phaseindex) {
		print $gnh ", \\\n";
		my $filename="phase_" . $ph . ".dat";
		if (exists $cmdline{"-1f"}) {
		    $filename="-";
		}
		print $gnh "'" . $filename . "' u (\$6+0.5*\$7):(0.866*\$7) w l t '" . escape_bad($ph) . "'";
	    }
	}
	else {
	    if (exists $cmdline{"-vtk"}) {
		my $zscale=($tempmax-$tempmin);
		my $allscaling="-zmin=$tempmin -zmax=$tempmax -zscale=$zscale";
		system("${mkaxes} -tri -col=0 $allscaling > axes.vtk");
		system("${mkaxes} -tri -col=0 $allscaling -vl=" . join(',',@elem) . "> labels.vtk");
		if ( ! exists $cmdline{"-nt"} ) {
		    system("awk '{print \$4,\$5,\$6,\$7}' $tiefile | ${mkaxes} -trans -tri $allscaling | ${simplexize} -t=- -l -col=0  > tielines.vtk");
		}
		my $color=1;
		my $cmd="";
		for my $ph (keys %phaseindex) {
		    my $filename="phase_" . $ph . ".dat";
		    my $vtkfilename="phase_" . $ph . ".vtk";
		    $cmd = $cmd . "awk '{print \$4,\$5,\$6,\$7}' $filename | ${mkaxes} -trans -tri $allscaling | ${simplexize} -t=- -rmin=$pttooclose -rmax=$pttoofar -col=$color > $vtkfilename " . ( $nbproc > 1 ?  "& " : "; ");
		    $color++;
		}
		if ($nbproc>1) {$cmd = $cmd . "wait";}
		print $cmd,"\n";
		system($cmd);
	    }
	    else {
		my $tempmid=($tempmin+$tempmax)/2;
		print $gnh "unset border\nunset tics\n";
		print $gnh "set size square 1,1\n";
		print $gnh "set label '$elem[0]' at -0.05,0,$tempmid\n";
		print $gnh "set label '$elem[1]' at 1.05,0,$tempmid\n";
		print $gnh "set label '$elem[2]' at 0.5,0.866+0.05,$tempmid\n";
		print $gnh "set arrow from 0,0,$tempmin to 1,0,$tempmin nohead\n";
		print $gnh "set arrow from 1,0,$tempmin to 0.5,0.866,$tempmin nohead\n";
		print $gnh "set arrow from 0.5,0.866,$tempmin to 0,0,$tempmin nohead\n";
		print $gnh "set arrow from 0,0,$tempmax to 1,0,$tempmax nohead\n";
		print $gnh "set arrow from 1,0,$tempmax to 0.5,0.866,$tempmax nohead\n";
		print $gnh "set arrow from 0.5,0.866,$tempmax to 0,0,$tempmax nohead\n";
		print $gnh "set arrow from 0,0,$tempmin to 0,0,$tempmax nohead\n";
		print $gnh "set arrow from 1,0,$tempmin to 1,0,$tempmax nohead\n";
		print $gnh "set arrow from 0.5,0.866,$tempmin to 0.5,0.866,$tempmax nohead\n";
		print $gnh "set label '$tempmin K' at -0.15,0,$tempmin\n";
		print $gnh "set label '$tempmax K' at -0.15,0,$tempmax\n";
		print $gnh "splot [0:1] [0:0.866] [$tempmin:$tempmax]";
		print $gnh "'" . $tiefile . "' u (\$6+0.5*\$7):(0.866*\$7):(\$4) w l $tielinecolor notitle";
		for my $ph (keys %phaseindex) {
		    print $gnh ", \\\n";
		    my $filename="phase_" . $ph . ".dat";
		    if (exists $cmdline{"-1f"}) {
			$filename="-";
		    }
		    print $gnh "'" . $filename . "' u (\$6+0.5*\$7):(0.866*\$7):(\$4) w p t '" . escape_bad($ph) . "'";
		}
	    }
	}
    }
    elsif (@elem==4 && ($tempmin == $tempmax)  && (! exists $cmdline{"-nd"}) ) {
	if (exists $cmdline{"-vtk"}) {
	    system("${mkaxes} -tetra -col=0 > axes.vtk");
	    system("${mkaxes} -tetra -col=0 -vl=" . join(',',@elem) . "> labels.vtk");
	    if ( ! exists $cmdline{"-nt"} ) {
		system("awk '{print \$4,\$5,\$6,\$7,\$8}' $tiefile | ${mkaxes} -trans -tetra | ${simplexize} -t=- -l -col=0  > tielines.vtk");
	    }
	    my $color=1;
	    my $cmd="";
	    for my $ph (keys %phaseindex) {
		my $filename="phase_" . $ph . ".dat";
		my $vtkfilename="phase_" . $ph . ".vtk";
		$cmd = $cmd . "awk '{print \$4,\$5,\$6,\$7,\$8}' $filename | ${mkaxes} -trans -tetra | ${simplexize} -t=- -rmin=$pttooclose -rmax=$pttoofar -col=$color > $vtkfilename" . ( $nbproc > 1 ?  "& " : "; ");
		$color++;
	    }
	    if ($nbproc>1) {$cmd = $cmd . "wait";}
	    system($cmd);
	}
	else {
	    print $gnh "unset border\nunset tics\n";
	    print $gnh "set size square 1,1\n";
	    print $gnh "set arrow from 0,0,0 to 1,0,0 nohead\n";
	    print $gnh "set arrow from 1,0,0 to 0.5,0.866,0 nohead\n";
	    print $gnh "set arrow from 0.5,0.866,0 to 0,0,0 nohead\n";
	    print $gnh "set arrow from 0,0,0 to 0.5,0.289,0.816 nohead\n";
	    print $gnh "set arrow from 1,0,0 to 0.5,0.289,0.816 nohead\n";
	    print $gnh "set arrow from 0.5,0.866,0 to 0.5,0.289,0.816 nohead\n";
	    print $gnh "set label '$elem[0]' at -0.05,-0.05,0\n";
	    print $gnh "set label '$elem[1]' at 1.05,-0.05,0\n";
	    print $gnh "set label '$elem[2]' at 0.5,0.866+0.05,0\n";
	    print $gnh "set label '$elem[3]' at 0.5,0.289,0.816+0.05\n";
	    print $gnh "splot [0:1] [0:0.866] [0:0.816] ";
	    print $gnh "'" . $tiefile . "' u (\$6+0.5*\$7+0.5*\$8):(0.866*\$7+0.289*\$8):(0.816*\$8) w l $tielinecolor notitle";
	    for my $ph (keys %phaseindex) {
		print $gnh ", \\\n";
		my $filename="phase_" . $ph . ".dat";
		if (exists $cmdline{"-1f"}) {
		    $filename="-";
		}
		print $gnh "'" . $filename . "' u (\$6+0.5*\$7+0.5*\$8):(0.866*\$7+0.289*\$8):(0.816*\$8) w p t '" . escape_bad($ph) . "'";
	    }
	}
    }
    else {
	if ( ! exists $cmdline{"-nd"} ) {
	    if (@elem == 4) {print "Try setting T0=T1.\n";}
	    die "Sorry, case not yet covered for plotting. The files phase_xxxx.dat may still be useful... or try experimental -nd option";
	}
	else {
	    my $usempi=(exists $cmdline{"-mpi"});
	    my $dim;
	    my $axistype;
	    my $allscaling;
	    if ($tempmin == $tempmax) {
		$dim=@elem-1;
		$axistype="-tetra";
		$allscaling="";
	    }
	    else {
		$dim=@elem;
		$axistype="-tri";
		my $zscale=($tempmax-$tempmin);
		$allscaling="-zmin=$tempmin -zmax=$tempmax -zscale=$zscale";
	    }
	    system("${mkaxes} ${axistype} -ef=bi -rgb=255,255,255,50,1 -hidim=${dim} $allscaling > axes.nd");
	    system("${mkaxes} ${axistype} -mat -hidim=${dim} $allscaling -vl=" . join(',',@elem) . " > ndaxes.in");
	    if ( ! exists $cmdline{"-nt"} ) {
		system("awk '{nc=split(\$0,a,\" \"); for (i=4; i<=nc; i++) {printf a[i] \" \";} print \"\";}' $tiefile | ${mkaxes} -trans ${axistype} ${allscaling} -hidim=${dim} | ${simplexize} -t=- -tie -rgb=33,33,33 -ef=bi > tielines.nd");
	    }

	    my $icolor=0;
	    my $ncolor=scalar keys %phaseindex;
	    my $cmd="";
	    for my $ph (keys %phaseindex) {
		$theta=2*$pi*$icolor/$ncolor;
		$rgb=int(100+100*cos($theta)) . "," . int(100+100*cos($theta+2*$pi/3)) . "," . int(100+100*cos($theta+4*$pi/3)) . ",255";
		my $filename="phase_" . $ph . ".dat";
		my $ndfilename="phase_" . $ph . ".nd";
		if ($usempi) {
		    my $tmpfilename="phase_" . $ph . ".tmp";
		    $cmd = $cmd . "awk '{nc=split(\$0,a,\" \"); for (i=4; i<=nc; i++) {printf a[i] \" \";} print \"\";}' $filename | " . 
			"${mkaxes} -trans ${axistype} ${allscaling} -hidim=${dim} > ${tmpfilename} ; " . 
			"${mpirun} mpi${simplexize} -t=${tmpfilename} -rmin=$pttooclose -rmax=$pttoofar -rgb=$rgb -ef=bi > $ndfilename ; ";
		    if (! exists $cmdline{"-debug"} ) {
			$cmd = $cmd . "rm ${tmpfilename} ; ";
		    }
		}
		else {
		    $cmd = $cmd . "awk '{nc=split(\$0,a,\" \"); for (i=4; i<=nc; i++) {printf a[i] \" \";} print \"\";}' $filename | " . 
			"${mkaxes} -trans ${axistype} ${allscaling} -hidim=${dim} | " . 
			"${simplexize} -t=- -rmin=$pttooclose -rmax=$pttoofar -rgb=$rgb -ef=bi > $ndfilename " . ( $nbproc > 1 ?  "& " : "; ");
		}
		$icolor++;
	    }
	    if ($nbproc>1 && ! $usempi) {$cmd = $cmd . "wait";}
	    print $cmd . "\n";
	    system($cmd);
	    #if ( ! exists $cmdline{"-nt"} ) {
	    #    system("awk '{print \$4,\$5,\$6,\$7,\$8}' $tiefile | ${mkaxes} -trans -tetra | ${simplexize} -t=- -l -col=0  > tielines.vtk");
	    #}
	}
    }
    if ( (! exists $cmdline{"-vtk"}) && (! exists $cmdline{"-nd"}) ) {
	print $gnh "\n";
	if (exists $cmdline{"-1f"}) {
	    copy_file_to_handle($gnh,"tielines.dat");
	    unlink("tielines.dat");
	    print $gnh "e\n";
	    for my $ph (keys %phaseindex) {
		my $filename="phase_" . $ph . ".dat";
		copy_file_to_handle($gnh,$filename);
		unlink($filename);
		print $gnh "e\n";
	    }
	}
	print $gnh "pause -1\n";
	close($gnh);
	print "You can now run:\n gnuplot ocplotpd.gnu\n";
    }
}

