Up 自転球体上移動方程式の計算プログラム 作成: 2022-09-15
更新: 2022-10-07


    「自転球体上の移動の方程式」で示した計算手順をプログラムにした (言語 : PHP):

    < ? /****** constants ***************************************************************/ /**** 自転球の半径と角速度 ************/ $R = 6370000; $Omega = 0.00007292; /**** 区分求積間隔 : 5分 (300秒) ************/ $Delta_t = 300; /****** variables *******************************************************************/ /* _fx, _fy, _fx は固定座標,_x, _y, _z は (P,v)座標 */ /**** 位置 (固定座標) ************/ $P_fx = 0; $P_fy = 0; $P_fz = 0; /**** 緯度・経度 (固定座標) ************/ $lattitude = 0; $longitude = 0; /**** 速度 (固定座標) ************/ $v_fx = 0; $v_fy = 0; $v_fz = 0; /**** (P,v)-座標系 ***********/ $Q_fx = 0; $Q_fy = 0; $alpha = 0; /**** 位置 (((P,v)座標) ************/ $P_x = 0; $P_y = 0; $P_z = 0; /**** 速度 (((P,v)座標) ************/ $v_x = 0; $v_y = 0; $v_z = 0; /**** 加速度 (((P,v)座標) ************/ $a_x = 0; $a_y = 0; $a_z = 0; /**** 加速度×時間 (((P,v)座標) ************/ $at_x = 0; $at_y = 0; $at_z = 0; /**** 元速度+加速分 (((P,v)座標) ************/ $v1_x = 0; $v1_y = 0; $v1_z = 0; /**** 元位置+移動距離 (((P,v)座標) ************/ $P1_x = 0; $P1_y = 0; $P1_z = 0; /**** 元位置+移動距離 (固定座標) ************/ $P1_fx = 0; $P1_fy = 0; $P1_fz = 0; /**** 元速度+加速分 (固定座標) ************/ $v1_fx = 0; $v1_fy = 0; $v1_fz = 0; /**** 表タイトル ************/ $subject = ""; /**** 追跡 flag ************/ $flag = ""; $counter = 0; /****** functions ******************************************************************/ /**** 緯度・経度 ************/ function lattitude_longitude() { global $R; global $P_fx, $P_fy, $P_fz; global $latitude, $longitude; $latitude = rad2deg( asin( $P_fz / $R ) ); $longitude = rad2deg( asin( $P_fy / $R ) ); } /**** 固定座標x軸と (P,v)座標軸xの開き (角度) ************/ function alpha() { global $R; global $P_fx, $P_fy, $P_fz; global $v_fx, $v_fy, $v_fz; global $alpha; global $Q_fy; if( $P_fz == 0 ){ $alpha = asin( $P_fy / $R ); $Q_fy = $P_fy; } else { $R2 = $R * $R; $A = $P_fx * $P_fx * $v_fy - $P_fx * $P_fy * $v_fx - $R2 * $v_fy; $B = $P_fx * $P_fy * $v_fy - $P_fy * $P_fy * $v_fx + $R2 * $v_fx; $C = $A / sqrt( $A * $A + $B * $B ); if( $P_fz > 0 ){ $alpha = asin( $C ); $Q_fy = $R * $C; } else { $alpha = - asin( $C ); $Q_fy = - $R * $C; } } } /**** 位置座標の変換 (固定座標 → ((P,v)座標) ************/ function coordinates_p() { global $R; global $P_fx, $P_fy, $P_fz; global $P_x, $P_y, $P_z; global $alpha; $P_x = $P_fx * cos( $alpha ) + $P_fy * sin( $alpha ); $P_y = - $P_fx * sin( $alpha ) + $P_fy * cos( $alpha ); $P_z = $P_fz; } /**** 速度座標の変換 (固定座標 → ((P,v)座標) ************/ function coordinates_v() { global $v_fx, $v_fy, $v_fz; global $v_x, $v_y, $v_z; global $alpha; $v_x = $v_fx * cos( $alpha ) + $v_fy * sin( $alpha ); $v_y = - $v_fx * sin( $alpha ) + $v_fy * cos( $alpha ); $v_z = $v_fz; } /**** 加速度 (((P,v)座標) ************/ function a() { global $R, $Omega; global $P_x, $P_y, $P_z; global $v_x, $v_y, $v_z; global $a_x, $a_y, $a_z; global $alpha; $R2 = $R * $R; $Omega2 = $Omega * $Omega; $v2 = $v_x * $v_x + $v_y * $v_y + $v_z * $v_z; $v = sqrt ( $v_x * $v_x + $v_y * $v_y + $v_z * $v_z ); /**** 直進大円が赤道 ****/ if( ( $P_y == 0 ) && ( $P_z == 0 ) && ( $v_x == 0 ) && ( $v_z == 0 ) ){ $a_x = - $v2 / $R + $R * $Omega2; $a_y = 0; $a_z = 0; } /**** 直進大円が経線 ****/ elseif( ( $P_y == 0 ) && ( $v_y == 0 ) ){ $a_x = - $v2 * $P_x / $R2 + $Omega2 * $P_x; $a_y = 0; $a_z = - $v2 * $P_z / $R2; } /**** 位置が赤道上 ****/ elseif( ( $P_y == 0 ) && ( $P_z == 0 ) ){ $a_x = - $v2 / $R + $R * $Omega2; $a_y = 0; $a_z = 0; } /**** それ以外 ****/ else { $a_x = - $v2 * $P_x / $R2 + $Omega2 * $P_x; $a_y = - $v2 * $P_y / $R2 + $Omega2 * $P_y; $a_z = - $v2 * $P_z / $R2; } } /**** 加速度×時間 (((P,v)座標) ************/ function at() { global $Delta_t; global $a_x, $a_y, $a_z; global $at_x, $at_y, $at_z; $at_x = $a_x * $Delta_t; $at_y = $a_y * $Delta_t; $at_z = $a_z * $Delta_t; } /**** 元位置+移動距離 (((P,v)座標) ************/ function p1() { global $R, $Delta_t; global $P_x, $P_y, $P_z; global $P1_x, $P1_y, $P1_z; global $v_x, $v_y, $v_z; $v2 = $v_x * $v_x + $v_y * $v_y + $v_z * $v_z; $v = sqrt ( $v2 ); $theta = $v * $Delta_t / $R; $cos_theta = cos( $theta ); $sin_theta = sin( $theta ); /**** 直進大円が赤道 ****/ if( ( $P_y == 0 ) && ( $P_z == 0 ) && ( $v_x == 0 ) && ( $v_z == 0 ) ){ $P1_x = $R * $cos_theta; $P1_y = $R * $sin_theta; $P1_z = 0; } /**** 直進大円が経線 ****/ elseif( ( $P_y == 0 ) && ( $v_y == 0 ) ){ $P1_x = $P_x * $cos_theta - $P_z * $sin_theta; $P1_y = 0; $P1_z = $P_z * $cos_theta + $P_x * $sin_theta; } /**** 位置が赤道上 ****/ elseif( ( $P_y == 0 ) && ( $P_z == 0 ) ){ $P1_x = $R * $cos_theta; $P1_y = $R * ( $v_y / $v ) * $sin_theta; $P1_z = $R * ( $v_z / $v ) * $sin_theta; } /**** それ以外 ****/ else { $A = sqrt( pow( $R, 2 ) - pow( $P_x, 2 ) ); $P1_x = $P_x * $cos_theta - $A * $sin_theta; $P1_y = $P_y * $cos_theta + $P_x * $P_y * $sin_theta / $A; $P1_z = $P_z * $cos_theta + $P_z * $P_x * $sin_theta / $A; } } /**** 元速度+加速分 (((P,v)座標) ************/ function v1() { global $R; global $v_x, $v_y, $v_z; global $v1_x, $v1_y, $v1_z; global $P1_x, $P1_y, $P1_z; global $at_x, $at_y, $at_z; $vx = $v_x + $at_x; $vy = $v_y + $at_y; $vz = $v_z + $at_z; /* $v = sqrt( pow( $vx, 2 ) + pow( $vy, 2 ) + pow( $vz, 2 ) ); $A = ( $P1_x * $vx + $P1_y * $vy + $P1_z * $vz ) / ( $R * $v ); $theta = asin( $A ); $cos_theta = cos( $theta ); */ $cos_theta = 1; $v1_x = $vx * $cos_theta; $v1_y = $vy * $cos_theta; $v1_z = $vz * $cos_theta; } /**** 位置座標の変換 (((P,v)座標 → 固定座標) ************/ function coordinates_p1f() { global $R, $Omega, $Delta_t; global $P1_x, $P1_y, $P1_z; global $P1_fx, $P1_fy, $P1_fz; global $alpha; $P1_fx = $P1_x * cos( - $alpha ) + $P1_y * sin( - $alpha ); $P1_fy = - $P1_x * sin( - $alpha ) + $P1_y * cos( - $alpha ); $P1_fz = $P1_z; } /**** 速度座標の変換 (((P,v)座標 → 固定座標) ************/ function coordinates_v1f() { global $R, $Omega, $Delta_t; global $v1_x, $v1_y, $v1_z; global $v1_fx, $v1_fy, $v1_fz; global $alpha; $v1_fx = $v1_x * cos( - $alpha ) + $v1_y * sin( - $alpha ); $v1_fy = - $v1_x * sin( - $alpha ) + $v1_y * cos( - $alpha ); $v1_fz = $v1_z; } /**** 表ヘッダ ************/ function table_head() { global $subject; print " ".$subject." <table cellpadding=2><tr><td></td></tr></table> <table border=1 cellpadding=2 cellspacing=0> <tr> <td align=center rowspan=3>経過時間</td> <td align=center colspan=8>固定座標</td> <td align=center rowspan=3 bgcolor=\"#ccccff\">速度 v</td> <td align=center colspan=2>(P,v)-座標枠</td> <td align=center colspan=15>(P,v)-座標</td> </tr> <tr> <td align=center rowspan=2 colspan=2 bgcolor=\"#eeee00\">緯度・経度</td> <td colspan=3 align=center bgcolor=\"#eeeeee\">位置 P</td> <td colspan=3 align=center>速度 vf</td> <td align=center rowspan=2 bgcolor=\"#ccffcc\">alpha</td> <td align=center bgcolor=\"#ccffcc\">Q</td> <td colspan=3 align=center bgcolor=\"#eeeeee\">位置 P</td> <td colspan=3 align=center>速度 v</td> <td colspan=3 align=center bgcolor=\"#ffcccc\">加速 at</td> <td colspan=3 align=center>元速度+加速 = v1</td> <td colspan=3 align=center bgcolor=\"#eeeeee\">元位置+移動= P1</td> </tr> <tr> <td align=center bgcolor=\"#eeeeee\">fx</td> <td align=center bgcolor=\"#eeeeee\">fy</td> <td align=center bgcolor=\"#eeeeee\">fz</td> <td align=center>fx</td> <td align=center>fy</td> <td align=center>fz</td> <td align=center bgcolor=\"#ccffcc\">fy</td> <td align=center bgcolor=\"#eeeeee\">x</td> <td align=center bgcolor=\"#eeeeee\">y</td> <td align=center bgcolor=\"#eeeeee\">z</td> <td align=center>x</td> <td align=center>y</td> <td align=center>z</td> <td align=center bgcolor=\"#ffcccc\">x</td> <td align=center bgcolor=\"#ffcccc\">y</td> <td align=center bgcolor=\"#ffcccc\">z</td> <td align=center>x</td> <td align=center>y</td> <td align=center>z</td> <td align=center bgcolor=\"#eeeeee\">x</td> <td align=center bgcolor=\"#eeeeee\">y</td> <td align=center bgcolor=\"#eeeeee\">z</td> </tr> "; } /**** 表末尾 ************/ function table_tail() { print " </table> <br> "; } /**** 表作成 *****************************************************************************************/ function table_make() { global $R; global $latitude, $longitude; global $P_fx, $P_fy, $P_fz; global $v_fx, $v_fy, $v_fz; global $alpha; global $Q_fy; global $P_x, $P_y, $P_z; global $v_x, $v_y, $v_z; global $at_x, $at_y, $at_z; global $P1_x, $P1_y, $P1_z; global $v1_x, $v1_y, $v1_z; global $P1_fx, $P1_fy, $P1_fz; global $v1_fx, $v1_fy, $v1_fz; global $counter; table_head(); for( $i = 0; $i < 500; $i++ ){ lattitude_longitude(); if( $latitude > 62 ){ break; } print " <tr> <td align=center>".$i."</td> <td align=left bgcolor=\"#eeee00\">".round( $latitude, 3 )."</td> <td align=left bgcolor=\"#eeee00\">".round( $longitude, 3 )."</td> <td align=left bgcolor=\"#eeeeee\">".round( $P_fx, 0 )."</td> <td align=left bgcolor=\"#eeeeee\">".round( $P_fy, 0 )."</td> <td align=left bgcolor=\"#eeeeee\">".round( $P_fz, 0 )."</td> <td align=left>".round( $v_fx, 2 )."</td> <td align=left>".round( $v_fy, 2 )."</td> <td align=left>".round( $v_fz, 2 )."</td> "; $v = sqrt( $v_fx * $v_fx + $v_fy * $v_fy + $v_fz * $v_fz ); print " <td align=left bgcolor=\"#ccccff\">".round( $v, 2 )."</td> "; alpha(); print " <td align=left bgcolor=\"#ccffcc\">".round( $alpha, 3 )."</td> <td align=left bgcolor=\"#ccffcc\">".round( $Q_fy, 0 )."</td> "; coordinates_p(); coordinates_v(); print " <td align=left bgcolor=\"#eeeeee\">".round( $P_x, 0 )."</td> <td align=left bgcolor=\"#eeeeee\">".round( $P_y, 0 )."</td> <td align=left bgcolor=\"#eeeeee\">".round( $P_z, 0 )."</td> <td align=left>".round( $v_x, 2 )."</td> <td align=left>".round( $v_y, 2 )."</td> <td align=left>".round( $v_z, 2 )."</td> "; a(); at(); print " <td align=left bgcolor=\"#ffcccc\">".round( $at_x, 3 )."</td> <td align=left bgcolor=\"#ffcccc\">".round( $at_y, 3 )."</td> <td align=left bgcolor=\"#ffcccc\">".round( $at_z, 3 )."</td> "; p1(); v1(); print " <td align=left>".round( $v1_x, 2 )."</td> <td align=left>".round( $v1_y, 2 )."</td> <td align=left>".round( $v1_z, 2 )."</td> <td align=left bgcolor=\"#eeeeee\">".round( $P1_x, 0 )."</td> <td align=left bgcolor=\"#eeeeee\">".round( $P1_y, 0 )."</td> <td align=left bgcolor=\"#eeeeee\">".round( $P1_z, 0 )."</td> "; /*** 北半球→南半球の接続が未解決なので,南半球に入ったら終了 ***/ if( $P1_z < 0 ){ break; } coordinates_p1f(); coordinates_v1f(); $P_fx = $P1_fx; $P_fy = $P1_fy; $P_fz = $P1_fz; $v_fx = $v1_fx; $v_fy = $v1_fy; $v_fz = $v1_fz; print " </tr>"; $counter++; } table_tail(); } /**************************************************************************************************************************************/ /**************************************************************************************************************************************/ /** 初期位置 **/ $P_fx = $R; $P_fy = 0; $P_fz = 0; /** 初速 : 50 m/s **/ $v_fx = 0; $v_fy = 0; $v_fz = 50; /** 表タイトル **/ $subject = "赤道から北へ──初速度 ; ".$v_fz." m/s"; table_make(); /********************************************************************/ /** 初期位置 **/ $lat = 50; $theta = pi() * $lat / 180;; $P_fx = $R * cos( $theta ); $P_fy = 0; $P_fz = $R * sin( $theta ); /** 初速 : 50 m/s **/ $v_fx = 0; $v_fy = 50; $v_fz = 0; /**** 区分求積間隔 : 5分 (300秒) ************/ $Delta_t = 300; /** 表タイトル **/ $subject =$subject0. "<br><br>北緯 ".$lat." 度で東へ──初速度 ; ".$v_fy." m/s"; table_make(); ?>