<?php
//================================================
//	[工事名称]
//		-名古屋 航跡波警報システム-
//	[概要]
//		引数で与えられた秒数sleep後AISデータを取得、波高の計算を行う
//		計算結果が波高80cm以上で施工箇所までの予想到達秒数が180秒の場合注意となり
//		300秒の場合警告とする
//		計算の結果は日付毎に../output/yyyyMMdd.txtに出力を行う
//	[備考]
//		http://penta-02.ptgw.net/ttest/nagoya/map/map.phpから計算部分を流用
//	[更新]
//		2014/06/30 KMK-Teraoka
//================================================

//-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
//	指定位置から、距離、方位角を指定して到達点を算出
//	関数名	:	vincentyDirect
//	引数	:	$lon 経度
//				$lat 緯度
//				$azimuth 方位角(rad型式)
//				$distance 距離
//				$imax 収束計算回数
//	戻り値	:	array(
//					lon 経度
//					lat 緯度
//					azimuth2 算出位置からみた指定位置への角度
//				)
//	備考	:	参考ページ - http://www.gammasoft.jp/direct/
//	更新	:	2014/06/30 KMK-Teraoka
//-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
function vincentyDirect($lon, $lat, $azimuth, $distance, $imax)
{
	$a = 6378137.0;
	$b = 6356752.31414036;
	$f = 1 / 298.257222101;
	
	$sinAz1 = sin($azimuth);
	$cosAz1 = cos($azimuth);
	$U1 = atan((1 - $f) * tan($lat * pi() / 180.0));
	$sigma1 = atan2((1 - $f) * tan($lat * pi() / 180.0), $cosAz1);
	$sinU1 = sin($U1);
	$cosU1 = cos($U1);
	$sinAlfa = $cosU1 * $sinAz1;
	$cosAlfa2 = 1 - pow($sinAlfa, 2);
	$uu = $cosAlfa2 * ($a * $a - $b * $b) / ($b * $b);
	$A = 1 + ($uu / 16384.0) * (4096 + $uu * (-768 + $uu * (320 - 175 * $uu)));
	$B = ($uu / 1024.0) * (256 + $uu * (-128 + $uu * (74 - 47 * $uu)));
	$sigma = $distance / ($b * $A);
	
	$sigma0 = 0;
	$sinSigma = 0;
	$cosSigma = 0;
	$cos2Sigmam = 0;
	$dSigma = 0;
	
	for ($i = 0; $i < $imax; $i++)
	{
		$cos2Sigmam = cos(2 * $sigma1 + $sigma);
		$sinSigma = sin($sigma);
		$cosSigma = cos($sigma);
		$dSigma = $B * $sinSigma * ($cos2Sigmam + ($B / 4.0) * ($cosSigma * (-1 + 2 * pow($cos2Sigmam, 2)) - ($B / 6.0) * $cos2Sigmam * (-3 + 4 * pow($sinSigma, 2)) * (-3 + 4 * pow($cos2Sigmam, 2))));
		$sigma0 = $sigma;
		$sigma = $distance / ($b * $A) + $dSigma;
		
		if (abs($sigma - $sigma0) <= 1.0e-12)
			break;
	}
	$cos2Sigmam = cos(2 * $sigma1 + $sigma);
	$sinSigma = sin($sigma);
	$cosSigma = cos($sigma);
	
	$phi2 = atan2($sinU1 * $cosSigma + $cosU1 * $sinSigma * $cosAz1, (1 - $f) * sqrt($sinAlfa * $sinAlfa + pow($sinU1 * $sinSigma - $cosU1 * $cosSigma * $cosAz1, 2)));
	$lamuda = atan2($sinSigma * $sinAz1, $cosU1 * $cosSigma - $sinU1 * $sinSigma * $cosAz1);
	$C = ($f / 16.0) * $cosAlfa2 * (4 + $f * (4 - 3 * $cosAlfa2));
	$omega = $lamuda - (1 - $C) * $f * $sinAlfa * ($sigma + $C * $sinSigma * ($cos2Sigmam + $C * $cosSigma * (-1 + 2 * pow($cos2Sigmam, 2))));
	$lamuda2 = $lon * pi() / 180.0 + $omega;
	$az2 = atan2($sinAlfa, -($sinU1) * $sinSigma + $cosU1 * $cosSigma * $cosAz1) + pi();
	if ($az2 < 0)
		$az2 = $az2 + pi() * 2.0;
	if ($i >= $imax)
		$i = $imax - 1;
	return array(
		"lon" => $lamuda2 * 180 / pi(),
		"lat" => $phi2 * 180 / pi(),
		"azimuth2" => $az2
	);
}

//-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
//	2点間の距離を求める
//	関数名	:	vincentyInverse
//	引数	:	$lon1 点1の経度
//				$lat1 点1の緯度
//				$lon2 点2の経度
//				$lat2 点2の緯度
//				$imax 収束計算回数
//	戻り値	:	array(
//					d 距離
//					azimuth1 点1から点2への方位角
//					azimuth2 点2から点1への方位角
//				)
//	備考	:	参考ページ - http://www.gammasoft.jp/distance/
//	更新	:	2014/06/30 KMK-Teraoka
//-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
function vincentyInverse($lon1, $lat1, $lon2, $lat2, $imax)
{
	$a = 6378137.0;
	$b = 6356752.31414036;
	$f = 1 / 298.257222101;
	$omega = ($lon2 - $lon1) * pi() / 180.0;
	$U1 = atan((1 - $f) * tan($lat1 * pi() / 180.0));
	$U2 = atan((1 - $f) * tan($lat2 * pi() / 180.0));
	$sinU1 = sin($U1);
	$cosU1 = cos($U1);
	$sinU2 = sin($U2);
	$cosU2 = cos($U2);
	$lamuda = $omega;
	$lamuda0 = 0;
	
	for ($i = 0; $i < $imax; $i++)
	{
		$sinSigma = sqrt(pow($cosU2 * sin($lamuda), 2) + pow($cosU1 * $sinU2 - $sinU1 * $cosU2 * cos($lamuda), 2));
		$cosSigma = $sinU1 * $sinU2 + $cosU1 * $cosU2 * cos($lamuda);
		$sigma = atan2($sinSigma, $cosSigma);
		
		$sinAlfa = $cosU1 * $cosU2 * sin($lamuda) / $sinSigma;
		$cosAlfa2 = 1 - pow($sinAlfa, 2);
		if ($cosAlfa2 == 0)
			return 0;
		
		$cos2Sigmam = $cosSigma - (2 * $sinU1 * $sinU2 / $cosAlfa2);
		$C = ($f / 16.0) * $cosAlfa2 * (4 + $f * (4 - 3 * $cosAlfa2));
		$lamuda0 = $lamuda;
		$lamuda = $omega + (1 - $C) * $f * $sinAlfa * ($sigma + $C * $sinSigma * ($cos2Sigmam + $C * $cosSigma * (-1 + 2 * pow($cos2Sigmam, 2))));
		
		if (abs($lamuda - $lamuda0) <= 1.0e-12)
			break;
	}
	
	$uu = $cosAlfa2 * ($a * $a - $b * $b) / ($b * $b);
	$A = 1 + ($uu / 16384.0) * (4096 + $uu * (-768 + $uu * (320 - 175 * $uu)));
	$B = ($uu / 1024.0) * (256 + $uu * (-128 + $uu * (74 - 47 * $uu)));
	$dSigma = $B * $sinSigma * ($cos2Sigmam + ($B / 4.0) * ($cosSigma * (-1 + 2 * pow($cos2Sigmam, 2)) - ($B / 6.0) * $cos2Sigmam * (-3 + 4 * pow($sinSigma, 2)) * (-3 + 4 * pow($cos2Sigmam, 2))));
	$s = $b * $A * ($sigma - $dSigma);
	$az1 = atan2(($cosU2 * sin($lamuda)), ($cosU1 * $sinU2 - $sinU1 * $cosU2 * cos($lamuda)));
	$az2 = atan2(($cosU1 * sin($lamuda)), (-($sinU1) * $cosU2 + $cosU1 * $sinU2 * cos($lamuda))) + pi();
	if ($az1 < 0)
		$az1 = $az1 + pi() * 2.0;
	if ($i >= $imax)
		$i = $imax - 1;
	
	return array(
		"d" => floor($s * 1000.0) / 1000.0,
		"azimuth1" => $az1,
		"azimuth2" => $az2
	);
}

function pre_dump($print)
{
	print "<pre>";
	var_dump($print);
	print "</pre>";
}

function getWaveH($ms, $length)
{
	if (is_nan($ms) || is_nan($length))
		return 0.0;
	$fn =  $ms / sqrt(9.8 * $length);
	return (0.9 * pow($fn, 3.5) * $length);
}

function distance2Sekou($lon, $lat)
{
	$SEKOU1 = array(136.825722, 34.993696);
	$SEKOU2 = array(136.819583, 34.997795);
	$SEKOU3 = array(136.801773, 35.009234);
	
	$ret = array();
	$tmp = vincentyInverse($SEKOU1[0], $SEKOU1[1], $lon, $lat, 10);
	$ret["0"] = $tmp["d"];
	
	$tmp = vincentyInverse($SEKOU2[0], $SEKOU2[1], $lon, $lat, 10);
	$ret["1"] = $tmp["d"];
	
	$tmp = vincentyInverse($SEKOU3[0], $SEKOU3[1], $lon, $lat, 10);
	$ret["2"] = $tmp["d"];
	return $ret;
}

//================================================
// main
//================================================
$slsec = $argv[1];
sleep($slsec);
//print "aaaaa";

//
// const
//
$DIR = "/www/ttest/nagoya/period/lst/";
$WSEC = array(300.0, 180.0);
$WARNING_WAVE_HEIGHT = 0.8;

//
// AIS情報を取得
//
$body = file_get_contents("http://www4.ptgw3.net/078_nagoya/term_out/ais_out_ict.php");
$body = explode("\\n", $body);

$ofile = "/www/ttest/nagoya/output/". date("Ymd"). ".csv";
if (file_exists($ofile))
{
	$ofp = fopen($ofile, "a");
	if ($ofp == NULL)
		return;
}
else
{
	$ofp = fopen($ofile, "w");
	if ($ofp == NULL)
		return;
	fputs($ofp, mb_convert_encoding("警報,mmsi,船名,緯度,経度,速度[kt],速度[m/s],方位角[°],波高[m],船の長さ[m]", "Shift-JIS", "UTF-8"). "\r\n");
}

//pre_dump($body);
$cnt = count($body);
for($i = 0; $i < $cnt; $i++)
{
	$body[$i] = trim($body[$i]);
	if (strcmp($body[$i], "") == 0)
		continue;
	$datas = explode(",", $body[$i]);
	
	//pre_dump($datas);
	$mmsi = $datas[0];
	$shipLength = $datas[12];
	$shipWidth = $datas[13];
	$sog = $datas[4];
	// cogではなくhdgを使用
	//$azimuth = $datas[5];
	$azimuth = $datas[6];
	$lon = $datas[8];
	$lat = $datas[9];
	$shipType = $datas[11];
	$shipName = $datas[10];
	
	// 一つ前のデータを読込
	$lsfile = $DIR . $mmsi. ".txt";
	if (file_exists($lsfile))
	{
		$tmp = file($lsfile, FILE_SKIP_EMPTY_LINES);
		$tmp = explode(",", $tmp[0]);
		//pre_dump($tmp);
		$lsdatas = array();
		$lsdatas["lon"] = $tmp[8];
		$lsdatas["lat"] = $tmp[9];
		
		if ($lon == $lsdatas["lon"] && $lat == $lsdatas["lat"])
		{
			continue;
		}
		else
		{
			$fp = fopen($lsfile, "w");
			if ($fp == NULL)
				continue;
			fwrite($fp, $body[$i]);
			fclose($fp);
		}
	}
	else
	{
		$lsdatas = NULL;
		$fp = fopen($lsfile, "w");
		if ($fp == NULL)
			continue;
		fwrite($fp, $body[$i]);
		fclose($fp);
	}
	
	// 航跡波の到達箇所を計算
	$retR = vincentyDirect($lon, $lat, deg2rad($azimuth + 180 - 19.46), 500, 50);
	
	if ($sog >= 0.1 && $sog < 100)
	{
		$ms = $sog * 1852 / 3600;
		$waveH = getWaveH($ms, $shipLength);
		if (!is_nan($waveH))
		{
			if ($waveH >= $WARNING_WAVE_HEIGHT)
			{
				$sekouDis["ps"] = distance2Sekou($lon, $lat);
				// 施工箇所までの距離を計算
				$sekouDis["nw"] = distance2Sekou($retR["lon"], $retR["lat"]);
				$clwave = 0;
				if ($lsdatas == NULL)
				{
					// 前回のデータなし
					$clwave = 1;
				}
				else
				{
					// 前回データの施工箇所までの距離を計算
					//$sekouDis["ls"] = distance2Sekou($lsdatas["lon"], $lsdatas["lat"]);
					//$lsL = vincentyDirect($lsdatas["lon"], $lsdatas["lat"], deg2rad($azimuth + 180 + 19.46), 500, 50);
					$lsR = vincentyDirect($lsdatas["lon"], $lsdatas["lat"], deg2rad($azimuth + 180 - 19.46), 500, 50);
					$sekouDis["ls"] = distance2Sekou($lsR["lon"], $lsR["lat"]);
					
					// 前回の施工箇所の距離と比較
					if (($sekouDis["ls"]["0"] - $sekouDis["nw"]["0"]) > 0 || ($sekouDis["ls"]["1"] - $sekouDis["nw"]["1"]) > 0 || ($sekouDis["ls"]["2"] - $sekouDis["nw"]["2"]) > 0)
					{
						$clwave = 1;
					}
				}
				if ($clwave == 1)
				{
					$arsec = array();
					$arsec[0] = $sekouDis["ps"]["0"] / $ms;
					$arsec[1] = $sekouDis["ps"]["1"] / $ms;
					$arsec[2] = $sekouDis["ps"]["2"] / $ms;
					
					//print "距離<br>";
					//pre_dump($sekouDis);
					//pre_dump($ms);
					//print "秒数<br>";
					//pre_dump($arsec);
					
					if ($arsec[0] <= $WSEC[0] || $arsec[1] <= $WSEC[0] || $arsec[2] <= $WSEC[0]) // 到達までの時間[s]は300秒以下
					{
						if ($arsec[0] <= $WSEC[1] || $arsec[1] <= $WSEC[1] || $arsec[2] <= $WSEC[1]) // 到達までの時間[s]は180秒以下
						{
							print $shipName. ",red<br>";
							$op = "2,". $mmsi. ",". $shipName. ",". $lat. ",". $lon. ",". $sog. ",". $ms. ",". $azimuth. ",". $waveH. ",". $shipLength;
							fputs($ofp, $op. "\r\n");
						}
						else
						{
							// 警告黄色
							print $shipName. ",yellow<br>";
							$op = "1,". $mmsi. ",". $shipName. ",". $lat. ",". $lon. ",". $sog. ",". $ms. ",". $azimuth. ",". $waveH. ",". $shipLength;
							fputs($ofp, $op. "\r\n");
						}
					}
				}
			}
		}
	}
}
fclose($ofp);

?>