#!/bin/perl -w

#
# wconv8 - take 8-byte UKMO unified model ancillary fields and convert them to pp-fields
#
# Use: wconv8.pl infile outfile
#
# method:
#	1. read in all of file (first argument or stdin)
#	2. diagnostic info to wconv.info
#	3. interpret header (mostly for info + to make sure we know what we're doing)
#	4. read pp-headers from lookup tables and data from data block and write to stdout
#
# kludges:
#   You can call as wconv8.pl arbitrary_skip=nnn (for example, 480) to add nnn*field-n
#
# reference:
#	UMDP f3 (format of files), R. S. Bell, v 22 27/1/95

# 0. Setup
$Verbose=1;
$arbitrary_skip=0;
# Command-line setup
eval "\$$1=\$2" while $ARGV[0] =~ /^(\w+)=(.*)/ && shift;
if (scalar @ARGV != 2) { die "use: wconv8.pl infile outfile" };

# 1. readin
#
$INFILE=shift;
$OUTFILE=shift;
undef $/;
open IN,$INFILE or die "Can't open $INFILE"; $IN=<IN>; close IN;
$NL=length($IN);

# 2. Diags to stdout, output to OUTFILE
#
if (!$OUTFILE) { die "No output specified" };
open OUT,"> $OUTFILE" or die "Can't open $OUTFILE";

# 3. Info
#
@H1=unpack("i512",substr($IN,0,256*8));
for ($I=0; $I<256; $I++) { $H[$I]=$H1[$I*2]*(2**32)+$H1[$I*2+1] };
if ($Verbose > 1) { for ($I=0; $I<256; $I++) { print sprintf("%4d%15d",$I+1,$H[$I]); if (($I+1) % 4 == 0) { print "\n" }; }; };
if ($H[0] == 15 or $H[0] == -32768) { print "H0=$H[0] (= 15 or -32768); OK\n" } else { warn "H0 <> 15 or -32768: is: $H[0]\n" };
print "Vn: ",$H[11]/100.,"\n";
print "(t0) Y/M/D:H/M/S:D $H[20]/$H[21]/$H[22]:$H[23]/$H[24]/$H[25]:$H[26]\n";
print "(t1) Y/M/D:H/M/S:D $H[27]/$H[28]/$H[29]:$H[30]/$H[31]/$H[32]:$H[33]\n";
print "(t2) Y/M/D:H/M/S:D $H[34]/$H[35]/$H[36]:$H[37]/$H[38]/$H[39]:$H[40]\n";
print "Global (0): $H[3], Type (4=ancil): $H[4]\n";
print "Start, len of ints: $H[99],$H[100]\n";
@INTS1=unpack("i*",substr($IN,$H[99]*8-8,$H[100]*8));
for ($I=0; $I<$H[100]; $I++) { $INTS[$I]=$INTS1[$I*2]*(2**32)+$INTS1[$I*2+1] };
if ($Verbose > 1) { print "  ints: ",join(",",@INTS),"\n" };
$NX=$INTS[5]; $NY=$INTS[6];		# Ancil specific
$NF=$INTS[14]*$H[151];				# Ancil specific
# Probably should be *H[151]... see wunconv.pl
print "  ${NX}x${NY}x$NF\n";

print "Start (StartI+LenI), len of real: $H[104] (",$H[99]+$H[100],"),$H[105]\n";
print "Start (StartR+lenR), len of Lookup: $H[149] (",$H[104]+$H[105],"), \[$H[150],$H[151]\]\n";
print "Start (StartR+lenR), len of data: $H[159] (",$H[149]+$H[150]*$H[151],"),$H[160]\n";
if ($H[160] == $NX*$NY*$NF) { 
  print "Length of data matches NX*NY*NF; good\n" 
} else { 
  warn "Length of data ($H[160]) != NX*NY*NF (",$NX*$NY*$NF,")\n" 
};

# 4. Output fields
#
$Tot=0;
for ($F=0; $F<$H[151]; $F++) {
  print "F$F";
  print " Head: ",$H[149]+$F*$H[150],"-",$H[149]+(1+$F)*$H[150]-1;
  print " Data: ",$H[159]+$F*$NX*$NY,"-",$H[159]+(1+$F)*$NX*$NY-1;
  $h0=($H[149]+$F*$H[150]-1)*8;
#  print " Lbuser(1-7): ", join ",",unpack("l7",substr($IN,$h0+38*8,7*8))," (".$F*$NX*$NY.")\n";
  @lb=unpack("i14",substr($IN,$h0+38*8,7*8));
  for ($I=0; $I<7; $I++) { $lb1[$I]=$lb[$I*2]*(2**32)+$lb[$I*2+1] };
  print " Lbuser(1-7): ", join ",",@lb1," (".$F*$NX*$NY.")\n";
  print OUT pack "i", $H[150]*4 or die "Failed to write Header-1"; $Tot+=4;

  if ($H[150] != 64) { die "Header length is ne 64. Oh no! I just can't cope. I'm going to top myself... goodbye...\n" };

# Extract the 8-byte integer portion of the pp-field header and
# pack as 4-byte. "l" doesn't seem to work.
  @IH=unpack("i*",substr($IN,($H[149]+$F*$H[150]-1)*8,45*4*2));
  for ($I=0; $I<45; $I++) { $IH1[$I]=$IH[$I*2]*(2**32)+$IH[$I*2+1] };
  $IH2=pack "i*", @IH1;
  print OUT $IH2 or die "Failed to write Integer Header"; $Tot+=45*4;

# Now deal with the real portion. This is easier, as "f" and "d" work.
  $RH=pack "f*", unpack "d*", substr($IN,($H[149]+$F*$H[150]-1+45)*8,(64-45)*4*2);
  print OUT $RH or die "Failed to write Real Header"; $Tot+=(64-45)*4;

  print OUT pack "i", $H[150]*4 or die "Failed to write Header-2"; $Tot+=4;
  print OUT pack "i", $NX*$NY*4 or die "Failed to write data header"; $Tot+=4;
  $p=($H[159]+$F*($NX*$NY+$arbitrary_skip)-1)*8;
  $l=$NX*$NY*8;
  @data=unpack "d*",substr($IN,$p,$l);
  $data1=pack "f*",@data;
  print "data range: ".makerange(\@data,unpack("d",substr($IN,$h0+(46+16)*8,8))).", [0,0]: $data[0], data from: $p to ",$p+$l,"\n";
  print OUT $data1 or die "Failed to write data"; $Tot+=$NX*$NY*4;
  if (length($IN) < $p+$l-1) {
    print "\n***Warning! Not enough data available...\n***You asked for $p to ".($p+$l)." and I only had up to ".length($IN)."\n\n";
  };
  print OUT pack "i", $NX*$NY*4 or die "Failed to write tail"; $Tot+=4;
};
print "(Length (words) of file read in: $NL (",$NL/8,")). Written: $Tot bytes in $H[151] fields\n";


sub makerange {

  ($rdata,$bmdi)=@_;

  $begun=0;
  for (@$rdata) {
    if (abs($_-$bmdi) > abs($bmdi)/100) {
      if (!$begun) { $begun=1; $mn=$_; $mx=$_ };
      if ($_>$mx) { $mx=$_ };
      if ($_<$mn) { $mn=$_ }
    }
  };
  return "[".sprintf("%.2f",$mn).",".sprintf("%.2f",$mx)."]";

}