#!/usr/bin/perl
#version 1.30, July 7, 2000.  Author: Brian Fiedler bfiedler@ou.edu
use Getopt::Long;
use Cwd;
GetOptions("b","px","pxx","jflip","k:s","dk:s","u:i","v:i","s:i","o:s","p:s","t:s","dots:s");
$path_to_here=$0;
$path_to_here=~s/\/conplot$//;
$path_to_desc=$path_to_here;
$path_to_desc=~s/conplot$/conplot_desc/;
$cwd=cwd();
$datafile=$ARGV[0];
$descfile=$path_to_here.'/'.'default.desc';
$descfile='default.desc' if -e 'default.desc';
$descfile=$path_to_desc.'/'.$opt_dk if defined($opt_dk);
$descfile=$opt_k if defined($opt_k);
$outeps=$ARGV[1];
$outeps='temp.eps' unless defined($outeps);
$outeps=~m/\.eps$/ || die ("$outeps not .eps");
$pixplot=1 if $opt_px or $opt_pxx;

while ($opt_o=~s/(\$.*?=.*?;)//){
	print $1,"\n";
	eval($1);
}


if ($descfile){
	open (DESC,"<$descfile") || die ("cannot open data file $descfile");
	foreach $line (<DESC>) {
		next if $line=~m/^\s*$/ or $line=~m/^\#.*/;
		if ($line=~m/^%/) { push @descom,$line}
		elsif ($line=~m/^\s*\@/){eval($line)} 
		elsif ($line=~m/^\s*\$/){eval($line)} 
	}
	close DESC;
}

open (DATA,"<$datafile") || die ("cannot open data file $datafile");
$inn=0;
$jnn=0;
$ispc=0;
$font='Geneva' unless defined $font;
$x_col=1 unless defined $x_col;
$y_col=2 unless defined $y_col;
$s_col=3 unless defined $s_col;
$s_col=$opt_s if defined $opt_s;
$u_col=$opt_u if defined $opt_u;
$v_col=$opt_v if defined $opt_v;
$splitter='\|' unless defined $splitter;
print "using \$s_col= $s_col  \$u_col=$u_col \$v_col=$v_col\n";

foreach $line (<DATA>) {
	$ispace=$ispc;
	$ispc=($line=~m/^\s*$/);
	next if $line=~m/^\s*$/;
	if ($line=~m/^\#\-/) {
		$tline=$line
	}
	elsif ($line=~m/^\#.*/) {
		$temp=$line;
		while ($temp=~s/(\$.*?=.*?;)//){
			print $1,"\n";
			eval($1);
		}
		$line='%'.$line;
		push @header,$line;
	}
	else {
		$inn=0 if $ispace;
		++$jnn if $inn==0;
		++$inn;
		chomp($line);
		$line=~s/^\s+//;
		@tmp=split /\s+/,$line;
		push @xa,$tmp[$x_col-1];
		push @ya,$tmp[$y_col-1];
		push @data,$tmp[$s_col-1];
		push @ua,$tmp[$u_col-1] if defined $u_col;
		push @va,$tmp[$v_col-1] if defined $v_col;
	}
}
chomp($tline);
$tline=~s/\#\-//;
$tline=~s/\s+$//;
$tline=~s/^\s+//;
@tita=split /$splitter/,$tline;

while ($opt_p=~s/(\$.*?=.*?;)//){
	print $1,"\n";
	eval($1);
}

$in=$inn unless defined($in);
$jn=$jnn unless defined($jn);
print " using $in $jn\n";
$info="";
die ("\$in not defined") if (!defined($in)); 
die ("\$jn not defined") if (!defined($jn)); 
$colortable="default" if (!defined($colortable));
$iplot=5 if (!defined($iplot));
$tickx_scale=1 if (!defined($tickx_scale));
$ticky_scale=1 if (!defined($ticky_scale));
$ismooth=0 if (!defined($ismooth));

open (PSOUT,">$outeps") || die ("cannot open outputfile $outeps");

$im=$in-1;
$jm=$jn-1;
#	@data=@all[0..($in*$jn-1)];
#	@xa=@all[$in*$jn..$in*$jn+$in*$jn-1];
#	@ya=@all[$in*$jn+$in*$jn..$in*$jn+2*$in*$jn-1];
($xmax,$xmin)=&maxmin(\@xa);
($ymax,$ymin)=&maxmin(\@ya);
($umax,$umin)=&maxmin(\@ua);
($vmax,$vmin)=&maxmin(\@va);
$info.="%data range of x:  $xmin $xmax\n";
$info.="%data range of y: $ymin $ymax\n";
if ($xmx-$xmn==0) {($xmx,$xmn)=($xmax,$xmin)}
if ($ymx-$ymn==0) {($ymx,$ymn)=($ymax,$ymin)}
$info.="%plot range of x:  $xmn $xmx\n";
$info.="%plot range of y: $ymn $ymx\n";
@all=();
($mean)=&subtract_mean(\@data) if $subtract_mean;
&rescale_lin(\@data,$rescale) if defined ($rescale);
$title=$opt_t or $title=$datafile.' :';
$title.='  '.@tita[$s_col-1];
$stats.='mean='.sprintf('%10.3e',$mean) if $subtract_mean;
($maxdata,$mindata)=&maxmin(\@data);
$stats.='  min='.sprintf('%10.3e',$mindata).'  max='.sprintf('%10.3e',$maxdata);
for ($j==0; $j<=$jn-1; $j++){
push @dataa,[@data[$j*$in..($j+1)*$in-1]];
push @xaa,[@xa[$j*$in..($j+1)*$in-1]];
push @yaa,[@ya[$j*$in..($j+1)*$in-1]];
push @uaa,[@ua[$j*$in..($j+1)*$in-1]];
push @vaa,[@va[$j*$in..($j+1)*$in-1]];
}
@xa=();
@ya=();
@za=();
@ua=();
@va=();

$info.="% range of data values: $mindata $maxdata \n";
if ($smin==0 && $smax==0) {
$smin=$mindata;
$smax=$maxdata;
if ($smin < 0 && $smax > 0){
	$smin=-$smax if -$smax<$smin;
	$smax=-$smin if -$smin>$smax;
}
}
print "\n color table will span: $smin $smax \n";
$nocon=1 if $smax==$smin;
$norm=1./($smax-$smin) unless $nocon;

unless ($nocon){
if (scalar(@contu)){
	$ncontours=@contu;
	$icmax=$ncontours+1;
}
else {
	$icmax=$ncontours+1;
	$dconval=($smax-$smin)/($icmax-1);
	for ($n=1; $n<=$icmax-1; $n++){
		$conval=$smin+($n-.5)*$dconval;
		push @contu, $conval;
	}
}
print "\n contour levels will be:\n";
foreach (@contu){ print "$_\n"}

&make_colors;
}
if ($pixplot==1){
&print_first_pix;
#&print_image_pix;
&print_image_pix_big;
&print_last_pix;
}
else{
if ($xs==0 or $ys==0){ 
    $xs=$xmx-$xmn;
    $ys=$ymx-$ymn;
	&reset_size;
}
&check_size;
$ibb=int(512*$xs+.9999);
$jbb=int(512*$ys+.9999);
if ($noaxes){
	$xf=$xs*12800/($xmx-$xmn);
	$yf=$ys*12800/($ymx-$ymn);
	$xoff=0;
	$yoff=0;
}
else{
	$xf=$xs*9600/($xmx-$xmn);
	$yf=$ys*9600/($ymx-$ymn);
	$xoff=2560;
	$yoff=1920;
}
&begin_regions;
print PSOUT "%begin colortable\n","% 0. $icmax $colortable 1.\n",
	 $colordefs,"%end colortable\n";
print PSOUT "%transformation: $xf*(x-$xmn), $yf*(y-$ymn)\n";
print PSOUT "%%EndProlog\n";
&arrow_header if scalar(@uaa);
print "please be patient...\n";
&axes unless $noaxes;
#sometimes it is better to put these next three calls before &axes
&contours unless $nocon;
&gridlines unless $nogrid;
&boundarylines;
&arrows if (scalar(@uaa) && !($noarrows));
&dots if defined($opt_dots) and -e $opt_dots;
&labels;
}
print PSOUT "\n","%%%%#######%%%%%%######\n";
print PSOUT "%COMMENTS IN DESC FILE:\n",@descom,"%END DESC COMMENTS\n";
print PSOUT "%COMMENTS IN DATA FILE:\n",@header,"%END DATA COMMENTS\n";
print PSOUT "%COLUMN TITLES WERE:\n";
$i=0;
for $q (@tita){
$i++;
print PSOUT "% $i  $q\n";
}
print PSOUT $info;
$tim = `date +%d-%m-%Y/%I.%M%p`;
print PSOUT "% conplot used in directory = $cwd\n";
print PSOUT "% used data file = $datafile\n";
print PSOUT "% x_col=$x_col  y_col=$y_col  s_col=$s_col u_col=$u_col v_col=$v_col\n";
print PSOUT "% used desc file = $descfile\n";
print PSOUT "% plotted by conplot on $tim";
close PSOUT;
print "all done!\n";
system ("gv -noantialias $outeps") if !($opt_b);
unlink $outeps if $outeps eq 'temp_kon.eps';

################################
######subroutines follow########
################################
sub make_colors{
$ctab="(\$red,\$green,\$blue)=&".$colortable."_ctab"."(\$v)";
#print $ctab,"\n";
if (! eval $ctab){
	 print " no colortable $colortable, change to default\n";
	$ctab="(\$red,\$green,\$blue)=&default_ctab"."(\$v)";
}
eval $ctab or die "no default color table";
local ($red,$green,$blue,$norm,$icdef);
@colors=();
$colordefs="";
$icdef=0;
for ($iv=0; $iv<$icmax; $iv++){
	$v=$iv/($icmax-1);
	eval $ctab;
	$temp='-9.0e+30';
	$temp=$contu[$iv-1] if $iv>0;
	$colordefs.=&make_color_def(\$icdef,$red,$green,$blue,$temp);
}
#print $colordefs;
}
###################################
sub make_color_def{
	my($icdef,$red,$green,$blue,$valu)=@_;
	my ($aline,$adummy);
	$aline=sprintf("%6.4f %6.4f %6.4f",$red,$green,$blue);
	$adummy="/C".sprintf("%3.3u",$$icdef);
	$rgb=sprintf("%2.2x%2.2x%2.2x",255*$red,255*$green,255*$blue);
	$hexcolor[$$icdef]=$rgb;
	$$icdef++;
	return "$adummy {  $aline C} bind def % $valu  $rgb\n";
}

################################
sub axes{
local($ibl,$jbl,$ibr,$jbr,$itr,$jtr,$itl,$jtl);
$astr.=&mv($xmn,$ymn);
$astr.=&ln($xmn,$ymx);
$astr.=&ln($xmx,$ymx);
$astr.=&ln($xmx,$ymn);
$astr.=&ln($xmn,$ymn);
$astr.="AX\n";
if (!defined($tickx)) {
	$tickx=int(($xmx-$xmn)/$tickx_scale/10+1)
}
if (!defined($ticky)) {
	$ticky=int(($ymx-$ymn)/$ticky_scale/10+1)
}
$xtic=0.;
while ($xtic > $xmn/$tickx_scale){
	$xtic-=$tickx;
}
while ($xtic < $xmn/$tickx_scale){
	$xtic+=$tickx;
}
while ($xtic <= $xmx/$tickx_scale){
	$astr.=&mv($xtic*$tickx_scale,$ymn);
	$astr.="0 -100 V \n";
	$astr.="0 -300 R \n";
	$astr.="($xtic) Cshow\n";
	$xtic+=$tickx;
}
$ytic=0.;
while ($ytic > $ymn/$ticky_scale){
	$ytic-=$ticky;
}
while ($ytic < $ymn/$ticky_scale){
	$ytic+=$ticky;
}
while ($ytic <= $ymx/$ticky_scale){
	$astr.=&mv($xmn,$ytic*$ticky_scale);
	$astr.="-100 0 V \n";
	$astr.="-300 0 R \n";
	$astr.="($ytic) Rshow\n";
	$ytic+=$ticky;
}

print PSOUT "/$font 480 selectfont\n";

if (defined($title)) {
	print PSOUT "
6400 12360 M
($title) Cshow
"
}
if (defined($stats)) {
print PSOUT "/$font 400 selectfont\n";
	print PSOUT "
6400 11800 M
($stats) Cshow
";
print PSOUT "/$font 480 selectfont\n";
}
if (defined($xlabel)) {
	print PSOUT "
7680 640 M
($xlabel) Cshow
"
}
if (defined($ylabel)) {
	print PSOUT "
1280 6400 M
90 rotate
($ylabel) show
-90 rotate

"
}

print PSOUT "%start_axes\n";
print PSOUT "ACOLOR\n";
print PSOUT "/$font 400 selectfont\n";

print PSOUT $astr;
$x=$xmn;
$i1=&ix;
$y=$ymn;
$j1=&jy;
$x=$xmx;
$i2=&ix-$i1;
$y=$ymx;
$j2=&jy-$j1;
print PSOUT  "[$i1 $j1 $i2 $j2] rectclip\n"; 
}
#####################################
sub labels{
print PSOUT 
"%start_labels
25 25  scale
/$font 12 selectfont
\%uncomment next four lines for example labels
\%105 105 BX
\%(this) show
\%305 205 BX
\%(that) show
"
}
################################
sub gridlines{
local($ibl,$jbl,$ibr,$jbr,$itr,$jtr,$itl,$jtl);
$cstr="";

for ($j=0; $j<=$jm; $j++){
	$i=0;
	$cstr.=&mv($xaa[$j][$i],$yaa[$j][$i]);
	for ($i=1; $i<=$im; $i++){
		$cstr.=&ln($xaa[$j][$i],$yaa[$j][$i]);
	}
	$cstr.="G\n";
}

for ($i=0; $i<=$im; $i++){
	$j=0;
	$cstr.=&mv($xaa[$j][$i],$yaa[$j][$i]);
	for ($j=1; $j<=$jm; $j++){
		$cstr.=&ln($xaa[$j][$i],$yaa[$j][$i]);
	}
	$cstr.="G\n";
}

print PSOUT "%start_grid\n";
print PSOUT "GCOLOR\n";
print PSOUT $cstr;
}
#################################################
sub boundarylines{
local($ibl,$jbl,$ibr,$jbr,$itr,$jtr,$itl,$jtl);
$cstr="";

for ($j=0; $j<=$jm; $j+=$jm){
	$i=0;
	$cstr.=&mv($xaa[$j][$i],$yaa[$j][$i]);
	for ($i=1; $i<=$im; $i++){
		$cstr.=&ln($xaa[$j][$i],$yaa[$j][$i]);
	}
	$cstr.="BG\n";
}

for ($i=0; $i<=$im; $i+=$im){
	$j=0;
	$cstr.=&mv($xaa[$j][$i],$yaa[$j][$i]);
	for ($j=1; $j<=$jm; $j++){
		$cstr.=&ln($xaa[$j][$i],$yaa[$j][$i]);
	}
	$cstr.="BG\n";
}

print PSOUT "%start_boundary\n";
print PSOUT "GCOLOR\n";
print PSOUT $cstr;
}
#########################
#########################
sub arrow_header{
my ($f,$temp);
$info.="%umin=$umin\n";
$info.="%umax=$umax\n";
$info.="%vmin=$vmin\n";
$info.="%vmax=$vmax\n";
if ($vec_scale==0){
	$vmax=-$vmin if -$vmin > $vmax;
	$umax=-$umin if -$umin > $umax;
	$f=1;
	$temp=1;
#	print "$xmax $xmin $in $umax\n";
	$f=($xmax-$xmin)/$in/$umax if $umax > 0;
	$temp=($ymax-$ymin)/$jn/$vmax if $vmax > 0;
	$f=$temp if $temp < $f;
	$vec_scale=$f;
}
$info.="%vec_scale=$vec_scale\n";
}
###############################
sub arrows{
my($x,$y,$u,$v,$ft,$fh,$vec_max);
$arrow_cen=.5 unless defined($arrow_cen);
$ft=(1-$arrow_cen)*$vec_scale;
$fh=$arrow_cen*$vec_scale;
$vstr="";
$vec_max=0;
for ($j=1; $j<=$jn; $j++){
	for ($i=1; $i<=$in; $i++){
		$x=$xaa[$j-1][$i-1];
		$y=$yaa[$j-1][$i-1];
		$u=$uaa[$j-1][$i-1];
		$v=$vaa[$j-1][$i-1];
		$temp=sqrt($u*$u+$v*$v);
		$vec_max=$temp if $temp > $vec_max;
		$vstr.=&arrow($x-$ft*$u,$y-$ft*$v,$x+$fh*$u,$y+$fh*$v);
	}
}
print PSOUT "%start_arrows, vec_max=$vec_max\n";
print PSOUT "VCOLOR\n";
print PSOUT $vstr;
$vstr=0;
}
################################
sub contours{
print PSOUT "%start_contour_fills\n";
print PSOUT "FLW\n";
$rstr="";
$cstr="";
for ($j=1; $j<=$jm; $j++){
	for ($i=1; $i<=$im; $i++){

		$xbl=$xaa[$j-1][$i-1];
		$ybl=$yaa[$j-1][$i-1];
		$zbl=$dataa[$j-1][$i-1];

		$xbr=$xaa[$j-1][$i];
		$ybr=$yaa[$j-1][$i];
		$zbr=$dataa[$j-1][$i];

		$xtl=$xaa[$j][$i-1];
		$ytl=$yaa[$j][$i-1];
		$ztl=$dataa[$j][$i-1];

		$xtr=$xaa[$j][$i];
		$ytr=$yaa[$j][$i];
		$ztr=$dataa[$j][$i];

		@zcr=($zbl,$zbr,$ztr,$ztl);
		($zcmx,$zcmn)=&maxmin(\@zcr);
		$thecolor=&onecolor($zcmx,$zcmn); #check if square is all one color
	 	if ($thecolor > 0){
				$rstr.=&cl($thecolor).&mv($xbl ,$ybl).&ln($xbr, $ybr).&ln($xtr, $ytr).&fl($xtl ,$ytl);
		}
		else{  #divide square into four triangles and plot the contours
		$xx=.25*($xbl+$xbr+$xtl+$xtr);
		$yy=.25*($ybl+$ybr+$ytl+$ytr);
		if ($j==1 or $j==$jm or $i==im or $i==1 or $ismooth==0 ){  
				$zz=.25*($zbl+$zbr+$ztl+$ztr);
		} 
		else{ # this does some funky interpolation to get the value at the center of the square 
				$zz=.25*($zbl+$zbr+$ztl+$ztr)
				+.03125*(
					$dataa[$j-2][$i-1]+$dataa[$j-2][$i]+
					$dataa[$j+1][$i-1]+$dataa[$j+1][$i]+
					$dataa[$j  ][$i-2]+$dataa[$j-1][$i-2]+
					$dataa[$j  ][$i+1]+$dataa[$j-1][$i+1])
				-.0625*(
					$dataa[$j-2][$i-2]+
					$dataa[$j+1][$i-2]+
					$dataa[$j-2][$i+1]+
					$dataa[$j+1][$i+1])
		}
		
  &tri($xbl,$xx,$xbr,$ybl,$yy,$ybr,$zbl,$zz,$zbr);
  &tri($xtr,$xx,$xbr,$ytr,$yy,$ybr,$ztr,$zz,$zbr);
  &tri($xtr,$xx,$xtl,$ytr,$yy,$ytl,$ztr,$zz,$ztl);
  &tri($xbl,$xx,$xtl,$ybl,$yy,$ytl,$zbl,$zz,$ztl);
	}
	}
}
print PSOUT $rstr;
print PSOUT "%start_contour_lines\n";
print PSOUT "SLW\n";
print PSOUT $cstr;
print PSOUT "[] 0 setdash\n";
}

####################################################
sub onecolor{
my($zcmx,$zcmn)=@_;
$cindx=1;
$clo=1;
$chi=1;
foreach $c (@contu){
 $cindx++;
	$clo=$cindx if ($zcmn > $c );
	$chi=$cindx if ($zcmx > $c );
}
if ($clo == $chi) {
		return ($clo)
}
else {
	return (-1)
}		
}


####################################################

sub tri{
my($x1,$x2,$x3,$y1,$y2,$y3,$z1,$z2,$z3)=@_;
my($xa,$xb,$ya,$yb,$xc,$yc,$xd,$yd,$m,$n,$c,$cindx,$cindxu);
($x2,$x1,$y2,$y1,$z2,$z1)=($x1,$x2,$y1,$y2,$z1,$z2) if ($z1 > $z2);
($x3,$x2,$y3,$y2,$z3,$z2)=($x2,$x3,$y2,$y3,$z2,$z3) if ($z2 > $z3);
($x2,$x1,$y2,$y1,$z2,$z1)=($x1,$x2,$y1,$y2,$z1,$z2) if ($z1 > $z2);
$xc=$x1;
$yc=$y1;
$xd=$x1;
$yd=$y1;
$m=0;
$n=0;
$cindx=0;
$cindxu=0;
foreach $c (@contu){
$cindx++;
	next if ($c > $z3);
	$cindxu=$cindx;
	 next if ($c < $z1 );
	 next if ($z3 == $z1 );
	$t13=($c-$z1)/($z3-$z1);
	$xa=$t13*($x3-$x1)+$x1;
	$ya=$t13*($y3-$y1)+$y1;
	$colcon="NCOLOR\n"; $colcon="PCOLOR\n" if ($c>0);
	if ($c < $z2) {
		next if ($z1==$z2);
	 	$m++;
		$t12=($c-$z1)/($z2-$z1);
		$xb=$t12*($x2-$x1)+$x1;
		$yb=$t12*($y2-$y1)+$y1;
	}
	else{
		next if ($z3==$z2);
		$n++;
		$t23=($c-$z2)/($z3-$z2);
		$xb=$t23*($x3-$x2)+$x2;
		$yb=$t23*($y3-$y2)+$y2;
	}
		$cstr.=$colcon.&mv($xa ,$ya).&sl($xb, $yb);
	if ($n==1){
		$rstr.=&cl($cindx).&mv($xa ,$ya).&ln($xb, $yb).&ln($x2, $y2).&ln($xc ,$yc).&fl($xd, $yd);
	}
	else{
		$rstr.=&cl($cindx).&mv($xa ,$ya).&ln($xb, $yb).&ln($xc ,$yc).&fl($xd, $yd);
	}
	$xc=$xb;
	$yc=$yb;
	$xd=$xa;
	$yd=$ya;
}
$cindx=$cindxu+1;
if ($m + $n ==0){
	$rstr.=&cl($cindx).&mv($x1, $y1).&ln($x2, $y2).&fl($x3, $y3);
}
elsif ($n==0){
	$rstr.=&cl($cindx).&mv($x3, $y3).&ln($x2, $y2).&ln($xc, $yc).&fl($xd, $yd);
}
else{
	$rstr.=&cl($cindx).&mv($x3, $y3).&ln($xc, $yc).&fl($xd, $yd);
}
}
####################################################


sub cl{
my($cindx)=@_;
my($cdef,$cadj);
$cadj=$cindx-1;
$cadj=$icmax-1 if $cadj>$icmax-1;
$cdef="C".sprintf("%3.3u",$cadj);
"$cdef N
";
}

sub ix{
#print "$x $xmn $xf $xoff\n";
int(($x-$xmn)*$xf+$xoff)
}

sub jy{
int(($y-$ymn)*$yf+$yoff)
}

sub mv{
local($x,$y)=@_;
my($i,$j);
$i=&ix;
$j=&jy;
"$i $j M
";
}

sub ln{
local($x,$y)=@_;
my($i,$j);
$i=&ix;
$j=&jy;
"$i $j L
";
}

sub fl{
local($x,$y)=@_;
my($i,$j);
$i=&ix;
$j=&jy;
"$i $j L P
";
}

sub sl{
local($x,$y)=@_;
my($i,$j);
$i=&ix;
$j=&jy;
"$i $j L S
";
}

sub arrow{
local($x1,$y1,$x2,$y2)=@_;
my($i1,$j1,$i2,$j2,$ia,$ja,$ib,$jb,$a);
$a=-.1; # determines arrowhead size
$x=$x1;
$y=$y1;
$i1=&ix;
$j1=&jy;
$x=$x2;
$y=$y2;
$i2=&ix;
$j2=&jy;
$ia=int($i2+$a*($i2-$i1)+$a*($j2-$j1));
$ja=int($j2-$a*($i2-$i1)+$a*($j2-$j1));
$ib=int($i2+$a*($i2-$i1)-$a*($j2-$j1));
$jb=int($j2+$a*($i2-$i1)+$a*($j2-$j1));
"$i1 $j1 M
$i2 $j2 L VS 
$i2 $j2 M
$ia $ja L VS
$i2 $j2 M
$ib $jb L VS
";
}
####################################################

sub begin_regions{
$conyes='';
$gridyes='';
$arrowyes='';
$conyes='%' if ($iplot==3 or $iplot==5);
$gridyes='%' if ($iplot>=4);
$gridyes='%' if $nocon;
$arrowyes='%' unless $noarrow;
$gridyes='' if $nogrid;
print PSOUT
"%!PS-Adobe-2.0 EPSF-2.0
%%Creator: conplot 
%%DocumentFonts: $font 
%%BoundingBox: 0 0 $ibb $jbb 
%%EndComments
/vshift -46 def
/M {moveto} bind def
/L {lineto} bind def
/R {rmoveto} bind def
/V {rlineto} bind def
/C {setrgbcolor} bind def
/CMYK {setcmykcolor} bind def
/N {newpath} bind def
/S {stroke} bind def
$conyes/S {newpath} bind def %uncomment to NOT stroke contours
/VS {stroke} bind def
$arrowyes/VS {newpath} bind def %uncomment to NOT stroke arrows 
/P {closepath  fill newpath} bind def
%/P {closepath gsave stroke grestore fill newpath} bind def % works better sometimes
%/P {newpath} bind def  %uncomment to NOT fill regions
/G {stroke} bind def
$gridyes/G {newpath} bind def %uncomment to NOT stroke grid
/BG {stroke} bind def 
%/BG {newpath} bind def %uncomment to NOT stroke boundary 
/AX {stroke} bind def
%/AX {newpath} bind def %uncomment to NOT stroke axes
/Ashow {show} bind def
%/Ashow {} bind def
/BX {M save 1 1 1 setrgbcolor 24 0 V 0 12 V -24 0 V P restore
 0 1 R 0 0 0 setrgbcolor} bind def %used for labels, see end of file
/LA {show}
/Lshow { currentpoint AX M
  0 vshift R Ashow } def
/Rshow { currentpoint AX M
  dup stringwidth pop neg vshift R Ashow } def
/Cshow { currentpoint AX M
  dup stringwidth pop -2 div vshift R Ashow } def
/GCOLOR {  0.0000 0.0000 0.0000 C [] 0 setdash} bind def %grid color
/ACOLOR {  0.0000 0.0000 0.0000 C} bind def % axis color
/VCOLOR {  0.0000 0.5000 0.0000 C} bind def % arrow color
/PCOLOR {  1.0000 0.0000 0.0000 C [] 0 setdash} bind def % positive contour color
/NCOLOR {  0.0000 0.0000 1.0000 C [100 100] 0 setdash} bind def % negative contour color
/SLW {6 setlinewidth} bind def
/FLW {6 setlinewidth} bind def % set to 25 for some ghostscript conversions
SLW
0 0 translate
0.040 0.040 scale
0 setgray
"
}

####################################################
sub reset_size{
    if ($ys/$xs>0.3 and $ys/$xs<3.0) {
        if ($ys>$xs){ 
            $xs=$xs/$ys;
            $ys=1
		}
        else{
            $ys=$ys/$xs;
            $xs=1
		}
	}
    else{
        $ys=1;
        $xs=1
	}
}

####################################################
sub check_size{
if ($xs<=0 or $xs>=1) {$xs=1};
if ($ys<=0 or $ys>=1) {$ys=1};
#print "\n plot dimensions will be: $xs $ys \n";
}
####################################################

################################
sub grey_ctab{
my($v)=@_;
my($red,$green,$blue);
$red=$v;
$green=$v;
$blue=$v;
return ($red,$green,$blue);
}
################################
sub green_ctab{
my($v)=@_;
my($red,$green,$blue);
$red=(-6*($v*$v*$v/3-$v*$v/2));
$green=$v*$v*$v-3*$v*$v+3*$v;
$blue=.3+.65*$v;
return ($red,$green,$blue);
}
################################
sub luciano_ctab{
my($v)=@_;
my($red,$green,$blue);
$blue=sqrt($v);
$red=1;
$green=sqrt($v);
return ($red,$green,$blue);
}
################################
sub red_ctab{
my($v)=@_;
my($red,$green,$blue);
$blue=(-6*($v*$v*$v/3-$v*$v/2));
$red=$v*$v*$v-3*$v*$v+3*$v;
$green=.3+.65*$v;
return ($red,$green,$blue);
}
################################
sub blue_ctab{
my($v)=@_;
my($red,$green,$blue);
$green=(-6*($v*$v*$v/3-$v*$v/2));
$blue=$v*$v*$v-3*$v*$v+3*$v;
$red=.3+.65*$v;
return ($red,$green,$blue);
}
################################
sub invgreen_ctab{
my($v)=@_;
my($red,$green,$blue);
$v=1-$v;
$red=(-6*($v*$v*$v/3-$v*$v/2));
$green=$v*$v*$v-3*$v*$v+3*$v;
$blue=.3+.65*$v;
return ($red,$green,$blue);
}
################################
sub invred_ctab{
my($v)=@_;
my($red,$green,$blue);
$v=1-$v;
$blue=(-6*($v*$v*$v/3-$v*$v/2));
$red=$v*$v*$v-3*$v*$v+3*$v;
$green=.3+.65*$v;
return ($red,$green,$blue);
}
################################
sub invblue_ctab{
my($v)=@_;
my($red,$green,$blue);
$v=1-$v;
$green=(-6*($v*$v*$v/3-$v*$v/2));
$blue=$v*$v*$v-3*$v*$v+3*$v;
$red=.3+.65*$v;
return ($red,$green,$blue);
}
################################
################################
sub default_ctab{
my($v)=@_;
my ($q,$ex,$s);
my($red,$green,$blue);
if ($v < .5 ){
$q=(2*$v);
$ex=$q+(1-$q)*.5;
$q=$q**$ex;
$red=12*($q*$q*$q/3-$q*$q/2+$q/4);
$blue=1.;
$green=3*$q*$q*(1-2*$q/3);
}
else {
$s=(2-2*$v);
$ex=$s+(1-$s)*.5;
$s=$s**$ex;
$blue=12*($s*$s*$s/3-$s*$s/2+$s/4);
$red=1.;
$green=3*$s*$s*(1-2*$s/3);
}
return ($red,$green,$blue);
}
################################
sub blue_red_ctab{
my($v)=@_;
my($red,$green,$blue);
$red=$v;
$green=0;
$blue=1-$v;
return ($red,$green,$blue);
}
################################
sub blue_white_red_ctab{
my($v)=@_;
my($red,$green,$blue);
$green=&min(1.9*$v,1.9*(1-$v));
$red=.95;
$blue=.95;
$blue=$green if ($v>0.5);
$red=$green if ($v<0.5);
return ($red,$green,$blue);
}
################################
sub blue_grey_red_ctab{
my($v)=@_;
my($red,$green,$blue);
$red=$v**2;
$green=$v*(1-$v);
$blue=(1-$v)**2;
return ($red,$green,$blue);
}
################################
sub grey_short_ctab{
my($v)=@_;
my($red,$green,$blue);
$red=.3+.6*$v;
$green=$red;
$blue=$red; 
return ($red,$green,$blue);
}
################################
sub cyan_white_yellow_ctab{
my($v)=@_;
my($red,$green,$blue);
$red=.95;
$blue=.95;
$green=min(1.9*$v,1.9*(1-$v));
$blue=$green if ($v > .5);
$red=$green if ($v <.5);
$green=max($green,1-1.9*(1-$v));
$green=max($green,1-1.9*$v);
return ($red,$green,$blue);
}
###########################
sub min{
my (@dummy);
@dummy=sort @_;
return $dummy[0];
}
###########################
sub max{
my (@dummy);
@dummy=sort @_;
return $dummy[-1];
}
###################################
sub maxmin{
my ($da)=@_;
my ($max,$min);
$max=@$da[0];
$min=$max;
foreach $one (@$da){
if ($one > $max) {$max=$one};
if ($one < $min) {$min=$one};
}
chomp($max);
chomp($min);
return ($max,$min);
}
###################################
sub subtract_mean{
my ($da)=@_;
my ($mean);
foreach $one (@$da){
$mean+=$one;
}
$mean=$mean/(scalar(@$da));
foreach $one (@$da){
$one-=$mean;
}
return ($mean);
}
###################################
sub rescale_lin{
my ($da,$fac)=@_;
foreach $one (@$da){
$one=$one*$fac;
}
}

################################
sub print_first_pix{
$isize=int $xs*512/$in;
$jsize=int $ys*512/$jn;
$isize=1 if defined($opt_pxx);
$jsize=1 if defined($opt_pxx);
#$ib=$in*$isize-1;
#$jb=$jn*$jsize-1;
$ib=$in*$isize;
$jb=$jn*$jsize;
print PSOUT
"%!PS-Adobe-2.0 EPSF-2.0
%%Creator: conplot 
%%DocumentFonts: $font 
%%BoundingBox: 0 0 $ib $jb %size close to 512 512, but consistent with scale
";
print PSOUT "%begin colortable\n","% 0. $icmax $colortable 1.\n",
	 $colordefs,"%end colortable\n";
}

################################
sub print_last_pix{
$rin=1./$in;
$rjn=1./$jn;
print PSOUT "
%uncomment the next four lines if you want a frame
%0 setlinewidth
%0 0
%$in $jn 
%rectstroke\n
%uncomment the next four lines if you want a comment
%$rin $rjn scale
%/$font 12 selectfont
%100 100 moveto
%(a comment) show\n  
";
}

################################
sub print_image_pix_big{
print PSOUT "
$isize $jsize scale  %these have been set to integers, the number of pts per pixel
%will write strips of colors:\n
";
for ($j=0; $j<=$jm; $j++){
$jj=$j;
$jj=$jm-$j if $opt_jflip;
print PSOUT "0 1 translate\n" if $j>=1;
print PSOUT "<\n";
for ($i=0; $i<=$im; $i++){
$val=$dataa[$jj][$i];
$cindx=0;
foreach $c (@contu){
	last if ($val < $c);
	$cindx++;
}
print PSOUT $hexcolor[$cindx],"\n";
}
print PSOUT ">
$in 1 
8 
[1 0 0 1 0 0]  
{} 
false 
3 
colorimage

"; 
}
}
################################
sub print_image_pix{
#obsolete version, works only for small arrays (less than 100x100)
#print_image_pix_big writes in strips
print PSOUT "
$isize $jsize scale  %these have been set to integers, the number of pts per pixel
<\n";
for ($j=0; $j<=$jm; $j++){
for ($i=0; $i<=$im; $i++){
$val=$dataa[$j][$i];
$cindx=0;
foreach $c (@contu){
	last if ($val < $c);
	$cindx++;
}
print PSOUT $hexcolor[$cindx],"\n";
}
}
print PSOUT ">
$in $jn %describes size of the array of pixels
 8 %8 bit colors, requiring two hexadecimal symbols, e.g. ff, a9, etc.  
[1 0 0 1 0 0]  %defines orientation? ignore this...
{} %ignore this too!
false %and this!
3 % use three colors, red green blue, e.g., ffa9d1
colorimage % plot the image\n"
}
sub dots{
print PSOUT "%start_dots
/hpt_ 100 def
/vpt_ 100 def
/hpt hpt_ def
/vpt vpt_ def
/vpt2 vpt 2 mul def
/hpt2 hpt 2 mul def
/Pnt { stroke [] 0 setdash
   gsave 1 setlinecap M 0 0 V stroke grestore } def
/Dia { stroke [] 0 setdash 2 copy vpt add M
  hpt neg vpt neg V hpt vpt neg V
  hpt vpt V hpt neg vpt V closepath stroke
  Pnt } def
/Pls { stroke [] 0 setdash vpt sub M 0 vpt2 V
  currentpoint stroke M
  hpt neg vpt neg R hpt2 0 V stroke
  } def
/Box { stroke [] 0 setdash 2 copy exch hpt sub exch vpt add M
  0 vpt2 neg V hpt2 0 V 0 vpt2 V
  hpt2 neg 0 V closepath stroke
  Pnt } def
/Crs { stroke [] 0 setdash exch hpt sub exch vpt add M
  hpt2 vpt2 neg V currentpoint stroke M
  hpt2 neg 0 R hpt2 vpt2 V stroke } def
/TriU { stroke [] 0 setdash 2 copy vpt 1.12 mul add M
  hpt neg vpt -1.62 mul V
  hpt 2 mul 0 V
  hpt neg vpt 1.62 mul V closepath stroke
  Pnt  } def
/Star { 2 copy Pls Crs } def
/BoxF { stroke [] 0 setdash exch hpt sub exch vpt add M
  0 vpt2 neg V  hpt2 0 V  0 vpt2 V
  hpt2 neg 0 V  closepath fill } def
/TriUF { stroke [] 0 setdash vpt 1.12 mul add M
  hpt neg vpt -1.62 mul V
  hpt 2 mul 0 V
  hpt neg vpt 1.62 mul V closepath fill } def
/TriD { stroke [] 0 setdash 2 copy vpt 1.12 mul sub M
  hpt neg vpt 1.62 mul V
  hpt 2 mul 0 V
  hpt neg vpt -1.62 mul V closepath stroke
  Pnt  } def
/TriDF { stroke [] 0 setdash vpt 1.12 mul sub M
  hpt neg vpt 1.62 mul V
  hpt 2 mul 0 V
  hpt neg vpt -1.62 mul V closepath fill} def
/DiaF { stroke [] 0 setdash vpt add M
  hpt neg vpt neg V hpt vpt neg V
  hpt vpt V hpt neg vpt V closepath fill } def
/Pent { stroke [] 0 setdash 2 copy gsave
  translate 0 hpt M 4 {72 rotate 0 hpt L} repeat
  closepath stroke grestore Pnt } def
/PentF { stroke [] 0 setdash gsave
  translate 0 hpt M 4 {72 rotate 0 hpt L} repeat
  closepath fill grestore } def
/Circle { stroke [] 0 setdash 2 copy
  hpt 0 360 arc stroke Pnt } def
/CircleF { stroke [] 0 setdash hpt 0 360 arc fill } def
/Dot1 {0 0 0 C Pnt} def
/Dot2 {0 1 0 C DiaF} def
/Dot3 {1 0 0 C Pls} def
/Dot4 {0 0 1 C TriDF} def
/Dot5 {1 1 0 C PentF} def
/Dot6 {0 1 1 C BoxF} def
/Dot7 {1 0 1 C CircleF} def
/Dot8 {0 0 0 C TriUF} def
";
open DOTS,"<$opt_dots" or die "no dot file $opt_dots\n";
print "plotting dots for $opt_dots:\n";
foreach $line (<DOTS>) {
		print $line;
		chomp($line);
		$line=~s/^\s+//;
		($x,$y,$doti)=split /\s+/,$line;
		$dotline=&ix.' '.&jy.' Dot'.$doti;
		print PSOUT $dotline,"\n";  
}
}