#!/bin/perl -w

# Module to decode GRIB.
#
# Copyright (C) W. M. Connolley 2002
# wmc@bas.ac.uk / http://www.antarctica.ac.uk/met/wmc
#
# $Author: wmc $:$Date: 2002/07/26 08:43:09 $:$Revision: 1.3 $
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# http://www.gnu.org/copyleft/gpl.html
#
# Useful methods:
#
# new - create new GRIB object
# set_message(ref-to-string) - setup the message to decode
# decode - decode
# as_text - output as a text string
#
# Thanks to for the PDB stuff beyond octet 23...

package GRIB;

require Exporter;
@ISA = qw(Exporter);
@EXPORT = qw();
@EXPORT_OK = ();

use strict;

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

sub new {

  my $self = { 
    Error_message => "",
    Warn_message => "",
    Error => 0,
    Warn => 0,
    Skipped => 0,
    DEBUG => 0,
    data => [[]],
    data_max => -1e30,
    data_min => 1e30,
    @_
  };
  bless $self;

  return $self;

};

# ----------------------------------------------
 
sub as_grid {

  my $self = shift;

  bless $self,"GRID";

#
# Conversion
#

# Nlat and lon
  $self->{nana} = $self->{nlat};
  $self->{ngng} = $self->{nlon};

# Data: 2d to 1d
  $self->{data2} = [[]];
  my $ind=0;
  for (my $i=0; $i<$self->{nlat}; $i++) {
    $self->{data2}->[$i]=[];
    for (my $j=0; $j<$self->{ngng}; $j++) {
      $self->{data2}->[$i]->[$j] = $self->{data}->[$ind++]
    }
  };

# Ref val
  $self->{rrrrrrr} = $self->{data_min};
# Plausible guess at n digits to use...
  $self->{n1} = int($self->{bits_per_value} / 3 +0.99);
  $self->{n2} = 0;
# Sign: since we subtract ref val, all are +ve
  $self->{is} = 1;
# We don't code a pole row, since we don't need to save space
  $self->{pole_row_index} = -1;
# Lon origin
  $self->{iaiaia} = $self->{zlon};

# Segment index and number (both always 1 for us)
  $self->{ntnt} = "01";
  $self->{nn} = "01";

# Center and Catalog number (hardcode 999)
  $self->{NNN} = "999";
  $self->{F1F2} = sprintf("%2.2i",$self->{centre});
# Grid geometry... not really sure about this either...
  $self->{grgr} = "99";

# Variable types
  my $var = sprintf("%3.3i",$self->{var});
  $self->{a1a1} = substr($var,-2); 
  $self->{a2a2} = "00";

# ref lev
  $self->{h1h1h1h1} = sprintf("%4.4i",$self->{level_value});

# Date stuff
  $self->{yy} = sprintf("%2.2i",$self->{year} % 100);
  $self->{mm} = sprintf("%2.2i",$self->{month});
  $self->{dd} = sprintf("%2.2i",$self->{day});
  $self->{hh} = sprintf("%2.2i",$self->{hour});
  $self->{mm1} =sprintf("%2.2i",$self->{min});

# Other stuff
  $self->{q1} = 2;			# 2 = no spaces between data
  $self->{q2} = 0;			# 0 = all groups included?
  $self->{np} = 1;			# grid pts per data group? Not used?

# Note: this is now a call in the GRID package
  return $self->as_grid()

}

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

sub as_grib {

  my $self = shift;

# All the encode routines default their values from "self"

  my $pdb = $self->encode_pdb();
  my $gdb = $self->encode_gdb() or return 0;
  my $db = $self->encode_db();
  my $end = $self->encode_end();

  return "GRIB".$pdb.$gdb.$db.$end

}

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

sub encode_gdb {

# Expect all arguments passed in as "name" => value pairs
# or defaulted from self.

  my $self = shift;

  my %opts = (
               grid_type => 0
             , res_flag => 0x80
             , @_
             );

# Populate opts with self where undefined... 
  for (keys %$self) { if (!$opts{$_}) { $opts{$_} = $$self{$_} } };

  if ($self->{DEBUG}) { print " --- begin encode GDB --- \n" };

  if ($opts{grid_type} != 0) { $self->error("Grid type must be lat-lon (0)") };

  my $gdb = n2three(32); 			# length
  $gdb .= "\0\0";				# 2 bytes of junk
  $gdb .= pack("C",$opts{grid_type});
  $gdb .= n2two($opts{nlon});
  $gdb .= n2two($opts{nlat});
  $gdb .= sn2three($opts{zlat}*1000);
  $gdb .= sn2three($opts{zlon}*1000);
  $gdb .= pack("C",$opts{res_flag});
  my $xlat = $opts{zlat} + $opts{dlat}*($opts{nlat}-1);
  my $xlon = $opts{zlon} + $opts{dlon}*($opts{nlon}-1); if ($xlon > 360) { $xlon -= 360 };
  $gdb .= sn2three($xlat*1000);
  $gdb .= sn2three($xlon*1000);
  my $scan_mode = 0;
  if ($opts{dlat} > 0) { $opts{dlat} *= -1; $scan_mode |= 0x40 };
  if ($opts{dlon} < 0) { $opts{dlon} *= -1; $scan_mode |= 0x80 };
  $gdb .= n2two(abs($opts{dlat}*1000));
  $gdb .= n2two(abs($opts{dlon}*1000));
  $gdb .= pack("C",$scan_mode);
  $gdb .= "\0\0\0\0";

  return $gdb

};

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

sub encode_pdb {

  my $self = shift;

  my %opts = (
               hour => 0
             , min => 0
             , @_
             );

# Populate opts with self...
  for (keys %$self) { if (!$opts{$_}) { $opts{$_} = $$self{$_} } };

  my $pdb;

  if ($self->{DEBUG}) { print " --- begin encode PDB ---\n" };

# Octets:
#   1-3: length1
#   4:   edition
#        if edition = 0,			if edition = 1,
#        length1 = of pdb			length1 = of whole message 
#   5:   centre...				5-7: length (of pdb)
#                                               8: don't know
#						5+4: centre...

  $pdb = n2three($self->{pdb_n});
  $pdb .= pack("C",1);				# We write edition 1 GRIB
  $pdb .= n2three(28);
  $pdb .= "\0";
  $pdb .= pack("C",$self->{centre});
  $pdb .= pack("C",$self->{modelid});
  $pdb .= pack("C",$self->{grid});
  $pdb .= pack("C",$self->{flag});
#  if (!($self->{flag} & 128)) { return $self->error("encode_pdb: there is no gdb ($self->{flag})") };
#  if ($self->{flag} & 64) { return $self->error("encode_pdb: there is a bdb (bit-map-block) ($self->{flag})") };
  $pdb .= pack("C",$opts{var});
  $pdb .= pack("C",$opts{level_type});
  $pdb .= n2two($opts{level_value});
  $pdb .= pack("C5",$opts{year}, $opts{month}, $opts{day}, $opts{hour}, $opts{min});
  $pdb .= pack("C",$self->{tr});
  $pdb .= pack("C2",$self->{time1}, $self->{time2});
  $pdb .= pack("C",$self->{tri});
  $pdb .= n2two($self->{navg});
  $pdb .= "\0";

  if ($self->{DEBUG}) {
    print "So: pdb length is: ",length($pdb),"\n"
  };

  return $pdb."\0\0\0\0"

};

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

sub encode_db {
  my $self = shift;

  my $db;
  my $x;
  my $data = $self->{data};

  if ($self->{DEBUG}) { print " --- Begin encode DataB --- \n" };

  $db = n2three($self->{db_n});					# 1-3 : length
  $db .= pack("C",$self->{db_flag});				# 4   : first 4 bits: flag; last 4: unused
  my $unused = $self->{db_flag} & 0x0F;
  $db .= sn2two($self->{sf1});					# 5-6 : scale factor (E)
  $db .= write_ref_val($self->{ref_value});			# 7-10: ref val
  $db .= pack("C",$self->{bits_per_value});			# 11  : bits-per-value

# Calculate expected length
  my $n_points = $self->{nlat} * $self->{nlon};
  my $n_data_bytes = int ( $self->{bits_per_value} * $n_points / 8 +0.999 );
  my $unused1 = $n_data_bytes*8 - $n_points*$self->{bits_per_value} ;

  if ($self->{DEBUG}) { print "unused: $unused $unused1\n" };

# Decode and output data
  my $dbb; my $dbbb;
  for (my $i= 0; $i<$n_points; $i++) { 
    $x=($$data[$i] - $self->{ref_value})/$self->{sf};
    $dbbb = substr(unpack ( "B32", pack("N",$x)),-$self->{bits_per_value});
    $dbb .= $dbbb;
if ($i > $n_points - 20) { print "..$dbbb ($x).. " };
  };
  $db .= pack("B*",$dbb);

print "\n -- " . unpack("B*",substr($db,-$self->{bits_per_value})) . " -- \n";

  return $db

};

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

sub encode_end {

  return "7777"

};

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

sub as_text {

  my $self = shift;

  my %opts = (
               data => "full"
             , @_
             );

  my $text;

  $text = "";

  $text.= "Grid: $self->{grid_type} $self->{zlat} to $self->{xlat} by $self->{dlat} ($self->{nlat}); $self->{zlon} to $self->{xlon} by $self->{dlon} ($self->{nlon}) $self->{res_flag} $self->{scan_mode}\n";

  $text.= "Ed: $self->{edition} Cen: $self->{centre} Mod: $self->{modelid} Grid: $self->{grid} Var: $self->{var} LevT: $self->{level_type} LevV: $self->{level_value}\n";

  $text.="$self->{year}/$self->{month}/$self->{day}:$self->{hour}-$self->{min} Tr: $self->{tr} T1: $self->{time1} T2: $self->{time2} Tri: $self->{tri} Navg: $self->{navg}\n";

  $text.="$self->{nlat} $self->{nlon}\n";
  if ($opts{data} eq "full") {
    my $rdata=$self->{data}; for (@$rdata) { $text.="$_\n" }
  } elsif ($opts{data} eq "range") {
    $text .= "Data range: [$self->{data_min},$self->{data_max}]\n" 
  };
  return $text;

};

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

sub Warn {

  my $self = shift;
  my $message = shift;

  $self->{Warn_message}=$message;
  $self->{Warn}=1;

  return(0)

};

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

sub error {

  my $self = shift;
  my $message = shift;

  $self->{Error_message}=$message;
  $self->{Error}=1;

  return(0)

};

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

sub set_message {

  my $self = shift;
  
  $self->{message} = shift;

};

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

sub check5 {

  my $self = shift;
  my $message = shift;
  my $group;

  while($group = shift) {
    if (length($group) != 5) { 
      return $self->error("[check5] $message group <$group>")
    }
  };

  return 1;

};

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

sub decode {

  my $self = shift;

  my $rmessage = $self->{message};
  my $message = $$rmessage;

  my $extra;				# Extra stuff (which we don't expect)
  my $junk;				# Guess
  my $line="";                  	# One line of message
  my $in;
  my $ind = 0;

  if ($self->{DEBUG}) { print "starting GRIB decode (".length($message).")\n" };

# We may have to skip padding
  while (defined($in=substr($message,$ind,4)) and $in ne "GRIB") { $ind++ };
  $self->{Skipped}=$ind;

# Adjust $ind to point to just after GRIB
  $ind += 4;

  if ($in ne "GRIB") { return $self->error("No GRIB found") };

  if ($self->{DEBUG}) { print "found GRIB after skipping $self->{Skipped} bytes\n" };

  $self->decode_pdb($message,\$ind) or return 0;
  $self->decode_gdb($message,\$ind) or return 0;
  $self->decode_db($message,\$ind) or return 0;
  $self->decode_end($message,\$ind) or return 0;

};

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

sub two2n {

  my $self = shift; 		
  my $in;
  if (ref($self)) {
    $in = shift
  } else {
    $in = $self
  };

  return vec($in,0,8)*256+vec($in,1,8);

};

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

sub stwo2n {

  my $self = shift; 		
  my $in;
  if (ref($self)) {
    $in = shift
  } else {
    $in = $self
  };

  my $sign = (vec($in,0,8) & 128) ? -1 : 1;

  return $sign*((vec($in,0,8)&0x7F)*256+vec($in,1,8));

};

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

sub n2two {

  my $self = shift; 		
  my $n;
  if (ref($self)) {
    $n = shift
  } else {
    $n = $self
  };

#  if ($n > 256*256) { print "too big!\n" };

  my $n1 = int($n/256);
  my $n2 = $n - $n1*256;

  return pack("C2",$n1,$n2)

};

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

sub sn2two {

  my $self = shift; 		
  my ($n,$sgn);
  if (ref($self)) {
    $n = shift
  } else {
    $n = $self
  };

#  if (abs($n) > 128*256) { print "too big!\n" };

  if ($n < 0) { $sgn = 1; $n = -$n } else { $sgn = 0 };
  my $n1 = int($n/256);
  my $n2 = $n - $n1*256;

  return pack("C2",$n1 | 0x80*$sgn,$n2)

};

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

sub sn2three {

  my $self = shift; 		
  my ($n,$sgn);
  if (ref($self)) {
    $n = shift
  } else {
    $n = $self
  };

#  if (abs($n) > 128*256*256) { print "too big!\n" };

  if ($n < 0) { $sgn = 1; $n = -$n } else { $sgn = 0 };
  my $n1 = int(($n/256)/256);
  my $n2 = int(($n-$n1*256*256)/256);
  my $n3 = $n - $n1*256*256 - $n2*256;

  return pack("C3",$n1 | 0x80*$sgn,$n2,$n3)

};

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

sub n2three {

  my $self = shift; 		
  my $n;
  if (ref($self)) {
    $n = shift
  } else {
    $n = $self
  };

#  if ($n > 256*256*256) { print "too big!\n" };

  my $n1 = int(($n/256)/256);
  my $n2 = int(($n-$n1*256*256)/256);
  my $n3 = $n - $n1*256*256 - $n2*256;

  return pack("C3",$n1,$n2,$n3)

};

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

sub three2n {

  my $self = shift; 		
  my $in;
  if (ref($self)) {
    $in = shift
  } else {
    $in = $self
  };

  return (vec($in,0,8)*256+vec($in,1,8))*256+vec($in,2,8);

};

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

sub sthree2n {

  my $self = shift; 		
  my $in;
  if (ref($self)) {
    $in = shift
  } else {
    $in = $self
  };

  my $sign = (vec($in,0,8) & 128) ? -1 : 1;

  return $sign*(((vec($in,0,8)&0x7F)*256+vec($in,1,8))*256+vec($in,2,8));

};

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

sub write_ref_val {

  my $self = shift; 		
  my $in;
  if (ref($self)) {
    $in = shift
  } else {
    $in = $self
  };

  my $sgn;
  if ($in < 0) { $in = - $in; $sgn = 1 } else { $sgn = 0 };
  my $exp = int(64 + log_n(16,$in) + 0.99999);
  $in = $in / 2**-24 / 16**($exp-64);

# For now, fake it...
  return pack("C",$exp | 0x80*$sgn).n2three($in)

};

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

sub get_ref_val {

  my $self = shift; 		
  my $in;
  if (ref($self)) {
    $in = shift
  } else {
    $in = $self
  };

  my $ref_frac=three2n(substr($in,1,3));
  my $ref_mant=vec($in,0,8);
  my $ref_sign = ($ref_mant & 0x80) ? -1 : 1;
  $ref_mant=($ref_mant & 0x7F) - 64;

  my $ref_val = $ref_sign*$ref_frac*2**(-24)*16**($ref_mant);
  # print "(".display_raw($in).") ref_sign: $ref_sign; ref_frac: $ref_frac; ref_mant: $ref_mant [$ref_val]\n" };

  return $ref_val 

};

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

sub decode_gdb {

  my $self = shift;
  my $message = shift;			# The message to decode
  my $rind = shift;			# Pointer into the message
  my $ind = $$rind;

  my $in;

# Save ind. We shall set it to this, plus block length, at end
  my $ind_save=$ind;

  $self->{gdb_n} = three2n(substr($message,$ind,3)); 

  if ($self->{DEBUG}) { 
    print " --- begin GDB --- ($ind: ".display_raw(substr($message,$ind,$self->{gdb_n}))." [" 
                                      .display_raw(substr($message,$ind+$self->{gdb_n},3))."])\n" 
  };
  $ind += 3;

# Skip 2 bytes of junk
  $ind +=2;

#
# Interpret GDB info
#
$self->{grid_type}=vec(substr($message,$ind++,1),0,8);
if ($self->{grid_type} == 0) { 
#
# Regular lat-lon grid
#
  $self->{nlon}=two2n(substr($message,$ind,2)); $ind+=2;
  $self->{nlat}=two2n(substr($message,$ind,2)); $ind+=2;
  my $totalpoints=$self->{nlat}*$self->{nlon};
  $self->{zlat}=sthree2n(substr($message,$ind,3))/1000; $ind+=3;
  $self->{zlon}=sthree2n(substr($message,$ind,3))/1000; $ind+=3;
  $self->{res_flag}=vec(substr($message,$ind++,1),0,8);
  $self->{xlat}=sthree2n(substr($message,$ind,3))/1000; $ind+=3;
  $self->{xlon}=sthree2n(substr($message,$ind,3))/1000; $ind+=3;
  $self->{dlat}=two2n(substr($message,$ind,2))/1000; $ind+=2;
  $self->{dlon}=two2n(substr($message,$ind,2))/1000; $ind+=2;
  $self->{scan_mode}=vec(substr($message,$ind++,1),0,8);
  if ($self->{scan_mode} & 128)    { $self->{dlon} *= -1 };
  if (!($self->{scan_mode} &  64)) { $self->{dlat} *= -1 };
  if ($self->{DEBUG}) { 
    print "Scan mode: $self->{scan_mode}\n";
    print "Grid type is regular lat/regular long: OK\n";
    print "Nlat: $self->{nlat}, Nlon: $self->{nlon}\n";
    print "Origin (lat, lon) [may not be coded]: $self->{zlat}, $self->{zlon}\n";
    print "Xtreme (lat, lon) [may not be coded]: $self->{xlat}, $self->{xlon}\n";
    print "Delta (lat, lon) [may not be coded]: $self->{dlat}, $self->{dlon}\n"; 
  };

} else {

  die "Unknown grid type: $self->{grid_type}"

};

  $$rind = $ind_save+$self->{gdb_n};

};

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

sub decode_pdb {

  my $self = shift;
  my $message = shift;			# The message to decode
  my $rind = shift;			# Pointer into the message
  my $ind = $$rind;

  my ($in,$n);

  if ($self->{DEBUG}) { print " --- begin PDB --- ($ind: ".display_raw(substr($message,$ind,28)).")\n" };

#
# Read PDB block length. 
#
# Octets:
#   1-3: length1
#   4:   edition
#        if edition = 0,			if edition = 1,
#        length1 = of pdb			length1 = of whole message 
#   5:   centre...				5-7: length (of pdb)
#                                               8: don't know
#						5+4: centre...
#
  $self->{pdb_n} = $self->three2n(substr($message,$ind,3)); $ind += 3;
  $self->{edition} = vec($message,$ind++,8);
  if ($self->{edition} == 1) { 
    $n = three2n(substr($message,$ind,3));
    if ($self->{DEBUG}) { print "(edition=1) pdb length: $n\n" };
    if ($n != 28) { return $self->error("pdb length a mystery: ($n)") };
    $ind +=4
  } elsif ($self->{edition} == 0) {
    if ($self->{pdb_n} != 24) { return $self->error("pdb length a mystery: ($n)") };
    $n = $self->{pdb_n};
  } else {
    return $self->error("decode_pdb: edition not 0 or 1: <$self->{edition}>")
  };
# Note: if edition=1, the pdb length seems to be considereed as excluding bytes 1-4. Hmm...
  my $ind_save=$ind-4;

# Centre - octet 5
  $self->{centre}=vec($message,$ind++,8);
# Model id - octet 6
  $self->{modelid}=vec($message,$ind++,8);
# Grid - octet 7
  $self->{grid}=vec($message,$ind++,8);
# Flag - octet 8
  $self->{flag}=vec($message,$ind++,8);
  if (!($self->{flag} & 128)) { return $self->error("decode_pdb: there is no gdb ($self->{flag})") };
  if ($self->{flag} & 64) { return $self->error("decode_pdb: there is a bdb (bit-map-block) ($self->{flag})") };
# Variable type - octet 9
  $self->{var}=vec($message,$ind++,8);
# Level type - octet 10
  $self->{level_type}=vec($message,$ind++,8);
# Level value - octets 11,12
  $self->{level_value}=$self->two2n(substr($message,$ind,2)); $ind += 2;

# Year, Month, Day, Hour, Min 
  $self->{year}=vec($message,$ind++,8);
  $self->{month}=vec($message,$ind++,8);
  $self->{day}=vec($message,$ind++,8);
  $self->{hour}=vec($message,$ind++,8);
  $self->{min}=vec($message,$ind++,8);

# Time range
  $self->{tr}=vec($message,$ind++,8);
# Time 1
  $self->{time1}=vec($message,$ind++,8);
# Time 2
  $self->{time2}=vec($message,$ind++,8);
# Time range indicator
  $self->{tri}=vec($message,$ind++,8);
# N averaged
  $self->{navg}=$self->two2n(substr($message,$ind,2)); $ind += 2;
# Last byte missing
# octet 24 Number missing from averages or accumulations
  $self->{nmaa}=vec($message,$ind++,8);
# octet 25 Century of reference time of data
  $self->{crtd}=vec($message,$ind++,8);
# octet 26 Sub-center identification
  $self->{cci}=vec($message,$ind++,8);
# octet 27-28 Units decimal scale factor (D)
  $self->{udsf}=$self->two2n(substr($message,$ind,2)); $ind += 2;
  $self->{udsf1}=10**$self->{udsf};

  if ($self->{DEBUG}) {
    if ($self->{time1} == 0) { print "Analysis\n" };
    print "Var, Level type, value, date, udsf: " .
    "$self->{var}, $self->{level_type}, $self->{level_value} " .
    "$self->{year}/$self->{month}/$self->{day}:$self->{hour}.$self->{min}, $self->{udsf1} ($self->{udsf})\n";
  };

  $$rind = $ind_save+$n;

};

sub decode_db {
  my $self = shift;
  my $message = shift;			# The message to decode
  my $rind = shift;			# Pointer into the message
  my $ind = $$rind;

  my $x;

#
# Read DataB
#
if ($self->{DEBUG}) { print " --- Begin DataB --- ($ind: ".display_raw(substr($message,$ind,21)).")\n" };
my $ind_save=$ind;

$self->{db_n}=three2n(substr($message,$ind,3)); $ind+=3;

# Decode DataB

# Octet 4 - flag (code table 11) and unused bits
$self->{db_flag}=vec(substr($message,$ind++,1),0,8);
my $unused = $self->{db_flag} & 0x0F;

# Octets 5-6 - scale factor (E)
$self->{sf1}=stwo2n(substr($message,$ind,2));
$self->{sf}=2**($self->{sf1}); $ind+=2;

# Octets 7-10 - ref value
$self->{ref_value}=get_ref_val(substr($message,$ind,4));  $ind+=4;

# Octet 11 - bits per packed value
$self->{bits_per_value}=vec(substr($message,$ind++,1),0,8);

if ($self->{DEBUG}) {
  print "Length of DataB block: $self->{db_n}\n";
  print "Scale factor: $self->{sf} ($self->{sf1})\n";
  print "Ref value: $self->{ref_value}\n";
  print "Bits per packed value: $self->{bits_per_value}\n"
};

# Calculate expected length
my $n_points = $self->{nlat} * $self->{nlon};
my $n_data_bytes = $self->{bits_per_value} * $n_points / 8;
my $n1 = $n_data_bytes + $unused/8 + 11;
if ($self->{db_n} == $n1) {
  if ($self->{DEBUG}) { print "Length of DataB is OK (= $n_data_bytes + $unused/8 + 11)\n" }
} else {
  print "DataB length ($self->{db_n}) not consistent with data it is supposed to contain ($n1) ...Going with DataB version\n" 
};

# Decode and output data
my $unpacked=unpack("B*",substr($message,$ind,int($n_data_bytes+0.999)));
my @data;
for (my $i=0; $i<$n_points; $i++) { 
  $x=unpack("N",pack("B32",substr("0"x32 . substr($unpacked,$i*$self->{bits_per_value},$self->{bits_per_value}),-32)));
  $x=$x*$self->{sf}+$self->{ref_value}; 
# Additionally, scale by the udsf
  $x/=$self->{udsf1};
  $data[$i]=$x;
  if ($i == 0) {
    $self->{data_min} = $x;
    $self->{data_max} = $x
  } else {
    if ($self->{data_min} > $x) { $self->{data_min} = $x };
    if ($self->{data_max} < $x) { $self->{data_max} = $x }
  }
};
$self->{data}=\@data;

$$rind = $ind_save + $self->{db_n};

};

sub decode_end {

  my $self = shift;
  my $message = shift;			# The message to decode
  my $rind = shift;			# Pointer into the message

#
# Read "End"
#
  if ($self->{DEBUG}) { print " --- begin End --- ($$rind: ".display_raw(substr($message,$$rind,4)).")\n" };

  if (substr($message,$$rind,4) eq "7777") { 
    if ($self->{DEBUG}) { print "Read End (7777) OK\n" }
  } else { 
    print "Didn't find 7777 at End\n"
  };

  $$rind+=4;

};

# Utility routine to output the next few values of the message, rather
# helpful if we're lost...

sub display_raw {

  my $self = shift;
  my $message;
  if (ref($self)) { $message = shift } else { $message = $self };
  my $result = "";

  for (unpack("C*",$message)) { $result .= "$_ " }
  chop $result;

  return $result;

};

# z = log_n(base,x)

sub log_n {

  my $a = log(shift);
  return log(shift)/$a

};

1;