#!perl -w
#
# GPSy拡張形式のトラックログをMPS形式に
# Nowral
# 00/9/21
#
use Time::Local;
$tint = 3600; # 表示間隔 [秒]

# 定数
$pi  = 4 * atan2(1,1); # 円周率
$rd  = $pi / 180;      # [ラジアン/度]
$g = 9.80665;   # 重力加速度 [m/s^2]
# データム諸元 (Tokyo)
$a = 6377397.155; # 赤道半径
$f = 1 / 299.152813; # 扁平率
$e2 = 2*$f - $f*$f; # 第1離心率

# 出力ファイル準備
$0 =~ /(.+:)/;
$newargv = "$1UserLine.mps";
open(OUT, ">$newargv");
MacPerl::SetFileInfo('pa4M', 'TEXT', $newargv);

# ヘッダ出力
print OUT <<END_OF_HEADER;
// MapServer Script : generated with 'GPSy2MPS2'
mapserver script ver="1.0" savable="yes";
theme
  str="";
message
  str="";
END_OF_HEADER

# データ読み込み
undef @ts;
undef @bs;
undef @ls;
undef @hs;
undef @strts;
$oldargv = "";
$n = 0;
while (<>) {
  next unless /^T\t/;
  @ed = split("\t");
  if ($ARGV ne $oldargv || $ed[1] ne $tid) {
    push(@strts, $n);
    $oldargv = $ARGV;
    $tid = $ed[1];
  }
  if($ed[9] ne "DMS") { die "$ed[9] Not DMS !"; }
#  if($ed[10] ne "Tokyo") { die "$ed[10] Not Tokyo Datum !"; }

  $ed[4] =~ /(\d+)\/(\d+)\/(\d+)\s+(\d+):(\d+):(\d+)/; # 時刻
  push(@ts, timelocal($6,$5,$4,$2,$1-1,$3-1900));
  $ed[5] =~ /(\d+)｡(\d+)'(\d+\.?\d*)"/; # 緯度
  push(@bs, $1 + $2/60 + $3/3600);
  $ed[6] =~ /(\d+)｡(\d+)'(\d+\.?\d*)"/; # 経度
  push(@ls, $1 + $2/60 + $3/3600);
  push(@hs, $ed[8]); # 高さ <>楕円体高
  ++$n;
}
push(@strts, $n);
print "# 全データ数： $n\n";

# グループループ
$ldt = 0;
foreach $ig (0 .. $#strts-1) {
  $strt = $strts[$ig];
  $end = $strts[$ig+1] - 1;

  print OUT <<END_OF_HEADER2;
group
  gid = "$ig"
  gname = ""
  textrgb = "(0,0,0)"
  textstyle = "plain"
  textsize = "12"
  textoff = "100000";
END_OF_HEADER2

# トラックループ
  undef @srs;
  undef @vr2s;
  undef @ar2s;
  ($t2, $b2, $l2, $h2) = ($ts[$strt], $bs[$strt], $ls[$strt], $hs[$strt]);
  ($x2, $y2, $z2) = &llh2xyz($b2, $l2, $h2, $a, $e2);
  $t0 = $t2;
  ($sr, $vmax, $amax, $dt1) = (0, 0, 0, 0);
  $td2 = 0;
  push(@srs, $sr);
  foreach $i ($strt+1 .. $end){
    ($t3, $b3, $l3, $h3) = ($ts[$i], $bs[$i], $ls[$i], $hs[$i]);
    ($x3, $y3, $z3) = &llh2xyz($b3, $l3, $h3, $a, $e2);
# 距離
    ($dt2, $dx2, $dy2, $dz2) = ($t3-$t2, $x3-$x2, $y3-$y2, $z3-$z2);
    $dr2 = sqrt($dx2*$dx2 + $dy2*$dy2 + $dz2*$dz2);
    $sr += $dr2;
    push(@srs, $sr);
# 速度、方位
    if($dt2) {
      ($vx2, $vy2, $vz2)  = ($dx2/$dt2, $dy2/$dt2, $dz2/$dt2);
      $vr2 = $dr2 / $dt2;
    }
    else { $vr2 = 0; }
    if($vr2>$vmax) { $vmax = $vr2; }
    push(@vr2s, $vr2);
# 加速度
    if($dt1) { # 以降、3点必要な量
      $dtm = ($dt1+$dt2) / 2;
      ($ax2, $ay2, $az2) = (($vx2-$vx1)/$dtm, ($vy2-$vy1)/$dtm, ($vz2-$vz1)/$dtm);
      $ar2 = sqrt($ax2*$ax2 + $ay2*$ay2 + $az2*$az2);
    }
    else { $ar2 = 0; }
    if($ar2>$amax) { $amax = $ar2; }
    push(@ar2s, $ar2);

# 次段準備
    ($t2, $x2, $y2, $z2) = ($t3, $x3, $y3, $z3);
    ($b2, $l2) = ($b3, $l3);
    ($vx1, $vy1, $vz1) = ($vx2, $vy2, $vz2);
    $dt1 = $dt2;
  }
  $st = $t3 - $t0;

#集計
  print "\n";
  printf "# 開始時刻  ： %s\n", &formd($t0);
  printf "# 終了時刻  ： %s\n", &formd($t3);
  printf "# データ数  ： %d\n", $end-$strt+1;
  printf "# 総時間    ： %d:%02d:%02d\n", int($st/3600), ($st/60)%60, $st%60;
  printf "# 移動距離  ： %.1f [km]\n", $sr/1000;
  printf "# 平均速度  ： %.1f [km/h]\n", $sr/$st*3.6;
  printf "# 最高速度  ： %.1f [km/h]\n", $vmax*3.6;
  printf "# 最大加速度： %.2f [G]\n", $amax/$g;

# トラックループ
  ($t2, $b2, $l2) = ($ts[$strt], $bs[$strt], $ls[$strt]);
  $td2 = 0;
  foreach $i ($strt+1 .. $end){
    $j = $i - $strt - 1;
    ($t3, $b3, $l3) = ($ts[$i], $bs[$i], $ls[$i]);
    ($sr3, $vr2, $ar2) = ($srs[$j], $vr2s[$j], $ar2s[$j]);

# 線巾
    $lw = $ar2 * 10 + 2;
    $lw = 12 if $lw > 12;

# カラー
    $yy = - 64 * cos(int($t3/3600)%24/24*$rd);
    $th = sqrt($vr2) * 15;
#    $th = 90 if $th > 90;
    $th *= $rd;
    $rr = $yy + 256 * sin($th);
    $gg = $yy + 256 * cos($th);
    $bb = $yy;
    $rr = 255 if $rr > 255;
    $rr = 0 if $rr < 0;
    $gg = 255 if $gg > 255;
    $gg = 0 if $gg < 0;
    $bb = 255 if $bb > 255;
    $bb = 0 if $bb < 0;
    $clr = sprintf("%d,%d,%d", $rr, $gg, $bb);

# 線素片表示
    print OUT "line\n\tgid=\"$ig\"\n\tpos=\"",
     &deg2dms($b2), ",", &deg2dms($l2), ",",
     &deg2dms($b3), ",", &deg2dms($l3), "\"";
    printf OUT "\n\tlinewidth=\"%d\"", $lw;
    print OUT "\n\tlinergb=\"$clr\"";
    if($t2>=$td2 || ($ar2<0.1 && $vr2<0.3)) { # 時刻表示のタイミング?
      if($t2>$ldt+$tint/2) {
        $ldt = $t2;
        print OUT "\n\tstr=\"", &formd($t2), "\"\n\ttextrgb = \"$clr\"";
      }
      $td2 = int($t2/$tint + 1) * $tint;
    }
    print OUT ";\n";

# 円表示
    if($ar2*3.6>1000) {
    print OUT "circle\n\tgid = \"$ig\"\n\tpos=\"",
    &deg2dms($b2), ",", &deg2dms($l2), "\"";
    printf OUT "\n\tr=\"%.1fm\"", $ar2/3.6;
    print OUT "\n\tbodrrgb=\"$clr\"";
    print OUT "\n\tbodrstyle = \"line\#0\"";
    print OUT "\n\tbodrwidth= \"0\"";
    print OUT "\n\tareargb=\"$clr\"";
    print OUT "\n\tareastyle=\"fill\"";
    print OUT ";\n";
    }

# 次段準備
    ($t2, $b2, $l2) = ($t3, $b3, $l3);
  }
  print OUT "point\n\tgid=\"$ig\"\n\tpos=\"",
   &deg2dms($b2), ",", &deg2dms($l2), "\"\n";
  print OUT "\tstr=\"", &formd($t2), "\";\n";

# 終了
  print " ... done\n";
}

# フッタ出力
print OUT <<END_OF_FOOTER;
disp all;
mapserver end;
END_OF_FOOTER
close(OUT);

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



sub formd { # 日付、時刻の表示形式
  my($t) = @_;
  my($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst) = localtime($t);
  return sprintf("%d:%02d", $hour,$min);
#  return sprintf("%04d/%d/%d %d:%02d:%02d", $year+1900,$mon+1,$mday,$hour,$min,$sec);
}

sub deg2dms { # 経緯度、DD.DDDD形式 -> MPS形式
  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;
  return sprintf("%d/%02d/%02d.%d", $d, $m, $s, $sf);
}

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