File "wave_calc_20140708.php"
Full path: /www/ttest/nagoya/period/Backup/wave_calc_20140708.php
File size: 9.93 KiB (10172 bytes)
MIME-type: text/x-php
Charset: utf-8
<?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);
}
if ($sog >= 0.1 && $sog < 100)
{
$ms = $sog * 1852 / 3600;
$waveH = getWaveH($ms, $shipLength);
if (!is_nan($waveH))
{
if ($waveH >= $WARNING_WAVE_HEIGHT)
{
// 施工箇所までの距離を計算
$sekouDis["nw"] = distance2Sekou($lon, $lat);
$clwave = 0;
if ($lsdatas == NULL)
{
// 前回のデータなし
$clwave = 1;
}
else
{
// 前回データの施工箇所までの距離を計算
$sekouDis["ls"] = distance2Sekou($lsdatas["lon"], $lsdatas["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["nw"]["0"] / $ms;
$arsec[1] = $sekouDis["nw"]["1"] / $ms;
$arsec[2] = $sekouDis["nw"]["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);
?>