#!/usr/bin/perl -w
#
# japan1km.filled.grd -> TEM
# Nowral
# 00/4/5
#
# 内部の経緯度は全て秒単位
$0 =~ /^(.*):(.*)$/;
$cwd = $1;
&UTM::utminit();

# japan1km.filled.grdデータ諸元 > グローバル変数
$hds = 0x043C;
$npix = 2080;
$nlin = 2573;
$lat0 =  46 * 3600;
$lon0 = 122 * 3600;
$lat2 =  24 * 3600;
$lon2 = 148 * 3600;
$dlat = 30; # 0.00833333333333333 * 3600;
$dlon = 45; # 0.0125 * 3600;

$lat2 = $lat0 - $nlin*$dlat;
$lon2 = $lon0 + $npix*$dlon;
$hlat = $dlat / 2;
$hlon = $dlon / 2;

# 入力ファイルのループ
@files = sort grep {-f} @ARGV;
foreach $oldargv (@files) {
$oldargv =~ /^(.*):(.*)$/;
print "\n$2\n";

# ファイル読み込み >経緯度から標高ファイル決定
undef @ofss;
undef @hofs;
undef @hs;
$i = 0;
open(IN, $oldargv);
while(<IN>){
  next unless /^T\t/;
  @ed = split("\t");
  if($ed[9] ne "DMS") { die "Not DMS !"; }
#  die unless($ed[10]=~/WGS\s*84/i || $ed[10]=~/<NO TRANSLATION>/);

  $ed[5] =~ /(\d+)｡(\d+)'(\d+\.?\d*)"/; # 緯度
  $b = $1*3600 + $2*60 + $3;
  $ed[6] =~ /(\d+)｡(\d+)'(\d+\.?\d*)"/; # 経度
  $l = $1*3600 + $2*60 + $3;
# 近接三点を選択 < ロジック甘い?
  $ofs1 = &ll2ofs($b, $l);
  push(@ofss, $ofs1);
  $ofsd = &ll2ofs($b, $l-$hlon); # 左
  push(@ofss, $ofsd) unless($ofsd == $ofs1);
  $ofsd = &ll2ofs($b-$hlat, $l); # 下
  push(@ofss, $ofsd) unless($ofsd == $ofs1);
  $ofsd = &ll2ofs($b, $l+$hlon); # 右
  push(@ofss, $ofsd) unless($ofsd == $ofs1);
  $ofsd = &ll2ofs($b+$hlat, $l); # 上
  push(@ofss, $ofsd) unless($ofsd == $ofs1);
  die unless($#ofss == $i+2);

  @hs[$i .. $i+2] = ("-", "-", "-");
  $hofs{$ofss[$i]} .= " $i";
  $i++;
  $hofs{$ofss[$i]} .= " $i";
  $i++;
  $hofs{$ofss[$i]} .= " $i";
  $i++;
}
$nr = $i / 3;

# 水深読み込み <まとめて
open(IN, "$cwd:japan1km.filled.grd") || die;
foreach $ofs (sort keys %hofs) {
  seek(IN, $ofs, 0);
  read(IN, $_, 4);
  $h = unpack("f", $_);
  foreach $i (split(" ", $hofs{$ofs})) { $hs[$i] = $h; }
}

# 出力ファイル準備
$newargv = $oldargv . '.w/d';
open(OUT, ">$newargv");
MacPerl::SetFileInfo('GPSy', 'TEXT', $newargv);

# 標高補足
$i = 0;
open(IN, $oldargv);
while(<IN>){
  if(/^T\t/) {
    @ed = split("\t");
    $ed[5] =~ /(\d+)｡(\d+)'(\d+\.?\d*)"/; # 緯度
    $b = $1*3600 + $2*60 + $3;
    $ed[6] =~ /(\d+)｡(\d+)'(\d+\.?\d*)"/; # 経度
    $l = $1*3600 + $2*60 + $3;
# 平面で補間
    $zn = ($l/3600+180) / 6 % 60 + 1; # ゾーン < 西境界含む
    ($x, $y) = &UTM::sec2utm($b, $l, $zn);
    ($h1, $h2, $h3) = @hs[$i .. $i+2];
    ($x1, $y1) = &UTM::sec2utm(&ofs2cll($ofss[$i]), $zn);
    $i++;
    ($x2, $y2) = &UTM::sec2utm(&ofs2cll($ofss[$i]), $zn);
    $i++;
    ($x3, $y3) = &UTM::sec2utm(&ofs2cll($ofss[$i]), $zn);
    $i++;

    $d = $x1*$y2*$h3 + $x2*$y3*$h1 + $x3*$y1*$h2 - $x1*$y3*$h2 - $x2*$y1*$h3 - $x3*$y2*$h1;
    $u = $y2*$h3 + $y3*$h1 + $y1*$h2 - $y3*$h2 - $y1*$h3 - $y2*$h1;
    $v = $x1*$h3 + $x2*$h1 + $x3*$h2 - $x1*$h2 - $x2*$h3 - $x3*$h1;
    $w = $x1*$y2 + $x2*$y3 + $x3*$y1 - $x1*$y3 - $x2*$y1 - $x3*$y2;
    $h = ($d - $u*$x - $v*$y) / $w;
    $ed[8] = sprintf("%.1f", $h);
    $_ = join("\t", @ed);
  }
  else{ s/^(#\tRecords:\t)\d+/$1$nr/; } # データ数は正確に

  print OUT; # 書き出し
}
close(IN);
close(OUT);

# 1ファイル終了
print " ... done\n";
}

MacPerl::Quit(2);
die "The Unhappy End";



sub ll2ofs { # 経緯度 -> 先頭からのオフセット
  my($b, $l) =@_;
  my($x, $y);
  $x = int(($l - $lon0)/$dlon);
  $y = int(($lat0 - $b)/$dlat);
  return ($y*$npix + $x)*4 + $hds;
}

sub ofs2cll { # 先頭からのオフセット -> メッシュ中心の経緯度
  my($ofs) =@_;
  my($b, $l);
  my($x, $y);
  $ofs = ($ofs - $hds) / 4;
  $x = $ofs % $npix;
  $y = ($ofs - $x) / $npix;
  $b = $lat0 - $y * $dlat;
  $l = $lon0 + $x * $dlon;
  return ($b - $hlat, $l + $hlon);
}

package UTM; # グローバル変数隔離のため
# $pi, $rd, $a, $f, $e2, $ed2, $cad, $cbd, $ccd, $cdd, $ced

sub utminit { # グローバル変数準備
#定数
  $pi  = 4 * atan2(1, 1); # 円周率
  $rd  = $pi / 180;       # radian / degree

# 日本測地系
  $a   = 6377397.155;    # 赤道半径 (Tokyo)
  $f   = 1 / 299.152813; # 扁平率 (Tokyo)
  $e2  = 2*$f - $f*$f;   # 第1離心率
  $ed2 = $e2 / (1-$e2);  # 第2離心率

# 子午線弧長計算式係数
  my($e4, $e6, $e8);
  $e4 = $e2 * $e2;
  $e6 = $e4 * $e2;
  $e8 = $e6 * $e2;

  $cad = $a * (1-$e2) * ( 1 + 3/4*$e2 + 45/64*$e4 + 175/256*$e6 + 11025/16384*$e8 );
  $cbd = $a * (1-$e2) * (     3/4*$e2 + 15/16*$e4 + 525/512*$e6 +   2205/2048*$e8 );
  $ccd = $a * (1-$e2) * (               15/64*$e4 + 105/256*$e6 +   2205/4096*$e8 );
  $cdd = $a * (1-$e2) * (                            35/512*$e6 +    315/2048*$e8 );
  $ced = $a * (1-$e2) * (                                           315/16384*$e8 );

  return 1;
}

sub sec2utm { # 経緯度 -> UTM座標
  my($b, $l, $zn) =@_;
  my($x, $y) = &ll2xy($b/3600*$rd, ($l/3600-$zn*6+183)*$rd);
  return ($x * 0.9996, $y * 0.9996 + 500000);
}

sub ll2xy { # Gauss-Klueger
  my($b, $dl) = @_;
  my($dl2, $dl4, $sb, $cb, $cb3, $cb5, $tb2, $tb4, $eta2, $rn, $x, $y);

  $dl2 = $dl * $dl;
  $dl4 = $dl2 * $dl2;
  
  $sb = sin($b);
  $cb = cos($b);
  $cb3 = $cb * $cb * $cb;
  $cb5 = $cb3 * $cb * $cb;
  $tb2 = $sb * $sb / $cb / $cb;
  $tb4 = $tb2 * $tb2;
  
  $eta2 = $ed2 * $cb * $cb;
  $rn = $a / (1-$f) / sqrt(1+$eta2); # 卯酉線曲率半径
  
  $x = $cad*$b - $cbd/2*sin(2*$b) + $ccd/4*sin(4*$b)
    - $cdd/6*sin(6*$b) + $ced/8*sin(8*$b); # 赤道から緯度Bまでの子午線弧長
  $x += $rn/2*$sb*$cb*$dl2;
  $x += $rn/24*$sb*$cb3*(5-$tb2+9*$eta2+4*$eta2*$eta2)*$dl4;
  $x += $rn/720*$sb*$cb5*(61-58*$tb2+$tb4)*$dl4*$dl2;

  $y = $rn*$cb*$dl;
  $y += $rn/6*$cb3*(1-$tb2+$eta2)*$dl*$dl2;
  $y += $rn/120*$cb5*(5-18*$tb2+$tb4)*$dl*$dl4;

  return ($x, $y);
}