#!/bin/perl -w

# wunconv8 - re-insert modified pp-fields into unified model ancillary fields, 8-byte version
#
# use:	  - wunconv.pl ancil-file pp-file output-ancil-file
#
# kludge: can use "arbitrary_skip" in way that mirrors wconv.pl
#
# inputs:
#		ancil-file: the original ancil file, probably the one you used wconv8 on (8-byte)
#               pp-file: source for the modified fields (4-byte)
#
# output:
#               output-ancil-file: ancil-file, with the data segments replaced from pp-file (8-byte)
#
# diags to: stdout
#
# method:
#	1. read in all of ancillary file (first argument)
#	2. read in pp-field(s) (second arg)
#	3. diagnostic info to wunconv.info
#	4. interpret header (mostly for info + to make sure we know what we're doing)
#	5. write pp-headers into lookup tables and data into data block 
#          (check that the headers are the same; warn if not)
#	6. write out result to stdout
#
# reference:
#	UMDP f3 (format of files), R. S. Bell, v 22 27/1/95
#       Code sourced from "wunconv2.pl" (~wmc/bin) 2/11/2001

$N_ins=0;
$MAX_warn=5;
$Verbose=1;
$arbitrary_skip=0;
%FIXHD10=qw(0 single  1 timeseries  2 repeating);

# Command-line setup
eval "\$$1=\$2" while $ARGV[0] =~ /^(\w+)=(.*)/ && shift;

if (scalar @ARGV != 3) {
  die "use: wunconv.pl ancil-file pp-file output-ancil-file"
};

# 1. readin ancil
#
undef $/;
$ANCIL=shift;
open(IN,$ANCIL) or die "Failed to open ancillary file <$ANCIL> (first arg)\n"; $IN=<IN>; close(IN);
$NL=length($IN);

# 2. readin pp(s)
#
$PPF=shift;
open(IN,$PPF) or die "Failed to open pp file <$PPF> (second arg)\n"; $PIN=<IN>; close(IN);
$NLP=length($PIN);

# 3. Output
#
$OUTFILE=shift;
if (!$OUTFILE) { die "No output file specified" };
open(OUT,"> $OUTFILE") or die "Can't open output to $OUTFILE";

# 4. 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++) { 
  if ($H[$I] == -32768) {
    print sprintf("%4d%7s",$I+1,"null")
  } else {
    print sprintf("%4d%7d",$I+1,$H[$I])
  }; 
  if (($I+1) % 8 == 0) { print "\n" }; 
} };
if (($H[0] == -32768) or ($H[0] == 15)) { print "H0=15 or null; OK\n" } else { warn "H0 <> 15: $H[0]\n" };
print "Vn: ",$H[11]/100.,"\n";
print "t0 (start) Y/M/D:H/M/S:D $H[20]/$H[21]/$H[22]:$H[23]/$H[24]/$H[25]:$H[26]\n";
print "t1 (end)   Y/M/D:H/M/S:D $H[27]/$H[28]/$H[29]:$H[30]/$H[31]/$H[32]:$H[33]\n";
print "t2 (skip)  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
print "We seem to have ${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";
if ($NF == $H[151]) { print "N Field types matches 2nd dim of lookup; good ($NF)\n" 
} else { warn "N Field types ($NF) does not match 2nd dim of lookup ($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" };

# Check expected pp-length
$LREC=$NX*$NY;
$PPL=4*(4+64+$LREC);
$LREC1=unpack("i",substr($PIN,15*4,4));
$NY1=unpack("i",substr($PIN,18*4,4));
$NX1=unpack("i",substr($PIN,19*4,4));
$PPL1=($LREC1 +64 +4)*4;
if ($PPL != $PPL1) { print "  these PP-fields are not the same as the ancil. Ah well. $PPL. $PPL1\n" };
$PPL=$PPL1; $LREC=$LREC1;
$NX=$NX1; $NY=$NY1;
$NPP=$NLP/$PPL;
print "Expected pp-field length is: $PPL; that would mean we read in ",$NPP," fields\n";

# 5. Insert fields
#
if ($N_ins) { $N=$N_ins } else { $N=$H[151] };
for ($F=0; $F<$N; $F++) {
  if ($F >= $NPP) {
    print "Run out of pp-fields to insert; stopping\n"; last
  };
  print "Inserting field number $F / Data from  : ",$H[159]+$F*$LREC," to ",$H[159]+(1+$F)*$LREC-1,"\n";

# For this version, I'm a bit confused about 8-byte integers, so
# we're not going to allow modification of the headers, just the data fields

# This is the current pp-field
  $PP=substr($PIN,$PPL*$F,$PPL);
# Get header and data portions
  $PPH=substr($PP,4,64*4); $PPD=substr($PP,4*(3+64),4*$LREC);
# Unpack the 4-byte reals...
  @data1=unpack "f*",$PPD;
# And repack them as 8-byte
  $data=pack "d*",@data1;
# Insert data into ancillary
  substr($IN,($H[159]+$F*$LREC-1+$F*$arbitrary_skip)*8,$LREC*8)=$data;

};
print "(Length (words) of file read in: $NL (",$NL/4,"))\n";

# 6. Write out result
print OUT $IN;



sub wwarn {

  $MAX_warn--;
  if ($MAX_warn >= 0) {
    print @_
  };
  if ($MAX_warn == 0) { print "...and this is the last warning! (set MAX_warn=nn on the command line if you wnat more)\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)."] ($bmdi)";

}