#!/usr/bin/perl
#
# 差分
# Nowral
# 99/9/6
#
#時間間隔
use Time::Local;

# 定数
$pi = 4 * atan2(1,1); # 円周率
$rd = $pi / 180;      # radian / degree

# 前段
open(A, $ARGV[0]) || die "Open A Failed :$!";
($ta0, $xa0, $ya0, $ha0) = &getg(*A);
($ta1, $xa1, $ya1, $ha1) = &getg(*A);
open(B, $ARGV[1]) || die "Open B Failed :$!";
($tb0, $xb0, $yb0, $hb0) = &getg(*B);
($tb1, $xb1, $yb1, $hb1) = &getg(*B);

$newargv = $ARGV[0] . ".dif";
open(OUT, ">$newargv");
&MacPerl'SetFileInfo('ttxt', 'TEXT', $newargv);
$fh = select(OUT); # 前の出力先を記憶

print "#時刻\t緯度[m]\t経度[m]\t高度[m]\t距離[m]\n";
undef($n, $sdx, $sdy, $sdh, $sr2, $t0);
undef($rmax, $rmin, $amax, $amin);

# ループ
while( ! (eof(A) || eof(B)) ) {
  undef $u;
  if($ta0>=$tb0 && $ta0<$tb1) {
    $ts = &formd($ta0);
    if(! defined $t0) { $t0 = $ta0; }
    $t1 = $ta0 - $t0;
    $u = ($tb1-$ta0) / ($tb1-$tb0);
    $dx = $xa0 - $u*$xb0 - (1-$u)*$xb1;
    $dy = $ya0 - $u*$yb0 - (1-$u)*$yb1;
    $dh = $ha0 - $u*$hb0 - (1-$u)*$hb1;
  }
  elsif($tb0>$ta0 && $tb0<$ta1) {
    $ts = &formd($tb0);
    if(! defined $t0) { $t0 = $tb0; }
    $t1 = $tb0 - $t0;
    $u = ($ta1-$tb0) / ($ta1-$ta0);
    $dx = $u*$xa0 + (1-$u)*$xa1 - $xb0;
    $dy = $u*$ya0 + (1-$u)*$ya1 - $yb0;
    $dh = $u*$ha0 + (1-$u)*$ha1 - $hb0;
  }
  if(defined $u) {
    ++$n;
#    printf "%4d\t", $n;
# 距離
    $r = sqrt($dx*$dx + $dy*$dy); # 高度抜き
#    $r = sqrt($dx*$dx + $dy*$dy + $dh*$dh); # 高度入り
    if($r>$rmax || ! defined $rmax) { $rmax = $r; }
    if($r<$rmin || ! defined $rmin) { $rmin = $r; }
    $sdx += $dx;
    $sdy += $dy;
    $sdh += $dh;
    $sr2 += $r * $r;
    printf "%s\t%.0f\t%.0f\t%.0f\t%.1f", $ts, $dx, $dy, $dh, $r;
#    printf "%5d\t%.0f\t%.0f\t%.0f\t%.2f", $t1, $dx, $dy, $dh, $r;
# 方位
    if($r) {
      $a = 450 - atan2($dx,$dy)/$rd; # 北から東回り
      if($a>=360) { $a -= 360; } # modを使うと整数にされてしまう
      if($a>$amax || ! defined $amax) { $amax = $a; }
      if($a<$amin || ! defined $amin) { $amin = $a; }
      printf "\t%.0f", $a;
    }
    print "\n";
  }

  if($ta1<=$tb1) { # スライド
    ($ta0, $xa0, $ya0, $ha0) = ($ta1, $xa1, $ya1, $ha1);
    ($ta1, $xa1, $ya1, $ha1) = &getg(*A);
  }
  else {
    ($tb0, $xb0, $yb0, $hb0) = ($tb1, $xb1, $yb1, $hb1);
    ($tb1, $xb1, $yb1, $hb1) = &getg(*B);
  }
}
$rav = sqrt($sdx*$sdx + $sdy*$sdy + $sdh*$sdh) / $n;
$rer = sqrt( ($sr2 - $n*$rav*$rav) / ($n-1) );
$aav = 450 - atan2($sdx,$sdy)/$rd; # 北から東回り
if($aav>=360) { $aav -= 360; }
#$aer = atan2($rer, sqrt($rav*$rav-$rer*$rer)) / $rd;
$aer = atan2($rer, $rav) / $rd;

close(A);
close(B);
close(OUT);
select($fh);

# 集計
print "\n";
printf "# データ数 ： %d\n", $n;
printf "# 平均間隔 ： %.1f ± %.2f [m]\n", $rav, $rer/sqrt($n);
printf "# 1点の誤差： ± %.1f [m]\n", $rer;
printf "# 最大間隔 ： %.1f [m]\n", $rmax;
printf "# 最小間隔 ： %.1f [m]\n", $rmin;
print "\n";
printf "# 平均角度 ： %.1f [度]\n", $aav;
printf "# 1点の誤差： ± %.1f [度]\n", $aer;
printf "# 最大角度 ： %.0f [度]\n", $amax;
printf "# 最小角度 ： %.0f [度]\n", $amin;

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



sub getg { # GPSy拡張形式の読み込み
  local(*IN) = @_;
  undef $_;
  $_ = <IN> until /^T\t/;
  my(@ed) = split("\t");
  if($ed[9] ne "UTM") { die "$ed[9] Not UTM!\n"; }
  $ed[4] =~ /(\d+)\/(\d+)\/19(\d+)\s+(\d+):(\d+):(\d+)/;
  (timelocal($6,$5,$4,$2,$1-1,$3), $ed[5], $ed[6], $ed[8]);
}

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

sub avcalc { # 平均、誤差
  my($s,$s2,$n) =@_;
#  ( $s/$n , sqrt(($s2-$s*$s/$n)/($n-1)) ); # 1点の誤差
  ( $s/$n , sqrt(($s2-$s*$s/$n)/$n/($n-1)) ); # 平均の誤差
}