#!/usr/bin/perl -w
#
# この座標の水深は?
# 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;

# 経緯度入力
MacPerl::Ask("緯度? (北緯：入力例  43/38/27.9)", "43/38/27.9") =~ m/(\d+)[^\d](\d+)[^\d](\d+\.?\d*)/ || die;
$b = $1*3600 + $2*60 + $3;
MacPerl::Ask("経度? (東経：入力例  144/23/40.8)", "144/23/40.8") =~ m/(\d+)[^\d](\d+)[^\d](\d+\.?\d*)/ || die;
$l = $1*3600 + $2*60 + $3;
die unless($b<$lat0-$hlat && $b>$lat2+$hlat && $l>$lon0+$hlon && $l<$lon2-$hlon);

undef @ofss;
undef @hofs;
undef @hs;

# 近接三点を選択 < ロジック甘い?
$i = 0;
$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++;
#print join("\n", @ofss), "\n";

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

# 平面で補間
$i = 0;
$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;

# 結果出力
print "緯度：", &sec2dms($b), "\n";
print "経度：", &sec2dms($l), "\n";
printf("水深： %.02f [m]\n", $h);
MacPerl::Ask("水深 [m]", sprintf("%.02f", $h));

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

sub sec2dms { # 秒から度分秒への変換、精度指定あり
  my($s) = @_;
  my($dfe) = 3; # 小数点以下必要な桁数
  my($df) = 10 ** $dfe;
  my($m, $d, $sf);
  $sf = int($s*$df + 0.5);
  $s  = $sf / $df % 60;
  $m  = $sf / $df / 60 % 60;
  $d  = int($sf / $df / 3600);
  $sf %= $df;
#  ($d, $m, $s, $sf);
  sprintf("%d｡%02d\'%02d\.%0$dfe\d\"", $d, $m, $s, $sf);
}

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