###############################################################
#
#   Converting a frame data file to gnuplot input file for 
#   colorflow plot of a scalar field
#
#   Janine Glänzel, Roman Unger ... 07/2014 
#
###############################################################
#
#   Copyright (C) 2014  Janine Glänzel, Roman Unger
#
#   janine.glaenzel@mathematik.tu-chemnitz.de
#   roman.unger@mathematik.tu-chemnitz.de
#
#   This program is free software: you can redistribute it 
#   and/or modify it under the terms of the GNU General Public 
#   License as published by the Free Software Foundation, 
#   either version 3 of the License, or (at your option) any 
#   later version.
#
#   This program is distributed in the hope that it will be 
#   useful, but WITHOUT ANY WARRANTY; without even the implied 
#   warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
#   PURPOSE. 
#   See the GNU General Public License for more details.
#
#   You should have received a copy of the GNU General Public 
#   License along with this program.  
#   If not, see <http://www.gnu.org/licenses/>.
#
###############################################################

use strict;


# input data file
my $f ="frame1.txt";

# gnuplot command file to create
my $f2="../example_v3.txt";

# tex input name 
my $f3="example_v3-input.tex";


# read all quads 
open(EIN,"<$f") || die "Cant read $f \n";

my @quads;
while (<EIN>)
{
  my $l=$_;
  chomp $l;
  if ( $l =~ /^\s*#/) {next;}
  if ( $l =~ /^\s*$/) {next;}
  push (@quads,$l);

}
close EIN;

my $n=@quads;

print "Found $n quads \n";

# get min max
my $xmin= 1.0e200;
my $xmax=-1.0e200;
my $ymin= 1.0e200;
my $ymax=-1.0e200;
my $zmin= 1.0e200;
my $zmax=-1.0e200;

for (my $i=0;$i<$n;$i++)
{

 (my $x1, my $y1, my $z1, my $x2, my $y2, my $z2, my $x3, my $y3, my $z3, my $x4, my $y4, my $z4 )=split(" ",$quads[$i]);

  if ( $x1 < $xmin ) { $xmin=$x1;}
  if ( $y1 < $ymin ) { $ymin=$y1;}
  if ( $z1 < $zmin ) { $zmin=$z1;}

  if ( $x1 > $xmax ) { $xmax=$x1;}
  if ( $y1 > $ymax ) { $ymax=$y1;}
  if ( $z1 > $zmax ) { $zmax=$z1;}

  if ( $x2 < $xmin ) { $xmin=$x2;}
  if ( $y2 < $ymin ) { $ymin=$y2;}
  if ( $z2 < $zmin ) { $zmin=$z2;}

  if ( $x2 > $xmax ) { $xmax=$x2;}
  if ( $y2 > $ymax ) { $ymax=$y2;}
  if ( $z2 > $zmax ) { $zmax=$z2;}

  if ( $x3 < $xmin ) { $xmin=$x1;}
  if ( $y3 < $ymin ) { $ymin=$y1;}
  if ( $z3 < $zmin ) { $zmin=$z1;}

  if ( $x3 > $xmax ) { $xmax=$x3;}
  if ( $y3 > $ymax ) { $ymax=$y3;}
  if ( $z3 > $zmax ) { $zmax=$z3;}

  if ( $x4 < $xmin ) { $xmin=$x4;}
  if ( $y4 < $ymin ) { $ymin=$y4;}
  if ( $z4 < $zmin ) { $zmin=$z4;}

  if ( $x4 > $xmax ) { $xmax=$x4;}
  if ( $y4 > $ymax ) { $ymax=$y4;}
  if ( $z4 > $zmax ) { $zmax=$z4;}


}



print "zmin= $zmin  zmax=$zmax \n";


open (AUS,">$f2") || die "Cant write $f2\n";

print AUS "########################################################\n";
print AUS "#\n";
print AUS "#  This file is generated by $0 \n";
print AUS "#\n";
print AUS "#  DO NOT EDIT !!! \n";
print AUS "#\n";
print AUS "#  EDIT $0 and run it again !!! \n";
print AUS "#\n";
print AUS "########################################################\n";

print AUS "reset\n";

print AUS "set terminal epslatex size 16cm,10cm color colortext dashed \n";
print AUS "set output \"$f3\" \n";
print AUS "set size  0.9,0.55 \n";


print AUS "set multiplot\n";


# plot the colorbar
print AUS "set size 0.23,0.5 \n";
print AUS "set origin 0.78,0.0\n";
print AUS "unset xtics\n";
print AUS "unset ytics\n";
print AUS "set y2tics\n";
print AUS "set y2label '\$|\$u\$|\$[mm]' rotate by 0 offset -7.5,6.5 \n";


my $ncols=150;
print AUS "plot [:] [$zmin:$zmax] \\\n";
for (my $i=0;$i<$ncols;$i++)
{
 
(my $red, my $green, my $blue)=colmap99($i / ($ncols-1));


 my $col=sprintf("#%02x%02x%02x",$red*255,$green*255,$blue*255);
 print AUS " \"-\" title \"\" with filledcurve lc rgb \"$col\" fillstyle transparent solid 1.000000  ";
 if ($i<$ncols-1) 
   {print AUS "\\\n,";}
 else 
   {print AUS ";\n";}
}


for (my $i=0;$i<$ncols;$i++)
{
 my $y1=$zmin+$i*($zmax-$zmin)/$ncols;
 my $y2=$y1+($zmax-$zmin)/$ncols;

 my $x2=0.1*($zmax-$zmin);
 print AUS "0.0 $y1\n";
 print AUS "$x2 $y1\n";
 print AUS "$x2 $y2\n";
 print AUS "0.0 $y2\n";
 print AUS "e\n";
}



# and plot the mesh
print AUS "set size 0.83,0.54 \n";
print AUS "set origin 0.0,0.0\n";

print AUS "set xtics\n";
print AUS "set ytics\n";
print AUS "unset y2tics\n";
print AUS "unset y2label\n";
print AUS "set xlabel 'x\$_{\\mbox{\\scriptsize 1}}\$[mm]' offset 0,0.75\n";
print AUS "set ylabel 'x\$_{\\mbox{\\scriptsize 2}}\$[mm]' offset 1.6,0\n";

my $xrange=sprintf("[%f:%f]",$xmin-0.03*($xmax-$xmin),$xmax+0.03*($xmax-$xmin));
my $yrange=sprintf("[%f:%f]",$ymin-0.1*($ymax-$ymin),$ymax+0.1*($ymax-$ymin));

print AUS "plot $xrange $yrange \\\n";
for (my $i=0;$i<$n;$i++)
{

 (my $x1, my $y1, my $z1, my $x2, my $y2, my $z2, my $x3, my $y3, my $z3,my $x4, my $y4, my $z4 )=split(" ",$quads[$i]);

 my $u = ($z1+$z2+$z3+$z4) / 4;

 # map u to [0,1]
 $u = ($u-$zmin) / ($zmax-$zmin);

 # now u can mapped to a color,
 (my $red, my $green, my $blue)=colmap99($u);
 

 # and make the color string
 my $col=sprintf("#%02x%02x%02x",$red*255,$green*255,$blue*255);

 print AUS " \"-\" title \"\" with filledcurve lc rgb \"$col\" fillstyle transparent solid 1.000000  ";
 if ($i<$n-1) 
   {print AUS "\\\n,";}
 else 
   {print AUS ";\n";}
}


for (my $i=0;$i<$n;$i++)
{

 (my $x1, my $y1, my $z1, my $x2, my $y2, my $z2, my $x3, my $y3, my $z3,my $x4, my $y4, my $z4)=split(" ",$quads[$i]);
 print AUS "$x1 $y1\n";
 print AUS "$x2 $y2\n";
 print AUS "$x3 $y3\n";
 print AUS "$x4 $y4\n";
 print AUS "e\n";

}


print AUS "unset multiplot\n";

close AUS;



# colormap like  "99 colors" in workadag
sub colmap99
{
  my $u=shift;
  my $red;
  my $green;
  my $blue;
  
  if($u<0){ $red=1;$green=0;$blue=1;return ($red,$green,$blue);}
  if($u>1){ $red=1;$green=1;$blue=1;return ($red,$green,$blue);}


  if  ( $u < 0.2 ) 
      {
      $red    = 1.0-5*$u;
      $green  = 0;
      $blue   = 1; 
      return ($red,$green,$blue);
      }

    if  ( $u < 0.4 ) 
      {
      $red   = 0;
      $green = 5*($u-0.2);
      $blue  = -5.0 * ($u-0.2) + 1.0;
      return ($red,$green,$blue);
      }
    
    if  ( $u < 0.6  ) 
      {
      $red = 5.0 * ($u-0.4);
      $green = 1;
      $blue  = 0;
      return ($red,$green,$blue);
      }

    if  ( $u < 0.8  ) 
      {
      $red = 1;
      $green = -5.0 * ($u-0.6) + 1.0;
      $blue  = 0;
      return ($red,$green,$blue);
      }

    # case u \in [ 0.8 , 1 ]
    $red   = 1;
    $green = 5.0 * ($u-0.8);
    $blue  = 5.0 * ($u-0.8);
    return ($red,$green,$blue);

}