|
|
!--------------------------------------------------------------------------------
!日の出、日の入りと、均時差。
!そのままであったが、今回は、それに習い、独自にやってみる。
!
!プログラムは、2010年と、2011年の２ヵ年間について、
!全４０ポイントで、暦象年表(2010年) 理科年表(2011年)の同項 と比較する。
!--------------------------------------------------------------------------------
!近日点通過 2010-1- 3(日) JST 08:59 J=2455200.374306（平均近点角0°)
! 春分 2010-3-21(日) JST 02:32 J=2455277.105556（太陽黄経0°)
!遠日点通過 2010-7- 6(火) JST 20:14 J=2455384.843056（平均近点角180°）
! 秋分 2010-9-23(木) JST 12:09 J=2455463.50625 （太陽黄経180°)
!--------------------------------------------------------
OPTION ANGLE DEGREES
LET e_=0.0167 !惑星地球の離心率
LET ε=-23.44 !地心黄道 → 赤道傾斜角ε
LET μ=0.84 !(太陽の視半径0.267°)＋地平線軸 屈折角0.57°≒0.84°
!
LET Jm0= CJDnum(2010,1, 3)+JFrac( 8,59,0) !近日点 JST 2010.1.03_08:59 (平均近点角 M=0)
LET Js0= CJDnum(2010,3,21)+JFrac( 2,32,0) !春分点 JST 2010.3.21_02:32 (太陽の黄経λ=0)
LET Ms= (Js0-Jm0)/D_Y(2010)*360 !春分点の、平均近点角 Ms
LET Ts= T_( EnM2( Ms )) !春分点の、真近点角 Ts
!
! E135°00' N34°38' 明石市相生町（現在の天文町2丁目）
! E139°44' 28.8759" N35°39' 29.1572" 東京（日本経緯度原点）旧東京天文台
!
LET Ls=135 !当地の中央標準時の経度（明石市）
LET Lw=139 +44/60 +28.8759/3600 ! 当地の経度（東京）
LET φ= 35 +39/60 +29.1572/3600 ! 当地の緯度（東京）
PRINT USING "東京都の、緯度＝##.###°経度＝###.###°" :φ,Lw
!
!均時差＝ 視太陽時－平均太陽時 ・・明石市でのみ、JST 12:00－均時差 に南中する。
!
!東京都の 均時差 出 南中 入
DATA 2010, 1, 1, " -3m18.2s 6:51 11:44:23 16:38"
DATA 2010, 1, 3, " -4m14.2s 6:51 11:45:19 16:40" !≒近日点(2010-1- 3 JST 08:59)
DATA 2010, 1, 7, " -6m01.6s 6:51 11:47:07 16:43"
DATA 2010, 2, 7, "-14m04.6s 6:36 11:55:07 17:14"
DATA 2010, 3, 7, "-11m08.7s 6:04 11:52:09 17:41"
DATA 2010, 3,21, "- 7m21.0s 5:44 11:48:21 17:53" !≒春分点(2010-3-21 JST 02:32)
DATA 2010, 4, 7, " -2m17.1s 5:20 11:43:17 18:07"
DATA 2010, 4,16, " +0m04.0s 5:08 11:40:57 18:15" !均時差≒０
DATA 2010, 5, 7, " +3m25.7s 4:44 11:37:36 18:32"
DATA 2010, 6, 7, " +1m15.6s 4:25 11:39:48 18:55"
DATA 2010, 6,13, " +0m03.9s 4:25 11:41:00 18:58" !均時差≒０
DATA 2010, 7, 7, " -4m49.6s 4:31 11:45:53 19:00"
DATA 2010, 8, 7, " -5m50.5s 4:53 11:46:52 18:40"
DATA 2010, 9, 2, " +0m08.1s 5:13 11:40:52 18:08" !均時差≒０
DATA 2010, 9, 7, " +1m47.0s 5:17 11:39:13 18:01"
DATA 2010, 9,23, " +7m25.5s 5:29 11:33:34 17:37" !≒秋分点(2010-9-23 JST 12:09)
DATA 2010,10, 7, "+12m00.0s 5:40 11:29:00 17:17"
DATA 2010,11, 7, "+16m20.5s 6:08 11:24:42 16:41"
DATA 2010,12, 7, " +8m47.8s 6:37 11:32:17 16:28"
DATA 2010,12,25, " +0m14.0s 6:48 11:40:51 16:33" !均時差≒０
DATA 2011, 1, 1, " -3m10.9s 6:50 11:44:16 16:38"
DATA 2011, 1, 4, " -4m35.1s 6:51 11:45:40 16:41" !≒近日点(2011-1-4 JST 03:32)
DATA 2011, 1, 7, " -5m55.9s 6:51 11:47:01 16:43"
DATA 2011, 2, 7, "-14m04.5s 6:37 11:55:07 17:14"
DATA 2011, 3, 7, "-11m12.9s 6:04 11:52:13 17:41"
DATA 2011, 3,21, "- 7m25.3s 5:44 11:48:25 17:53" !≒春分点(2011-3-21 JST 08:21)
DATA 2011, 4, 7, " -2m22.5s 5:20 11:43:23 18:07"
DATA 2011, 4,16, " +0m00.3s 5:08 11:41:01 18:14" !均時差≒０
DATA 2011, 5, 7, " +3m23.0s 4:44 11:37:39 18:32"
DATA 2011, 6, 7, " +1m16.4s 4:25 11:39:47 18:55"
DATA 2011, 6,13, " +0m05.7s 4:25 11:40:58 18:57" !均時差≒０
DATA 2011, 7, 7, " -4m48.7s 4:31 11:45:52 19:00"
DATA 2011, 8, 7, " -5m52.9s 4:53 11:46:54 18:40"
DATA 2011, 9, 2, " +0m03.3s 5:13 11:40:57 18:08" !均時差≒０
DATA 2011, 9, 7, " +1m42.1s 5:17 11:39:18 18:01"
DATA 2011, 9,23, " +7m21.5s 5:29 11:33:38 17:38" !≒秋分点(2011-9-23 JST 18:05)
DATA 2011,10, 7, "+11m56.1s 5:40 11:29:04 17:18"
DATA 2011,11, 7, "+16m21.6s 6:08 11:24:41 16:41"
DATA 2011,12, 7, " +8m54.1s 6:37 11:32:11 16:28"
DATA 2011,12,26, " -0m08.7s 6:49 11:41:14 16:34" !均時差≒０
DATA 2012, 1, 1, " 比較データーなし。"
!----------------------------------
! 日の出、日の入りと、均時差の計算
!----------------------------------
SUB sun( YY,MM,DD )
!
! Jx ＝ 当地の視太陽時が、12:00(南中) になる時の当地の中央標準時
! M ＝ その時の、平均近点角
!
LET Jx=CJDnum(YY,MM,DD)+.5 -((Lw-Ls)/360) !CJD.年月日(12:00) -(当地～当地の中央標準時との時差)
!
LET M= MOD( (Jx-Jm0)/D_Y(2010)*360, 360) !(太陽の)M= 当日の平均近点角、Jm0= 近日点 M=0°のCJD
LET Tm= MOD( T_( EnM2( M )), 360) ! " " Tm= 当日の真近点角
LET λ= MOD( Tm-Ts, 360) ! " " λ= 黄経、Ts= 春分点の真近点角 (λ黄経=0)
LET α= α_(λ,0) ! " " α= 赤経 (λ黄経,β黄緯=0)
!
LET ⊿e= ( M -Tm ) /360 !⊿e 当日の差角「平均近点角－真近点角」
LET ⊿s= ( λ-α ) /360 !⊿s 当日の差角「 黄経( 真近点角)－赤経(その投影)」
LET ⊿= ⊿e+⊿s !⊿均時差(合計)
LET Jtransit= Jx -⊿ !南中の、平均太陽時= (CJD 視太陽時)-⊿均時差
!
LET δ= δ_(λ,0) !太陽のδ赤緯 ←(λ黄経,β黄緯=0)、φ当地の緯度
LET W0= ACOS( -TAN(φ)*TAN(δ) -SIN(μ)/(COS(φ)*COS(δ)) )/360 !南中～SunSetの時角( 日単位)
LET Jset= Jtransit +W0 !日の入
LET Jrise= Jtransit -W0 !日の出
END SUB
!-------------------------
! MAIN 表示
!-------------------------
SET WINDOW -17,360,-15,+20
SET COLOR MIX(15) .4, .4, .4
DRAW grid(360/12, 5)
SET TEXT background "opaque"
ASK TEXT WIDTH(" ") tw
!----
PLOT TEXT,AT 9*tw,-13.2:"縦軸：±分(time)均時差"
PLOT TEXT,AT 9*tw,-14.5:"横軸：0~360°平均近点角( 近日点～近日点 )365.25日"
!----
PLOT TEXT,AT 4.5*tw,18.2:"視太陽時－平均太陽時"
PLOT TEXT,AT 4.5*tw,17 :"太陽の"
PLOT TEXT,AT 19*tw,17 :"＋"
SET TEXT COLOR "red"
PLOT TEXT,AT 0,18.2:"均時差＝"
PLOT TEXT,AT 3*tw,17:"＝"
SET TEXT COLOR "blue"
PLOT TEXT,AT 8*tw,17:"{平均近点角－真近点角}"
SET TEXT COLOR 10!"green"
PLOT TEXT,AT 20*tw,17:"{黄経－赤経}"
!----
PRINT
DO
READ IF MISSING THEN EXIT DO: YY,MM,DD,w$
!----
PRINT USING "####/##/## ":YY,MM,DD;
PRINT TAB(18);"均時差"; TAB(27);"日の出 南中 日の入"
PRINT "暦象年表の東京→ ";w$
!----
CALL sun( YY,MM,DD ) !均時差、日の出、日の入 の計算
!----
PRINT TAB(9);"計算値→";TAB(18);ms$(⊿); TAB(27);hm$(Jrise); hms$(Jtransit); hm$(Jset);
PRINT USING " Jx=#######.#####":Jx;
PRINT USING " Jrise=#######.##### Jtransit=#######.##### Jset=#######.#####":Jrise,Jtransit,Jset;
PRINT USING " 時角ω=###.####°":W0*360;
!----
PRINT USING " 平均近点角 M=###.####°黄経λ=###.####°太陽の赤緯δ=###.####°太陽の赤経α=###.####°=": M,λ,δ,α;
PRINT hms$(α/360);
PRINT
!----
IF M< Mw THEN LET Mw=Mw-360
PRINT
!---- 均時差のグラフ
LET ⊿e=⊿e*24*60
LET ⊿s=⊿s*24*60
LET ⊿=⊿*24*60 !均時差±分(time) ← ±日(time)
PLOT POINTS: M,⊿ !M= 平均近点角 0~360°365.25日(近日点～近日点)
PLOT POINTS: M,⊿e
PLOT POINTS: M,⊿s
IF Mw<>0 THEN
SET LINE COLOR "blue"
PLOT LINES: Mw,⊿ew; M,⊿e !{平均近点角－真近点角}の成分
SET LINE COLOR "green"
PLOT LINES: Mw,⊿sw; M,⊿s !{黄経－赤経}の成分
SET LINE COLOR "red"
PLOT LINES: Mw,⊿w ; M,⊿ !合計の均時差±分(time)
END IF
LET Mw=M
LET ⊿w=⊿
LET ⊿ew=⊿e
LET ⊿sw=⊿s
LOOP
!------------------------------
! 分、秒 "±00m00.0s" の文字列
!------------------------------
FUNCTION ms$(J) !0~±1 日
LET mm=1440*FP( ROUND( 86400*ABS(J),1) /86400) !0.1 秒未満、四捨五入
LET ss=60*FP(mm)
IF J< 0 THEN LET s$="-" ELSE LET s$="+"
LET ms$=s$& USING$("%%",IP(mm))& "m"& USING$("%%.#",ss)& "s"
END FUNCTION
!--------------------------
! 時、分 "-00:00" の文字列
!--------------------------
FUNCTION hm$(J) !0~±1 日
LET hh=24*FP( ROUND( 1440*ABS(J)) /1440) !1 分未満、四捨五入
LET mm=60*FP(hh)
IF J< 0 THEN LET s$="-" ELSE LET s$=" "
LET hm$=s$& USING$("##",IP(hh))& ":"& USING$("%%",mm)
END FUNCTION
!---------------------------------
! 時、分、秒 "-00:00:00" の文字列
!---------------------------------
FUNCTION hms$(J) !0~±1 日
LET hh=24*FP( ROUND( 86400*ABS(J)) /86400) !1 秒未満、四捨五入
LET mm=60*FP(hh)
LET ss=60*FP(mm)
IF J< 0 THEN LET s$="-" ELSE LET s$=" "
LET hms$=s$& USING$("##",IP(hh))& ":"& USING$("%%",IP(mm))& ":"& USING$("%%",ss)
END FUNCTION
!-------------------------------------------------------------
! Chronological Julian Day(CJD) 1900.3.1~ 2100.2.28 の範囲
!
! 小数下の端数を、時間に対応させるため、( 0~ 1 → 00:00~24:00)
! JD=0 の基点( 12:00 January 1, 4713 BC, Monday)を、0.5日過去へ下げたもの。
! CJD=0 の基点( 00:00 January 1, 4713 BC, Monday)
!-------------------------------------------------------------
!西暦年月日 YY.MM.DD から､ユリウス日 CJD の計算
FUNCTION CJDnum( YY,MM,DD)
IF 2< MM THEN
LET WM=MM+1
LET WY=YY
ELSE
LET WM=MM+13
LET WY=YY-1
END IF
LET CJDnum=INT(365.25*WY)+INT(30.6001*WM)+DD+1720982
END FUNCTION
!-------------------------
! 時、分、秒 → 日( 0~1)
!-------------------------
DEF JFrac(hh,mm,ss)= hh/24+mm/1440+ss/86400
!-------------------
! 離心近点角 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 EnM2( M ) !入出力は、degree 単位
LET M2=M *PI/180
LET En1=M2 !初期値= M2
DO
LET En=En1
LET En1=En-(En-e_*SIN(En *180/PI)-M2)/(1-e_*COS(En *180/PI))
LOOP UNTIL ABS(En1-En)< 1e-10 !7 !収束まで 3～4回
LET EnM2=En1 *180/PI
END FUNCTION
!-------------------
!その他の計算
!-------------------
! DEF T_(E)=ACOS( (COS(E)-e_)/(1-e_*COS(E)) ) 真近点角 T( 0~π) ← 離心近点角 E
!-----------------------------------------
!cos(T) = (cos(E)-e)/(1-e*cos(E)) 入出力角を T(-π~π) に、拡張する。
!cos(T/2+T/2)= (cos(E)-e)/(1-e*cos(E))
!1-2*sin(T/2)^2= (cos(E)-e)/(1-e*cos(E)) → 2*sin(T/2)^2=1-(cos(E)-e)/(1-e*cos(E)) …(1)
!2*cos(T/2)^2-1= (cos(E)-e)/(1-e*cos(E)) → 2*cos(T/2)^2=1+(cos(E)-e)/(1-e*cos(E)) …(2)
!tan(T/2)^2= (1-e*cos(E)-cos(E)+e)/(1-e*cos(E)+cos(E)-e) …(1)/(2)
! = (1+e-(1+e)*cos(E)) /(1-e+(1-e)*cos(E))
! = { (1+e)*(1-cos(E)) } /{ (1-e)*(1+cos(E)) }
! = (1+e)/(1-e) *(1-cos(E))/(1+cos(E))
! = (1+e)/(1-e) *tan(E/2)^2
!tan(T/2) = √{(1+e)/(1-e)} *tan(E/2)
DEF T_(E)=2*ATN( SQR((1+e_)/(1-e_))*TAN(E/2+EPS(E)) ) !真近点角 T(-π~π) ← 離心近点角 E
!-----------------------------------------
! 赤道座標 黄道座標
!
! |緯度δ経度αの方向余弦ｘ| |1 0 0 | |緯度β経度λの方向余弦ｘ|
! |緯度δ経度αの方向余弦ｙ|＝|0 cos(ε) sin(ε)| |緯度β経度λの方向余弦ｙ|
! |緯度δ経度αの方向余弦ｚ| |0 -sin(ε) cos(ε)| |緯度β経度λの方向余弦ｚ|
!
! |cos(δ)*cos(α)| |1 0 0 | |cos(β)*cos(λ)|
! |cos(δ)*sin(α)|＝|0 cos(ε) sin(ε)| |cos(β)*sin(λ)|
! |sin(δ) | |0 -sin(ε) cos(ε)| |sin(β) |
!
!緯度δ経度αの方向余弦ｘ cos(δ)*cos(α)= cos(β)*cos(λ)
!緯度δ経度αの方向余弦ｙ cos(δ)*sin(α)= cos(β)*sin(λ)*cos(ε) +sin(β)*sin(ε)
!緯度δ経度αの方向余弦ｚ sin(δ)= -cos(β)*sin(λ)*sin(ε) +sin(β)*cos(ε) ※1
!
!経度α=mod( angle(ｘ,ｙ), 360) ※2 (output 0~360°)
!経度α=mod( angle(ｘ/cos(β),ｙ/cos(β)), 360) ※2'
!-----------------------------
DEF δ_(λ,β)= ASIN( -COS(β)*SIN(λ)*SIN(ε) +SIN(β)*COS(ε) ) !緯度δ ※1
DEF α_(λ,β)= MOD( ANGLE( COS(λ), SIN(λ)*COS(ε)+TAN(β)*SIN(ε) ), 360) !経度α ※2'
!太陽年（年毎の、春分点～春分点の日数）
DEF D_Y(yy)= 365.24219878-(6.14e-8)*(yy-1900)
!--------------------------------------------------------
END
|
|