#!/bin/perl -w

# wunconv - re-insert modified pp-fields into unified model ancillary fields
#
# use:	  - wunconv.pl ancil-file pp-file output-ancil-file
#
# 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

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

# 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
#
@H=unpack("i256",substr($IN,0,256*4));
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";
@INTS=unpack("i*",substr($IN,$H[99]*4-4,$H[100]*4));
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] };
if ($H[151] != $NPP or $H[151] != $N_ins) { 
  print "***Warning: I think I read in $NPP fields but I'm inserting $N and the original ancil had $H[151]\n";
  print "***I'm going to change H[152(151)] to $N (or try to...)\n";  
  substr($IN,151*4,4)=pack "i",$N;
#  print "I think I'll change NX and NY too... (to ${NX}x${NY})\n";
  substr($IN,$H[99]*4-4+5*4,4)=pack "i",$NY;
  substr($IN,$H[99]*4-4+6*4,4)=pack "i",$NX;
  print "I think I'll change 160,161 too... (to: ".($H[149]+$N*$H[150]).", ".($N*$LREC).")\n";
  substr($IN,159*4,4)=pack "i",$H[149]+$N*$H[150];
  substr($IN,160*4,4)=pack "i",$N*$LREC;
# Redo this with the new values.
  @H=unpack("i256",substr($IN,0,256*4));
};
if ($FIXHD10 != -1) {
  print "Changing Fixhd(10 (ftn indexing)) from $H[9] ($FIXHD10{$H[9]}) to $FIXHD10 ($FIXHD10{$FIXHD10})\n";
  substr($IN,9*4,4)=pack "i",$FIXHD10;
# Redo this with the new values.
  @H=unpack("i256",substr($IN,0,256*4));
};
if ($DATE != -1) {
  print "changing Fixhd(21 (ftn indexing)) from $H[20] to $DATE \n";
  substr($IN,20*4,4)=pack "i",$DATE;
# Redo this with the new values.
  @H=unpack("i256",substr($IN,0,256*4));  
};
if ($ENDDATE != -1) {
  print "change Fixhd(28 (ftn indexing)) from $H[27] to $ENDDATE \n";
  substr($IN,27*4,4)=pack "i",$ENDDATE;
# Redo this with the new values.
  @H=unpack("i256",substr($IN,0,256*4));  
};  
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 / Header from: ",$H[149]+$F*$H[150]," to ",$H[149]+(1+$F)*$H[150]-1, " / Data from  : ",$H[159]+$F*$LREC," to ",$H[159]+(1+$F)*$LREC-1,"\n";
  $PP=substr($PIN,$PPL*$F,$PPL);
  $PPH=substr($PP,4,64*4); $PPD=substr($PP,4*(3+64),4*$LREC);
  if ($PPH ne substr($IN,($H[149]+$F*$H[150]-1)*4,$H[150]*4)) {
    wwarn("*** Warning: header inserting for field $F does not match the one in the ancil\n ")
  } else {
    print "  header to insert matches that in ancil; good\n"
  };
  substr($IN,($H[149]+$F*$H[150]-1)*4,$H[150]*4)=$PPH;
  substr($IN,($H[159]+$F*$LREC-1)*4,$LREC*4)=$PPD;
  @data=unpack "f*",$PPD;
  print "data range: ".makerange(\@data,unpack("f",substr($PPH,(46+16)*4,4))).", [0,0]: $data[0]\n";

};
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)";

}