#!/opt/perl/bin/perl 
use warnings;
use strict;
use File::Basename;
use Getopt::Long;
use Cwd;
$ENV{GDFONTPATH}="/opt/lib/fonts/TrueType";
#$ENV{LD_LIBRARY_PATH}="$ENV{LD_LIBRARY_PATH}:/opt/lib/gd-2.0.23";
#print "$ENV{LD_LIBRARY_PATH}\n";
#exit;
#get this program name
my $progname = basename ($0);

#define usage
my $usage= qq {
$progname is a wraper of gnuplot to plot timeseries 
obtained with the program dqrfit.
The output of $progname  is to the screen ( default ) or 
to a file (postscript or jpg or both, see options).
$progname takes advantage of the information in the header 
provided by dqrfit to allow you to custumize the plot
using the functions and the parameters founded in the fitting.

$progname plots standard timeseries files in 3 column format:
	T Y SIGY
T: time (decimal year) ; Y: data ; SIGY: Y error.
comments are allowed with \# at the begining of the line.

Files must have the same name and extensions lat,lon or rad 
in order to allow $progname to automatically look 
and plot  the 3 components.

You only need to provide the path to 1 of the files 
to plot the 3 components.
Usage:
$progname -i 	FILENAME.[lat][lon][rad] [OPTIONS..]

OPTIONS:
    -xrange 		XMIN:XMAX 
    -yrange 		YMIN:YMAX
    -xunit		UNITSTRING (ex. "Years" , "DOY") 
    -yunit		UNITSTRING (ex. "mm" , "cm")
    -vgrid      	Turn on vertical grid
    -hgrid		Turn on horizontal grid
    -width		WIDTH in inches
    -heigth		HEIGTH in inches
    -vscl       	SCL vertical scale. Ex 1.5 is 1.5 bigger.
    -hscl       	SCL horizontal scale. Ex 1.5 is 1.5 bigger.
    -font   		FONTNAME   (default: Helvetica)	
    -fontsize 		SIZE     (default: 8 );
    -llr		[lat,lon,rad] select component to plot.
   				Ex. "-llr  lat,rad" (plot lat and rad only) 
     				Otherwise it will plot all of them by default.
    -ffun			plot fitted function
    -pf			[func_name1,func_name2,...] plot fitted functions 
    			(if they are available)
    -line		plot line
    -season 		plot season
    -det		detrend	 ( removes trend )
    -des		de-season (removes seasonal variations )
    -deo		de-offset (removes offsets )
    -rmf		[func_name1,func_name2,...] remove fitted functions from data (if they are available)
    -offsets 		mark offsets with a vertical line
    -resd		plot residuals *
    -jumps		plot sum of jumps
    -eb			plot with errorbars
    -ps			POINT_SIZE ( ex. 0.5 is half size of the default) 
    -lw			LINE_WIDTH ( ex. 0.5 is half size of the default )
    -jpg		make jpg file
    -png                make png file
    -term  		[ps|x11] output device
    -psfile 		PS_FILENAME
    -orient 		[portrait|landscape]  (default:portrait)
    -rotjpg		rotate jpg 90 degrees
    -show		[title,header,footer,logo,timestamp,ps]
    			title 		: SITE name
    			header		: statistics
    			footer		: units
    			logo		: agency (Geodesy Lab & PANGA,Central Washington University)
    			timestamp	: logo + time stamp
    			ps		: show postscript file using gv
    -[help|h|H]
	
* : "-resd" Actually show the difference between the data values 
    and the fitted function,the error bars still belong
    to the original data.
    To plot the true uncertainties of the residuals 
    you should plot the file with the residuals without -resd or 
    -de*  flags.

Example: 
$progname -i ALBH.lon  -xrange 1998:2004 -yrange -20:20 -term ps -jpg -ffun
 plots the 3 components lat,lon, and rad and fitted function between 1998 and  2004 , 
 and between -20 and 20 and makes the jpg file and the ps file.



};
#default settings
my($ymin,$ymax,
   $tmin,$tmax,
   $nymin,$nymax,
   $eymin,$eymax,
   $uymin,$uymax,
   $iymin,$iymax) = (0)x12;
my $llr='';
my $site='';
my $SITE='';
my $dtype='';
my $show='';
my $showps='';
my $showheader=1;
my $showtitle=1;
my $showtype=0;
my $showfooter=1;
my $nrange='';
my $erange='';
my $urange='';
my $yrange='';
my $xrange='';
my $mxtics=2;
my $mytics=2;
my $xtics=1;
my $ytics=5;
my $xunit='Years';
my $yunit='';
my $setgrid=''; 
my $xgrid=0;
my $ygrid=0;
my $help=0;
my $det=0;
my $des=0;
my $def=0;
my $rf='';
my @rmfun=();
my $dump= 0;
my $scriptfile = '';
my $line=0;
my $season=0;
my $jumps=0;
my $offsets=0;
my $event="step";
my $deo=0;
my $ffun=0;
my $pfun='';
my @plotfun=();
my $resd=0;
my $data=0;
my $eb=0;
my $pt="pt 7";
my $ps=0.2;
my $lw=2;
my $with='';
my $nsmpl=1;
my $dir='';
my $outdir='';
my $file='';
my $outfile='';
my $jpg=0;
my $png=0;
my $pfile='';
my $par;
# Default size of a landscape PostScript plot is 10 inches wide and 7 inches high
my $orient='portrait';
my $W=7.0;
my $H=10.0;
my $AR = $H/$W; 
my $width=800;
$width=500;
my $heigth=0;
my $hscl=1.0;
my $vscl=1.0;
my $dim='';
my $hsize=1;
my $vsize=1;
my $media = "letter";
my $rotjpg=0;
my $tstamp=0;
my $logo=1;
my $PTITLE='';
my $FOOTER='';
my $LOGO = '';
my $fontname="Helvetica";
my $fontsize=12;
my $here=cwd();
my @ftype=();
my $term='x11';
my $termstr='x11';
my $nargs=@ARGV;
my $ylabel = '';
my @comps=( "lat","lon","rad");
my %NEV=( "lat" => "N",
	  "lon" => "E",
	  "rad" => "U",
	  "x"	=> "X",
	  "y"	=> "Y",
	  "z"	=> "Z",
	  "l"	=> "L",
	  "t"	=> "T",
	  "u"	=> "U"
	);

my %NEU=( "lat" => "North",
	  "lon" => "East",
	  "rad" => "Up", 
	  "x"	=> "X",
	  "y"	=> "Y",
	  "z"	=> "Z",
	  "l"	=> "Length",
	  "t"	=> "Traverse",
	  "u"	=> "Up"
	 );

my %DTYPE=( "raw" => "RAW",
	    "ldet" => "DETRENDED",
	    "data" => "DATA",
	    "resd" => "RESIDUALS",
	    "crsd" => "CLEANED-RESD",
 	    "ltu"  => "BASELINE"	
	);

#Get command line options
GetOptions( "i=s"    		=> \$file ,
            "term=s"    	=> \$term ,
            "psfile=s"  	=> \$outfile , #ps file
            "outdir=s"  	=> \$outdir , 
            "orient=s" 		=> \$orient ,  #portrait | landscape
	    "rotjpg"    	=> \$rotjpg,
	    "font=s"		=> \$fontname,
	    "fontsize=i" 	=> \$fontsize,
	    "llr=s"		=> \$llr,
	    "det"		=> \$det,
	    "des"  		=> \$des,
	    "deo"		=> \$deo,	    
    	    "rmf=s"		=> \$rf,
	    "dump!"		=> \$dump,
	    "script=s"		=> \$scriptfile,		
	    "par"		=> \$par,
	    "resd"		=> \$resd,
	    "ffun"		=> \$ffun,
	    "pf=s"		=> \$pfun,		
	    "line"		=> \$line,
	    "season"		=> \$season,
	    "jumps"		=> \$jumps,
 	    "offsets"		=> \$offsets,
 	    "events=s"		=> \$event,
	    "xrange=s"  	=> \$xrange ,
	    "yrange=s"  	=> \$yrange ,
	    "nrange=s"      => \$nrange ,
	    "erange=s"      => \$erange ,
	    "urange=s"      => \$urange ,
	    "xunit=s"   	=> \$xunit ,
	    "yunit=s"   	=> \$yunit ,	    
	    "ylabel=s"   	=> \$ylabel,	    
	    "width=f"    	=> \$width,
	    "heigth=f"    	=> \$heigth,
	    "dim=s"		=> \$dim,
	    "hscl=f"    	=> \$hscl,
	    "vscl=f"    	=> \$vscl,    
	    "media=s"    	=> \$media,	    
	    "vgrid"     	=> \$xgrid ,
	    "hgrid"     	=> \$ygrid ,
	    "eb"		=> \$eb,
	    "ps=f"		=> \$ps,
	    "lw=f"		=> \$lw,
        "jpg"		=> \$jpg,
	"png"           => \$png,
	    "show=s"		=> \$show,
	    #"showtstamp"	=> \$tstamp,
	    #"showlogo!"	=> \$logo,
	    #"showheader!"   	=> \$showheader ,
	    #"showtitle!"    	=> \$showtitle ,
	    #"showfooter!"    	=> \$showfooter ,
	    #"showps"   		=> \$showps ,
            "help"      	=> \$help
           ) || die "error in options";
if ( ! $nargs || $help ) { print "$usage\n"; exit; };
# validate arguments
if ( ! -e $file ) { die "$usage\n" }; 

$dir=dirname($file);
$outdir = $dir if ! $outdir;
$file=basename($file);



if ( $xgrid ) {
	if ( ! $setgrid ) { $setgrid="set grid"};
	$setgrid="$setgrid xtics mxtics";
};
if ( $ygrid ) {
	if ( ! $setgrid ) { $setgrid="set grid"};
	$setgrid="$setgrid ytics mytics";
};
my @c=();
if ( $llr ) { 
	my @fields=split(/,/,join(',',$llr));
	foreach (@fields ) {
		$_=lc($_);
		push (@c,$_) if /(lat|lon|rad|x|y|z|l|t|u)/ ;
	};
	#if ( @c ) { @comps=@c };
};
if ( $show ) {
	$showps=0;
	$showheader=0;
	$showtitle=0;
	$showtype=0;
	$showfooter=0;
	$logo=0;
	$tstamp=0;
	my @fields = split (/,/,join(',',$show));
	foreach my $choice ( @fields ) {
		print "$choice\n";
		$showps=1     if $choice eq 'ps';
		$showheader=1 if $choice =~ /header/;
		$showtitle=1  if $choice =~ /title/;
		$showtype=1  if $choice =~ /dtype/;
		$showfooter=1 if $choice =~ /footer/;
		$tstamp =1    if $choice =~ /(time|stamp)/;
		$logo = 1 if $choice =~ /logo/;
	};
	print "$showtitle $showheader $showfooter\n";
}


if ( $det )  { $rf .= "line" };
if ( $des )  { $rf .= ",season" };
if ( $deo )  { $rf .= ",jumps" }; 
if ( $resd ) { $rf = "f" };
if ( $rf ) {
	$rf = lc($rf);
	@rmfun=split(/,/,join(',',$rf));
	$det=1 if $rf =~ /line/;
	$deo=1 if $rf =~ /jumps/;
#	push(@rmfun,@f);
};
if ( $ffun ) { $pfun = "f" };
if ( $pfun ) {
	$pfun = lc($pfun);
	@plotfun=split(/,/,join(',',$pfun));
};  
$orient = lc($orient);
$orient = "landscape" if $orient =~ /^l/;
$orient = "portrait" if $orient =~ /^p/;
if ( $dim  ) {
	if ( $dim =~ /\d+x\d+/ ) {
		my ($w,$h)=split(/x/,$dim);
		if ( $w < 0 || $h < 0 ) {  
			warn "invalid dimensions\n";
		} else { $width = $w; $heigth=$h };
	} else {
		warn "invalid format try 800x600\n";
	}; 
};
my $sx=$width*$hscl;
my $sy=$sx*$AR*$vscl;
$sy = $heigth*$vscl if $heigth;
if ( $orient =~ /landscape/ ) { $sx=$sy; $sy=$width };
#print @rmfun ;
#exit;
$term=lc($term);
if ( $jpg ) { $term = 'ps' };
if ( $png ) { $term = 'ps' };
if ( $term !~ /(ps|pdf|png|jpeg|x11|x)/ ) { die "invalid terminal $term [term = ps pdf x11 ]\n" };
if ( $term eq 'ps' ) {
	$termstr="postscript $orient enhanced color solid \"$fontname\" $fontsize"; 
}elsif ( $term eq 'pdf' ) {
	 $termstr="pdf  $fontname $fontsize enhanced";
	 # set terminal pdf {fname "<font>"} {fsize <fontsize>}
         #               {{no}enhanced}
} elsif (  $term eq 'jpeg' ) {
	my $sx=$width*$hscl;
	my $sy=$sx*$AR*$vscl;
	$sy = $heigth*$vscl if $heigth;
	if ( $orient =~ /landscape/ ) { $sx=$sy; $sy=$width };
	$termstr="jpeg nointerlace font $fontname $fontsize size $sx,$sy crop enhanced";
	# set terminal jpeg ",set terminal jpeg large font arial size 800,600
       #                 {{no}interlace}",
       #                 {tiny | small | medium | large | giant}",
       #                 {font <face> {<pointsize>}}",
       #                 {size <x>,<y>} {{no}crop}",
       #                 {{no}enhanced}",
       #                 {<color0> <color1> <color2> ...}",
	#Five basic fonts are supported directly by the gd library. These are",
	#" `tiny` (5x8 pixels), `small` (6x12 pixels), `medium`, (7x13 Bold), ",
	#" `large` (8x16) or `giant` (9x15 pixels). These fonts cannot be scaled",
	#" or rotated (pure horizontal or vertical text only).",
} elsif (  $term eq 'png' ) {
	my $sx=$width*$hscl;
	my $sy=$sx*$AR*$vscl;
	$sy = $heigth*$vscl if $heigth;
	if ( $orient =~ /landscape/ ) { $sx=$sy; $sy=$width };
	$termstr="png nointerlace font $fontname $fontsize size $sx,$sy crop enhanced";
}else { 
	$fontname = "Helvetica";
	$fontsize = 8;
	$termstr="x11 enhanced font \"Helvetica,8\" ";
	$pt=" ";
};
if ( $xrange ) {
	  die "$usage \n" unless $xrange =~ /(\d+|\*)\:(\d+|\*)/ ; 
	  ($tmin,$tmax)=split(/\:/,$xrange);
	  $xrange="xrange [$xrange]\n";
	  print "$xrange\n";
};

#if ( $det ) { $yrange="set yrange [-30:30]\n"; };

if ( $yrange ) {
	  die "$usage \n" unless $yrange =~ /(\d+|\*)\:(\d+|\*)/ ; 
	  ($ymin,$ymax)=split(/\:/,$yrange);
	  $yrange="yrange [$yrange]\n";
  	  print "$yrange\n";
};
if ( $nrange ) {
	  die "$usage \n" unless $nrange =~ /(\d+|\*)\:(\d+|\*)/ ; 
	  ($nymin,$nymax)=split(/\:/,$nrange);
	  $nrange="yrange [$nrange]\n";
  	  print "$yrange\n";
};
if ( $erange ) {
	  die "$usage \n" unless $erange =~ /(\d+|\*)\:(\d+|\*)/ ; 
	  ($eymin,$eymax)=split(/\:/,$erange);
	  $erange="yrange [$erange]\n";
  	  print "$erange\n";
};
if ( $urange ) {
	  die "$usage \n" unless $urange =~ /(\d+|\*)\:(\d+|\*)/ ; 
	  ($uymin,$uymax)=split(/\:/,$urange);
	  $urange="yrange [$urange]\n";
  	  print "$urange\n";
};


my @events;
if ( $event ) {
	@events = split(/,/,join(',',$event ));
};
my @YLABELS=();
if ( $ylabel ) {
	@YLABELS=split(/,/,join(',',$ylabel));
};
# try to determine file type 
@ftype=split(/\./,$file);
#print "Data type @ftype\n";
if ( $ftype[0] !~ /\b\w{4}\b/ ) { warn "WARNING!: $site it is not a 4 character site name \n"; };
if ( $ftype[$#ftype] !~ /(lat|lon|rad|x|y|z)/ ) { warn "WARNING !: $site has no lat, lon or rad extension\n" };

$site = shift @ftype;
$SITE=uc($site);

#next field is a component
if ( $ftype[$#ftype] =~ /(x|y|z)/ ) { @comps = qw( x y z ) };
#if ( $ftype[$#ftype] =~ /(l|t|u)/ ) { @comps = qw( l t u ) };
if ( @c ) { @comps=@c };

if ( $ftype[0] =~ /(lat|lon|rad|x|y|z)/ ) { $dtype="raw" }
else { $dtype=shift @ftype };

#print "Plotting $SITE , Data type: $dtype , comp: $ftype[$#ftype]\n";
if ( $resd || $dtype eq "resd" ) {
	$det=0;
	$des=0;
	$deo=0;
};

#make file names for 3 components
my %FILE=();
my %HEADER=();
my %PAR=();
my %FUNC=();
my %OFFSETS;
my (%JUMPS,%EVENTS);
my %TITLE=();
my %MINMAX=();
my $f = '';
my $t='';
@c=@comps;
@comps=();

my $header= "$SITE $DTYPE{$dtype}";

print "$header\n";
my ($Ymin,$Ymax);
foreach my $comp ( @c ) {
	if ( $dtype eq  "raw" and -e "$dir/$site.$comp" ) {
		$f="$dir/$site.$comp";
	}elsif ( -e "$dir/$site.$dtype.$comp" ) {
		$f="$dir/$site.$dtype.$comp";
	}else { $f='' };
        #print "dir $dir $site $dtype $comp $f\n";
	next if ! $f ;
	push(@comps,$comp);
	if ( $ylabel ) {
		if ( @YLABELS == 1 ) { $NEU{$comp} = $YLABELS[0] }
		else { $NEU{$comp} = shift @YLABELS};
	};
	$FILE{$comp} = $f;
	push( @{$MINMAX{$comp}},fminmax($f,$det,$deo,$comp));
	($Ymin,$Ymax,$tmin ,$tmax)= @{$MINMAX{$comp}} if ! $xrange;
	if ( $par ) {
		push(@{$HEADER{$comp}},getheader("$here/fit/$site.fit.par",$comp));
	} else {		
		push(@{$HEADER{$comp}},getheader($f,$comp));
	};
	if ( @{$HEADER{$comp}} ) {
		push(@{$PAR{$comp}},getpar(@{$HEADER{$comp}}));
		push(@{$FUNC{$comp}},getfunc(@{$HEADER{$comp}}));
		push(@{$OFFSETS{$comp}},getoffsets('OFFSET-',@{$HEADER{$comp}}))  if $offsets;
		push(@{$JUMPS{$comp}},getoffsets('\sOFFSET\s',@{$HEADER{$comp}})) if $jumps;
		#print (map{"$_\n"} @{$OFFSETS{$comp}});
		my ($vel, $velsig, $rms, $chis_dof, $dof, $tspan)=split (/,/,join(",",${$PAR{$comp}}[0])); 
		#print "$vel, $velsig, $rms, $chis_dof, $dof, $tspan\n";
		#$t = "'{/$fontname=$fontsize V=$vel {/Symbol \261} $velsig mm/yr  WRMS=$rms mm  Chi^2/DOF=$chis_dof  DOF=$dof TSPAN=$tspan}'";
		$t = "'{/$fontname=$fontsize V=$vel {/Symbol \261} $velsig mm/yr  WRMS=$rms mm  Chi^2/DOF=$chis_dof  TSPAN=$tspan}'";
		push(@{$TITLE{$comp}},$t);
	} else {
		print "Plotting without parameters values\n";	
		push(@{$FUNC{$comp}},"\n");
	};
	print "$f\n";
	#print " @{$FUNC{$comp}}\n";
};	

$hsize=1*$hscl;
$vsize=1*$vscl;
#print "$hsize $vsize\n";
my %PLOT=();
my $nfiles = keys %FILE;
my $k=0;
my $cols=1;
my $rows=$nfiles/$cols;
my $ys=sprintf("%3.2f",$vsize/$rows);
my $xsize=sprintf("%3.2f",$hsize/$cols-0.05);
my $ysize=sprintf("%3.2f",$ys-0.03);
my $xoff=0.01;
my $ydelta=$ys-0.015;
my $yoff=1-$ydelta;
my %XORIG = ( );
#print "$ydelta\n";
my %YORIG = ( );

#set plot margins
foreach my $comp ( @comps ) {
	push(@{$XORIG{$comp}},$xoff);
	
	push(@{$YORIG{$comp}},$yoff);
	
	$yoff-=$ydelta;
};
my $gnpath ="/opt/gnuplot/bin";
$gnpath = dirname `which gnuplot`;
chomp $gnpath;
#my $gnpath ="";
if ( $dump ) {
	$scriptfile = "commands.gnp" if ! $scriptfile;
	open(GP,"> $scriptfile") or die "couldn't open commands.gnp for writing\n";
} else {
	open(GP,"| ${gnpath}/gnuplot -persist") or die "couldn't open gnuplot for writing\n";
};

#set terminal
print GP qq{ set terminal $termstr };
if ( $term ne 'x11' ) { 
	if ( ! $outfile ) {
		$outfile="'$outdir/$site.$dtype.$term'";
	} else { $outfile="'$outfile'";};
	print "output file: $outfile\n";
	print GP qq {
		set output $outfile 
	}; 

};


# line styles for the plot
# linestyle <index> ; it is selected with ls <index> in the plot
# <index > = 1,2,3 ...
# lt (linetype) : -2  thin black; 
# colors: 
#	1 red 
#	2 green 
#	3 blue
#	4 magenta
#	5 cyan
#	6 brown
#	7 orange
#	8 pink
#	-8 small black dots
# lw (linewidth ): 0.5 times default width
# pt (pointtype) : 7 filled circle; 
# ps (pointsize) : 0.5 half default size

# "set linestyle 1 lt -2 lw 0.5" 	# for black errorbars
# "set linestyle 2 lt  1 pt 7 ps 0.5" 	# for red filled points
# "set linestyle 3 lt 3 lw 2" 		# for dashed line fitted function 
# set linestyle 4 lt 3 lw 2		# for offsets (blue)

my $timestamp='';
my $x1=0.65*$hsize;
my $x2=0.1*$hsize;
my $ytstamp=1-$vsize+0.02;
my $cwu="{/Symbol \323} {/Times-Roman-bold=9  Central Washington University }";
my $panga="{/Helvetica-bold=9 Geodesy Lab. \\\\& PANGA}";
if ( $tstamp ) {
	$LOGO=qq{
	set label \"\`date "+ \%Y-\%m-\%d; \%H:\%M; DOY: \%j"\`\" at screen $x1,$ytstamp 
	set label \"$cwu\" at screen $x2,$ytstamp
	};
}elsif (  $logo ) {	
	#my $x1=0.76*$hsize;
	#my $x2=0.1*$hsize;
	$LOGO=qq{
		set label \"$cwu\" at screen $x1,$ytstamp
		set label \"$panga\" at screen $x2,$ytstamp
	 };
};

#($ymin,$ymax,$tmin ,$tmax)=fminmax($FILE{$comp},$det,$deo);# @{$MINMAX{$comp}};

my $dy = $Ymax - $Ymin;
print "dy: $dy \n";
if ( $dy > 1  ) { 
	$ytics = 1; 
	if ( $dy  <= 5 ) { 
	   $mytics=1;
	} else {
	   $mytics=2 
	}  
} 
elsif ( $dy > 1/10 ) { 
	$ytics = 1/10;
	$mytics=4 
} else {
	$ytics=1/100; 
	$mytics= 9 ;
}


#print "$tmin $tmax\n";
my $dx= $tmax -$tmin;
if ( $dx > 1  ) { 
	$xtics = 1; 
	$xunit="Years";
	if ( $dx  <= 5 ) { 
		$mxtics=12;
		$xunit="Years-months"; 
	} else { $mxtics=2 };  
} elsif ( $dx > 1/12 ) { 
	$xtics = 1/12;
	$xunit="Months-weeks";
	$mxtics=4 
} else {
	$xtics=1/(12*4); 
	$mxtics=7 ;
	$xunit="Weeks-days";
};

	
my $xheader =0.48*$hsize;
my $yheader =1-0.015;
my $xxunit  =0.48*$hsize;
my $yxunit  =1+0.015-$vsize;
my $multiplot = " ";
$header = "{/Helvetica-bold=10 $SITE }";
if ( $showtype ) { $header="$header $DTYPE{$dtype}" };		
if ( $showtitle )  { $PTITLE = qq{  set label  \"$header\" at screen $xheader,$yheader \n} };
if ( $showfooter )  { $FOOTER = qq{ set label \"$xunit\" at screen $xxunit,$yxunit \n } };
if ( $#comps > 0 ) { $multiplot = "set multiplot" };   
# set plot frame
#set linestyle 1 lt -1 lw 0.5#
#set linestyle 2 lt 1 pt 7 ps $ps
#set linestyle 3 lt 5 lw $lw pt 7 ps 0.2
#set linestyle 4 lt 3 lw 2 pt 7 ps 0.1

print GP qq{
set dummy t
set nokey
set size 1,1
set origin 0,0 
set mytics $mytics
set mxtics $mxtics
set xtics $xtics
set ytics $ytics 
$setgrid
$PTITLE
$FOOTER
$LOGO
$multiplot
set size $xsize,$ysize
set samples 3000
set function style lines
};

foreach my $comp ( @comps ) {
	my $ylabel="$NEU{$comp} \[$yunit\]";
    my $YL=$ylabel;
	my $rfun= 0;
	my $fitfun='';
	$fitfun = "f(t)" if grep { /^f/ } @plotfun;
	my @title=();
	@title=@{$TITLE{$comp}} if defined @{$TITLE{$comp}} && $showheader;
	
    if ( defined $FUNC{$comp} ) {
    	print "Available Functions in $FILE{$comp}:\n";
	foreach my $f ( @{$FUNC{$comp}} ) {
		
		my @lhs = grep { /\(t\)/ } split(/=/,$f);
		
		next if ! defined $lhs[0] || $lhs[0] =~ /f\(t\)/ || $lhs[0] =~ /heav\(t\)/;
		print "@lhs\n";
		foreach my $rf ( @rmfun ) {
			if ( $f =~ /$rf/ )  { $rfun = "$rfun + $lhs[0] " };
		};
		
		foreach my $pf ( @plotfun ) {
			if ( $f =~ /$pf/)  { $fitfun = "$fitfun + $lhs[0] " };
		};
		
	};
	$rfun =~ s/^(0\s+\+\s+)// if $rfun =~ /0/;
	$fitfun =~ s/^(\s+\+\s+)//;
	print "Removing: $rfun\n";
	print "Plotting: $fitfun\n";
	$det=1 if $rfun =~ /line/;
	$deo=1 if $rfun =~ /(jumps|offsets)/;
#	print "@{$FUNC{$comp}}\n";
#	exit;
   };
	
	
	
	my $XRANGE="set xrange [\*:\*]\n";
	my $YRANGE="set yrange [\*:\*]\n";
    	
    #$iymax=$ymax;
    #$iymin=$ymin;
	if (  $xrange ) {	
		$XRANGE="set xrange [$tmin:$tmax]\n";
	}
    
    if ( $comp =~ /(lat|x)/  && $nrange ) {
      $YRANGE="set yrange [$nymin:$nymax]\n";
      $iymax=$nymax;
      $iymin=$nymin;
    }
    elsif ( $comp =~ /(lon|y)/  && $erange ) {
      $YRANGE=" set yrange [$eymin:$eymax]\n";
      $iymax=$eymax;
      $iymin=$eymin;
    }
    elsif ( $comp =~ /(rad|z)/  && $urange ) {
      $YRANGE="set yrange [$uymin:$uymax]\n";
      $iymax=$uymax;
      $iymin=$uymin;
    }
	elsif (  $yrange  ) {
		$YRANGE="set yrange [$ymin:$ymax]\n";
        $iymax=$ymax;
        $iymin=$ymin;
	} else {
		($ymin,$ymax,undef ,undef)= fminmax($FILE{$comp},$det,$deo,$comp);# @{$MINMAX{$comp}};
		$YRANGE="set yrange [$ymin:$ymax]\n";
        $iymax=$ymax;
        $iymin=$ymin;
	};
	
	my $dy= $iymax -$iymin;
	#print "$dy $ymax $ymin\n";	
    my $dy_exp = log($dy*1.5)/log(10);
	if ( $dy > 50 )  { $ytics = 10 };
	if ( $dy > 100 ) { $ytics = 25 };
	if ( $dy > 500 && $dy_exp > 2 ) { $ytics = 10**($dy_exp-1) };
	if ( $dy <= 50 ) { $ytics = 10 };
	if ( $dy <= 40 ) { $ytics =  5 };
	if ( $dy <= 10 ) { $ytics =  1 };
	if ( $dy <= 1  )  { $ytics =  1/10; $mytics=2; $YL="$ylabel 10^-1";  };
	if ( $dy <= 0.1 )  { $ytics= 1/100; $mytics=9; $YL="$ylabel 10^-2"; };
		
#	};	
	
	my $file = $FILE{$comp};	
	my $every="every $nsmpl";
	my $usingf=" using 1:($fitfun) ";
	my $using3=" using 1:(\$2 - ( $rfun )):3 ";
	my $using2=" using 1:(\$2 - ( $rfun )) ";
	$usingf =~ s/(\(t\))/\(\$1\)/g;
	$using2 =~ s/(\(t\))/\(\$1\)/g;
	$using3 =~ s/(\(t\))/\(\$1\)/g;
	my $linestyle="ls 1";
	
	my $plotstr='';
	
	if ( $eb ) {  # plot errorbars
		$plotstr="'$file'  $every $using3 with errorbars lt -1 lw 0.5,";
	#	$plotstr="'$file'  $every $using3 with errorbars lt 15 lw 0.5,";
	};

	#$plotstr="$plotstr '$file' $every $using2 with points lt 1 $pt ps $ps";
	#red
    $plotstr .= " '$file' $every $using2 with points lt 1 $pt ps $ps";
    # blue
	#$plotstr .= " '$file' $every $using2 with points lt 3 $pt ps $ps";

	if ( $fitfun ) { #plot fit function
		#$plotstr="$plotstr , '$file' $every $usingf with points lt 5 lw $lw $pt ps 0.2";
		$plotstr .= " , '$file' $every $usingf with points lt 5 lw $lw $pt ps 0.2";
	};
	
	my $offsetstr='';
	my @off=();
#	if (  defined @{$EVENTS{$comp}} ) { # plot vertical line at antenna jump
	if ( $event ) {
		my $evcolor = 2; #1=red, 2=green, 3=blue
		foreach my $event ( @events ) {
			$event = lc($event);
			$evcolor = 1 if $event eq "eqk";
			$evcolor = 2 if $event eq "slow";
			$evcolor = 3 if $event eq "step";
			my @imp;
			my @epochs = get_events($event,$site) ;
			next if ! $epochs[0];
			foreach my $e ( get_events($event,$site) ) {
				push(@imp,"$e\n") if  $e > $tmin && $e < $tmax ;
			};
			if ( @imp ) {
				push(@off,"@imp e \n @imp e\n");
				my $vline= ",'-' using 1:($iymin*1) with impulses lt $evcolor lw 2 ,'-' using 1:($iymax*1) with impulses lt $evcolor lw 2";
				$offsetstr .= " $vline";
			};
		};
	};


	if (  defined @{$JUMPS{$comp}} ) { # plot vertical line at antenna jump
		my @imp;
		foreach  (  @{$JUMPS{$comp}} ) {
			my ($e,$o)=split;
			push(@imp,"$e\n");
		};
		push(@off,"@imp e \n @imp e\n");
		my $vline= ",'-' using 1:($ymin*1) with impulses lt 3 lw 2 ,'-' using 1:($ymax*1) with impulses lt 3 lw 2";
		$offsetstr .= " $vline";
	};

	if (  defined  @{$OFFSETS{$comp}} ) { # plot vertical line at offset
		my @imp=();
		foreach  (  @{$OFFSETS{$comp}} ) {
			my ($e,$o)=split;
			push(@imp,"$e\n");
		};
		push(@off,"@imp e \n @imp e\n");
		my $vline= ",'-' using 1:($ymin*1) with impulses lt 2 lw 2 ,'-' using 1:($ymax*1) with impulses lt 2 lw 2";
		#my $vline= ",'-'  using 1 with impulses lt 2 lw 2 ,'-' using 1 with impulses lt 2 lw 2";
		$offsetstr .= " $vline";
	};
	$plotstr .= " $offsetstr \n @off";
#	print "$plotstr\n";
	print  GP  qq{ 
@{$FUNC{$comp}}
set mxtics $mxtics 
set mytics $mytics 
set xtics $xtics 
set ytics $ytics
set noautoscale y
#set autoscale yfix
$XRANGE 
$YRANGE
set origin @{$XORIG{$comp}},@{$YORIG{$comp}} 
set title @title 
set ylabel '$ylabel'
set dummy t 
plot $plotstr 
	};		

};
print "comps ", $#comps+1,"\n";
print GP "set nomultiplot \n " if $#comps > 0 ;
print GP "quit\n";
close(GP);


if ( $outfile && $jpg  && ! $dump ) {
    my $jpgfile = $outfile;
    $jpgfile =~ s/(\.ps)/\.jpeg/g;
   my $rotate = "";
   $rotate = "-rotate 90" if $rotjpg;   
   print "creating $jpgfile from $outfile\n";
   system("convert -size $sx,$sy $outfile -quality 100 -trim  $rotate $jpgfile");
	   
# an alternative to convert is with gs, although doesn't work that well.
#	gs -r100 -q -g900x1100 -sDEVICE=jpeg -dBATCH -dNOPAUSE -sOUTPUTFILE=$jpg $outfile

};

if ( $outfile && $png && ! $dump ) {
   my $pngfile = $outfile;
    $pngfile =~ s/(\.ps)/\.png/g;
   my $rotate = "";
   $rotate = "-rotate 90" if $rotjpg;   
   print "creating $pngfile from $outfile\n";
   system("convert -size $sx,$sy $outfile -alpha off -quality 100 -trim  $rotate $pngfile");
};
#show ps if it is required
system("gv $outfile &") if ( $term ne "x11" and $showps );

0;
#==================SUBROUTINES=============================================
sub fminmax {
	my $file=shift;
	my $det=shift || 0;
	my $deo=shift || 0;
	my $comp = shift;
	my (@T,@S);
	my ($t,$b,$s);
	my $k=0;
	my ($bmin,$bmax);
	my ($tmin,$tmax);
	my ($smin,$smax);
	my @header=getheader($file,$comp);
	my @jumps=getoffsets('OFFSET ',@header);
	my ($m,$h)=getline(@header);
	my $line=0;
	my $offs=0;
	open ( FILE , "<","$file") || warn "couldn't open file $file\n";
	while ( <FILE>) {
  		chomp;                  #no newline
  		next if /^\s*$/;        #skip blank lines
  		next if /^\#/;
  		($t,$b,$s)=split;
		$tmin = $t if $k == 0;
		$offs=0;
		$line=0;
		if ( $deo && @jumps ) {
			   foreach ( @jumps ) {
			  	my ( $e,$a ) = split;
			  	$offs=$offs+$a*heav($t - $e);
			   };
			   $line=$offs;
		};
		if ( $det ) {
			
			$line=$h + $m * $t + $line; 
			#print "$line $offs\n";
		};
		#print "$b - $line = ";
		$b=$b - $line;
		#print "$b \n";
		if ( $k == 0 ) {
			$tmin=$t;
			$bmin=$b;
			$bmax=$b;
			$smin=$s;
			$smax=$s;
			$k++;
		};
		if ( $b > $bmax ) { 
			$bmax = $b;
			$smax = $s;
		}elsif ( $b < $bmin ) { 
			$bmin = $b;
			$smin = $s;
		}; 
		
	};
	#print "minmax $bmin $bmax $smin $smax\n";
	$tmax=$t;
	my $dy=abs($bmax - $bmin);

        #print "$bmin+",sign($bmin)*3,"*$smin,$bmax+",sign($bmax)*3,"*$smax\n";
	#return($bmin-2*$smin,$bmax+2*$smax,$tmin,$tmax);
	return($bmin-(2*$smin+0.05*$dy),$bmax+(2*$smax+0.05*$dy),$tmin,$tmax);

};
sub sign{
	my $x=shift;
	my $s=0;
	if ( $x < 0) {
		$s=-1;
	} elsif ( $x > 0 ) {
		$s=1;
	};
	return ( $s );
}

sub heav {
	my $x=shift;
	return( ( $x >= 0 ? 1:0 ))
}
sub minmax {
	my @B=@_;
	my $bmin=$B[0];
	my $bmax=$B[0];
	
	foreach my $b (@B) {
	
		if ( $b > $bmax ) { 
			$bmax = $b;
		}elsif ( $b < $bmin ) { 
			$bmin = $b;
		};

	};	
	return($bmin,$bmax);
};


sub getheader{
	my $file=shift;
	my $comp = shift;
	my $flag=0;
	my $par = ( $file =~ /\.par$/ );
	my $start = 0;
	my @header;
	my @fields;
	#print "$start $par $file\n";
	open ( FILE , "<","$file") || warn "couldn't open file $file\n";
	while ( <FILE>) {
  		chomp;                  #no newline
		s/^\s+//;  #eliminate spaces at the begining of the line
  		next if /^\s*$/;        #skip blank lines
  		next if /^(##)/;
		#if ( /INPUT FILE/ && /\.$comp$/ ) { $flag=1 };
		if ( /INPUT FILE/ && ! $par ) { $flag=1 }
		elsif ( /INPUT FILE/ && /\.$comp$/ && $par ) { $flag =1 }; 
		next if ! $flag;
		last if ! /^\#/  && ! $par ;        #
		#last if $start > 1;
  		push(@header,$_);  
		#print "$_\n";
	};	
	
	return(@header);
};

sub get_par {
	my $fname = shift;
	
	open ( FILE, "< $fname") || die "Can't open $fname $! \n";
	my @lines = <FILE>;
	chomp(@lines);
	my @header = grep { "$_ =~ /^\#/" } @lines;
	my %P;

	foreach ( @header ) {
		s/^\s+//;  # no spaces at the beggining
		s/^\s*$//; # no blank lines
		s/\s+/ /g;
	    #last if ! /^\#/;
		s/^\# //;
		s/P /P/;
		#push(@header,$_);
		my @f = split;	
	
		if ( /f\(t\)=/ )	{ push(@{$P{F_T}},$f[0]); };
		if ( /OUTPUT/  )   	{ push(@{$P{DTYPE}},$f[1]) };
		if ( /MTOT/ )      	{ push(@{$P{MTOT}},$f[3]) };
		if ( /INTERVAL/ )  	{ push(@{$P{I_INTERVAL}},( $f[3],$f[5] )) };
		if ( /PREFILTER M=/ )   { push(@{$P{M}},$f[5]) };
		if ( /TIME-SPAN/) { 
				push(@{$P{F_INTERVAL}},($f[1],$f[3]));
				push(@{$P{TIME_SPAN}},$f[5]); 
		};
		if ( /USED M=/ )  { 
				push(@{$P{M_USED}},$f[3]);
				push(@{$P{M_REJECTED}},$f[6]);
				push(@{$P{M_PORC}},$f[8]); 
		};
		if ( /N=/  )  { push(@{$P{N}},$f[1]) };
		if ( /CHI/ )  { 
				push(@{$P{CHI2}},$f[1]);
				push(@{$P{DOF}},$f[3]);
				push(@{$P{CHI_DOF}},$f[5]);   		 
		};
		if ( /RMS/ ) {
				push(@{$P{RMS}},$f[1]);
				push(@{$P{SDEV}},$f[3]);
				push(@{$P{MEAN}},$f[5]);   				
		};
		if ( /^P\d+/ ) {
	    		if ( /OFFSET/ ) {
	    	    #(undef,$p,undef,$psig,undef,$name,$epoch)=split;	
		     		push(@{$P{$f[5]}},($f[1],$f[3],$f[6]));
	    		} else {
 	    	    #(undef,$p,undef,$psig,undef,$name)=split ;
				 push(@{$P{$f[5]}},($f[1],$f[3]));
		  		# push(@{$P{$name}},($p,$psig));
	    		};		
		};
	};
	return(%P);
}

sub getline {
	my @input=@_;
	my $vel=0;
	my $yzero=0;
	foreach (@input) { 
		chomp;
		s/^\s+//;  #eliminate spaces at the begining of the line
		s/^\#+//; 
		#next if ! /^P/ ; #next line if it is not a parameter line
		s/(P\s+)/P/g; 	# replace space between P's and numbers
		s/\s+/ /g ; #replace multiple white spaces with only one
		#print "$_\n";
		(undef,$vel,undef,undef,undef)=split if  / SLOPE / ;
		(undef,$yzero,undef,undef,undef)=split if  /\bYZERO\b/ ;
		
	};
	#print "$vel $yzero\n";
	return($vel,$yzero);
};

sub getpar { 
	my @input=@_;
	my $output='';
	my $vel=0;
	my $velsig=0;
	my $yzero=0;
	my $yzerosig=0;
	my $chis=0;
	my $chis_dof=0;
	my $sdev=0;
	my $dof=0;
	my $rms=0;
	my $mean=0;
	my $tspan=0;
	foreach (@input) { 
		chomp;
		s/^\s+//;  #eliminate spaces at the begining of the line
		s/^\#+//; 
		#next if ! /^P/ ; #next line if it is not a parameter line
		s/(P\s+)/P/g; 	# replace space between P's and numbers
		s/\s+/ /g ; #replace multiple white spaces with only one
		#print "$_\n";
		(undef,undef,$tspan)=split(/=/,$_) if  /TIME-SPAN=/ ; 
		(undef,undef,$vel,undef,$velsig,undef)=split if  /\sSLOPE\s/ ;
		(undef,undef,$yzero,undef,$yzerosig,undef)=split if  /\bYZERO\b/ ;
		(undef,$chis,undef, $dof,undef, $chis_dof)=split if  /CHI/ ; 
		(undef,$rms,undef, $sdev,undef,$mean)=split if  /RMS/g ; 
	};
	$output=sprintf("%5.2f,%4.2f,%4.2f,%5.2f,%d,%6.3f",$vel,$velsig,$rms,$chis_dof,$dof,$tspan);
	return($output);

};
sub getoffsets {
	my $key = shift; 
	my @input=@_;
	my @output=();
	my @offsets=();
	my $offset=0;
	my $offsetsig=0;
	my $epoch=0;
	foreach (@input) { 
		chomp;
		s/^\s+//;  #eliminate spaces at the begining of the line
		next if /\#\#/; 
		#next if ! /^P/ ; #next line if it is not a parameter line
		s/^\#//;
		s/(P\s+)/P/g;
#		s/^P//;  #erase the P
		
		s/\s+/ /g ; #replace multiple white spaces with only one
		if ( /$key/ ) {
			(undef,undef,$offset,undef,$offsetsig,undef,undef,$epoch)=split;
			push(@offsets,"$epoch $offset");
		}; 
	};
	push(@output,@offsets);
	@output=sort(@output);
	#print ( map {"$_\n"} @output );
	return(@output);
};
sub get_events {
	my $key = shift || 'step';
	my $site = shift; 
	my @output=();
    my $here = cwd;
	my $dir = -d "$here/step" ? "$here/step" : "/panga/officialresults/$key/all_$key/";
	#$dir = "/home/marcelo/processing/post_proc/tswork/panga_work/mmslow/" if $key eq 'slow';
	#$dir = "$ENV{GPS_HOME_DIR}/tswork/official/$key/";
	if ( -e "$dir/$site.$key" ) {
		open ( FE , "< $dir/$site.$key");
		my @input = <FE>;
		close (FE);
		my $epoch;
		foreach (@input) { 
			chomp;
			s/^\s+//;  #eliminate spaces at the begining of the line
			next if /^\#/;
			($epoch,undef)=split;
			push(@output,$epoch);
		};
	} else { push ( @output , 0 ) };
	#print ( map {"events: $key => $_\n"} @output );
	return(@output);
};

sub getfunc {
	my @input=@_;
	#print "@fin\n";
	my @fout=();
	my @funcs=();
	foreach ( @input ){
		chomp;
		s/\s+//g;      
		next if /^\s*$/;        #skip blank lines
		next if /\#\#/;
		#last if ! /^\#/;        #use with the file
		s/^\#//;
		s/(P\s+)/P/g; 	# replace space between P's and numbers
		#get functions
		if ( /\(t\)=/ ) {
			s/(=\+)/=/;	# remove '+' after '='
	#		s/\(t/\(x/g;	#replace t by x
	#		s/\*t/\*x/g;	#replace t by x
			push(@funcs, $_); 
		};
		#print "$_\n";
		#get parameters
		if ( /(P\d+\=)/ ) {
			my @fields=split(/\+\-/,$_ );
			#print FOUT "$fields[0]\n";
			push(@fout, "$fields[0]\n");
		};
	#	push(@funcs, $_) if /\(t\)=/
	};
	foreach my $f ( reverse @funcs ) {
#		print FOUT "$f\n";
		push(@fout,"$f\n");
	};
	return( @fout );

}
