|
|
! 惑星の軌道と、近点角
!------------------
!楕円軌道上の惑星は、楕円焦点(＝両者の重心 ≒主星) と結ぶ線の描く扇形の面積の、
!その時間的変化が一定になるように回る。（ケプラーの法則 Kepler's laws ）
! その回転角を、（真近点角 True anomaly 、又は、真近点離角）
!一様でない角速度は扱いにくく、上の扇形の面積に比例した角度を、抽象的に考える。
!これは一定速なので、比例計算に便利。 （平均近点角 Mean anomaly 、又は、平均近点離角）
!平均近点角 → 真近点角への仲介の角度。（離心近点角 Eccentric anomaly 、又は、離心近点離角）
!
!このプログラムは、上の３つの近点角を、アニメーションにしたものです。
!（近点角、又は、近点離角）とは、主星に最も近い軌道位置(近点)からの、左回り角度。
!------------------
!離心率 ＝ √(長半径^2－短半径^2) /長半径
!
LET e_=0.6 !離心率
LET a=1 !長半径
LET b=SQR(a^2*(1-e_^2)) !短半径
!
LET wh=1.2
SET WINDOW -wh,wh,-wh,wh
DRAW grid(.2,.2)
!
DRAW circle WITH SCALE(a)
DRAW disk WITH SCALE(.03)*SHIFT(a*e_,0)
FOR i=0 TO 2*PI STEP PI/64
PLOT LINES: a*COS(i),b*SIN(i);
NEXT i
PLOT LINES
!
PLOT TEXT,AT .86*wh,-.03*wh :"近点"
SET TEXT COLOR "blue"
PLOT TEXT,AT .6*wh,.76*wh :"平均近点角 Ｍ"
SET TEXT COLOR 64 !dark green
PLOT TEXT,AT .6*wh,.68*wh :"離心近点角 Ｅ"
SET TEXT COLOR "red"
PLOT TEXT,AT .6*wh,.60*wh :" 真近点角 Ｔ"
!
SET LINE width 2
SET DRAW MODE NOTXOR !２度書きで消す。<---> OVERWRITE
!
FOR i=-180 TO 360+15.1
LET M=i*PI/180 !平均近点角 M
!----
CALL calculate
CALL plotter(MM,EE,xx,yy) !前回の描画を消す。初回目 xx=yy=0 は、無動作
CALL plotter(M,E,x,y) !今回の描画
WAIT DELAY 0.01
NEXT i
!
FOR i=30 TO 360+31 STEP 15
LET M=i*PI/180 !平均近点角 M
!----
WAIT DELAY 0.5
CALL calculate
CALL plotter(MM,EE,xx,yy) !前回の描画を消す
CALL plotter(M,E,x,y) !今回の描画
NEXT i
WAIT DELAY 1
CALL plotter(MM,EE,xx,yy) !前回の描画を消す
!
SET DRAW MODE OVERWRITE !描画を通常へ戻す
!
PRINT USING "軌道長半径 a=#.#### 離心率 e=#.####":a,e_
FOR i=0 TO 360-1 STEP 15
LET M=i*PI/180 !平均近点角 M
!----
CALL calculate
CALL plotter(M,E,x,y) !重ね 描画
PRINT USING "平均近点角 M=####.#°離心近点角 E=####.##°真近点角 T=####.##°":M*180/PI,E*180/PI,T*180/PI;
PRINT USING "軌道半径 r=#.#### 惑星座標(x,y)=##.#### ##.####":r,x,y
NEXT i
!-------------------
SUB plotter(M,E,x,y)
IF x=0 AND y=0 THEN EXIT SUB
!----
SET LINE COLOR "blue"
PLOT LINES : 0,0; a*COS(M),a*SIN(M) !平均近点角 M
SET LINE COLOR 64 !dark green
PLOT LINES :0,0; x,a*SIN(E) !離心近点角 E
SET LINE width 1
PLOT LINES :x,0; x,a*SIN(E) !離心近点角の補助線
SET LINE width 2
SET LINE COLOR "red"
PLOT LINES :a*e_,0; x,y !真近点角 T 軌道半径 r
DRAW disk WITH SCALE(.02)*SHIFT(x,y) !惑星 (x,y)
!----
LET MM=M
LET EE=E
LET xx=x
LET yy=y
END SUB
!-------------------
SUB calculate
LET E=EnM(M) !離心近点角 E !LET E=E_M(M)
LET T=T_(E) !真近点角 T
LET r=r_E(E) !軌道半径 r !LET r=r_T(T)
LET x=a*COS(E) !惑星座標 (x,y)
LET y=b*SIN(E) ! !LET y=r*SIN(T)
END SUB
!-------------------
! 離心近点角 E ← 平均近点角 M
!
! M= E - e_*SIN(E) ・・・ケプラーの方程式( 第２法則「面積速度は一定」から導出)
!
! 解けないので、 E= M + e_*SIN(E)
! 漸化式 En+1= M + e_*SIN(En) ・・として、En+1==En まで収束させる。
!
! E0=M, E1=M+e_*sin(M), E2=M+e_*sin(M+e_*sin(M)),,,
!-------------------
FUNCTION E_M( M )
LET En1=M !初期値= M
DO
LET En=En1
LET En1=M+e_*SIN(En)
LOOP UNTIL ABS(En1-En)< 1e-7 !収束まで 15～34回
LET E_M=En1
END FUNCTION
!-------------------
! 離心近点角 E ← 平均近点角 M （ニュートン近似法）※この方が速い。
!
! M= E - e_*SIN(E) ・・・ケプラーの方程式
!
! f (E)= E - e_*SIN(E) - M ・・として、f(E)=0 となる E を求める。
! f'(E)= 1 - e_*COS(E)
! En+1= En- f(En)/f'(En)
!--------------------
FUNCTION EnM( M )
LET En1=M !初期値= M
DO
LET En=En1
LET En1=En-(En-e_*SIN(En)-M)/(1-e_*COS(En))
LOOP UNTIL ABS(En1-En)< 1e-7 !収束まで 3～4回
LET EnM=En1
END FUNCTION
!-------------------
!その他の計算
!-------------------
!DEF T_(E)=ACOS( (COS(E)-e_)/(1-e_*COS(E)) ) !真近点角 T(範囲狭い0~π) ← 離心近点角 E
DEF T_(E)=2*ATN( SQR((1+e_)/(1-e_))*TAN(E/2) ) !真近点角 T(-π~π) ← 離心近点角 E
DEF r_E( E)=a*(1-e_*COS(E)) !軌道半径 r ← 離心近点角 E
DEF r_T( T )=a*(1-e_^2)/(1+e_*COS(T)) !軌道半径 r ← 真近点角 T
END
|
|