fargnuplot2D1D

#!/usr/bin/perl -s
# =====================================================
# FARGNUPLOT 
# =====================================================
# Plotting Tool for FARGO (Frederic Masset). Fargo is a
# simple and fast 2D polar hydrodynamic code dedicated
# to planet-disk tidal interactions. This perl script
# generates 2D plots, to quickly check the results of a
# run.
#
# Run './fargoplot -h' to see the syntax and examples.
#
# Note that this script needs a gnuplot version >= 4.2
#
# -----------------------------------------------------
# Hanno Rein, Cambridge University, 2008-2009
# http://damtp.cam.ac.uk/user/hr260/
# Thanks to Aurelien Crida and Tobias Mueller.
# -----------------------------------------------------

### CONSTANTS & INCLUDES########################################################

use constant GNUPLOTCMD => "gnuplot";
use constant PI2 => "6.28318530717959";

# needed for random hash
use Digest::MD5;
use File::Basename;

### SUB ROUTINES ###############################################################

# prints command line options supported
sub printusage
{
	print <<EUS;
FarGnuPlot 0.53 of 21 March 2012. Maintained by Hanno Rein.
Send comments to hr260\@damtp.cam.ac.uk

Usage: fargnuplot [options] [in/fargo.par | directory]

  -h          print this message
  -d          do not delete temporary files, usefull 
              for further processing
  -p          persistent X11 window
  -x          exit gnuplot after plotting (implies -p)
  -num=NUM    plot file with number NUM in given
              directory instead of newest file
  -file=FILE  plot FILE, need to specify nr and ns as well
  -nr         number of radial grid cells
  -ns         number of azimuthal grid cells
  -nr1D       number of radial grid cells in the 1D grid

  -temp=DIR   directory for creating temporary files
              (ideally a tmpfs filesystem)

  -g          plot gas density (default)  gasdens*.dat
  -vr         plot radial velocity        gasvrad*.dat
  -vt         plot azimuthal velocity     gasvtheta*.dat
  -t          plot gas temperature        gasTemperature*.dat
  -e          plot gas energy             gasenergy*.dat
  -q          plot toomre parameter	  gasToomre*.dat

  -l          use logarithmic scale
  -lcb        use logarithmic scale for cb axis
  -scale=FACTOR	scale by factor
  -r=RANGE    range for cb axis (e.g. 0:1)
  -xr=RANGE   range for x axis (e.g. 0:1)
  -yr=RANGE   range for y axis (e.g. 0:1)
  -palette=PALETTE set gnuplot palette
              (e.g. "( 0 \"blue\", 1 \"red\" )" )
  -c          circular plot (default: rectangular plot)
  -axis       change axis in circular plot to have positive x axis in "right" direction
  -clean      plots only the data, no borders, no axis
  -curve      plot average values as a function of radius
  -all        plot also the 1D grid (for FARGO-2D1D)
  -s          square plot
  -png        output png file and exit gnuplot
  -eps        output eps file and exit gnuplot
  -title=TITLE title
  -xlabel=XLABEL x-axis label
  -ylabel=YLABEL y-axis label
  -cblabel=CBLABEL cb-axis label

Examples:
  fargnuplot -nr=384 -ns=384 out5       => plot the lates
                     density file in the directory out5
  fargnuplot in/stockholm_MJ_visc.par   => plot the latest
                     density file in the output directoty
                     specified in the .par file
  fargnuplot -nr=128 -ns=384 -file=out5/gasdens125.dat
                     => plot the specified file
EUS
	if ($_[0]){
		print "\n*** Error: ";
		print $_[0]."\n";
	}
	exit(-1);
}

# cleanUp deletes all temporary files
sub cleanUp {
	unless ($d) {
		print ("Deleting temporary data files.\n");
		unlink ("$tempdir.fargnuplot.$random_hash.plt");
		unlink ("$tempdir.fargnuplot.$random_hash.dat");
		unlink ("$tempdir.fargnuplot1D.$random_hash.dat");
	}
}

# intHandler is called when SIGINT is received, e.g. if Ctrl+C is pressed
sub intHandler {
	print "SIGINT received! Aborting... \n";
	&cleanUp();
	# return error code != 0
	exit(1);
}

$SIG{'INT'} = 'intHandler';

### MAIN ROUTINE ###############################################################

## search for file in commandline arguments
print @ARGV;
$argfile = -1;
for($argvn=0;$argvn<=$#ARGV;$argvn++){
	print $argvn."  .  ".$ARGV[$argvn].  "\n";
	if ((-e $ARGV[$argvn])){ 
		$argfile = $argvn;
	}
}

# generate random hash
my $digest = new Digest::MD5;
$digest->add(rand());
$random_hash = $digest->hexdigest();

# check if a special directory for temp files was specified
if ($temp) {
	$tempdir = "$temp/";
} else {
	$tempdir = "./";
}

if ($vr){
	$tpprefix = "gasvrad";
	$zlabel = "v_rad"
} elsif($vt) {
	$tpprefix = "gasvtheta";
	$zlabel = "v_theta"
} elsif($t) {
	$tpprefix = "gasTemperature";
	$zlabel = "Temperature"
} elsif($e) {
	$tpprefix = "gasenergy";
	$zlabel = "Energy"
} elsif($q) {
	$tpprefix = "gasToomre";
	$zlabel = "Toomre"
} else {
	# this is the default
	$tpprefix = "gasdens";
	$zlabel = "Gas Density"
}

if ($cblabel) {
	$zlabel = $cblabel;
}

$dy = 1;
$rmin =0.;
$dx = 1;

if ($file) {
	# user specified file
	if (!(-e $file)){
		&printusage("File ".$file." does not exist.");
	}
	if (-d $file){
		&printusage("Specified file ".$file." is a directory.");
	}
	$tp = $file;
	if ($nr==0 || $ns ==0){
		&printusage("Number of gridcells not specified.");
	}
	$dx = PI2 / ($ns-1);
} else {
	# automatically find newest file
	if (!(-e $ARGV[$argfile])) {
		 &printusage("Directory or file '".$ARGV[$argfile]."' does not exist.");
	}

	if (!(-d $ARGV[$argfile])) {
	# parse fargo .par file 
		open(FARGOPAR, $ARGV[$argfile]);
		while (<FARGOPAR>) {
			chomp;                  # no newline
			s/#.*//;                # no comments
			s/^\s+//;               # no leading white
			s/\s+$//;               # no trailing white
			next unless length;     # anything left?
			my ($var, $value) = split(/[\ \t]/, $_, 2);
			if ($var =~ m/\s*nrad\s*/i) {
				$_ = $value;
				chomp;
				m/(\S+)[\ \t]*/;
				$nr=$1;
			}
			if ($var =~ m/\s*nsec\s*/i) {
				$_ = $value;
				chomp;
				m/(\S+)[\ \t]*/;
				$ns=$1;
			}
			if ($var =~ m/\s*outputdir\s*/i) {
				$_ = $value;
				chomp;
				m/(\S+)[\ \t]*/;
				$ARGV[$argfile]=$1;
				$dir=$1;
			}
			if ($var =~ m/\s*rmin\s*/i) {
				$_ = $value;
				chomp;
				m/(\S+)[\ \t]*/;
				$rmin=$1;
			}
			if ($var =~ m/\s*rmax\s*/i) {
				$_ = $value;
				chomp;
				m/(\S+)[\ \t]*/;
				$rmax=$1;
			}
		} 
		close(FARGOPAR);
		if ((!$dir)||(!$nr)||(!$ns)) {
			&printusage("File ".$ARGV[$argfile]." is not a valid FARGO .par file.");
		}
		if ($rmin!=$rmax) {
			$dy = $rmax-$rmin;
			$dy /= $nr;
		}
		$dx = PI2 /($ns-1);
		if ($all) {
		    print "Using directory ".$ARGV[$argfile]. ", nr=".$nr. ", ns=".$ns. ", nr1D=".$nr1D."\n";
		} else {
		    print "Using directory ".$ARGV[$argfile]. ", nr=".$nr. ", ns=".$ns."\n";
	        }
	}
	if ($num =~/^\d+$/) {
		## user specified file number
		if (!(-e $ARGV[$argfile]."/".$tpprefix.$num.".dat")){
			&printusage("File ".$ARGV[$argfile]."/".$tpprefix.$num.".dat does not exist.");
		}
		$tp = $ARGV[$argfile]."/".$tpprefix.$num.".dat";
		if ($all) {
			if (!(-e $ARGV[$argfile]."/".$tpprefix."1D".$num.".dat")) {
				&printusage("File ".$ARGV[$argfile]."/".$tpprefix."1D".$num.".dat does not exist.");
		   	}
			$tp1D = $ARGV[$argfile]."/".$tpprefix."1D".$num.".dat";
		}
	} else {
		if (-d $ARGV[$argfile]) {
		## find newest file in directory	
			@sortedfiles =`ls -lt $ARGV[$argfile]/$tpprefix*.dat | awk '{print \$9;}'| grep -v 1D | grep -v ascii`;
			$tp = $sortedfiles[0];
			$tp =~ s/\n//;
			$num = $tp;
			$num =~ s/.*\/\/[^0-9]*(\d*)\.dat/\1/;
			print "Newest file: $tp (num=$num)\n";
			if (!(-e $tp)) {
				&printusage("Directory does not contain FARGO files.");
			}
			if ($all) {
				@sortedfiles =`ls -lt $ARGV[$argfile]/$tpprefix*.dat | awk '{print \$9;}'| grep 1D`;
				$tp1D = $sortedfiles[0];
				$tp1D =~ s/\n//;
				print " ".$tp1D;
				if (!(-e $tp1D)) {
					&printusage("Directory does not contain FARGO 1D files.");
				}
			}
		}
	}

	# search for used_rad.dat and read it if available
	if (-e $ARGV[$argfile]."/used_rad.dat") {
		open(USED_RAD, $ARGV[$argfile]."/used_rad.dat");
		@radii = <USED_RAD>;
		close(USED_RAD);
	} else {
		print "Couldn't find any information about spacing (used_rad.dat). Using linear spacing\n";
	}

	# search for used_rad1D.dat and read it if available
	if ($all) {
		if (-e $ARGV[$argfile]."/used_rad1D.dat" ) {
			open(USED_RAD, $ARGV[$argfile]."/used_rad1D.dat");
			@radii1D = <USED_RAD>;
			close(USED_RAD);
		} else {
			print "Couldn't find any information about spacing (used_rad1D.dat). Using linear spacing\n";
		}
	}
}

print "Plotting ".$tp."\n";

# logarithmic scaling?
if ($l){
	# logarithmic scaling
	$plot3 = "(log10(\$3))";
	$zlabel = $zlabel ." (log10)";
} else {
	# linear scaling
	$plot3 = "3";
}

# circular plot?
if ($c) {
	# circular plot
	if ($axis) {
		$plot1 = "(\$2*sin(\$1))";
		$plot2 = "(\$2*cos(\$1))";
	} else {
		$plot1 = "(\$2*cos(\$1))";
		$plot2 = "(\$2*sin(\$1))";
	}
	if (!($xlabel)) {
		$xlabel = "radius";
	}
} else {
	# rectangular plot
	$plot1 = "1";
	$plot2 = "2";
	if (!($xlabel)) {
		$xlabel = "theta";
	}
	$xtics='("0" 0, "{/Symbol p}/2" pi/2.0, "{/Symbol p}" pi, "3{/Symbol p}/2" 3.0/2.0*pi, "2{/Symbol p}" 2*pi)';
}

# ylabel
if (!($ylabel)) {
	$ylabel = "radius"
}

$ns++;
open(GNUPLOTW, ">$tempdir.fargnuplot.$random_hash.plt");

# enable EPS/PNG output
if ($eps) {
	if ($file) {
		my ($filename, $directories, $suffix) = fileparse($file, qr/\.[^.]*/);
		$outputfile = $directories.$filename;
	} else {
		$outputfile = $ARGV[$argfile]."/".$tpprefix.$num;
	}
	print("EPS mode. Writing to ".$outputfile.".eps\n");
	# set terminal
	print GNUPLOTW "\tset terminal postscript color enhanced\n";
	# set output filename
	print GNUPLOTW "\tset output \"$outputfile.eps\"\n";
	$x=1;
} elsif ($png) {
	if ($file) {
		my ($filename, $directories, $suffix) = fileparse($file, qr/\.[^.]*/);
		$outputfile = $directories.$filename;
	} else {
		$outputfile = $ARGV[$argfile]."/".$tpprefix.$num;
	}
	print("PNG mode. Writing to ".$outputfile.".png\n");
	# set terminal
	print GNUPLOTW "\tset terminal pngcairo enhanced truecolor notransparent size 1280,1024\n";
	# set output filename
	print GNUPLOTW "\tset output \"$outputfile.png\"\n";
	$x=1;
}

# square plot
if ($s) {
	print GNUPLOTW "\tset size square\n";
}

print GNUPLOTW "\tset tics out\n";
if ($title) {
	print GNUPLOTW "\tset title \"$title\"\n";
} else {
	print GNUPLOTW "\tset title \"FARGO file: $tp\"\n";
}

if (!$scale) {
	$scale = 1.0;
}

if ($curve) {
	# average values plot
	$plot3 =~ s/3/2/g; 
	print GNUPLOTW "\tset xlabel \"radius\"\n";
	print GNUPLOTW "\tset autoscale xfix\n";
	print GNUPLOTW "\tset ylabel \"$zlabel\"\n";
	print GNUPLOTW "\tset style data lines\n";
	print GNUPLOTW "\tplot \"$tempdir.fargnuplot.$random_hash.dat\" binary record=${nr} format=\"%double\" u 1:$plot3 notitle \n";
} else {
	# normal plot
	print GNUPLOTW "\tset pm3d map\n";
	print GNUPLOTW "\tset autoscale fix\n";
	print GNUPLOTW "\tset xlabel \"$xlabel\"\n";
	if ($xtics) {
	print GNUPLOTW "\tset xtics $xtics\n";
	}
	print GNUPLOTW "\tset ylabel \"$ylabel\"\n";
	print GNUPLOTW "\tset cblabel \"$zlabel\"\n";

	# logarithmic scale for cb axis
	if ($lcb) {
		print GNUPLOTW "\tset logscale cb\n"
	}
	
	# range for cb axis
	if ($r) {
		print GNUPLOTW "\tset cbrange [$r]\n"
	}

	# palette for cb axis
	if ($palette) {
		print GNUPLOTW "\tset palette defined $palette\n";
	}

	# range for x axis
	if ($xr) {
		print GNUPLOTW "\tset xrange [$xr]\n"
	}

	# range for y axis
	if (yr) {
		print GNUPLOTW "\tset yrange [$yr]\n"
	}

	if ($clean) {
		print GNUPLOTW "\tunset cbtics\n";
		print GNUPLOTW "\tunset ytics\n";
		print GNUPLOTW "\tunset xtics\n";
		print GNUPLOTW "\tunset tics\n";
		print GNUPLOTW "\tunset title\n";
		print GNUPLOTW "\tunset colorbox\n";
		print GNUPLOTW "\tunset xlabel\n";
		print GNUPLOTW "\tunset ylabel\n";
		print GNUPLOTW "\tset size square\n";
		print GNUPLOTW "\tset border 0\n";
		print GNUPLOTW "\tset lmargin 0\n";
		print GNUPLOTW "\tset rmargin 0\n";
		print GNUPLOTW "\tset tmargin 0\n";
		print GNUPLOTW "\tset bmargin 0\n";
		print GNUPLOTW "\tset size ratio 1\n";
	}
	
	# finally plot
	if ($all) {
		# with 1D
		print GNUPLOTW "\tsplot \"$tempdir.fargnuplot1D.$random_hash.dat\" binary record=${ns}x${nr1D} format=\"%double\" u $plot1:$plot2:$plot3 notitle, \".fargnuplot.$random_hash.dat\" binary record=${ns}x${nr} format=\"%double\" u $plot1:$plot2:$plot3 notitle\n";
	} else {
		# without 1D
		print "\tsplot \"$tempdir.fargnuplot.$random_hash.dat\" binary record=${ns}x${nr} format=\"%double\" u $plot1:$plot2:(\$$plot3*$scale) notitle\n";

		print GNUPLOTW "\tsplot \"$tempdir.fargnuplot.$random_hash.dat\" binary record=${ns}x${nr} format=\"%double\" u $plot1:$plot2:(\$$plot3*$scale) notitle\n";
	}
}
$ns--;
close(GNUPLOTW);

# read data file and write temporary .dat file
open(GNUPLOTW, ">$tempdir.fargnuplot.$random_hash.dat"); binmode(GNUPLOTW);
open(INP,$tp); binmode(INP);
$i=0;
$j=$rmin;
$counti=0;
$countj=0;
$sum=0;
while (read(INP, $buff, 8 )) {
	# closing the gap by printing first row twice
	if ($i==0) {
		$printagain=$buff;
	}
	if (@radii) {
		$rad=$radii[$countj];
	} else {
		$rad=$j;
	}
	if (!$curve) {
		$coord=pack("dd",$i,$rad);
		print GNUPLOTW $coord;
		print GNUPLOTW $buff;
	} else {
		$sum += unpack("d",$buff)/$ns;
	}
	$counti++;
	$i=$dx*$counti;
	if ($i>=$ns*$dx){
		if (!$curve){
			# convert to binary format
			$coord=pack("dd",$i,$rad);
			print GNUPLOTW $coord;
			print GNUPLOTW $printagain;
		}else{
			print GNUPLOTW pack("dd", $rad, $sum);
			$sum=0;
		}
		$j+=$dy;
		$countj++;
		$i=0;
		$counti=0;
	}
}
close GNUPLOTW;

# read 1D data file and write 1D temporary .dat file
if($all) {
	open (GNUPLOTW, ">$tempdir.fargnuplot1D.$random_hash.dat"); binmode(GNUPLOTW);
	open(INP,$tp1D); binmode(INP);
	$i=0;
	$j=$rmin;
	$countj=0;
	$counti=0;
	while (read(INP, $buff, 8 ) && ($countj<=$nr1D)) {
		if (@radii1D) {
			$rad=$radii1D[$countj];
		} else {
			$rad=$j;
		}
		if (!$curve ){
			while ($counti<($ns+1)) {
				$i=$dx*$counti;
				$coord=pack("dd",$i,$rad);
				print GNUPLOTW $coord;
				print GNUPLOTW $buff;
				$counti++;
			}
		} else {
			print GNUPLOTW pack("dd", $rad, $buff);
		}
		$j+=$dy;
		$countj++;
		$i=0;
		$counti=0;
	}
	close GNUPLOTW;
}

# run gnuplot with persistent window?
if ($p||$x){
	$gnuplotoptions = " -persist ";
}else{
	$gnuplotoptions = " ";
}

# close gnuplot after plotting?
if ($x) {
	$gnuplotoptions .= " $tempdir.fargnuplot.$random_hash.plt";
}else{
	$gnuplotoptions .= " $tempdir.fargnuplot.$random_hash.plt -";
}

# Print info how to quit gnuplot only if gnuplot is not quit automatically
if (!($x)) {
	print "All data files written. Starting gnuplot. Enter 'q' to quit.\n";
}

# Run gunplot and clean up if some errors occur
if (system(GNUPLOTCMD.$gnuplotoptions)!=0) {
	print "Failed running $1\n";
	&cleanUp();
	# return error code != 0
	exit(1);
}

# Finally clean up
&cleanUp();