目次
サンプル・プログラムの実行例
YOLPコンテンツジオコーダAPI（Yahoo!JAPAN）は、世界のキーワードから緯度・経度を求める機能が弱い。たとえば「サンティアゴ」でキーワード検索すると、Google Geocoding API チリの首都を第一候補として返すが、YOLPコンテンツジオコーダAPIはドミニカ共和国の県都しか返さない。
そこで、地図上のマーカーをドラッグし、計算ボタンを押すことで射程を再描画する機能を採り入れた。これは Googleマップ でも利用できる。
サンプル・プログラム
| ballistic.php | サンプル・プログラム本体 |
| pahooGeoCode.php | 住所・緯度・経度に関わるクラス pahooGeoCode。 使い方は「PHPで住所・ランドマークから最寄り駅を求める」「PHPで住所・ランドマークから緯度・経度を求める」などを参照。include_path が通ったディレクトリに配置すること。 |
弾道方程式
末尾に、参考サイトや参考書籍を掲げた。ぜひご覧いただきたい。
発射点をa、着弾点をbとすると、地表面に沿った円弧acbを射程（range）と呼ぶ。acbの左側は地中に潜っていると考えてほしい。
射程を、地球中心から見た角度φを半射程角（semi-range angle）と呼ぶ。射程をS、地球の半径をR（6,378,136m）とすると、φは次の式によって求められる。
\[ \phi = \frac{S}{2R} \]
ミサイルの最終段のロケットモーターが燃え切ったときの速度比をvoとすると、次の式によって求められる。
\[
\displaystyle
\begin{align}
v_o &= \sqrt{\frac{-B+\sqrt{B^2-AC}}{A}} \\
A &= \sin^2\gamma \ (\cos^2\phi - \sin^2\gamma) \\
B &= \sin^2\gamma \ (1 - \cos^2\phi) \\
B &= \cos^2\phi - 1
\end{align}
\]
楕円の離心率e、長径a、短径bは、次の式によって求められる。
\displaystyle
\begin{align}
e &= \sqrt{1-{v_o}^2 (2 -{v_o}^2) \sin^2\gamma} \\
a &= \frac{R}{2 - {v_o}^2} \\
b &= a \sqrt{1 - e^2}
\end{align}
\]
地上から測った軌道の最高点の高さ（軌道高）hmは、次の式によって求められる。
軌道上の速度は変化する。地球中心に対する軌道上の角度をθ、地上の引力係数をμとすると、θにおける動径rと、動径上の速度vは、次の式によって求められる。
\displaystyle
\begin{align}
\mu &= 3.986044\times10^{14}m^3s^{-2} \\
r &= \frac{a (1 - e^2)}{1 - e \cdot \cos \theta} \\
v &= \sqrt{\mu (\frac{2}{r} - \frac{1}{a}})
\end{align}
\]
準備：pahooGeoCode クラス
pahooGeoCode.php
37: class pahooGeoCode {
38: var $items; // 検索結果格納用
39: var $error; // エラー・フラグ
40: var $errmsg; // エラー・メッセージ
41: var $hits; // 検索ヒット件数
42: var $webapi; // 直前に呼び出したWebAPI URL
43:
44: // Google Cloud Platform APIキー
45: // https://cloud.google.com/maps-platform/
46: // ※Google Maps APIを利用しないのなら登録不要
47: var $GOOGLE_API_KEY_1 = '**************************'; // HTTPリファラ用
48: var $GOOGLE_API_KEY_2 = '**************************'; // IP制限用
49: var $GOOGLE_MAP_ID = '*************************'; // GoogleMaps ID
50:
51: // Yahoo! JAPAN Webサービス アプリケーションID
52: // https://e.developer.yahoo.co.jp/register
53: // ※Yahoo! JAPAN Webサービスを利用しないのなら登録不要
54: var $YAHOO_APPLICATION_ID = '*****************************';
クラスについては「PHPでクラスを使ってテキストの読みやすさを調べる」を参照されたい。
地図や住所検索として Google を利用するのであれば、Google Cloud Platform APIキー が必要で、その入手方法は「Google Cloud Platform - WebAPIの登録方法」を、Yahoo!JAPAN を利用するのであれば、Yahoo! JAPAN Webサービス アプリケーションIDが必要で、その入手方法は「Yahoo!JAPAN デベロッパーネットワーク - WebAPIの登録方法」を、それぞれ参照されたい。
準備：地図サービスの選択
ballistic.php
34: //地図描画サービスの選択
35: // 0：Google
36: // 2：地理院地図・OSM
37: define('MAPSERVICE', 0);
38:
39: //住所検索サービスの選択
40: // 0：Google
41: // 1：Yahoo!JAPAN
42: // 11：HeartRails Geo API
43: // 12：OSM Nominatim Search API
44: define('GEOSERVICE', 11);
45:
46: //逆ジオコーディングサービスの選択
47: // 0：Google
48: // 1：Yahoo!JAPAN
49: // 11：HeartRails Geo API
50: // 21：簡易ジオコーディングサービス
51: define('REVGEOSERVICE', 1);
住所検索サービスは、Google、Yahoo!JAPAN、HeartRails Geo API、OSM Nominatim Search API から選べる。あらかじめ、定数 GEOSERVICE に値を設定すること。
逆ジオコーディングサービスは、Google、Yahoo!JAPAN、HeartRails Geo API、簡易ジオコーディングサービスから選べる。あらかじめ、定数 REVGEOSERVICE に値を設定すること。
解説：初期値
ballistic.php
53: //マップの表示サイズ（単位：ピクセル）
54: define('MAP_WIDTH', 600);
55: define('MAP_HEIGHT', 400);
56: //マップID
57: define('MAPID', 'map_id');
58:
59: //初期値
60: define('DEF_LONGITUDE0', 165.0); //中心座標（経度）
61: define('DEF_LATITUDE0', 30.5); // （緯度）
62: define('DEF_TYPE', 'roadmap'); //マップタイプ
63: define('DEF_ZOOM', 4); //ズーム
64: define('DEF_CATEGORY', 'world'); //カテゴリ
65:
66: //発射点
67: define('DEF_FROM', '東倉里ミサイル発射場'); //地名
68: define('DEF_CAT_FROM', 'world'); //カテゴリ
69: define('DEF_LATITUDE1', 39.660075); //緯度
70: define('DEF_LONGITUDE1', 124.705294); //経度
71:
72: //着弾点
73: define('DEF_TO', 'ホノルル'); //地名
74: define('DEF_CAT_TO', 'world'); //カテゴリ
75: define('DEF_LATITUDE2', 21.308889); //経度
76: define('DEF_LONGITUDE2', -157.826111); //経度
77:
78: define('DEF_HM', 1500); //弾道高（km）
初回起動時に表示する発射点は、北朝鮮の東倉里ミサイル発射場、着弾地点はハワイ・ホノルルとしてある。弾道高は1500kmである。
解説：マーカー設定など
解説：射程，弾道高から軌道要素を求める
ballistic.php
502: /**
503: * 射程，射角から軌道要素を求める
504: * @param double $dt 射程（メートル）
505: * @param double $aof 射角（度）
506: * @param array $elements 軌道要素を格納する配列
507: * @return なし
508: */
509: function getOrbitalElementSub($dt, $aof, &$elements) {
510: static $r = 6378136; //地球半径（メートル）
511:
512: //ラジアン変換
513: $aof = deg2rad(90 - $aof);
514:
515: $sra = $dt / (2 * $r); //半射程角
516:
517: $ax = pow(sin($aof), 2) * (pow(cos($sra), 2) - pow(sin($aof), 2));
518: $bx = pow(sin($aof), 2) * (1 - pow(cos($sra), 2));
519: $cx = pow(cos($sra), 2) - 1;
520:
521: $vo = sqrt((-$bx + sqrt(pow($bx, 2) - $ax * $cx)) / $ax); //燃え切り速度比
522: $e = sqrt(1 - pow($vo, 2) * (2 - pow($vo, 2)) * pow(sin($aof), 2)); //離心率
523: $a = $r / (2 - pow($vo, 2)); //長径
524: $b = $a * sqrt(1 - pow($e, 2)); //短径
525: $hmm = $a * (1 + $e) - $r; //弾道高
526:
527: $elements['sra'] = $sra;
528: $elements['vo'] = $vo;
529: $elements['e'] = $e;
530: $elements['a'] = $a;
531: $elements['b'] = $b;
532: $elements['hmm'] = $hmm;
533: }
ballistic.php
535: /**
536: * 射程，弾道高から軌道要素を求める
537: * @param double $dt 射程（メートル）
538: * @param double $hmm 弾道高（メートル）
539: * @param array $elements 軌道要素を格納する配列
540: * @return なし
541: */
542: function getOrbitalElement($dt, $hmm, &$elements) {
543: $flag = TRUE;
544: $aof = 0.0; //射角（度）
545:
546: //誤差0.1度で射角を求める
547: while ($flag && ($aof < 90.0)) {
548: $aof += 1.0;
549: getOrbitalElementSub($dt, $aof, $elements);
550: if ($elements['hmm'] == $hmm) {
551: $flag = FALSE;
552: break;
553: } else if ($elements['hmm'] > $hmm) {
554: $aof -= 1.0;
555: while ($flag && ($aof < 90.0)) {
556: $aof += 0.1;
557: getOrbitalElementSub($dt, $aof, $elements);
558: if ($elements['hmm'] >= $hmm) {
559: $flag = FALSE;
560: break;
561: }
562: }
563: }
564: }
565: $elements['dt'] = $dt;
566: $elements['aof'] = $aof;
567: }
そこで本プログラムでは、射角の代わりに弾道高を入力して、軌道要素を求める関数 getOrbitalElement を用意した。
弾道高から射角を逆算することは難しいので、射角を0.1度ずつ増やし、ユーザー関数 getOrbitalElementSub を呼び出し、与えられた弾道高にほぼ近づいたら計算を打ち切るという方法で軌道要素を求めている。
計算回数を減らすため、初段は1度ずつ射角を増やし、近づいてきたら0.1度ずつ増やすようにしている。
解説：軌道要素から、飛翔経路長、速度を求める
ballistic.php
470: /**
471: * 軌道要素から、飛翔経路長、速度を求める
472: * @param array $elements 軌道要素
473: * @return なし
474: */
475: function getBallistic(&$elements) {
476: static $u = 3.98602e+14; //地心引力係数
477: $n = 100; //分割数
478:
479: $dw = $elements['sra'] / $n;
480: $ds = 0.0; //飛翔経路長
481: $vs = 0.0;
482: $vmax = 0.0; //最大速度
483: for ($i = 0; $i <= 100; $i++) {
484: $r = $elements['a'] * (1.0 - pow($elements['e'], 2)) / (1.0 - $elements['e'] * cos($dw * $i));
485: $v = sqrt($u * ((2.0 / $r) - (1.0 / $elements['a'])));
486: if ($v > $vmax) $vmax = $v;
487: $vs += $v;
488: if ($i == 0) $tmp1 = $r;
489: else $tmp2 = $r;
490: if ($i > 1) {
491: $dx = $tmp2 * sin($dw);
492: $dy = $tmp1 - $tmp2 * cos($dw);
493: $ds += sqrt(pow($dx, 2) + pow($dy, 2));
494: $tmp1 = $tmp2;
495: }
496: }
497: $elements['ds'] = $ds * 2;
498: $elements['vm'] = $vs / ($n + 1);
499: $elements['vmax'] = $vmax;
500: }
そこで、「弾道方程式」で紹介した動径を使って、動径角を100分割して、それを足していくことで飛翔経路長を近似している。
同時に速度を100回求め、平均速度と最高速度を計算している。
解説：軌道要素の計算
ballistic.php
865: //軌道要素の計算
866: if (($hm > 0) && ($items['longitude1'] != '') && ($items['latitude1'] != '') && ($items['longitude2'] != '') && ($items['latitude2'] != '')) {
867: $dt = $pgc1->distance($items['longitude1'], $items['latitude1'], $items['longitude2'], $items['latitude2']);
868: $items['hmm'] = $hm * 1000;
869: getOrbitalElement($dt, $hm * 1000, $items);
870: getBallistic($items);
871: }
872:
873: //地図描画サービス
874: switch (MAPSERVICE) {
875: //Google
876: case 0:
877: $call = jsDragMarker_Gmap($items2[$i]['latitude'], $items2[$i]['longitude']);
878: $call .= jsDrawLine_Gmap(1, $items1, 'rgb(255, 0, 0)', 1);
879: $call .= jsDrawLine_Gmap(2, $items2, 'rgb(0, 0, 0)', 1);
880: $zd = 1;
881: break;
882: //Yahoo!JAPAN
883: case 1:
884: $call = jsDragMarker_YOLPmap($items2[$i]['latitude'], $items2[$i]['longitude']);
885: $call .= jsDrawLine_YOLPmap(1, $items1, 'FF0000', 1);
886: $call .= jsDrawLine_YOLPmap(2, $items2, '000000', 1);
887: $zd = 0;
888: break;
889: //地理院地図・OSM
890: case 2:
891: $call = jsDragMarker_Leaflet($items2[$i]['latitude'], $items2[$i]['longitude']);
892: $call .= jsDrawLine_Leaflet(1, $items1, 'FF0000', 1);
893: $call .= jsDrawLine_Leaflet(2, $items2, '000000', 1);
894: $zd = 0;
895: break;
896: default:
897: $call = '';
898: $zd = 0;
899: break;
900: }
901: $lat = $items2[$i]['latitude'];
902: $lng = $items2[$i]['longitude'];
903:
904: //エラー発生時
905: } else {
結語
参考サイト
- PHPで大圏航路を描く：ぱふぅ家のホームページ
- PHPで2地点間の直線距離を求める：ぱふぅ家のホームページ
参考書籍
|
惑星探査機の軌道計算入門 | ||
| 著者 | 半揚稔雄 | ||
| 出版社 | 日本評論社 | ||
| サイズ | 単行本 |
|
|
| 発売日 | 2017年09月19日頃 | ||
| 価格 | 2,420円（税込） | ||
| ISBN | 9784535788459 | ||
| 序章 宇宙飛翔力学とは 第1章 円錐曲線の幾何学 第2章 ケプラー運動 第3章 軌道遷移 第4章 軌道決定 第5章 ロケットの性能 第6章 惑星間飛行 付録 | |||
|
天文計算入門 | ||
| 著者 | 長谷川一郎 | ||
| 出版社 | 恒星社厚生閣 | ||
| サイズ | 単行本 |
|
|
| 発売日 | 1996年01月25日頃 | ||
| 価格 | 2,750円（税込） | ||
| ISBN | 9784769908180 | ||
そこで今回は、マップ上で発射地点と着弾地点を指定し、弾道の最高点の高さを入力するだけで、軌道要素と弾道を計算するPHPプログラムをつくってみることにする。
ユーザーインターフェースは、「PHPで大圏航路を描く」を流用した。
なお、作者の力不足から、コリオリの力、多段ロケットによる加速、弾頭の形状・材質、地球大気による減速など、物理的な諸条件は含めていない。大まかな検証をすることが目的だと考えていただきたい。
（2021年3月28日）PHP8対応，リファラ・チェック改良。
（2020年3月29日）地理院地図・OSMに対応した。OSM Nominatim Search APIを利用できるようにした。