#
#  commacont.gnu
#
#  - a Gnuplot script for making postscript contour plot of 2-dim. function 
#  tabulated on rectangular grid (e.g., charge density as exported from WIEN2k).
#  Written by Andrei Postnikov (andrei.postnikov@univ-lorraine.fr)
#  Universit\'e de Lorraine, Metz, February 2015
#
#  Check some settings : 
#
# 1) Input (data) and output (postscript) files:
    datafile = "rho1"; set output "rhoplot.ps" 
#   For the datafile format, see the Gnuplot documentation. Each line contains either a (X,Y,Z) triple,
#   or just a Z value, in datablock structure (blocks separated by blank lines): the datablock number 
#   will be used for Y, and the index of the data point in the datablock will be used for X.
#
# 2) Contour from .. to, linear or logarithmic, emphasize every NNth contour : 
    cont_mode = 2 #  1: for linear contours spacing, 2: for logarithmic spacing
    z_min = 0.01; z_max = 5.12; n_z = 10    # - define n_z contours from z_min to z_max 
    m_z = 5       #   every m_z'th contour will be drawn differently (set > n_z if not needed)
#
# 3) Monochrome or color ? 
    color_mode = 2 # 1: monochrome, 2: color
#
# 4) Set line properties: linewidth (=lw) and color (=lc) ;
#    "show colornames" for the list of 112 predefined color names
    lw_1 = 2; lc_1 = "forest-green"  # - for "regular" contours
    lw_2 = 4; lc_2 = "coral"         # - for "selected" contours and labels
    lw_3 = 2; lc_3 = "navy"          # - for the plot frame
    lw_t = 2                         # - line width global multiplicator in terminal settinig
#
# 5) Canvas size / plot dimensions (in units of canvas size):
    x_canvas = 7.86885; y_canvas = 7.86885  # - canvas size in inch (can be left as is)
    x_size = 0.40; y_size =0.90             # - may need fine-tuning to correct proportions
#
# 6) Define the places as which labels will be placed at "selected" contours.
#    These places are intersections of the contours with isolines flab(x,y)=0,
#    where the both arguments of flab() vary 0:1 throughout the plot window.
#    Leave one definition (or none); examples are:
#    flab(x,y) = y-0.5  #  labels at half-height of the plot
#    flab(x,y) = (3*x-1)*(3*x-2)  #  labels on two vertical lines, x=1/3 and x=2/3
#    flab(x,y) = y-x*x      #  on a parabola from (0,0) to (1,1)
#    flab(x,y) = 4*y-2*x-1  #  a diagonal cut from (0,1/4) to (1,3/4)
     flab(x,y) = 4*x-2*y-1  #  a diagonal cut from (1/4,0) to (3/4,1)
#    flab(x,y) = (8*y-2*x-1)*(8*y-2*x-5)  #  two diagonal cuts: (0,1/8)-(1,3/8) and (0,5/8)-(1,7/8)

# 7) Attributes of labels to place at contours.
     lab_font = "Helvetica";   lab_size = 18    # - font and size of labels
#    lab_font = "Times-Roman"; lab_size = 20    # - font and size of labels
#    lab_font = "Courier";     lab_size = 18    # - font and size of labels
     lab_form = "%4.2f"    # - label format
     aa = 0.12; bb = 0.06  # - dimension of (white) ellipses drawn beneath the labels,
#    to mask part of contours, (both) expressed in units of the X-axis span.
#    May need adjustment depending on font size.

# 8) Debugging
     debug = 0  #  =0 just does the job, otherwise verbose output and drawing, retaining temp files
#
# ----------- From here on, hopefully nothing needs to be touched: ---------
#
  if ( GPVAL_VERSION < 4.6 ) { print " This script needs gnuplot 4.6 installed "; exit }
# Explore initial data, extract min/max x/y values;
# for this we need to "splot" the original data once; do it into the 'dumb' terminal:
  set term dumb size 1,1; unset xtics; unset ytics
  splot datafile with points 
  if (debug!=0) {
  print " X min,max : ",GPVAL_X_MIN," ,  ",GPVAL_X_MAX, \
     " ;  Y min,max : ",GPVAL_Y_MIN," ,  ",GPVAL_Y_MAX  }
  x_min = GPVAL_X_MIN; x_max = GPVAL_X_MAX-1.0; y_min = GPVAL_Y_MIN; y_max = GPVAL_Y_MAX-1.0
# actual sizes (in units of X) of masking ellipses :
  a_axis=aa*(x_max-x_min); b_axis=bb*(x_max-x_min)
# set default color of masking ellipses :
  if (debug == 0) { ell_color = "white" } else { ell_color = "yellow" }
# if output monochrome force the labels to be black :
  if (color_mode == 1) { lc_2 = "black" }

  if ((color_mode != 1) && (color_mode != 2)) {
    print "color_mode = ", color_mode, " : not predefined ! - use 1 or 2 "; quit }
  if (color_mode == 1) {
    set terminal postscript eps monochrome solid lw lw_t  size x_canvas, y_canvas }
  if (color_mode == 2) {
  set terminal postscript eps color      solid lw lw_t  size x_canvas, y_canvas }
  set linetype 1  lw lw_1  lc rgb lc_1   # - for "regular" contours
  set linetype 2  lw lw_2  lc rgb lc_2   # - for "selected" contours and (color) labels
  set linetype 3  lw lw_3  lc rgb lc_3   # - for plot frame

  set noxtics; set noytics; set noztics         # - suppress all labels
  set nosurface; set hidden3d; set contour base # - show only contours (no surface) 
  set nokey; set view 0,0,1
  set cntrparam bspline; set cntrparam order 3  # - contour interpolation scheme
# set noclabel
# set lmargin 0; set bmargin 0; set tmargin 0; set rmargin 0

  if ((cont_mode != 1) && (cont_mode != 2)) {
    print "cont mode = ", cont_mode, " : not predefined ! - use 1 or 2 "; quit }
  if (cont_mode == 2) { 
    dz_log = (log(z_max)-log(z_min))/(log(2.0)*(n_z-1)) }  # - logarithmic scale (each interval doubles) 
  if (cont_mode == 1) { 
    dz_lin = (z_max - z_min)/(n_z-1) }                     # - linear scale (all intervals equal) 

  unset border  # - remove the frame, to begin with...
# - (if frame needed): draw frame (15 means closed box) by previously defined linetype
  set border 15 lt 3     

  unset clabel   # to draw all contours with the same style

# Define functions to save three variables
# (the values of X, Y and of the block number from the previous data point) :
   prev_x(x) = (x_p2=x_p1, x_p1=x, x_p2)  # the argument is shifted into x_p1 ant retrieved on next call
   prev_y(x) = (y_p2=y_p1, y_p1=x, y_p2)
   prev_b(x) = (b_p2=b_p1, b_p1=x, b_p2)
# Define function to set label (on contours) :
  set_label(x,y,cont_lab,angle) = sprintf( \
 "set label \"%s\"at first %f,%f center rotate by %f font \"%s,%f\" front textcolor rgb lc_2",\
  cont_lab,x,y,angle,lab_font,lab_size)

set multiplot  #  Now, the plotting part. Define dimensions :

set size x_size, y_size; set origin 0.05,0.05
set lmargin 0.1; set bmargin 0.05; # set tmargin 0; set rmargin 0
set xrange [x_min:x_max]; set yrange [y_min:y_max]

#  Loop over "regular" contours (each added by its own 'splot' ):
do for [t=0:n_z-1] {
  if (cont_mode == 2) { z_cont = z_min*exp(log(2.0)*t*dz_log) }   # - for logarithmic scale ! 
  if (cont_mode == 1) { z_cont = z_min+t*dz_lin }                 # - for linear scale !
  mod_z = t+1 - ((t+1)/m_z)*m_z    # - contour number modulo m_z ( =0 for every m_z'the contour)
  set cntrparam levels discrete z_cont
  if (mod_z != 0) { 
   set table 'contour.dat'; splot datafile with points; unset table
   plot 'contour.dat' using 1:2 with lines lt 1
  }  
}  #  ... do for [t=0:n_z-1]

# Debug: plot flab(x,y)=0
   if ( debug != 0 ) {
   set table 'flab.dat' 
   set cntrparam levels discrete 0.0
   splot flab(x/(x_max-x_min),y/(y_max-y_min)) with points
   unset table
   plot 'flab.dat' using 1:2 with lines lt 3
}

#  Loop over 'selected' contours (each added by its own 'splot' ):
do for [t=0:n_z-1] {
  if (cont_mode == 2) { z_cont = z_min*exp(log(2.0)*t*dz_log) }   # - for logarithmic scale ! 
  if (cont_mode == 1) { z_cont = z_min+t*dz_lin }                 # - for linear scale !
  mod_z = t+1 - ((t+1)/m_z)*m_z          # - contour number modulo m_z ( =0 for every m_z'the contour)
  set cntrparam levels discrete z_cont   # - set contour levels (only one here)
  cont_lab = sprintf(lab_form,z_cont)    # - transform contour level into formatted string
  if (debug != 0) { show contour }  # - prints out the contour values into the standard output
  if (mod_z == 0) {
    set table 'contour.dat'; splot datafile with lines lt 2; unset table
    plot 'contour.dat' using 1:2 with lines lt 2  # - draw "selected" contours

# The rest is about setting labels. We go along the contour line
# and check where the function flab(x,y) changes sign between two points,
# then we identify the intersection point and the slope of the contour.
# For this, we need access to actual and previous (X, Y, flab(X,Y)) triplets.
# The difficulty is, contours may be split into multiple fragments (blocks). When we switch
# to a new block we want to start the history of triplets anew. We use that
# the index of current data block is stored in $(-1)$ 
# (and the currect data line in $0, but this is not used). 
# Open 'block.dat' with the following data structure :
#  (1 if opens new block):(old_x):(old_y):(new_x):(new_y):(new_flab)
    x_p1 = 0.; y_p1 = 0.; b_p1 = -1
    set table 'block.dat'
    plot 'contour.dat' using (column(-1)-prev_b(column(-1))):(prev_x($1)):(prev_y($2)): \
      ($1):($2):(flab($1/(x_max-x_min),$2/(y_max-y_min))) with xyerrorbars 
    unset table
# Regroup the same data so as to have access, in each record, to old and new values :
#  (old_x):(old_y):(old_flab):(new_x):(new_y):(new_flab)
    set table 'trace.dat'
    plot 'block.dat' using ($1==1?$4:$2):($1==1?$5:$3): \
      ($1==1?$6:flab($2/(x_max-x_min),$3/(y_max-y_min))): \
      ($4):($5):($6) with xyerrorbars
    unset table
    if (debug==0) { system "rm ./block.dat" }

#  Note that for the first line in each block, three 'new' entries are identical
#  to three 'old' ones. If (old_flab) and (new_flab) have different sign,
#  we search for the interpolated (x,y) point where to place the label,
#  and for the slope (old)-(new) of the contour fragment, to rotate the label  :
    set table 'lbls0.dat'
    plot 'trace.dat' using (sgn($3*$6)<=0?(($6*$1-$3*$4)/($6-$3)):1/0): \
      (($6*$2-$3*$5)/($6-$3)):(atan(($5-$2)/($4-$1))*180/pi): \
      (($5-$2)/($4-$1)) with errorbars
    unset table
    if (debug==0) { system "rm ./trace.dat" }
#  Data structure in 'lbls0.dat' :
#  (x_center):(y_center):(rotation angle):(slope)

# This would be already OK for plotting; unfortunately 'plot with labels'
# does not support rotating the labels at arbitrary angle. Instead, we resort to
# 'set label' mechanism. For this, we prepare a compact list of labels
# (removing all junk from the 'lbls0.dat' file), format each line
# into a meaningful command and feed in the list with 'load'. The formatting
# is done via a sequence of 'sed' actions, split into several steps for clarity: 
#   1. lbls0.dat -> lbls1.dat :
# remove lines beginning with '#', then empty lines, and finally the "NaN" ones
# (those marked 'u' in the last column; as the rest in the line are numerics,
# we simply remove all the lines containing 'u'). 
    system "sed '/^#/d' ./lbls0.dat | sed '/^$/d' | sed '/u/d' > lbls1.dat"
#   2. lbls1.dat -> lbls2.dat :
# squeeze multiple spaces into a single one, insert ")" at the 4th space, delete till end of line:
    system "tr -s ' ' <  ./lbls1.dat | sed 's/ / )/4' | sed 's/).*/)/' > lbls2.dat"  
#   3. lbls2.dat -> lbls3.dat :
# add words and commas as needed to shape the command:
    system "sed 's/ /, cont_lab, /3' ./lbls2.dat | sed 's/ /, /2' | sed 's/^/ eval set_label(/1' > lbls3.dat"

    load "lbls3.dat"  # - load 'set label' commands from previously prepared file

# The ellipses (to mask part of the contours under the label) are set by 'plot with ellipses' :
    plot 'lbls1.dat' using ($1):($2):(a_axis):(b_axis):($3) with ellipses \
                    lc rgb ell_color fs solid
    if (debug==0) {system "rm lbls[0-3].dat"}
  }  # ... if (mod_z == 0) 
}  # ... do for [t=0:n_z-1]
unset multiplot