「自転球体上の移動の方程式」で示した計算手順をプログラムにした (言語 : 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."
経過時間 |
固定座標 |
速度 v |
(P,v)-座標枠 |
(P,v)-座標 |
緯度・経度 |
位置 P |
速度 vf |
alpha |
Q |
位置 P |
速度 v |
加速 at |
元速度+加速 = v1 |
元位置+移動= P1 |
fx |
fy |
fz |
fx |
fy |
fz |
fy |
x |
y |
z |
x |
y |
z |
x |
y |
z |
x |
y |
z |
x |
y |
z |
";
}
/**** 表末尾 ************/
function table_tail() {
print "
";
}
/**** 表作成 *****************************************************************************************/
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 "
|
".$i." |
".round( $latitude, 3 )." |
".round( $longitude, 3 )." |
".round( $P_fx, 0 )." |
".round( $P_fy, 0 )." |
".round( $P_fz, 0 )." |
".round( $v_fx, 2 )." |
".round( $v_fy, 2 )." |
".round( $v_fz, 2 )." |
";
$v = sqrt( $v_fx * $v_fx + $v_fy * $v_fy + $v_fz * $v_fz );
print "
".round( $v, 2 )." |
";
alpha();
print "
".round( $alpha, 3 )." |
".round( $Q_fy, 0 )." |
";
coordinates_p();
coordinates_v();
print "
".round( $P_x, 0 )." |
".round( $P_y, 0 )." |
".round( $P_z, 0 )." |
".round( $v_x, 2 )." |
".round( $v_y, 2 )." |
".round( $v_z, 2 )." |
";
a();
at();
print "
".round( $at_x, 3 )." |
".round( $at_y, 3 )." |
".round( $at_z, 3 )." |
";
p1();
v1();
print "
".round( $v1_x, 2 )." |
".round( $v1_y, 2 )." |
".round( $v1_z, 2 )." |
".round( $P1_x, 0 )." |
".round( $P1_y, 0 )." |
".round( $P1_z, 0 )." |
";
/*** 北半球→南半球の接続が未解決なので,南半球に入ったら終了 ***/
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 "
";
$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. "
北緯 ".$lat." 度で東へ──初速度 ; ".$v_fy." m/s";
table_make();
?>
|