pp_addpm({At=>Top},<<'EOD');
=head1 NAME
PDL::IO::Grib::Wgrib - shamelessly stolen code
=head1 SYNOPSIS
Helper routines for PDL::IO::Grib
=head1 DESCRIPTION
This module is for miscellaneous PP-defined utility routines for
the PDL::IO::Grib module.
=head1 AUTHOR
Copyright (C) 2000 James P. Edwards
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_addhdr('
void BDS_unpack(float *flt, unsigned char *bits, unsigned char *bitmap,
int n_bits, int n, double ref, double scale);
int flt2ibm(float x, unsigned char *ibm);');
pp_def(
'BDSunpack',
Pars => 'byte bits(n1);
byte bitmap(n2);
int bitcnt();
double ref();
double scale();
float [o] bdsdata(n3); ',
GenericTypes => [B,F],
Code => 'BDS_unpack($P(bdsdata), $P(bits), $P(bitmap), $bitcnt(), $SIZE(n3),
$ref(), $scale());'
);
pp_addpm('
sub get_ref_val{
my($d) = shift;
my $do = PDL->zeroes((new PDL::Type(float)),$d->nelem/4);
decode_ref_val($d,$do);
return($do);
}
');
pp_def(
'decode_ref_val',
Pars => 'byte datain(m);
float [o] dataout();',
RedoDimsCode => '$SIZE(m)=4;',
GenericTypes => [B,F],
Code => '
int positive, power;
unsigned int abspower;
long int mant;
double value, exp;
int m1,m2,m3,m4;
m1 = 0;
m2=m1+1;
m3=m2+1;
m4=m3+1;
positive = ($datain(m=>m1) & 0x80) == 0;
mant = ($datain(m=>m2) << 16) + ($datain(m=>m3) << 8) + $datain(m=>m4);
power = (int) ($datain(m=>m1) & 0x7f) - 64;
abspower = power > 0 ? power : -power;
/* calc exp */
exp = 16.0;
$dataout() = 1.0;
while (abspower) {
if (abspower & 1) {
$dataout() *= exp;
}
exp = exp * exp;
abspower >>= 1;
}
if (power < 0) $dataout() = 1.0 / $dataout();
$dataout() = $dataout() * mant / 16777216.0;
if (positive == 0) $dataout() = -$dataout();
'
);
pp_def(
'encode_ref_val',
Pars => 'float datain();
byte [o] dataout(n);',
GenericTypes => [B,F],
RedoDimsCode => '$SIZE(n)=4;',
Code => 'flt2ibm($datain(),$P(dataout));'
);
pp_addxs('','
int
flt2ibm(x,ibm)
float x
unsigned char *ibm
');
pp_done();