pp_addpm({At=>Top},<<'EOD');

=head1 NAME

PDL::IO::Misc - misc IO routines for PDL

=head1 DESCRIPTION

A mixture of basic I/O functionality

=head1 SYNOPSIS

 use PDL::IO::Misc;

=cut

EOD

pp_add_exported('',"rfits wfits rcols wcols rgrep rdsa");

############################## PM CODE ########################################


pp_addpm(<<'EOD');

use PDL::Primitive;
use PDL::Types;
use PDL::Options;
use Carp;

=head2 rfits()

=for ref

Simple piddle FITS reader.

=for example

  $pdl = rfits('file.fits');

Suffix magic:

  # Automatically uncompress via gunzip pipe
  $pdl = rfits('file.fits.gz');
  # Automatically uncompress via uncompress pipe 
  $pdl = rfits('file.fits.Z');  

FITS Headers stored in piddle and can be retrived with $a->gethdr.

This header is a reference to a hash where the hash keys are the 
keywords in the FITS header. It is important to note that for strings, 
the surrounding quotes are kept to ensure that strings that look like
numbers are kept as strings. This is also of importance if you create 
your own header information and you want the value to be printed out 
as a string.

Comments in headers are stored as $$h{COMMENT}{<Keyword>} where $h is
the header retrieved with $a->gethdr.
History entries in the header are stored as $$h{HISTORY}, which is an
anonymous array for each HISTORY entry in the header.

=cut

sub rfits {PDL->rfits(@_)}

sub PDL::rfits {
   my $class = shift;
   barf 'Usage: $a = rfits($file); $a = PDL->rfits("...")' if $#_!=0;
   my $file = shift; my $pdl  = $class->new;
   my($nbytes, $line, $name, $rest, $size, $i, $bscale, $bzero);

   $file = "gunzip -c $file |" if $file =~ /\.gz$/;    # Handle compression
   $file = "uncompress -c $file |" if $file =~ /\.Z$/;

   open(FITS, $file) || barf "FITS file $file not found";
   binmode FITS;
   $nbytes = 0; # Number of bytes read so far
   $line = "";

   my $foo={};       # To go in pdl
   $$foo{"BSCALE"}=1;
   $$foo{"BZERO"}=0;
   my @history=();
   while( !eof(FITS)) {
      read(FITS,$line,80);
      barf "file $file is not in FITS-format:\n$line\n"
                  if( $nbytes==0 && ($line !~ /^SIMPLE  = +T/));
      $nbytes += 80;

      $name = (split(' ',substr($line,0,8)))[0]; $rest=substr($line,8);
      if ($name =~ m/^HISTORY/) {
         push @history, $rest;
      } else {
         $$foo{$name} = "";
         undef $comment;
         $$foo{$name}=$1 if $rest =~ m|^= +(.*\S) *$| ;
         $$foo{$name}=$1 if $rest =~ m|^= +(.*\S) +/(.*)$| ;
         $$foo{$name}="'".$1."'" if $rest =~ m|^= '(.*)' *$| ;
         $$foo{$name}="'".$1."'" if $rest =~ m|^= '(.*)' +/(.*)$| ;
         $$foo{COMMENT}{$name} = $2 if defined($2);
         last if $name eq "END";
      }
   }
   $$foo{"HISTORY"} = \@history if $#history >= 0;
   $nbytes %= 2880;
   my $bar; read(FITS,$bar, 2880-$nbytes) if $nbytes!=0; # Skip to end of card

   # Setup piddle structure

   $pdl->set_datatype($PDL_B)    if $$foo{"BITPIX"} ==   8;
   $pdl->set_datatype($PDL_S)    if $$foo{"BITPIX"} ==  16;
   $pdl->set_datatype($PDL_L)    if $$foo{"BITPIX"} ==  32;
   $pdl->set_datatype($PDL_F)    if $$foo{"BITPIX"} == -32;
   $pdl->set_datatype($PDL_D)    if $$foo{"BITPIX"} == -64;

   my @dims; # Store the dimenions 1..N, compute total number of pixels
   $size = 1;  $i=1;
   while(defined( $$foo{"NAXIS$i"} )) {
     $size = $size*$$foo{"NAXIS$i"};
     push @dims, $$foo{"NAXIS$i"} ; $i++;
   }
   $pdl->setdims([@dims]);

   my $dref = $pdl->get_dataref();

   print "BITPIX = ",$$foo{"BITPIX"}," size = $size pixels \n"
         if $PDL::verbose;

   # Slurp the FITS binary data

   print "Reading ",$size*PDL::Core::howbig($pdl->get_datatype) , "bytes\n" if
   $PDL::verbose;

   read( FITS, $$dref, $size*PDL::Core::howbig($pdl->get_datatype) );

   close(FITS);
   $pdl->upd_data();

   if (!isbigendian() ) { # Need to byte swap on little endian machines
      bswap2($pdl) if $pdl->get_datatype == $PDL_S;
      bswap4($pdl) if $pdl->get_datatype == $PDL_L || $pdl->get_datatype ==
      $PDL_F;
      bswap8($pdl) if $pdl->get_datatype == $PDL_D;
   }

   $bscale = $$foo{"BSCALE"}; $bzero = $$foo{"BZERO"};
   print "BSCALE = $bscale &&  BZERO = $bzero\n" if $PDL::verbose;
   $bscale = 1 if $bscale eq "";
   $bzero  = 0 if $bzero  eq "";
   
   # Be clever and work out the final datatype before eating
   # memory
   
   my $tmp = $pdl->clump(-1)->slice("0:0");
   $tmp = $tmp*$bscale if $bscale != 1; # Dummy run on one element of $pdl
   $tmp = $tmp+$bzero  if $bzero  != 0;
  
   $pdl = $pdl->convert($tmp->type) if $tmp->get_datatype != $pdl->get_datatype;
   
   $pdl *= $bscale if $bscale != 1; 
   $pdl += $bzero  if $bzero  != 0;

   delete $$foo{"BSCALE"}; delete $$foo{"BZERO"};

   delete $$foo{"BSCALE"}; delete $$foo{"BZERO"};

   # Header

   $pdl->sethdr($foo);

   return $pdl;
}

=head2 wfits()

=for ref

Simple piddle FITS writer

=for example

  wfits $pdl, 'filename.fits', [$BITPIX];
  $pdl->wfits('foo.fits',-32);

Suffix magic:

  # Automatically compress through pipe to gzip
  wfits $pdl, 'filename.fits.gz';
  # Automatically compress through pipe to compress 
  wfits $pdl, 'filename.fits.Z';  

$BITPIX is optional and coerces the output format.

=cut

*wfits = \&PDL::wfits;

sub PDL::wfits { # Write a PDL to a FITS format file

   barf 'Usage: wfits($pdl,$file,[$BITPIX])' if $#_<1 || $#_>2;

   my ($pdl,$file,$BITPIX) = @_;
   my ($k, $buff, $off, $ndims, $sz);
   local($nbytes, %hdr);

   if ($file =~ /\.gz$/) {            # Handle compression
      $file = "|gzip -9 > $file";
   }
   elsif ($file =~ /\.Z$/) {
      $file = "|compress > $file";
   }
   else{
      $file = ">$file";
   }

   # Figure output type

   $BITPIX = "" unless defined $BITPIX;
   if ($BITPIX eq "") {
      $BITPIX =   8 if $pdl->get_datatype == $PDL_B;
      $BITPIX =  16 if $pdl->get_datatype == $PDL_S || $pdl->get_datatype == $PDL_US;
      $BITPIX =  32 if $pdl->get_datatype == $PDL_L;
      $BITPIX = -32 if $pdl->get_datatype == $PDL_F;
      $BITPIX = -64 if $pdl->get_datatype == $PDL_D;
   }
   my $convert = sub { return $_[0] };# Default - do nothing
   $convert = sub {byte($_[0])}   if $BITPIX ==   8;
   $convert = sub {short($_[0])}  if $BITPIX ==  16;
   $convert = sub {long($_[0])}   if $BITPIX ==  32;
   $convert = sub {float($_[0])}  if $BITPIX == -32;
   $convert = sub {double($_[0])} if $BITPIX == -64;

   # Automatically figure output scaling

   $bzero = 0; $bscale = 1;
   if ($BITPIX>0) {
      my $min = $pdl->min;
      my $max = $pdl->max;
      my ($dmin,$dmax) = (0, 2**8-1)     if $BITPIX == 8;
      ($dmin,$dmax) = (-2**15, 2**15-1)  if $BITPIX == 16;
      ($dmin,$dmax) = (-2**31, 2**31-1)  if $BITPIX == 32;

      if ($min<$dmin || $max>$dmax) {
         $bzero = $min;
         $max -= $bzero;
         $bscale = $max/$dmax if $max>$dmax;
       }
       print "BSCALE = $bscale &&  BZERO = $bzero\n" if $PDL::verbose;
   }

   open(FITS, "$file") || barf "Unable to create FITS file $file\n";
   binmode FITS;
   printf FITS "%-80s", "SIMPLE  =                    T ";

   $nbytes = 80; # Number of bytes written so far

   # Write FITS header

   %hdr = ();
   my $h = $pdl->gethdr;
   if (defined($h)) {
      for (keys %$h) { $hdr{$_} = $$h{$_} } # Copy
   }

   delete $hdr{SIMPLE}; delete $hdr{'END'};

   $hdr{BITPIX} =  $BITPIX;
   $hdr{BUNIT} = "Data Value" unless exists $hdr{BUNIT};
   wheader('BITPIX');

   $ndims = $pdl->getndims; # Dimensions of data array
   $hdr{NAXIS}  = $ndims;
   wheader('NAXIS');
   for $k (1..$ndims) { $hdr{"NAXIS$k"} = $pdl->getdim($k-1) }
   for $k (1..$ndims) { wheader("NAXIS$k") }

   if ($bscale != 1 || $bzero  != 0) {
      $hdr{BSCALE} =  $bscale;
      $hdr{BZERO}  =  $bzero;
      wheader('BSCALE');
      wheader('BZERO');
   }
   wheader('BUNIT');
   for $k (sort keys %hdr) { wheader($k) unless $k =~ m/HISTORY/}
   wheader('HISTORY'); # Make sure that HISTORY entries come last.
   printf FITS "%-80s", "END"; $nbytes += 80;
   $nbytes %= 2880;
   print FITS " "x(2880-$nbytes) if $nbytes != 0; # Fill up HDU

   # Decide how to byte swap - note does not quite work yet. Needs hack
   # to IO.xs

   my $bswap = sub {};     # Null routine
   if ( !isbigendian() ) { # Need to set a byte swap routine
      $bswap = \&bswap2 if $BITPIX==16;
      $bswap = \&bswap4 if $BITPIX==32 || $BITPIX==-32;
      $bswap = \&bswap8 if $BITPIX==-64;
   }

   # Write FITS data

   my $p1d = $pdl->clump(-1); # Data as 1D stream

   $off = 0;
   $sz  = PDL::Core::howbig(&$convert($p1d->slice('0:0'))->get_datatype);

   $nbytes = $p1d->getdim(0) * $sz;

   # Transfer data in blocks (because might need to byte swap)
   # Buffer is also type converted on the fly

   my $BUFFSZ = 360*2880; # = ~1Mb - must be multiple of 2880
   my $tmp;

   while ($nbytes - $off > $BUFFSZ) {

      # Data to be transferred

      $buff = &$convert( ($p1d->slice( ($off/$sz).":". (($off+$BUFFSZ)/$sz-1))
                         -$bzero)/$bscale );
      &$bswap($buff); print  FITS ${$buff->get_dataref};
      $off += $BUFFSZ;
   }
   $buff = &$convert( ($p1d->slice($off/$sz.":-1") - $bzero)/$bscale );
   &$bswap($buff); print  FITS ${$buff->get_dataref};
   print FITS " "x(($BUFFSZ - $buff->getdim(0) * $sz)%2880);  # Fill HDU

   close(FITS);

1;}

sub wheader {     # Local utility routine of wfits()
   my $k = shift;

   if ($k =~ m/HISTORY/) {
      return unless ref($hdr{$k}) eq 'ARRAY';
      foreach my $line (@{$hdr{$k}}) {
         printf FITS "HISTORY %-72s", substr($line,0,72);
       $nbytes += 80;
      }
      delete $hdr{$k};
   } else {
      # Check that we are dealing with a scalar value in the header
      # Need to make sure that the header does not include PDLs or
      # other structures. Return unless $hdr{$k} is a scalar

      return unless not ref($hdr{$k});

      if ($hdr{$k} eq "") {
         printf FITS "%-80s", substr($k,0,8);
      } else {
         printf FITS "%-8s= ", substr($k,0,8);
         if ($hdr{$k} =~ /^ *([+-]?)(?=\d|\.\d)\d*(\.\d*)?([Ee]([+-]?\d+))? *$/) { # Number?
         my $cl=60-((exists($hdr{COMMENT}{$k})) ? 2 : 0);
         my $end=' ' x $cl;
         $end =' /'.$hdr{COMMENT}{$k} if (exists($hdr{COMMENT}{$k}));
         printf FITS "%20s%-50s", substr($hdr{$k},0,20),
         substr($end, 0, 50);
        } elsif ($hdr{$k} eq 'F' or $hdr{$k} eq 'T') {
          # Logical flags ?
          printf FITS "%20s", $hdr{$k};
          my $end=' ' x 50;
          $end =' /'.$hdr{COMMENT}{$k} if (exists($hdr{COMMENT}{$k}));
          printf FITS "%-50s", $end;
        }
         else {
         $hdr{$k}=$1 if ($hdr{$k}=~m/'(.*)'/); # Take off surrounding ''
         my $sl=length($hdr{$k})+2;
         $sl=10 if $sl < 10;
         my $cl=68-$sl-(exists($hdr{COMMENT}{$k}) ? 2 : 0);
         printf FITS "'%-$ {sl}s'", substr($hdr{$k},0,$sl);
         if (exists($hdr{COMMENT}{$k})) {
               printf FITS " /%-$ {cl}s", substr($hdr{COMMENT}{$k}, 0, $cl);
           } else {
             printf FITS "%-$ {cl}s", ' ' x $cl;
            }
         }
      }
      $nbytes += 80; delete $hdr{$k};
   }
   delete $hdr{COMMENT}{$k};
1;}

# ***CROCK*** Internal routine to extend 1D PDL array by size $n - dirty hack
# - needs a proper extend function rather than this nasty recreation
# - changed to ensure input/output piddles have the same type

sub ext1D {
   my ($a,$n) = @_;
   my $nold   = $a->getdim(0);
   my $b      = zeroes($a->type,$nold+$n);  # New pdl
   my $bb     = $b->slice("0:".($nold-1));
   $bb       .= $a;
   $_[0]      = $b;
1;}

# taken outside of rcols() to avoid clutter
sub _handle_types ($$$) {
    my $ncols = shift;
    my $deftype = shift;
    my $types = shift;

    barf "Unknown PDL type given for DEFTYPE.\n"
	unless ref($deftype) eq "PDL::Type";

    my @cols = ref($types) eq "ARRAY" ? @$types : ();
	
    if ( $#cols > -1 ) {
  	# truncate if required
  	$#cols = $ncols if $#cols > $ncols;
	
  	# check input values are sensible
  	for ( 0 .. $#cols ) {
  	    barf "Unknown value '$cols[$_]' in TYPES array.\n" 
  		unless ref($cols[$_]) eq "PDL::Type";
  	}
    }

    # fill in any missing columns
    for ( ($#cols+1) .. $ncols ) { push @cols, $deftype; }

    return @cols;
} # sub: _handle_types

=head2 rcols()

=for ref

Read ASCII whitespaced cols from a file into piddles and perl arrays
(also see L</rgrep()>).

There are two calling conventions - the old version, where
a pattern can be specified after the filename/handle, and the new
version where options are given as as hash reference. This 
reference can be given as either the second or last argument.

The default behaviour is to ignore lines beginning with a # character
and lines that only consist of whitespace. Options exist to only read
from lines that match, or do not match, supplied patterns, and
to set the types of the created piddles.

Can take file name or *HANDLE, and if no columns are specified, 
all are assumed. For the allowed types, see
L<PDL::Core/Datatype_conversions>.

Options:

EXCLUDE or IGNORE
- ignore lines matching this pattern (default B<'/^#/'>).

INCLUDE or KEEP
- only use lines which match this pattern (default B<''>).

LINES   
- which line numbers to use. Line numbers start at 0 and the syntax 
is 'a:b:c' to use every c'th matching line between a and b
(default B<''>).

DEFTYPE
- default data type for stored data (if not specified, use the type 
stored in C<$PDL::IO::Misc::deftype>, which starts off as B<double>).

TYPES
- reference to an array of data types, one element for each column 
to be read in.
Any missing columns use the DEFTYPE value (default B<[]>).

PERLCOLS
- an array of column numbers which are to be read into perl arrays
rather than piddles. References to these arrays are returned after
the requested piddles (default B<undef>).

=for usage

  Usage:
    ($x,$y,...) = rcols( *HANDLE|"filename", { EXCLUDE => '/^!/' }, 
			 $col1, $col2, ... )
    ($x,$y,...) = rcols( *HANDLE|"filename", $col1, $col2, ..., 
			 { EXCLUDE => '/^!/' } )
    ($x,$y,...) = rcols( *HANDLE|"filename", "/foo/", $col1, $col2, ... )
 
e.g.,

=for example

  $x      = PDL->rcols 'file1';
  ($x,$y) = rcols *STDOUT;

  # read in lines containing the string foo, where the first
  # example also ignores lines that with a # character.
  ($x,$y,$z) = rcols 'file2', 0,4,5, { INCLUDE => '/foo/' };
  ($x,$y,$z) = rcols 'file2', 0,4,5, 
                 { INCLUDE => '/foo/', EXCLUDE => '' };

  # ignore the first 27 lines of the file, reading in as ushort's
  ($x,$y) = rcols 'file3', { LINES => '27:-1', DEFTYPE => ushort };
  ($x,$y) = rcols 'file3', 
              { LINES => '27:', TYPES => [ ushort, ushort ] };

  # read in the first column as a perl array and the next two as piddles
  ($x,$y,$name) = rcols 'file4', 1, 2, { PERLCOLS => [ 0 ] };
  printf "Number of names read in = %d\n", 1 + $#$name;

Notes:

1. Quotes are required on patterns.

2. Columns are separated by whitespace by default,
use C<$PDL::IO::Misc::colsep> to specify an alternate
separator.

3. For PDL-2.003, the meaning of the 'c' value in the LINES option has
changed: it now only counts matching lines rather than all lines as in
previous versions of PDL. 

4. LINES => '-1:0:3' may not work as you expect, since lines are skipped
when read in, then the whole array reversed.

=cut

$colsep = " ";      # Default column separator
$deftype = double;  # Default type for piddles

# NOTE: XXX
#  need to look at the line-selection code. For instance, if want
#   lines => '-1:0:3', 
#  read in all lines, reverse, then apply the step
#  -> fix point 4 above
# 

sub rcols{PDL->rcols(@_)}

sub PDL::rcols {
   my $class = shift;
   barf 'Usage ($x,$y,...) = rcols( *HANDLE|"filename", ["/pattern/" or \%options], $col1, $col2, ..., [ \%options] )' 
       if $#_<0;

   my $is_handle = ref(\$_[0]) eq "GLOB";
   my $fh = $is_handle ? $_[0] : "FILE";
   open $fh, $_[0] or die "File $_[0] not found\n" unless $is_handle;
   shift;

   # set up default options
   my $opt = new PDL::Options( {
       EXCLUDE => '/^#/',
       INCLUDE => '',
       LINES => '',
       DEFTYPE => $deftype,
       TYPES   => [],
       PERLCOLS => undef
       } );
   $opt->synonyms( { IGNORE => 'EXCLUDE', KEEP => 'INCLUDE' } );

   # has the user supplied any options
   if ( defined($_[0]) ) {
       # ensure the old-style behaviour by setting the exclude pattern to ''
       if ( $_[0] =~ m|^/.*/$| )      { $opt->options( { EXCLUDE => '', INCLUDE => shift } ); }
       elsif ( ref($_[0]) eq "HASH" ) { $opt->options( shift ); }
   }

   # maybe the last element is a hash array as well
   $opt->options( pop ) if defined($_[-1]) and ref($_[-1]) eq "HASH";

   # a reference to a hash array
   my $options = $opt->current();   

   # what are the patterns?
   foreach my $pattern ( qw( INCLUDE EXCLUDE ) ) {
       if ( $$options{$pattern} ne '' ) {
         barf "rcols() - unable to parse $pattern value.\n" if $$options{$pattern} !~ m|^/.*/$|;
         $$options{$pattern} =~ s|^/(.*)/$|$1|;
       }
   }

   # which columns are to be read into piddles and which into perl arrays?
   my @perl_cols = ();
   @perl_cols = @{ $$options{PERLCOLS} } if $$options{PERLCOLS};
   my ( @pdl_cols )  = @_;

   # work out which line numbers are required
   # - the regexp's are a bit over the top
   my ( $a, $b, $c );
   if ( $$options{LINES} ne '' ) {
       if ( $$options{LINES} =~ /^\s*([+-]?\d*)\s*:\s*([+-]?\d*)\s*$/ ) {
	   $a = $1; $b = $2;
       } elsif ( $$options{LINES} =~ /^\s*([+-]?\d*)\s*:\s*([+-]?\d*)\s*:\s*([+]?\d*)\s*$/ ) {
	   $a = $1; $b = $2; $c = $3;
       } else {
	   barf "rcols() - unable to parse LINES option.\n";
       }
   }

   # Since we do not know how many lines there are in advance, things get a bit messy
   my ( $index_start, $index_end ) = ( 0, -1 );
   $index_start  = $a if defined($a) and $a ne '';
   $index_end    = $b if defined($b) and $b ne '';
   my $line_step = $c || 1;

   # $line_rev = 0/1 for normal order/reversed
   # $line_start/_end refer to the first and last line numbers that we want
   # (the values of which we may not know until we've read in all the file)
   my ( $line_start, $line_end, $line_rev );
   if ( ($index_start >= 0 and $index_end < 0) ) {
       # eg 0:-1
       $line_rev = 0; $line_start = $index_start;
   } elsif ( $index_end >= 0 and $index_start < 0 ) {
       # eg -1:0
       $line_rev = 1; $line_start = $index_end; 
   } elsif ( $index_end >= $index_start and $index_start >= 0 ) {
       # eg 0:10
       $line_rev = 0; $line_start = $index_start; $line_end = $index_end;
   } elsif ( $index_start > $index_end and $index_end >= 0 ) {
       # eg 10:0
       $line_rev = 1; $line_start = $index_end; $line_end = $index_start;
   } elsif ( $index_start <= $index_end ) {
       # eg -5:-1
       $line_rev = 0;
   } else {
       # eg -1:-5
       $line_rev = 1;
   }
   
   my @ret;

   my (@v,$k); 

   my $line_num = -1;
   my $line_ctr = $line_step - 1;  # ensure first line is always included
   my $index    = -1;
   my $pdlsize  =  0;
   my $extend   = 10000;

   my $line_store;  # line numbers of saved data

   while(<$fh>) {

       $line_num++;

       # the order of these checks is important, particularly whether we
       # check for line_ctr before or after the pattern matching
       # Prior to PDL 2.003 the line checks were done BEFORE the
       # pattern matching
       last if     defined($line_end)   and $line_num > $line_end;
       next if     defined($line_start) and $line_num < $line_start;
       next if     $$options{EXCLUDE} ne '' and /$$options{EXCLUDE}/;
       next unless $$options{INCLUDE} eq '' or  /$$options{INCLUDE}/;
       next unless ++$line_ctr == $line_step;
       $line_ctr = 0;

       $index++;
       @v = $colsep eq ' ' ? split(' ') : split($colsep) ; 

       # if the first line, set up the output piddles
       # using all the columns if the user doesn't specify anything
       if ( $index == 0 ) {
         @pdl_cols = ( 0 .. $#v ) if $#pdl_cols < 0;

	 # sort out the types of the piddles
	 my @types = _handle_types( $#pdl_cols, $$options{DEFTYPE}, $$options{TYPES} );
	 if ( $PDL::verbose ) { # dbg aid
	     print "Reading data into piddles of type: [ ";
	     foreach my $t ( @types ) {
		 # not sure this is the official way of doing it
		 my $name = $PDL::Types::typehash{$PDL::Types::names[$t->[0]]}->{ctype};
		 $name =~ s/PDL_//;
		 print "$name ";
	     }
	     print "]\n";
	 }

         $k = 0;
         for (0..$#pdl_cols) { $ret[$_] = $class->zeroes($types[$_],1); $k++; }
         for (@perl_cols)    { $ret[$k++] = []; }
	 $line_store = $class->zeroes(long,1); # only need to store integers
       }

       # if necessary, extend PDL in buffered manner
       if ( $pdlsize < $index ) {
         for (0..$#pdl_cols) { ext1D( $ret[$_], $extend ); }
	   ext1D( $line_store, $extend );
	   $pdlsize += $extend;
       }

       # Set values - '1*' is split() bug workaround
       # - stick perl arrays onto end of $ret
       $k = 0;
       for (@pdl_cols)  { set $ret[$k++], $index, 1*$v[$_]; }
       for (@perl_cols) { push @{ $ret[$k++] }, $v[$_]; }

       # store the line number
       $line_store->set( $index, $line_num );
   }

   close($fh) unless $is_handle;

   # have we read anything in? if not, return empty piddles
   if ( $index == -1 ) {
       print "Warning: rcols() did not read in any data.\n" if $PDL::verbose;
       if ( wantarray ) {
         foreach ( 0 .. $#pdl_cols ) { $ret[$_] = PDL->null; }
         for ( ($#pdl_cols+1) .. ($#pdl_cols+1+$#perl_cols) ) { $ret[$_] = []; }
	   return ( @ret );
       } else { 
	   return PDL->null;
       }
   }

   # if the user has asked for lines => 0:-1 or 0:10 or 1:10 or 1:-1,
   # - ie not reversed and the last line number is known -
   # then we can skip the following nastiness
   if ( $line_rev == 0 and $index_start >= 0 and $index_end >= -1 ) {
       for ( 0 .. $#pdl_cols ) { $ret[$_] = $ret[$_]->slice("0:${index}"); };
       if ( $PDL::verbose ) {
	   if ( $#pdl_cols != -1 ) { print "Read in ", $ret[0]->nelem, " elements.\n"; }
	   else                    { print "Read in ", $#{ $ret[0] }, " elements.\n"; }
       }
       wantarray ? return(@ret) : return $ret[0];
   }

   # Work out which line numbers we want. First we clean up the piddle
   # containing the line numbers that have been read in
   $line_store = $line_store->slice("0:${index}");

   # work out the min/max line numbers required
   if ( $line_rev ) {
       if ( defined($line_start) and defined($line_end) ) {
	   my $dummy = $line_start;
	   $line_start = $line_end;
	   $line_end = $dummy;
       } elsif ( defined($line_start) ) {
	   $line_end = $line_start;
       } else {
	   $line_start = $line_end; 
       }
   }
   $line_start = $line_num + 1 + $index_start if $index_start < 0;
   $line_end   = $line_num + 1 + $index_end   if $index_end   < 0;

   my $indices;
   if ( $line_rev ) {
       $indices = which( $line_store >= $line_end & $line_store <= $line_start )->slice('-1:0');
   } else {
       $indices = which( $line_store >= $line_start & $line_store <= $line_end );
   }

   # truncate the piddles
   for ( 0 .. $#pdl_cols ) { $ret[$_] = $ret[$_]->index($indices); };

   # truncate/reverse/etc the perl arrays
   my @indices_array = list $indices;
   $k = $#pdl_cols + 1;
   foreach ( 0 .. $#perl_cols ) {
       my @temp = @{ $ret[$k] };
       $ret[$k] = [];
       foreach my $i ( @indices_array ) { push @{ $ret[$k] }, $temp[$i] };
       $k++;
   }

   if ( $PDL::verbose ) {
       if ( $#pdl_cols != -1 ) { print "Read in ", $ret[0]->nelem, " elements.\n"; }
       else                    { print "Read in ", $#{ $ret[0] }, " elements.\n"; }
   }
   wantarray ? return(@ret) : return $ret[0];
}


=head2 wcols()

=for ref

Write ASCII whitespaced cols into file from piddles efficiently.

If no columns are specified all are assumed.
Will optionally only process lines matching a pattern.
Can take file name or *HANDLE, and
if no file/filehandle is given defaults to STDOUT.

Options:

HEADER
- prints this string before the data. If the string
is not terminated by a newline, one is added
(default B<''>).

=for usage

 Usage: wcols $piddle1, $piddle2,..., *HANDLE|"outfile", [\%options];

e.g.,

=for example

  wcols $x, $y+2, 'foo.dat';
  wcols $x, $y+2, *STDERR;
  wcols $x, $y+2, '|wc';
  wcols $a,$b,$c; # Orthogonal version of 'print $a,$b,$c' :-)

  wcols "%10.3f", $a,$b; # Formatted
  wcols "%10.3f %10.5g", $a,$b; # Individual column formatting

  wcols $a,$b, { HEADER => "#   a   b" };

Note: columns are separated by whitespace by default,
use $PDL::IO::Misc::colsep to specify an alternate
separator.

=cut

*wcols = \&PDL::wcols;

sub PDL::wcols {
   barf 'Usage: wcols(@[$format_string], vectors,*HANDLE|"filename", [\%options])' if @_<1;

   my ($format_string, $step, $fh);
   if (ref(\$_[0]) eq "SCALAR") {
       $step = $format_string = shift; # 1st arg not piddle
       $step =~ s/(%%|[^%])//g;  # use step to count number of format items
       $step = length ($step);
   }
   # if last argument is a reference to a hash, parse the options
   my $header;
   if ( ref( $_[-1] ) eq "HASH" ) {
       my $opt = pop;
       foreach my $key ( keys %$opt ) {
	   if ( $key =~ /^H/i ) { $header = $opt->{$key}; }  # option: HEADER
	   else {
	       print "Warning: rcols does not understand option <$key>.\n"; 
	   }
       }
   }
   my $file = $_[-1];
   my $file_opened;
   if (ref(\$file) eq "GLOB") {  # file handle passed directly
       $fh = $file; pop;
   }
   else{
       if (ref(\$file) eq "SCALAR") {  # Must be a file name
          $fh = "FILE";
         if (!$is_handle) {
            $file = ">$file" unless $file =~ /^\|/ or $file =~ /^\>/;
             open $fh, $file or barf "File $file can not be opened for writing\n";
         }
          pop;
	  $file_opened = 1;
       }
       else{  # Not a filehandle or filename, assume something else
              # (probably piddle) and send to STDOUT
          $fh = *STDOUT;
       }
   }

   my @p = @_;
   my $n = $p[0]->nelem;
   for (@p) {
      barf "wcols: 1d args must have same number of elements\n"
         if $_->nelem != $n or $_->getndims!=1;
   }
   if ( defined $header ) {
       $header .= "\n" unless $header =~ m/\n$/;
       print $fh $header;
   }
   my $i;
   for ($i=0; $i<$n; $i++) {
       if ($format_string) {
	   my @d;
	   for (@p) {
	       push @d,$_->at($i);
	       if (@d == $step) {
		   printf $fh $format_string,@d;
		   printf $fh $colsep;
		   $#d = -1;
	       }
	   }
	   if (@d && !$i) {
	       my $str;
	       if ($#p>0) {
		   $str = ($#p+1).' columns don\'t';
	       } else {
		   $str = '1 column doesn\'t';
	       }
	       $str .= " fit in $step column format ".
	       '(even repeated) -- discarding surplus';
	       carp $str;
	       # printf $fh $format_string,@d;
	       # printf $fh $colsep;
	   }
       } else {
	   for (@p) {
	       print $fh $_->at($i),$colsep;
	   }
       }
       print $fh "\n";
   }
   close($fh) if $file_opened;
   return 1;
}



=head2 rgrep()

=for ref

Read columns into piddles using full regexp pattern matching.

Usage

=for usage

 ($x,$y,...) = rgrep(sub, *HANDLE|"filename")

e.g.

=for example

 ($a,$b) = rgrep {/Foo (.*) Bar (.*) Mumble/} $file;

i.e. the vectors C<$a> and C<$b> get the progressive values
of C<$1>, C<$2> etc.

=cut

  sub rgrep (&@) {
     barf 'Usage ($x,$y,...) = rgrep(sub, *HANDLE|"filename")'
         if $#_!=1;

     my (@ret,@v,$nret); my ($m,$n)=(-1,0); # Count/PDL size
     my $pattern = shift;

     my $is_handle = ref(\$_[0]) eq "GLOB";
     my $fh = $is_handle ? $_[0] : "FILE";
     open $fh, $_[0] or die "File $_[0] not found\n" unless $is_handle;

     if (ref($pattern) ne "CODE") {
         die "Got a ".ref($pattern)." for rgrep?!";
     }

     while(<$fh>) {
         next unless @v = &$pattern;

         $m++;  # Count got
         if ($m==0) {
	   $nret = $#v;   # Last index of values to return
	   for (0..$nret) {
	       $ret[$_] = double(pdl([0])); # Create PDLs
	   }
       } else { # perhaps should only carp once...
           carp "Non-rectangular rgrep" if $nret != $#v;
       }
       if ($n<$m) {
	   for (0..$nret) {
	       ext1D( $ret[$_], 10000 ); # Extend PDL in buffered manner
	   }
	   $n += 10000;
      }
       for(0..$nret) { set $ret[$_], $m, 1*$v[$_] } # Set values - '1*' is
   }                                                      # ensures numeric
   close($fh) unless $is_handle;
   for (@ret) { $_ = $_->slice("0:$m")->copy; }; # Truncate
   wantarray ? return(@ret) : return $ret[0];
}


=head2 rdsa()

=for ref

Read a FIGARO/NDF format file.

Requires non-PDL DSA module. Contact Frossie (frossie@jach.hawaii.edu)
Usage:

=for usage

 ([$xaxis],$data) = rdsa($file)

=for example

 $a = rdsa 'file.sdf'

Not yet tested with PDL-1.9X versions

=cut

sub rdsa{PDL->rdsa(@_)}

sub PDL::rdsa {
    my $class = shift;
    barf 'Usage: ([$xaxis],$data) = rdsa($file)' if $#_!=0;
    my $file = shift; my $pdl = $class->new; my $xpdl;
    eval 'use DSA' unless $dsa_loaded++;
    barf 'Cannot use DSA library' if $@ ne "";

    $status = 0;

    # Most of this stuff stolen from Frossie:

    dsa_open($status);
    dsa_named_input('IMAGE',$file,$status);
    goto skip if $status != 0;

    dsa_get_range('IMAGE',$vmin,$vmax,$status);
    dsa_data_size('IMAGE',5, $data_ndims, \@data_dims, $data_elements, $status);
    dsa_map_data('IMAGE','READ','FLOAT',$data_address,$data_slot,$status);

    @data_dims = @data_dims[0..$data_ndims-1];
    print "Dims of $file = @data_dims\n" if $PDL::verbose;
    $pdl->set_datatype($PDL_F);
    $pdl->setdims([@data_dims]);
    my $dref = $pdl->get_dataref;
    mem2string($data_address,4*$data_elements,$$dref);
    $pdl->upd_data();

    if (wantarray) { # Map X axis values
      dsa_axis_size('IMAGE',1,5, $axis_ndims, \@axis_dims,
                    $axis_elements, $status);
      dsa_map_axis_data('IMAGE',1,'READ','FLOAT',$axis_address,
                    $axis_slot,$status);
      @axis_dims = @axis_dims[0..$axis_ndims-1];
      $xpdl = $class->new;
      $xpdl->set_datatype($PDL_F);
      $xpdl->setdims([@axis_dims]);
      my $xref = $xpdl->get_dataref;
      mem2string($axis_address,4*$axis_elements,$$xref);
      $xpdl->upd_data();
    }

    skip: dsa_close($status);

    barf("rdsa: obtained DSA error") if $status != 0;

    return ($xpdl,$pdl);
}


################################ XS CODE ######################################
EOD

sub defpdl {
	pp_def(
		$_[0],
		Pars => $_[1],
		OtherPars => $_[2],
		Code => $_[3],
		Doc  => $_[4],
	);
}

pp_addpm(<<'EOD');

=head2 isbigendian()

=for ref

Determine endianness of machine - returns 0 or 1 accordingly

=cut
EOD

pp_addxs('','

MODULE = PDL::IO::Misc PACKAGE = PDL

int
isbigendian()
   CODE:
     unsigned short i;
     PDL_Byte *b;

     i = 42; b = (PDL_Byte*) (void*) &i;

     if (*b == 42)
        RETVAL = 0;
     else if (*(b+1) == 42)
        RETVAL = 1;
     else
         barf("Impossible - machine is neither big nor little endian!!\n");
     OUTPUT:
       RETVAL
');

pp_addpm("*isbigendian = \\&PDL::isbigendian;\n");
pp_add_exported("", "isbigendian");


###### Read ASCII Function ##########
pp_addhdr(<<'EOH');

#define SWALLOWLINE(fp) while ((s = PerlIO_getc(fp)) != '\n' && s != EOF)
#define    TRAILING_WHITESPACE_CHECK(s) \
   if (s!=' ' && s!='\t' && s!='\r' && s!='\n' && s!=',')  return -1

int getfloat(PerlIO *fp, PDL_Float *fz)
{
  PDL_Float f = 0;
  int nread = 0;
  int i, s = PerlIO_getc(fp);
  int afterp = 0, aftere=0;
  int expo = 0;
  PDL_Float sig = 1.0, esig = 1.0;
  PDL_Float div = 1.0;

  if (s == EOF) return 0;
  while (1) {
    if (s == EOF)
      return 0;   /* signal end of line */
    if (s == '#')
      SWALLOWLINE(fp);
    if ((s >='0' && s <='9') || s =='.' || s == 'e' || s == 'E' 
        || s == '+' || s == '-') break; 
    if (s!=' ' && s!='\t' && s!='\r' && s!='\n' && s!=',')
      return -1;  /* garbage */      
    s = PerlIO_getc(fp); /* else skip whitespace */
  }
  /* parse number */
  while (1) {
    switch (s) {
    case '0':
    case '1':
    case '2':
    case '3':
    case '4':
    case '5':
    case '6':
    case '7':
    case '8':
    case '9':
      if (aftere)
        expo = (expo*10) + (s - '0');
      else if (afterp) {
        div /= 10.0;
        f += div*(s - '0');
      } else
        f = (f*10) + (s - '0');
      break;
    case '+':
      /* ignore */
      break;
    case '-':
      if (aftere)
        esig = -1;
      else
        sig = -1;
      break;
    case 'e':
    case 'E':
      if (aftere)
        return -1;
      aftere = 1;
      break;
    case '.':
      if (afterp || aftere)
        return -1;
      afterp = 1;
      break;
    default:
      goto endread;
      break;
    }
    nread++;
    s = PerlIO_getc(fp);
  }
endread:
  f *= sig;
  for (i=0;i<expo; i++)
    f *= (esig > 0 ? 10.0 : 0.1);

  *fz = f;
  TRAILING_WHITESPACE_CHECK(s);
  return nread;
}

EOH

pp_add_exported('', 'rasc rcube');
pp_addpm(<<'EOPM');


=head2 rasc()

=for ref

Simple function to slurp in ASCII
numbers quite quickly, although error handling is marginal (to
nonexistent).

=for usage

  $pdl->rasc("filename" [,$noElements]);
      Where:
  	filename is the name of the ASCII file to read
	$noElements is the optional number of elements in the file to read.
	     (If not present, all of the file will be read to fill up $pdl)

=for example

  #  (test.num is an ascii file with 20 numbers. One number per line.)
  $in = PDL->null;
  $num = 20;
  $in->rasc('test.num',20);
  $imm = zeroes(float,20,2);
  $imm->rasc('test.num');


=cut

sub rasc {PDL->rasc(@_)}
sub PDL::rasc {
  my ($pdl, $file, $num) = @_;
  $num = -1 unless defined $num;
  open(FILE, "<$file") || barf "Can't open $file";
  $pdl->_rasc($num,'PDL::IO::Misc::FILE');
  close FILE;
  return $pdl;
}

# ----------------------------------------------------------

=head2 rcube

=for ref

Read list of files directly into a large data cube (for efficiency)

=for usage

 $cube = rcube \&reader_function, @files;

=for example

 $cube = rcube \&rfits, glob("*.fits");

This IO function allows direct reading of files into a large data cube,
Obviously one could use cat() but this is more memory efficient.

The reading function (e.g. rfits, readfraw) (passed as a reference)
and files are the arguments.

The cube is created as the same X,Y dims and datatype as the first
image specified. The Z dim is simply the number of images.

=cut

sub rcube {

    my $reader = shift;

    barf "Usage: blah" unless ref($reader) eq "CODE";

    my $k=0;
    my ($im,$cube,$tmp,$nx,$ny);
    my $nz = scalar(@_);

    for my $file (@_) {
       print "Slice ($k) - reading file $file...\n" if $PDL::verbose;
       $im = &$reader($file);
       ($nx, $ny) = dims $im;
       if ($k == 0) {
          print "Creating $nx x $ny x $nz cube...\n" if $PDL::verbose;
          $cube = $im->zeroes($im->type,$nx,$ny,$nz);
        }
        else {
          barf "Dimensions do not match for file $file!\n" if
             $im->getdim(0) != $nx or $im->getdim(1) != $ny ;

       }
       $tmp = $cube->slice(":,:,($k)");
       $tmp .= $im;
       $k++;
      }

      return $cube;
}


EOPM

# in the future this function should return a state indicating an error
# if appropriate
pp_def('_rasc',
       Pars => '[o] nums(n)',
       OtherPars => 'int num => n; char* fd',
       GenericTypes => [F],
       Code => q@ 
                 int ns, i;
                 PerlIO *fp;
                 IO *io;

                 io = GvIO(gv_fetchpv($COMP(fd),FALSE,SVt_PVIO));
                 if (!io || !(fp = IoIFP(io)))
                        croak("Can\'t figure out FP");
                 ns = $SIZE(n);
                 threadloop %{
                   for (i=0;i<ns; i++) {
                       if (getfloat(fp, &($nums(n=>i))) <= 0)
                              break;
                   }
                 %}
                @,
	Doc => 'Internal Function used by rasc. '
);



#pp_addpm(<<'EOD');

#=item bswap2( [o]x() )

#Swaps pairs of bytes in argument x()

#=cut
#EOD

defpdl(
	'bswap2',
	'x(); ',
	'',
	'
      int i;
      PDL_Short *aa; PDL_Short bb;
      PDL_Byte *a,*b;

      int n = sizeof($x()) / sizeof(PDL_Short);
      aa = (PDL_Short*) &$x();

      for(i=0;i<n; i++) {
         bb = aa[i]; a = (PDL_Byte*) (void*) (aa+i);
         b = (PDL_Byte*) &bb;
         *a = *(b+1);  *(a+1) = *b;
     }',
     "Swaps pairs of bytes in argument x()"
);


#pp_addpm(<<'EOD');
#
#=item bswap4( [o]x() )

#Swaps quads of bytes in argument x()

#=cut
#EOD

defpdl(
	'bswap4',
	'x(); ',
	'',
	'
      int i;
      PDL_Long *aa; PDL_Long bb;
      PDL_Byte *a,*b;

      int n = sizeof($x()) / sizeof(PDL_Long);
      aa = (PDL_Long*) &$x();

      for(i=0;i<n; i++) {
         bb = aa[i]; a = (PDL_Byte*) (void*) (aa+i);
         b = (PDL_Byte*) &bb;
         *a = *(b+3);  *(a+1) = *(b+2);  *(a+2) = *(b+1); *(a+3) = *b;
     }',
     "Swaps quads of bytes in argument x()"
);


#pp_addpm(<<'EOD');

#=item bswap8( [o]x() )

#Swaps octets of bytes in argument x()

#=cut
#EOD

defpdl(
	'bswap8',
	'x(); ',
	'',
	'
      int i;
      PDL_Double *aa; PDL_Double bb;
      PDL_Byte *a,*b;

      int n = sizeof($x()) / sizeof(PDL_Double);
      aa = (PDL_Double*) &$x();

      for(i=0;i<n; i++) {
         bb = aa[i]; a = (PDL_Byte*) (void*) (aa+i);
         b = (PDL_Byte*) &bb;
         *a     = *(b+7);  *(a+1) = *(b+6);  *(a+2) = *(b+5); *(a+3) = *(b+4);
         *(a+4) = *(b+3);  *(a+5) = *(b+2);  *(a+6) = *(b+1); *(a+7) = *b;
     }',
     "Swaps octets of bytes in argument x()"
);


pp_addpm({At=>Bot},<<'EOD');

=head1 AUTHOR

Copyright (C) Karl Glazebrook 1997.
All rights reserved. There is no warranty. You are allowed
to redistribute this software / documentation under certain
conditions. For details, see the file COPYING in the PDL
distribution. If this file is separated from the PDL distribution,
the copyright notice should be included in the file.


=cut

EOD

pp_done();