!
!satellite_orbital.f90
!
!system of unit 今回の計算で用いている、物理量の単位です。
!length 1AU (1AU = 1.495978707*10^{11} m)
!mass Msun = 1.0
!time 1yr
!Caltesian Coordinate
!The origin is sun. 太陽を原点としています。このため、太陽については軌道を解く必要がありません。
!
!
!
!
Moduleでは、program mainと関数ex2で使う変数や、初期軌道要素を定義しています。
Module constants_set
implicit none
real(kind=8),parameter :: pi=acos(-1.0d0)
real(kind=8),parameter :: GMsun = 4.0d0*pi*pi
!Mass & Orbital Elements
!Earth
real(kind=8), parameter:: M_E = 3.0D0*(10.0D0**(-6.0D0)) !Mass
!Mars
real(kind=8), parameter:: M_Ma = 3.0D0*(10.0D0**(-7.0D0)) !Mass
real(kind=8), parameter:: A_Ma_ini = 1.52368D0 !semi-major axis
real(kind=8), parameter:: e_Ma_ini = 0.09341233D0 !eccentricity
real(kind=8), parameter:: i_Ma_ini = 1.8497D0*pi/180.0D0 !inclination
real(kind=8), parameter:: omega_p_Ma_ini = 336.2075D0*pi/180.0D0 !pericentre
real(kind=8), parameter:: Omega_a_Ma_ini = 49.6198D0*pi/180.0D0 !asceding node
今回の計算では人工衛星は他に比べて圧倒的に軽いので、質量0としています。簡単のためというか、倍精度でも桁が足りないので。
end Module constants_set
program main
use constants_set Moduleを読んでいます
implicit none
integer(kind=4),parameter:: neq = 18 !number of equations, (x, y, z, vx, vy, vz) * 3 方程式の数は、位置速度の6つを地球と目標天体(火星)、人工衛星について解くので、18個です。
integer(kind=4),parameter:: nmax = 10000 !number of step 計算step数です。軌道を記述するのに必要な時間の刻み幅(=h)と、どれくらい長い時間(=T)の軌道を見たいかでT/hで決まりますが、今回は適当。
real(kind=8),parameter:: eps = 1.0d-10 ルンゲクッタの計算で使う微少量。
external ex2 今回は方程式をprogram mainの外で関数として書いている形式にしているので、ここでその関数ex2を宣言します。
integer(kind=4) ierr,i,n,iseed ルンゲクッタの計算で使う。
real(kind=8),dimension(neq,12):: work ルンゲクッタ。
real(kind=8),dimension(neq):: y0,yn !position & velocity 今回解く物理量は、このコードでは全て y で定義しています。なので、この y は3天体の位置と速度です。
real(kind=8) x0,xn,dx !time 今回は時間発展を解くので、xは時間です。
integer(4) counter, n_acc 何回人工衛星を加速させるかを表す量です。
real(8) a_stllt, e_stllt, cosf_stllt それぞれ、人工衛星の軌道長半径、離心率、近点からの角度(=真近点離角)です。人工衛星を太陽に対して一定の位置関係の時に加速するために用います。
counter = 0 以下、初期値の設定と時間刻み幅の設定です。
n_acc = 0
n = 0
x0 = 0.0d0 !x = time
dx = 0.001d0 !timestep 時間刻み幅は、見たい物理を十分見れるようにとります。ここでは衛星が地球を回るだいたい1ヶ月を、だいたい100回で記述するようにしています。
!Timesteps need enough short to describe satellite.
天体の初期条件です。
!initial_condition
!Caltesian coordinate
!Earth (simplified) 地球。今回は簡単化しています。円軌道(離心率0)です。
!position
y0(4) = 1.0D0
y0(5) = 0
y0(6) = 0
!velocity
y0(1) = 0
y0(2) = sqrt(GMsun) 円軌道のとき、1AUの惑星の速度はこうなります。重力と遠心力の釣り合いを考えると、すぐわかります。
y0(3) = 0
!2nd body 2つめの天体。多少まじめに解いています。ちょっと離心率と傾斜角があるだけでこんなに式が煩雑になります。
!x
y0(10) = A_Ma_ini* (1.0D0 - e_Ma_ini)* ( cos(omega_p_Ma_ini)* cos(Omega_a_Ma_ini) &
&- sin(omega_p_Ma_ini)* sin(Omega_a_Ma_ini)* cos(i_Ma_ini) )
y0(11) = A_Ma_ini* (1.0D0 - e_Ma_ini)* ( cos(omega_p_Ma_ini)* sin(Omega_a_Ma_ini) &
&+ sin(omega_p_Ma_ini)* cos(Omega_a_Ma_ini)* cos(i_Ma_ini) )
y0(12) = A_Ma_ini* (1.0D0 - e_Ma_ini)* ( sin(omega_p_Ma_ini)* sin(i_Ma_ini) )
!v
y0(7) = sqrt(GMsun)* (1.0D0 + e_Ma_ini)* (-cos(Omega_a_Ma_ini)*sin(omega_p_Ma_ini) &
&- sin(Omega_a_Ma_ini)* cos(omega_p_Ma_ini)* cos(i_Ma_ini) ) /sqrt(A_Ma_ini* (1.0D0 - e_Ma_ini* e_Ma_ini) )
y0(8) = sqrt(GMsun)* (1.0D0 + e_Ma_ini)* (-sin(Omega_a_Ma_ini)*sin(omega_p_Ma_ini) &
&+ cos(Omega_a_Ma_ini)* cos(omega_p_Ma_ini)* cos(i_Ma_ini) ) /sqrt(A_Ma_ini* (1.0D0 - e_Ma_ini* e_Ma_ini) )
y0(9) = sqrt(GMsun)* (1.0D0 + e_Ma_ini)* ( cos(omega_p_Ma_ini)* sin(i_Ma_ini) )/sqrt(A_Ma_ini* (1.0D0 - e_Ma_ini* e_Ma_ini) )
!satellite 衛星は、太陽の周りを円運動をしている地球の周りを円運動させています。初期位置があまり地球に近いと、計算コストがかかるのと、今回用いるルンゲクッタは解いてくれないので、多少離してあります。
!Initiallly on
!x
y0(16) = y0(4) + 0.0025D0
y0(17) = y0(5)
y0(18) = y0(6)
!v
y0(13) = 0
y0(14) = sqrt(GMsun/1.0D0) + sqrt((GMsun*M_E)/0.0025D0)
y0(15) = 0
衛星の軌道長半径、離心率を求めています。
a_stllt = 1.0D0/( (2.0D0/ sqrt( y0(16)*y0(16) + y0(17)*y0(17) &
&+ y0(18)*y0(18) ) ) - (( y0(13)*y0(13) + y0(14)*y0(14) + y0(15)*y0(15) )/GMsun ) )
e_stllt = sqrt(1.0D0 - ( y0(16)*y0(14)-y0(17)*y0(13) )**2.0D0/(GMsun*a_stllt) )
open(10,file='orbit.dat') データファイルの設定
write(10,'(10f20.10)') x0, y0(4), y0(5), y0(6), y0(10), y0(11), y0(12)&
&, y0(16), y0(17), y0(18) ファイルに初期値を書き込んでいます。順に時間、地球のx, y, z、火星のx, y, z、衛星のx, y, zです。
!Calculation
do n = 1,nmax nmax回計算を繰り返します。
人工衛星の加速条件の判断用に、真近点離角のcosをとったものを求めます。
if(e_stllt==0) then
cosf_stllt = 0
else
cosf_stllt = ( a_stllt* (1.0D0 - e_stllt* e_stllt)/sqrt(y0(16)*y0(16) + y0(17)*y0(17) &
&+ y0(18)*y0(18) ) - 1.0D0)/e_stllt
end if
!Kick 人工衛星の加速についてです。
if(counter < n_acc ) then ある回数だけ加速します。
if(cosf_stllt<1.000001D0 .and. cosf_stllt>0.999999D0) then 人工衛星が太陽に対してある位置関係にあるときのみ加速します。
y0(14) = y0(14) + 0.010D0*sqrt(GMsun/a_stllt) !Kick 加速。
counter = counter+1
end if
end if
xn = x0 + dx
call runge(neq,ex2,x0,xn,y0,yn,eps,ierr,work) !runge.f90 ルンゲクッタを呼び出します。
x0 = xn !time evolution ルンゲクッタの計算後、時間を1つ進めます。
do i = 1,neq
y0(i) = yn(i) !update y0 ルンゲクッタの計算後、天体の座標と速度を更新します。
end do
次の加速判断用に軌道長半径、離心率を求めておきます。
a_stllt = 1.0D0/( (2.0D0/ sqrt( y0(16)*y0(16) + y0(17)*y0(17) &
&+ y0(18)*y0(18) ) ) - (( y0(13)*y0(13) + y0(14)*y0(14) + y0(15)*y0(15) )/GMsun ) )
e_stllt = sqrt(1.0D0 - ( y0(16)*y0(14)-y0(17)*y0(13) )**2.0D0/(GMsun*a_stllt) )
if( mod(n, 10)==0 ) then 今見たい物理は衛星の太陽中心での軌道進化なので、1ヶ月を100回刻んだ分の結果は必要ではありません。そのため10回に1回(つまり0.1ヶ月~0.01年に1回)だけ書きだします。
write(10,'(10f20.10)') x0, y0(4), y0(5), y0(6), y0(10), y0(11), y0(12)&
&, y0(16), y0(17), y0(18) ファイルへの書きだしです。
end if
end do
close(10)
end program main
real(kind=8) function ex2(i,x,y)
use constants_set Moduleを使います。
implicit none
関数ex2で使う変数を定義します。
integer(kind=4),parameter:: neq=18
integer(kind=4) i
real(kind=8) x
real(kind=8),dimension(neq):: y
real(kind=8) indirect_x, indirect_y, indirect_z
real(kind=8) r_SE, r_S2, r_Ss
real(kind=8) r_E2, r_2s, r_Es
real(kind=8) M_2
M_2 = M_Ma !
天体の太陽からの距離です。
r_SE = sqrt( y(4)*y(4) + y(5)*y(5) + y(6)*y(6) )
r_S2 = sqrt( y(10)*y(10) + y(11)*y(11) + y(12)*y(12) )
r_Ss = sqrt( y(16)*y(16) + y(17)*y(17) + y(18)*y(18) )
天体同士の距離です。
r_E2 = sqrt( (y(4)-y(10))*(y(4)-y(10)) + (y(5)-y(11))*(y(5)-y(11)) + (y(6)-y(12))*(y(6)-y(12)) )
r_2s = sqrt( (y(10)-y(16))*(y(10)-y(16)) + (y(11)-y(17))*(y(11)-y(17)) + (y(12)-y(18))*(y(12)-y(18)) )
r_Es = sqrt( (y(4)-y(16))*(y(4)-y(16)) + (y(5)-y(17))*(y(5)-y(17)) + (y(6)-y(18))*(y(6)-y(18)) )
原点を中心星に固定しているため、惑星による中心星の摂動が惑星の軌道に及ぼす効果を計算する必要があります。この効果は中心星の揺れなので、全ての天体に同じように加わるので、先に計算しておきます。
indirect_x = M_E*y(4)/(r_SE**(3.0D0)) + M_2*y(10)/(r_S2**(3.0D0))
indirect_y = M_E*y(5)/(r_SE**(3.0D0)) + M_2*y(11)/(r_S2**(3.0D0))
indirect_z = M_E*y(6)/(r_SE**(3.0D0)) + M_2*y(12)/(r_S2**(3.0D0))
if (i==1) then
!Earth 地球についての方程式です。
!solar gravity + gravity of 2nd body + indirect tetrm 地球に加わる加速度は (太陽の重力)＋(火星の重力)＋(太陽の揺れ) です。
ex2= -GMsun* ( y(4)/(r_SE**(3.0D0)) + M_2*( y(4)-y(10) )/(r_E2**(3.0D0)) + indirect_x ) !dvx/dt =
else if (i==2) then
ex2= -GMsun* ( y(5)/(r_SE**(3.0D0)) + M_2*( y(5)-y(11) )/(r_E2**(3.0D0)) + indirect_y ) !dvy/dt =
else if (i==3) then
ex2= -GMsun* ( y(6)/(r_SE**(3.0D0)) + M_2*( y(6)-y(12) )/(r_E2**(3.0D0)) + indirect_z ) !dvz/dt =
else if (i==4) then
ex2= y(1) !dx/dt
else if (i==5) then
ex2= y(2) !dy/dt =
else if (i==6) then
ex2= y(3) !dz/dt =
!2 火星に加わる加速度は、地球と同様に求まります。
else if (i==7) then
ex2= -GMsun* ( y(10)/(r_S2**(3.0D0)) + M_E*( y(10)-y(4) )/(r_E2**(3.0D0)) + indirect_x )
else if (i==8) then
ex2= -GMsun* ( y(11)/(r_S2**(3.0D0)) + M_E*( y(11)-y(5) )/(r_E2**(3.0D0)) + indirect_y )
else if (i==9) then
ex2= -GMsun* ( y(12)/(r_S2**(3.0D0)) + M_E*( y(12)-y(6) )/(r_E2**(3.0D0)) + indirect_z )
else if (i==10) then
ex2= y(7)
else if (i==11) then
ex2= y(8)
else if (i==12) then
ex2= y(9)
!satellite 人工衛星には地球と火星の加速度が加わる点に注意。
else if (i==13) then
ex2= -GMsun* ( y(16)/(r_Ss**(3.0D0)) + M_E*( y(16)-y(4) )/(r_Es**3.0D0) + M_2*( y(16)-y(10) )/(r_2s**(3.0D0)) + indirect_x )
else if (i==14) then
ex2= -GMsun* ( y(17)/(r_Ss**(3.0D0)) + M_E*( y(17)-y(5) )/(r_Es**3.0D0) + M_2*( y(17)-y(11) )/(r_2s**(3.0D0)) + indirect_y )
else if (i==15) then
ex2= -GMsun* ( y(18)/(r_Ss**(3.0D0)) + M_E*( y(18)-y(6) )/(r_Es**3.0D0) + M_2*( y(18)-y(12) )/(r_2s**(3.0D0)) + indirect_z )
else if (i==16) then
ex2=y(13)
else if (i==17) then
ex2=y(14)
else
ex2=y(15)
end if
return
end function ex2
|