#!PERL_EXE

#	@(#)knet2xml 1.0	06/05/2008      vq

use File::Path;
use IO::File;
use FindBin;
use lib "$FindBin::Bin/../lib";
use ShakeConfig;
use lib "$shake_perl_lib";
use Shake::Die;
use Shake::Opt;
use XML::Writer;

use strict;

#######################################################################
# Global variables
#######################################################################

my $sac = '/home/shake/INSTALL/sac/bin/sac';
my $tar = '/bin/tar';
my $tmpdir;

my $alldata;
my $fh;			# The output (XML) filehandle
my $fin;		# The input (jb) filehandle
my $writer;		# The XML::Writer object
my $outfile;
my $filetag = 'knet';
my $dtdfile = "$shake_home/lib/xml/stationlist.dtd";
my $source = 'K-NET (NIED, Japan)';
my $netid = 'KNET';

my %BUTTER = (
  corner => 0.02,
  npoles => 4,
  passes => 2,
);

my %FILE_KEYS = (
  'Station Code' => 'code',
  'Station Lat.' => 'lat',
  'Station Long.' => 'lon',
  'Max. Acc. \(gal\)' => 'pga',
  'Scale Factor' => 'scale',
);

#######################################################################
# End global variables
#######################################################################

#######################################################################
# Stuff to handle command line options and program documentation:
#######################################################################

my $desc = "Generate XML file (following the stationlist.dtd spec) "
	 . "in '$shake_home/data/<EVENTID>/input' "
	 . "from ciim.cdi summary file.";

my $flgs = [{ FLAG => 'event',
              ARG  => 'event_id',
              TYPE => 's',
              REQ  => 'y',
              DESC => 'Specifies the id of the event to process'},
            { FLAG => 'file',
              ARG  => 'file',
              TYPE => 's',
              REQ  => 'y',
              DESC => 'Specifies the filename in the raw directory'},
            { FLAG => 'help',
              DESC => 'Print program documentation and quit.'}
           ];

my $options = Opt::setOptions($flgs) or Die "Error in setOptions";

if (defined $options->{'help'}) {
  Opt::printDoc($desc);
  exit 0;
}

defined $options->{'event'}
        or Die "Must specify an event with -event flag";

my $evid    = $options->{'event'};
my $file    = $options->{'file'};

$outfile = "${filetag}_dat.xml";

Die "Unknown argument(s): @ARGV" if (@ARGV);

#######################################################################
# End of command line option stuff
#######################################################################

main();
0;

sub main {

  my $ctime = time;
  #----------------------------------------------------------------------
  # Open the DTD file for reading
  #----------------------------------------------------------------------
  my $dtdfh = new IO::File($dtdfile) or Die "Couldn't open $dtdfile";

  #----------------------------------------------------------------------
  # Untar stations package
  #----------------------------------------------------------------------

  my $basepath = "$shake_home/data/$evid";
  my $rawpath = "$basepath/raw";

  $tmpdir = "$rawpath/tmp.knet";

  if (!$file) {
    die "Missing data tarfile" unless
      (glob("$tmpdir/*.NS") or glob("$tmpdir/*.EW"));
  }
  else {
    $file = "$rawpath/$file" unless ($file =~ /\//);
    Die "Couldn't find file '$file'" if (!-e $file);

    mkdir $tmpdir unless (-d $tmpdir);

    print "Extracting from tgz file to $tmpdir.\n";
    my $command = "$tar -xvzf $file -C $tmpdir . > /dev/null";
    print "$command\n"; system($command);
  }

  #----------------------------------------------------------------------
  # Extract data from station files
  #----------------------------------------------------------------------
  my @files;
  foreach my $ext qw(NS EW) {
    push @files, glob("$tmpdir/*.$ext");
  }
  die "Could not extract files in $tmpdir" unless (@files);
  print "Got ",(scalar @files)," files.\n";

  my $nf;
  foreach my $file (@files) {
    $nf++; print "File $nf.\n" unless ($nf % 100);
#    last if ($nf > 10);

    my $fin = new IO::File "< $file" or Die "Couldn't open $file";
    my $data = {};

    while (my $line = <$fin>) {
      last if ($line =~ /^Memo.\s+$/);
      chomp $line;
      my $key = substr($line,0,18); $key =~ s/\s+$//;
      my $val = substr($line,18);   $val =~ s/\s+$//;

      my $match;
      foreach my $testkey (keys %FILE_KEYS) {
        next unless ($key =~ /^$testkey/);

        $match = $FILE_KEYS{$testkey};
        last;
      }
      next unless ($match);

      $data->{$match} = $val;
    }
    my $scale = convert_scalefactor($data->{'scale'});
    my $rawfile = extract_waveform($file,$fin);

    my ($pga,$pgv,$sacfile) = correct_waveform($rawfile,$scale); 
    $fin->close;

    my $psa = response_spectra($sacfile,$scale);

    my $stacode = $data->{'code'};
    next unless ($stacode);

    my $dir = ($file =~ /NS$/) ? 'N' : 'E';

    unless (exists $alldata->{$stacode}) {
      $alldata->{$stacode} = {
        lat => $data->{'lat'},
        lon => $data->{'lon'},
        comp => { $dir => {} },
      };
    }
    my $ref = $alldata->{$stacode}{'comp'};

    unless (exists $ref->{$dir}) {
      my $stacomp = {};
      $ref->{$dir} = $stacomp;
    }
    $ref = $ref->{$dir};
    
    addpgm($ref,'pga',$pga);
    addpgm($ref,'pgv',$pgv);

    foreach my $comp qw( psa03 psa10 psa30 ) {
      next unless defined ($psa->{$comp});
      addpgm($ref,$comp,$psa->{$comp});
    }
    
#  use Data::Dumper; print Dumper($alldata); 
  } 
#  use Data::Dumper; print Dumper($alldata); die;

  #----------------------------------------------------------------------
  # Creae the XML::Writer object; the "NEWLINES" implementation is a 
  # disaster so we don't use it, and hence have to add our own all
  # over the place
  #----------------------------------------------------------------------
  $basepath = "$basepath/input";
  $fh = new IO::File "> $basepath/$outfile" 
    or Die "Couldn't open $basepath/$outfile";

  $writer = new XML::Writer(OUTPUT => \*$fh, NEWLINES => 0);

  #----------------------------------------------------------------------
  # Start the document with the XML declaration, and add the
  # document type declaration just in case anyone cares; also 
  # print the root element
  #----------------------------------------------------------------------
  $writer->xmlDecl("US-ASCII", "yes");
  print $fh '<!DOCTYPE stationlist [', "\n";
  while (<$dtdfh>) {
    print $fh $_;
  }
  $dtdfh->close;
  print $fh ']>', "\n";

  $writer->startTag("stationlist", "created" => $ctime);
  $writer->characters("\n");
    
  #----------------------------------------------------------------------
  # Write contents of alldata hash to XML
  #----------------------------------------------------------------------

  foreach my $stacode (sort keys %$alldata) {
    my $data = $alldata->{$stacode};
    new_sta($stacode,$stacode,$data->{'lat'},$data->{'lon'});

    my $ref = $data->{'comp'};
    foreach my $dir (sort keys %$ref) {
      new_comp("$dir");
      my $pga = fmt($ref->{$dir}{'pga'} / 9.81); # gals -> %g conversion 
      my $pgv = fmt($ref->{$dir}{'pgv'});
      write_amp('pga',$pga);
      write_amp('pgv',$pgv);
      foreach my $comp (qw( psa03 psa10 psa30)) {
        write_amp($comp,fmt($ref->{$dir}{$comp} / 9.81));
      }
      end_comp();
    } 
    end_sta();
  }

  $writer->characters("\n");
  $writer->endTag("stationlist");
  $writer->end();
  $fh->close;
  0;
}

sub new_sta {
 
  my $sta  = shift; 
  my $name = shift;
  my $lat  = shift;
  my $lon  = shift;
  my $cdi  = shift;

  $name = fix_name_encoding($name);

  $writer->startTag("station", "code"     => $sta, 
			       "name"     => $name, 
			       "insttype" => 'UNK', 
			       "lat"      => $lat, 
			       "lon"      => $lon, 
			       "source"   => $source,
		               "netid"    => $netid, 
			       "commtype" => 'UNK',
                               "intensity" => $cdi);
  $writer->characters("\n");
  return 1;
}

sub end_sta {

  $writer->endTag("station");
  $writer->characters("\n");
  return 0;
}

sub new_comp {

  my $comp = shift;

  $writer->startTag("comp", "name" => "$comp");
  $writer->characters("\n");
  return 0;
}

sub end_comp {

  $writer->endTag("comp");
  $writer->characters("\n");
  return 0;
}

sub write_amp {

  my ($type, $amp) = @_;

  if (not defined $type or not defined $amp) {
    Die "cdmg2xml::write_amp: (type, amp) must be specified";
    return 1;
  }
  $writer->emptyTag("$type", "value", "$amp");
  $writer->characters("\n");
  return 0;
}

sub fmt {

  return sprintf "%.4f", shift;
}

sub fix_name_encoding {
  my %conv = ( 236 => 'i', 239 => 'i', 249 => 'u' );
  my $name = shift;
  my $new = '';

  foreach my $c (split //,$name) {
    foreach my $asc (keys %conv) {
      $c = $conv{$asc} if (ord $c == $asc);
    }
    $new .= $c;
  }
  if ($name ne $new) {
    print "Changed $name to $new.\n";
    return $new;
  }
  return $name;
}

sub convert_scalefactor {
  my $string = shift;

  my ($num,$den) = split /\//,$string;
  $num =~ s/\(.*?\)//;
  my $scale = $num/$den;
#  $scale *= 10_000;

  return $scale;
}

sub extract_waveform {
  my ($file,$fin) = @_;
  my $c;
  my $fout;

  my $ascfile = "$file.asc";
  open $fout,">$ascfile";
  my $c = _decolumnize($fin,$fout);
  close $fout;

#  print "Got $c points for $ascfile.\n";
  return $ascfile;
}

sub _decolumnize {
  my ($fin,$fout) = @_;
  my $c=0;

  while (my $line = <$fin>) {
    $line =~ s/^\s+//;
    my @pts = split /\s+/,$line;
    $c += scalar @pts;

    foreach (@pts) {
      print $fout (sprintf "%.6f\n",$_);
    }
  }
  return $c;
}

sub correct_waveform {
  my ($file,$scale) = @_;
  my ($pga,$pgv);
  my $fh;
  my $sacfile = "$file.corr.sac";

  my $results = << `ENDCOMMAND`;
cd $tmpdir; $sac << ENDSAC
readalpha $file
rmean
chnhdr DELTA 0.01
hp butter c $BUTTER{'corner'} n $BUTTER{'npoles'} p $BUTTER{'passes'} 
taper
commit
lh DEPMIN DEPMAX DEPMEN
write alpha $sacfile
int
lh DEPMIN DEPMAX DEPMEN
quit
ENDSAC
ENDCOMMAND


  my ($minpga,$minpgv) = ($results =~ /DEPMIN = (.*?)\s/g);
  my ($maxpga,$maxpgv) = ($results =~ /DEPMAX = (.*?)\s/g);
  my ($mean) = ($results=~ /DEPMEN = (.*?)\s/);

  $pga = (-$minpga > $maxpga) ? -$minpga : $maxpga;
  $pgv = (-$minpgv > $maxpgv) ? -$minpgv : $maxpgv;
  $pga *= $scale; $pgv *= $scale;
  print "mean:$mean PGA:$pga PGV:$pgv scale:$scale\n";
  return ($pga,$pgv,$sacfile);
}

sub response_spectra {
  my ($sacfile,$scale) = @_;
  my $results;

  my ($fin,$fout);
  my $txtfile = $sacfile;
  $txtfile =~ s/\.sac$//;

  # Correct sacfile into 1 column ascii file
  open $fin,$sacfile;
  open $fout,">$txtfile";
  foreach (1..30) { <$fin>; }
  _decolumnize($fin,$fout);

  close $fin;
  close $fout;

  my $command = "../util/response_spectra $txtfile";
#  print "Command: $command\n";
  my @results = `$command`;

  while (@results) {
    my $check = shift @results;
    last if ($check =~ /^Sa =/);
  }

  foreach my $comp (qw(psa03 psa10 psa30)) {
    $results->{$comp} = (shift @results) * $scale;
  }

  return $results;
}

sub addpgm {
  my ($ref,$key,$val) = @_;

  unless (exists $ref->{$key} and $ref->{$key} > $val) {
    $ref->{$key}= $val;
  }
  return $ref->{$key};
}