#!/usr/bin/perl -w
#
# WGS 84 -> Tokyo
# Nowral
# 99/9/23
#
# 定数
$pi  = 4 * atan2(1,1); # 円周率
$rd  = $pi / 180;      # [ラジアン/度]

# データム諸元
# 変換元
# (WGS 84)
$a = 6378137.0; # 6378137; # 赤道半径
$f = 1 / 298.257223563; # 1 / 298.257223; # 扁平率
$e2 = 2*$f - $f*$f; # 第1離心率

# 変換先
# (Tokyo)
$a_ = 6378137.0 - 739.845; # 6377397.155;
$f_ = 1/298.257223563 - 0.000010037483; # 1 / 299.152813;
$e2_ = 2*$f_ - $f_*$f_;

# 並行移動量 [m]
# e.g. $x_ = $x + $dx etc.
$dx = +128; # -148;
$dy = -481; # +507;
$dz = -664; # +681;

@files = @ARGV;
foreach $oldargv (@files) {

$newargv = $oldargv . "(TKY)";
open(OUT, ">$newargv");
&MacPerl'SetFileInfo('GPSy', 'TEXT', $newargv);

open(IN, "$oldargv");
while(<IN>){
  if(/^T\t/) {
    @ed = split("\t");
    if($ed[9] ne "DMS") { die "Not DMS !"; }
#    unless($ed[10]=~/WGS\s*84/i || $ed[10]=~/<NO TRANSLATION>/) { die "$ed[10] Datum Not Supported"; }
    $ed[10] = "Tokyo";

    $ed[5] =~ /(\d+)｡(\d+)'(\d+\.\d+)"/; # 緯度
    $b = $1 + $2/60 + $3/3600;
    $ed[6] =~ /(\d+)｡(\d+)'(\d+\.\d+)"/; # 経度
    $l = $1 + $2/60 + $3/3600;
    $h = $ed[8]; # 高さ <>楕円体高
    ($x, $y, $z) = &llh2xyz($b, $l, $h, $a, $e2);
    ($b, $l, $h) = &xyz2llh($x+$dx, $y+$dy, $z+$dz, $a_, $e2_);
    $ed[5] = &deg2gdms($b);
    $ed[6] = &deg2gdms($l);
#    $ed[8] = sprintf("%.1f", $h);

    $_ = join("\t", @ed);
  }
  print OUT;
}
close(IN);
close(OUT);

}
#&MacPerl'Quit(2);
die "The Unhappy End";




sub llh2xyz { # 楕円体座標 -> 直交座標
  my($b, $l, $h, $a, $e2) = @_;
  my($sb, $cb, $rn, $x, $y, $z);

  $b *= $rd;
  $l *= $rd;
  $sb = sin($b);
  $cb = cos($b);
  $rn = $a / sqrt(1-$e2*$sb*$sb);

  $x = ($rn+$h) * $cb * cos($l);
  $y = ($rn+$h) * $cb * sin($l);
  $z = ($rn*(1-$e2)+$h) * $sb;

  ($x, $y, $z);
}

sub xyz2llh { # 直交座標 -> 楕円体座標
  my($x, $y, $z, $a, $e2) = @_;
  my($bda, $p, $t, $st, $ct, $b, $l, $sb, $rn, $h);
  $bda = sqrt(1-$e2); # b/a

  $p = sqrt($x*$x+$y*$y);
  $t = atan2($z, $p*$bda);
  $st = sin($t);
  $ct = cos($t);
  $b = atan2($z+$e2*$a/$bda*$st*$st*$st, $p-$e2*$a*$ct*$ct*$ct);
  $l = atan2($y, $x);

  $sb = sin($b);
  $rn = $a / sqrt(1-$e2*$sb*$sb);
  $h = $p/cos($b) - $rn;

  ($b/$rd, $l/$rd, $h);
}

sub deg2gdms { # GPSy形式の経緯度
  my($d) = @_;
  my($m, $s, $sf);
  $sf = int($d*36000 + 0.5);
  $s = int($sf/10) % 60;
  $m = int($sf/600) % 60;
  $d = int($sf/36000);
  $sf %= 10;
  sprintf("%d｡%02d\'%02d\.%d\"", $d, $m, $s, $sf);
}