#!/usr/bin/perl

use warnings;
use strict;

use Data::Dumper;
use Getopt::Long;
use PDL;
use PDL::Graphics::PGPLOT;
use PGPLOT;
#use Tie::IxHash;

# The RDB module is tested for and 'use'-ed in parse_args.

our ( $VERSION ) = '$Revision: 1.4 $' =~ /([\d.]+)/;

eval {
    main ();
};
if ( my $err = $@ ) {
    print STDERR "# $0: $_\n", foreach split /\n/, $err;
    exit 1 ;
}
exit 0;

sub main {
    my $r_pars = parse_args ();

    my $r_leg;
    my $r_points = {};
    #tie %$r_points, "Tie::IxHash";

    ( $r_points, $r_leg ) = read_points ( $r_pars, $r_points );
    draw_points ( $r_pars, $r_points, $r_leg );
}

sub read_points {
    my $r_pars   = shift () or die "no pars in read_points";
    my $r_points = shift () or die "no points in read_points";

    my @leg_arr = ( !defined $r_pars->{split} )
                      ? read_RDB  ( $r_pars, $r_points )
                      : read_file ( $r_pars, $r_points );

    print "read-in points: ", Dumper $r_points
        if $r_pars->{verbose} > 1;

    print "Legend array: ", Dumper [ map { $_->[1] } @leg_arr ]
        if $r_pars->{verbose} > 1;

    $r_points = loggify_data ( $r_points, $r_pars )
        if $r_pars->{xlog} || $r_pars->{ylog};

    return $r_points, [ map { $_->[1] } @leg_arr ];
}

sub read_RDB {
    my $r_pars   = shift () or die "no r_pars   in read_RDB";
    my $r_points = shift () or die "no r_points in read_RDB";

    print "reading with RDB.pm\n"
        if $r_pars->{verbose};

    my $rdb = new RDB;
    $rdb->open (
              $r_pars->{filename} eq 'stdin' 
                  ? \*STDIN 
                  : $r_pars->{filename} 
          ) or die "RDB open failed on file '$r_pars->{filename}' $!";

    my $r_line  = {};
    my @leg_arr = ();

    while ( $rdb->read( $r_line ) ) {

        my $brk = ( defined $r_pars->{break} )
                      ? $r_line->{$r_pars->{break}}
                      : $r_pars->{only};

        push @{$r_points->{$brk}{x}}, $r_line->{$r_pars->{x}};

        foreach ( 0 .. scalar @{$r_pars->{y}} - 1 ) {

            print "read: $_ $r_line->{$r_pars->{y}[$_]} $r_pars->{y}[$_]\n"
                if $r_pars->{verbose} > 3;

            push @{$r_points->{$brk}{y}[   $_]}, 
                 $r_line->{$r_pars->{y}[   $_]};
            push @{$r_points->{$brk}{yerr}[$_]}, 
                 $r_line->{$r_pars->{yerr}[$_]}
                if defined $r_pars->{yerr}[$_];

            my $leg_key = $r_pars->{y}[$_] . 
                          ( $brk ne $r_pars->{only} 
                                ? ", $brk:" 
                                : ""
                          );
            push @leg_arr, [ $brk, $leg_key ]
                unless grep { $_->[1] eq $leg_key } @leg_arr;
        }
    }
    return @leg_arr;
}

sub read_file {

    my $r_pars   = shift () or die "no r_pars   in read_RDB";
    my $r_points = shift () or die "no r_points in read_RDB";

    print "reading from regular file\n"
        if $r_pars->{verbose};

    # If a filename is given, read from it; otherwise, read from STDIN:
    #
    my $fh;
    if ( $r_pars->{filename} ne 'stdin' ) {
        open $fh, $r_pars->{filename}
    }
    else {
        $fh = 'STDIN';
    }

    my ( @col_names, %col_number, @leg_arr ) = ( (), (), () );
    while ( <$fh> ) {
        next if /^#/;          # skip comments
        next if /^[#SNM\t]+$/; # skip definition lines if you're 
                               #   kludge-reading a RDB file

        chomp ( my @vals = split /$r_pars->{split}/, $_ );

        die "Split returned only one item per line-- abort!\n"
            if 1 == scalar @vals;

        if ( 0 == scalar @col_names ) {
            @col_names = @vals;
            %col_number = map {$col_names[$_] => $_} 0..scalar @col_names-1;
            next;
        }
        my $brk =  ( defined $r_pars->{break} )
                      ? $vals[ $col_number{ $r_pars->{break} } ]
                      : $r_pars->{only};
        push @{$r_points->{$brk}{x}}, $vals[ $col_number{ $r_pars->{x} } ];

        foreach ( 0 .. scalar @{$r_pars->{y}} - 1 ) {
            push @{$r_points->{$brk}{y}[$_]}, 
                 $vals[ $col_number{$r_pars->{y}[$_]} ];

            push @{$r_points->{$brk}{yerr}[$_]}, 
                 $vals[ $col_number{$r_pars->{yerr}[$_]} ]
                if defined $r_pars->{yerr}[$_];

            my $leg_key = $r_pars->{y}[$_] . 
                          ( $brk ne $r_pars->{only} 
                                ? ", $brk:" 
                                : ""
                          );
            push @leg_arr, [ $brk, $leg_key ]
                unless grep { $_->[1] eq $leg_key } @leg_arr;
        }
    }
    close $fh
        if $r_pars->{filename} ne 'stdin';

    return @leg_arr;
}

sub draw_points {
    my ( $r_pars, $r_points, $r_leg ) = @_;

    my ( $r_x_data, $r_y_data, $r_y_errs, $r_res_data, $r_res_errs ) =
        make_pdls ( $r_points, $r_pars );

    my $win = setup_win ( $r_pars, $r_x_data, $r_y_data, $r_y_errs );

    my @lims1 = get_limits ( $r_pars, $r_x_data, $r_y_data, $r_y_errs, 0 );

    set_window ( $win, $r_pars, \@lims1, 0 );

    my $opts = 
        points_draw_loop ( $r_pars, $win, $r_x_data, $r_y_data, $r_y_errs );

    my $font_charsize_opts = {
        Font      => $r_pars->{font},
        HardFont  => $r_pars->{font},
        CharSize  => $r_pars->{char_size},
        HardCH    => $r_pars->{char_size},
    };
    $win->label_axes ( 
        ( defined $r_pars->{residuals} )
            ? ( "",      @{$r_pars}{'ylabel','title'}, $font_charsize_opts )
            : ( @{$r_pars}{'xlabel','ylabel','title'}, $font_charsize_opts )
    );
    pgmtxt ( 'T', 0.5, 0.5, 0.5, $r_pars->{subtitle} );

    $r_leg = $r_pars->{legend_text}
        if defined $r_pars->{legend_text};
    write_legend ( $r_pars, $win, \@lims1, $opts, $r_leg )
        if $r_pars->{legend};
    

    # Return if the residuals do not exist:
    #
    return 
        if not defined $r_pars->{residuals};


    # Draw the residuals:
    #
    my @lims2 = get_limits ( $r_pars, $r_x_data, $r_res_data, $r_res_errs, 1 );
    set_window ( $win, $r_pars, \@lims2, 1 );

    $opts = points_draw_loop ( $r_pars,$win,$r_x_data,$r_res_data,$r_res_errs, 1 );
    $win->label_axes (
        $r_pars->{xlabel},
        $r_pars->{delta_label}, 
        $font_charsize_opts
    );
}

sub points_draw_loop {
    my $r_pars       = shift () or die "no pars in points_draw_loop";
    my $win          = shift () or die "no win in points_draw_loop";
    my $r_x_loc      = shift () or die "no x data in points_draw_loop";
    my $r_y_loc      = shift () or die "no y data in points_draw_loop"; 
    my $r_e_loc      = shift () || undef;
    my $bRes         = shift () || 0; 

    my ( $opts, $opt_num )  = ( {}, 0 );
    my @syms   = split /,/, $r_pars->{symbols};
    my @colors = split /,/, $r_pars->{colors}; 

    foreach my $brk ( sort keys %$r_x_loc ) {
        foreach ( 0 .. scalar @{$r_y_loc->{$brk}} - 1 ) {

            $opts->{pt_col}[$opt_num] = $colors[ $opt_num % scalar @colors];
            $opts->{symbol}[$opt_num] = $syms[   $opt_num % scalar @syms  ];
            $opts->{ln_col}[$opt_num] = $colors[ $opt_num % scalar @colors];
            $opts->{ln_sty}[$opt_num] = 1 + $opt_num % 5;

            my $pointsopt = { color     => $opts->{pt_col}[$opt_num], 
                              symbol    => $opts->{symbol}[$opt_num],
                              linewidth => $r_pars->{line_width} };
            my $lineopt   = { color     => $opts->{ln_col}[$opt_num],
                              linestyle => $opts->{ln_sty}[$opt_num],
                              linewidth => $r_pars->{line_width} };

            my @data = ( $r_x_loc->{$brk}, $r_y_loc->{$brk}[$_]  );

            $win->line   ( @data, $lineopt   ) unless $r_pars->{noline};
            $win->points ( @data, $pointsopt ) unless $r_pars->{nopoints};
            $win->bin    ( hist $data[1] )     if     $r_pars->{histogram};

            my ( $err_lo, $err_hi );
            if ( defined $r_e_loc ) {
                $err_lo = $r_e_loc->{$brk}[$_];
                $err_hi = $r_e_loc->{$brk}[$_];

                # Data and err is already logged at this point:
                #
                if ( $r_pars->{ylog} ) {
                    my $log_kludge_lo = log10 ( 10**($r_y_loc->{$brk}[$_]) -
                                                10**($r_e_loc->{$brk}[$_]) );
                    my $log_kludge_hi = log10 ( 10**($r_y_loc->{$brk}[$_]) +
                                                10**($r_e_loc->{$brk}[$_]) );
                    $err_lo = $r_y_loc->{$brk}[$_] - $log_kludge_lo;
                    $err_hi = $log_kludge_hi - $r_y_loc->{$brk}[$_];
                }

                $win->errb (
                          @data, 
                          undef, undef, 
                          $err_lo, $err_hi,
                          $lineopt
                );
            }

            ++$opt_num;
        }
        last if $bRes;
    }
    return $opts;
}

sub make_pdls {
    my ( $r_points, $r_pars ) = @_;

    my $r_x_data;
    my $r_y_data;
    my $r_y_errs   if $r_pars->{yerr};
    my $r_res_data if $r_pars->{residuals};
    my $r_res_errs if $r_pars->{residuals} && $r_pars->{yerr};

    # Do x and y:
    #
    foreach my $brk ( sort keys %$r_points ) {

        $r_x_data->{$brk} = pdl @{ $r_points->{$brk}{x} };

        foreach ( 0 .. scalar @{$r_points->{$brk}{y}} - 1 ) {
            push @{$r_y_data->{$brk}}, pdl @{ $r_points->{$brk}{y}[   $_] };
            push @{$r_y_errs->{$brk}}, pdl @{ $r_points->{$brk}{yerr}[$_] } 
                if defined $r_pars->{yerr}[$_];
        }

    }
    my @ret = ( $r_x_data, $r_y_data );
    $ret[2] = ( $r_pars->{yerr}      ) ? $r_y_errs   : undef;
    return @ret
        unless $r_pars->{residuals};


    # Do residuals:
    #
    # Two column version-- just subract:
    #
    if ( scalar @{$r_pars->{y}} > 1 ) {
        foreach my $brk ( sort keys %$r_points ) {
            $r_res_data->{$brk} = [$r_y_data->{$brk}[0] - $r_y_data->{$brk}[1]];
            $r_res_errs->{$brk} = undef
                if $r_pars->{yerr};

            $r_res_errs->{$brk} = 
                [ sqrt ( $r_y_errs->{$brk}[0]**2 + $r_y_errs->{$brk}[1]**2 ) ]
                    if ( $r_pars->{residuals} && defined $r_pars->{yerr}[0]
                                              && defined $r_pars->{yerr}[1] )
        }
    }
    # Break column version-- need to interpolate set 1 to get to set 2:
    #
    elsif ( defined $r_pars->{break} )  {
        my @brk = ( sort keys %$r_points )[0..1]; 

        foreach ( 0 .. scalar @{$r_y_data->{$brk[0]}} - 1 ) {
            $r_res_data->{$brk[0]}[$_]  = 
                ( interpolate ( 
                      $r_x_data->{$brk[1]}, 
                      $r_x_data->{$brk[0]}, 
                      $r_y_data->{$brk[0]}[$_]
                  ) 
                )[0] - $r_y_data->{$brk[1]}[$_];
            $r_res_errs->{$brk[0]}[$_] = $r_y_errs->{$brk[0]}[$_];
        }
        $r_res_data->{$brk[1]} = $r_res_data->{$brk[0]};
    }
    else { 
        die "WTF in make_pdls.  trouble with making residuals";
    }


    $ret[3] = ( $r_pars->{residuals} ) ? $r_res_data : undef;
    $ret[4] = ( $r_pars->{residuals} && $r_pars->{yerr}) ? $r_res_errs : undef;

    return @ret;
}

sub setup_win {
    my $r_pars = shift () or die "no pars in setup_win";

    my $win = PDL::Graphics::PGPLOT::Window->new ( 
        Device     => $r_pars->{device}, 
        WindowName => "plot created by $0",
        AxisColor  => 'black',
        Color      => 'black',
        Axis       => $r_pars->{axis} ? 'axes' : 'normal',
        Font       => $r_pars->{font},
        HardFont   => $r_pars->{font},
        CharSize   => $r_pars->{char_size},
        HardCH     => $r_pars->{char_size},
    );

    return $win
}

sub set_window {
    my $win    = shift () or die "no window object in set_window";
    my $r_pars = shift () or die "no pars in set_window";
    my $r_lims = shift () or die "no limits in set_window";
    my $iWin   = shift () || 0;

    my @env_pars = (
        @$r_lims,
        {
            PlotPosition => $r_pars->{PlotPosition}[$iWin],
            Axis         => [ 'BCNST', 'BCMSTV' ],
        }
    );

    #$env_pars[-1]{Axis} = $r_pars->{WindowSetPars}{Axis}
    #    if defined $r_pars->{WindowSetPars}{Axis};

    $env_pars[-1]{Axis}[0] .= 'L' if $r_pars->{xlog};
    $env_pars[-1]{Axis}[1] .= 'L' if $r_pars->{ylog};


    if ( 0 == $iWin && $r_pars->{residuals} ) {
        $env_pars[-1]{Axis} = [ 
            "BCST"   . ( $r_pars->{xlog} ? 'L' : '' ),
            "BCSTNV" . ( $r_pars->{ylog} ? 'L' : '' ),
        ];
    }
    $env_pars[-1]{Axis}[1] =~ s/M/N/
        if !$r_pars->{delta_pos};

    $win->env ( @env_pars );
}

sub write_legend {
    my $r_pars = shift () or die "no pars in write_legend";
    my $win    = shift () or die "no win  in write_legend";
    my $lims   = shift () or die "no lims in write_legend";
    my $opts   = shift () or die "no opts in write_legend";
    my $r_leg  = shift () or die "no leg  in write_legend";

    my @loc = ( $lims->[0], $lims->[3] );

    my @deltas = ( $lims->[1] - $lims->[0], $lims->[3] - $lims->[2] );

    if ( $r_pars->{legend_location} ) {
        $loc[0] += $r_pars->{legend_location}[0] * $deltas[0];
        $loc[1] += $r_pars->{legend_location}[1] * $deltas[1];
    }
    else {
        $loc[0] +=  .1 * $deltas[0];
        $loc[1] += -.1 * $deltas[1];
    }

    # Legend Usage:
    #
    # [ names ],
    # x,y
    # { option hash }
    #
    $win->legend (
        $r_leg,
        #0.01, .7,
        @loc,
        {
            Symbol    => $opts->{symbol},
            LineStyle => $opts->{ln_sty},
            Color     => $opts->{ln_col},
            LineWidth => [ 50, 50 ],
            TextShift => 0,
            Font      => $r_pars->{font},
            HardFont  => $r_pars->{font},
            CharSize  => $r_pars->{char_size},
            HardCH    => $r_pars->{char_size},
            Fraction  => 0.5,
        }
    );
}

sub log_ten {
    my $num = shift () || 0.0;
    return ( $num <= 0.0 )
               ? undef
               : log ( $num ) / log ( 10.0 ); 
}

sub null_undef_points {
    my $r_pars   = shift () or die "no pars   in null_undef_points";
    my $r_points = shift () or die "no points in null_undef_points";

    # Determine which points are undefined:
    #
    my %null = ();
    foreach my $brk ( sort keys %$r_points ) {
        foreach ( 0 .. scalar @{$r_points->{$brk}{x}} - 1 ) {

            ++$null{$_}
                if not defined $r_points->{$brk}{x}[$_];
            
            foreach my $col ( 0 .. scalar @{$r_points->{$brk}{y}} - 1 ) {
                ++$null{$_} 
                    if  not defined $r_points->{$brk}{y}[$col][$_] ||
                       (not defined $r_points->{$brk}{yerr}[$col][$_] &&
                        defined $r_pars->{yerr}[$_]);
            }
        }
    }

    # And skip them in all columns:
    #
    my $ret_points;
    foreach my $brk ( sort keys %$r_points ) {
        foreach ( 0 .. scalar @{$r_points->{$brk}{x}} - 1 ) {
            next if exists $null{$_};

            push @{$ret_points->{$brk}{x}}, $r_points->{$brk}{x}[$_];

            foreach my $col ( 0 .. scalar @{$r_points->{$brk}{y}} - 1 ) {
                push @{$ret_points->{$brk}{y}[$col]},    
                    $r_points->{$brk}{y}[$col][$_];
                push @{$ret_points->{$brk}{yerr}[$col]}, 
                    $r_points->{$brk}{yerr}[$col][$_];
            }
        }
    }
    return $ret_points;
}

sub loggify_data {
    my ( $r_points, $r_pars ) = @_;

    # Deal with a logged x- or y-axis (or both).  Remove points
    #   that are <= 0.0:
    #
    my %log_axis = ();
    foreach my $brk ( sort keys %$r_points ) {
        foreach ( 0 .. scalar @{$r_points->{$brk}{x}} - 1 ) {

            if ( $r_pars->{xlog} ) {
                $log_axis{x} = 1;

                $r_points->{$brk}{x}[$_] =
                    log_ten ( $r_points->{$brk}{x}[$_] );
            }
            if ( $r_pars->{ylog} ) {
                $log_axis{y} = 1;

                $r_points->{$brk}{residuals}[$_] =
                    log_ten ( $r_points->{$brk}{x}[$_] );

                foreach my $col ( 0 .. scalar @{$r_points->{$brk}{y}} - 1 ) {

                    $r_points->{$brk}{y}[$col][$_] =
                        log_ten ( $r_points->{$brk}{y}[$col][$_] );

                    $r_points->{$brk}{yerr}[$col][$_] =
                        log_ten ( $r_points->{$brk}{yerr}[$col][$_] )
                            if defined $r_pars->{yerr}[$col];
                }
            }

        }
    }
    $r_points = null_undef_points ( $r_pars, $r_points );
    $r_pars->{WindowSetPars}{Axis} = 'log' . join ( "", sort keys %log_axis );

    return $r_points;
}

sub get_limits {
    my $r_pars   = shift () or die "no r_pars in get_limits";
    my $r_x_data = shift () or die "no x data in get_limits";
    my $r_y_data = shift () or die "no y data in get_limits";
    my $r_y_errs = shift () || undef;
    my $bIsResids = shift () || 0;

    # Limits:
    #            x_low  x_high, y_low, y_high
    my @lims = ( undef, undef,  undef, undef );

    # Get data extremes, and make them the limits:
    #
    my @brks = sort keys %$r_x_data;
    @brks = @brks[0,1] if $bIsResids && $r_pars->{break};

    foreach my $brk ( @brks ) {
        my ( $min, $max ) = ( $r_x_data->{$brk}->min, 
                              $r_x_data->{$brk}->max );
        $lims[0] = $min if ! defined $lims[0] || $min < $lims[0];
        $lims[1] = $max if ! defined $lims[1] || $max > $lims[1];

        foreach ( 0 .. scalar @{$r_y_data->{$brk}} - 1 ) {

            my ( $min, $max );
            if ( $r_pars->{ylog} ) {
                if ( defined $r_y_errs->{$brk}[$_] ) {
                    # Ignore negative err bars-- they can be huge:
                    $min = log10 ( ( 10**( $r_y_data->{$brk}[$_] )  )->min () );
                    $max = log10 ( ( 10**( $r_y_data->{$brk}[$_] ) +
                                     10**( $r_y_errs->{$brk}[$_] )  )->max () );
                }
                else {
                    $min = log10 ( ( 10**( $r_y_data->{$brk}[$_] )  )->min () );
                    $max = log10 ( ( 10**( $r_y_data->{$brk}[$_] )  )->max () );
                }
            }
            else {
                if ( defined $r_y_errs->{$brk}[$_] ) {
                    $min = ( $r_y_data->{$brk}[$_] - $r_y_errs->{$brk}[$_] )->min();
                    $max = ( $r_y_data->{$brk}[$_] + $r_y_errs->{$brk}[$_] )->max();
                }
                else {
                    $min = ( $r_y_data->{$brk}[$_] )->min();
                    $max = ( $r_y_data->{$brk}[$_] )->max();
                }
            }
            $lims[2] = $min if ! defined $lims[2] || $min < $lims[2]; 
            $lims[3] = $max if ! defined $lims[3] || $max > $lims[3];
        }
    }

    return pad_limits ( $r_pars, @lims );
}

sub pad_limits {
    my $r_pars = shift () or die "no pars in pad_limits";
    my @lims = @_;

    my $dx = $lims[1] - $lims[0];
    my $dy = $lims[3] - $lims[2];

    # Pad limits to create pleasing margins:
    #
    $lims[0] -= .1 * $dx;
    $lims[1] += .1 * $dx;
    $lims[2] -= .1 * $dy;
    $lims[3] += .1 * $dy;

    # Override if user-specified limits exist:
    #
    @lims[0,1] = split /,/, $r_pars->{xrange} if defined $r_pars->{xrange};
    @lims[2,3] = split /,/, $r_pars->{yrange} if defined $r_pars->{yrange};

    # Give a 2-unit range is either axis's range is zero:
    #
    my $epsilon = 1.0e-9;
    if ( not defined $r_pars->{xrange} && $lims[0] - $lims[1] < $epsilon ) {
        $lims[0] -= .1;
        $lims[1] += .1;
    }
    if ( not defined $r_pars->{yrange} && $lims[2] - $lims[3] < $epsilon ) {
        $lims[2] -= .1;
        $lims[3] += .1;
    }

    return @lims;
}

sub parse_args {

    # Do the getopt:
    #
    my @opts = (
        { name => 'title',           type => '=s', dval => '',         },
        { name => 'subtitle',        type => '=s', dval => '',         },
        { name => 'xlabel',          type => '=s', dval => '',         },
        { name => 'ylabel',          type => '=s', dval => '',         },
        { name => 'delta_label',     type => '=s', dval => '',         },
        { name => 'delta_pos',       type =>   '', dval =>  0,         },
        { name => 'legend',          type => '!',  dval => 1,          },
        { name => 'legend_location', type => '=s', dval => '.02,-.05', },
        { name => 'legend_text',     type => '=s', dval => undef,      },
        { name => 'x',               type => '=s', dval => undef,      },
        { name => 'y',               type => '=s', dval => undef,      },
        { name => 'yerr',            type => '=s', dval => undef,      },
        { name => 'break',           type => '=s', dval => undef,      },
        { name => 'residuals',       type => '',   dval => undef,      },

        { name => 'xrange',          type => '=s', dval => undef,      },
        { name => 'yrange',          type => '=s', dval => undef,      },
        { name => 'xlog',            type => '',   dval => 0,          },
        { name => 'ylog',            type => '',   dval => 0,          },
        { name => 'colors',          type => '=s', dval => 
                                'black,red,green,blue,cyan,magenta,gray' },
        { name => 'symbols',         type => '=s', dval => '999'       },

        { name => 'font',            type => '=f', dval => '1'         },
        { name => 'char_size',       type => '=f', dval => '1'         },
        { name => 'line_width',      type => '=f', dval => '2'         },


        { name => 'nopoints',        type => '',   dval => 0,          },
        { name => 'noline',          type => '',   dval => 0,          },
        { name => 'histogram',       type => '',   dval => 0,          },
                { name => 'axis',            type => '!',  dval => 0,          },

        { name => 'f',               type => '=s', dval => 'stdin',    },
        { name => 'filename',        type => '=s', dval => undef,      },

        { name => 's',               type => '=s', dval => undef,      },
        { name => 'split',           type => '=s', dval => undef,      },

        { name => 'device',          type => '=s', dval => '/xs',      },
        { name => 'only',            type => '=s', dval => 'only',     },
        { name => 'defaults',        type => '',   dval => 0,          },
        { name => 'verbose',         type => '=f', dval => 0,          },
        { name => 'help',            type => '',   dval => 0,          },
        { name => 'usage',           type => '',   dval => 0,          },
        { name => 'version',         type => '',   dval => 0,          },
    );
    my %pars = map { $opts[$_]{name} => $opts[$_]{dval} } 0..scalar @opts - 1;
    my @args = map { $opts[$_]{name} .  $opts[$_]{type} } 0..scalar @opts - 1;
    GetOptions ( \%pars, @args ) or die "bad GetOptions $!";

    die "$VERSION" if $pars{version};

    $pars{filename} = $pars{f} if !defined $pars{filename};
    $pars{split}    = $pars{s} if !defined $pars{split};

    $pars{symbols} = join ",", 3, 0, 5, 4, 6..99
        if $pars{symbols} eq '999';

    # Die if defaults are requested
    #
    die Data::Dumper->Dump( [ \%pars ], [ qw( pars ) ] ),
       "Printing defaults and exiting on user request"
            if ( $pars{defaults} );

    help ( 1 ) if $pars{help};
    help ( 2 ) if $pars{usage};

    # Determine if user wants to use RDB.pm. 
    #   Include it if it exists, otherwise,
    #   warn user, and try to soldier on.
    #
    if ( not defined $pars{split} ) {
        eval { require RDB; };
        if ( my $err = $@ ) {
            $pars{split} = "\t";
            printf STDERR "RDB.pm not found.  Setting --split='\t' and ",
                          "continuing with crossed fingers.  Here goes ",
                          "nothing.\n";
        }
    }

    # Massage arguments:
    #
    $pars{x}    = shift @ARGV unless $pars{x};
    $pars{y}    = shift @ARGV unless $pars{y};
    $pars{yerr} = shift @ARGV unless $pars{yerr};

    $pars{y}         = [ split /,/, $pars{y}         ] if defined $pars{y};
    $pars{yerr}      = [ split /,/, $pars{yerr}      ] if defined $pars{yerr};

    $pars{legend_location} = [ split /,/, $pars{legend_location} ] 
        if defined $pars{legend_location};
    $pars{legend_text} = [ split /,/, $pars{legend_text} ] 
        if defined $pars{legend_text};

    #die "Number of y data (" . scalar @{$pars{y}}    . 
    #    ") and error-bar ("  . scalar @{$pars{yerr}} . ") sets must be equal"
    #    if defined $pars{yerr} && scalar @{$pars{y}} != @{$pars{yerr}};

    die "Must have two y data columns or have a break column to use --residuals"
        if $pars{residuals} && scalar @{$pars{y}} < 2 && not defined $pars{break};

    deduce_x_y_names ( \%pars )
        if ! defined $pars{x} && ! defined $pars{y};

    $pars{PlotPosition} = 
        defined $pars{residuals}
            ? ( [ [ 0.1, 0.9, 0.25, 0.90 ], [ 0.1, 0.9, 0.1, 0.25 ] ] ) 
            : ( [ [ 0.1, 0.9, 0.1,  0.90 ] ] );

    $pars{xlabel}      = "$pars{x}"    if '' eq $pars{xlabel};
    $pars{ylabel}      = "@{$pars{y}}" if '' eq $pars{ylabel};
    $pars{delta_label} = "deltas"      if '' eq $pars{delta_label};

    my $tit_fmt   = $pars{break} ? "%s vs %s    grouped by %s" : "%s vs %s";
    my @print_arr = ( "@{$pars{y}}", "$pars{x}" );
    push @print_arr, $pars{break}
        if $pars{break};
    $pars{title} = sprintf $tit_fmt, @print_arr
        if '' eq $pars{title};

    print Dumper \%pars
        if $pars{verbose} > 3;

    return \%pars;
}

sub deduce_x_y_names {

    my $r_pars = shift () or die "no file in x_y_names";

    my $tmpfile = '.plot.hdr';

    if ( 'stdin' eq $r_pars->{filename} ) {
        open my $write_fh, '>', $tmpfile;
        print $write_fh $_ while <>;
        close $write_fh;

        $r_pars->{filename} = $tmpfile;
    }

    open my $fh, $r_pars->{filename};
    my $definition_line = '#';

    chomp ( $definition_line = <$fh> ) 
        while $definition_line =~ /^\s*#/;
    close $fh;

    # Non-RDB file case:
    #
    if ( defined $r_pars->{split} ) {
        $r_pars->{x} =   ( split /$r_pars->{split}/, $definition_line )[0];
        $r_pars->{y} = [ ( split /$r_pars->{split}/, $definition_line )[1] ];
    }
    # RDB file case:
    #
    else {
        $r_pars->{x} =   ( split /\t/, $definition_line )[0];
        $r_pars->{y} = [ ( split /\t/, $definition_line )[1] ];
    }
}


sub help {
    my ( $verbose ) = @_;
    require IO::Page; 
    require Pod::Usage;
    Pod::Usage::pod2usage ( { -exitval => 0, -verbose => $verbose } ); 
}

=pod

=head1 NAME

pdlplot- A (hopefully) easy to use plotter with a (relatively) simple interface.

=head1 SYNOPSIS

B<pdlplot> -f[ile] <datafile name> [-split '\t']
[-x] <x column name> [-y] <y column name> [ [-yerr] <y error bar column name>]
[I<options>]

=head1 DESCRIPTION

B<pdlplot> is a simple PDL-based plotter.  2-D plots can be easily 
made from a file containing columns of data.  The columns can be delimited
with any character(s) or regular expression.

=head1 OPTIONS

Options are specified using a getopt style interface.  Long names are available
when preceded with double hyphens C<-->, in which case only the minimal number
of unique characters are required.

=over 8

=item B<--file>=I<input data file>

The name of the file that contains the data.  If not specified, B<pdlplot> will
read from STDIN.  If B<--split> (see below) is not specified, the file is
assumed to be an RDB file (see
http://hea-www.harvard.edu/MST/simul/software/docs/rdb.html).  B<pdlplot> makes
an attempt to read the file if the RDB.pm module is not on the local system--
RDB comments (leading '#'s) are stripped, the column definition line is
ignored, and the body of the file is split into columns on tabs.

=item B<--split>=I<data separator for non-RDB input files>

Used for non-RDB files,  B<--split> specifies which character(s) (or regex) to
split the data from B<--file> on.  For a comma-delimitted file, --split ','
would be the correct usage, or --split '\|' for a pipe-delimitted file (the
pipe is a special char in perl regexes, so we must escape it).  If the file
specified by B<--file> is an RDB file, this switch should not be used.

=item B<--x>=I<X-axis column name> 

The name of the column in the rdb file that holds the independent variable.
This can also be given as the first unflagged argument, i.e., 'pdlplot -f 
test.rdb position velocity' would take 'position' as the X-axis column name.

=item B<--y>=I<Y-axis column name> 

The name of the column in the rdb file that holds the dependent variable.  To
plot multiple columns simultaneously, give a comma-separated list of
the names of the columns after the --y, e.g. B<--y position,velocity>.
This can also be given as the second unflagged argument, i.e., 'pdlplot -f 
test.rdb position velocity' would take 'velocity' as the Y-axis column name.

=item B<--yerr>=I<Y-data uncertainties column name>

The name of the column in the rdb file that holds the uncertainties for
the Y values (comma-separated list if the y-data is a comma-separated list).

=item B<--residuals>

A boolean that tells the program to plot the difference between the first two
y data columns.  This data will be plotted in a second, smaller
pane.  The smaller pane will have identical x-values to the main data pane.

For plots that pull two or more y-data sets from the same rows (i.e., no
break column), the residuals are the difference of the first two specified
y-data columns.  For plots that use a break column, the residuals are
the interpolated differences between the first and the second break-column
sets of y-data.


=item B<--xlog>

Plot the X-axis in log space.  Default is linear plotting.

=item B<--ylog>

Plot the Y-axis in log space.  Default is linear plotting.

This option does not currently work with --ylog or --xlog (see below).

=item B<--xrange>=I<Comma-separated X range>

Specify a comma-separated non-default range for the X values.  Example: -x -5,5
will plot the data from x=-5 to x=5.  If the I<xlog> flag is on, the xrange
values must be specified in powers of ten.  E.G. -xlog -x -1,2 will plot
the data on a logged X range from 0.1 to 100.

=item B<--yrange>=I<Comma-separated Y range>

Specify a comma-separated non-default range for the X values.  Example: -y -5,5
will plot the data from y=-5 to y=5.  If the I<ylog> flag is on, the yrange
values must be specified in powers of ten.  E.G. -ylog -y -1,2 will plot
the data on a logged Y range from 0.1 to 100.

=item B<--title>=I<Plot title>

A title to go at the top of the plot.  Default is the y-axis name vs the x-axis name.

=item B<--subtitle>=I<Plot subtitle>

A subtitle to go at the top of the plot.  No default.

=item B<--xlabel>=I<Label of X axis>

A title for the x-axis of the plot.  Default is the x-axis name.

=item B<--ylabel>=I<Label of Y axis>

A title for the y-axis of the plot.  Default is the y-axis name.

=item B<--delta_label>=I<Label of the residuals y-axis>

A title for the y-axis of the residuals plot.  Default is "deltas".

=item B<--delta_pos>

A flag to move the numbering on the residuals plot to the right side of the
pane.

=item B<--nopoints>

Suppresses data point drawing; does not suppress the line drawn between points.

=item B<--noline>

Suppresses data line drawing; does not suppress data point drawing.

=item B<--colors>

A comma-separated string that specifies the color set to use for drawing points
and lines.  For example, if 'red,blue,black' is the argument, the first set of
points will be drawn red, the second blue, the third black, and then the fourth
will be drawn red agan.  The default is 
'black,red,green,blue,yellow,cyan,magenta,gray'.

=item B<--symbols>

A comma-separated string that specifies the symbol set to use for drawing
points and lines.  For example, if '0,3,4' is the argument, the first set of
points will be drawn with symbol 0, the second with symbol 3, and the third
with symbol 4.  The fourth set of points will be drawn with symbol 0, and so
forth.  The default is '5,3,0,4" followed by the sequence 6...99.

=item B<--font>=I<integer>

A PGPlot font integer.  The default is 1, and the range is 1-4.

=item B<--char_size>=I<character size in multiples of the standard size>

The PGPlot character size.  The default is 1.

=item B<--line_width>=I<line width in multiples of the standard width>

The PGPlot line width.  The default is 2.

=item B<--histogram>

Experimental, incorrect, histogram flag.  DO NOT USE!

=item B<--nolegend>

Using the --nolegend flag will prevent a legend from being drawn.
The default is for a legend to be drawn.

=item B<--legend_location>=I<legend coordinates>

A comma-separated list that to specify a location for the plot's legend.
The default is .05,-.05.  The coordates are in the range [0-1] for x, and
[0,-1] for y, with the origin in the upper-left corner of the plot.

=item B<--legend_text>=I<legend text>

A comma-separated list, with one item to specify the text for each set of
dependent data.  The list must be given in in the same order as the data sets
are given.

=item B<--device>=I<plotting device, or output file name/device name>

Either the name of your plotting device, or a filename/device pair.
The default is '/xs', other likely options are '/xt', '/xw', and '/xw'.
For Windows, I '/gw' is rumored to work, although this has not been 
tested.

Example:  for a color postscript file to be written, use:
-device 'out.ps/cps'.


=item B<--verbose>=I<verbosity level>

Verbosity level 0 results in silent operation.  Verbosity level 1 prints out
some runtime information (the plot limits, at least), and verbosity level 2 
is pretty comprehensive, and level 3 is painfully chatty.

=item B<--help>

Print out brief help information.

=item B<--usage>

Print out complete help information.

=head1 EXAMPLES

Using an RDB file:

B<pdlplot> -f mydata.rdb -x time -y position 

B<pdlplot> -f mydata.rdb -x time -y position -yerr err 

B<pdlplot> -f mydata.rdb -x time -y position -xr -1,2 -xlog -yr 0,100 

B<pdlplot> -f mydata.rdb -x time -y position -xr -1,2 -xlog -yr 0,100 -row 'energy eq 0.277'  

B<pdlplot> -file mydata.rdb -x t -y p -title "My Plot" -xlab time -ylab pos -device 'out.ps/cps'

B<pdlplot> -file mydata.rdb -x time -y position,velocity

Using a comma-delimitted file:

B<pdlplot> -f mydata.dat -split ',' -x time -y position 

B<pdlplot> -f mydata.dat -split ',' -x time -y position -yerr err 


=head1 LICENSE

This software is released under the GNU General Public License.  You
may find a copy at 

   http://www.fsf.org/copyleft/gpl.html

=head1 AUTHOR

Kester Allen (callen@cfa.harvard.edu)

=cut