longitudeSun

ユリウス世紀(Julian Century)から太陽黄経を求めます。

構文
longitudeSun( JC )
引数
JC
ユリウス世紀
戻値
太陽黄経(λsun)

プログラム

//////////////////////////////////////////////////
// 【引数】
//   JC : ユリウス世紀 
// 【戻値】
//   太陽黄経(λsun) 
//////////////////////////////////////////////////
FUNCTION longitudeSun(JC)
	DIM A[18] = 0.0004, 0.0004, 0.0005, 0.0005, 0.0006, 0.0007, 0.0007, 0.0007, 0.0013, 0.0015, 0.0018, 0.0018, 0.0020, 0.0200, -0.0048*JC, 1.9147, 36000.7695*JC, 280.4659
	DIM k[18] = 31557.0, 29930.0, 2281.0, 155.0, 33718.0, 9038.0, 3035.0, 65929.0, 22519.0, 45038.0, 445267.0, 19.0, 32964.0, 71998.1, 35999.05, 35999.05, 0, 0
	DIM θ0[18] = 161.0, 48.0, 221.0, 118.0, 316.0, 64.0, 110.0, 45.0, 352.0, 254.0, 208.0, 159.0, 158.0, 265.1, 267.52, 267.52, 0, 0
	DIM λsun[18]
	FOR n = 0 TO 15
		DIM ang = normalizeAngle(k[n] * JC + θ0[n])
		λsun[n] = A[n] * COS(degToRad(ang))
	NEXT
	λsun[16] = normalizeAngle(A[16])
	λsun[17] = normalizeAngle(A[17])
	RESULT = normalizeAngle(CALCARRAY(λsun, CALC_ADD))
FEND

//////////////////////////////////////////////////
// 【引数】
//   deg : 度数法 
// 【戻値】
//   度数法から弧度法に変換した値 
//////////////////////////////////////////////////
FUNCTION degToRad(deg)
	DIM pi = 3.14159265358979
	RESULT = deg * pi / 180
FEND

//////////////////////////////////////////////////
// 【引数】
//   deg : 度数法 
// 【戻値】
//   度単位の角度を0~360度に正規化 
//////////////////////////////////////////////////
FUNCTION normalizeAngle(deg)
	SELECT TRUE
		CASE deg >= 360
			deg = deg - INT(deg / 360) * 360
		CASE deg < 0
			deg = deg + INT(ABS(deg / 360) + 1) * 360
	SELEND
	RESULT = deg
FEND

太陽黄経

\[\lambda_{\rm{sun}}=\sum_{n=1}^{18}A*\cos(k*t+\theta_{0})\]
\(A\)
振幅
\(k\)
角速度
\(t\)
時刻
\(\theta_{0}\)
初期位相角
\(n\) \(A\) \(k\) \(\theta_{0}\)
1 0.0004 31557 161
2 0.0004 29930 48
3 0.0005 2281 221
4 0.0005 155 118
5 0.0006 33718 316
6 0.0007 9038 64
7 0.0007 3035 110
8 0.0007 65929 45
9 0.0013 22519 352
10 0.0015 45038 254
11 0.0018 445267 208
12 0.0018 19 159
13 0.002 32964 158
14 0.02 71998.1 265.1
15 -0.0048t 35999.05 267.52
16 1.9147 35999.05 267.52
17 36000.7695t 0 0
18 280.4659 0 0

プログラム実行例

次の満月の日付を求める

GETTIME()
DIM year = G_TIME_YY
DIM month = G_TIME_MM
DIM day = G_TIME_DD
DIM JD = YMDToJD(year, month, day)

REPEAT
	DIM TD = JD - 9/24
	DIM JC = (TD + 0.5 - 2451545) / 36525
	DIM λsun = longitudeSun(JC)
	DIM λmoon = longitudeMoon(JC)
	DIM Δλ = ABS(λsun - λmoon)
	JD = JD + ABS(180 - Δλ) / 15	// 180°に遠いときは大きい値、近いときは小さい値を足す
UNTIL ABS(Δλ - 180) <= 1 // 誤差が1以下になるまで

d = JDToYMD(JD)

RESIZE(d, 2)
PRINT JOIN(d, "/")

//////////////////////////////////////////////////
// 【引数】
//   deg : 度数法 
// 【戻値】
//   度数法から弧度法に変換した値 
//////////////////////////////////////////////////
FUNCTION degToRad(deg)
	DIM pi = 3.14159265358979
	RESULT = deg * pi / 180
FEND

//////////////////////////////////////////////////
// 【引数】
//   JD : ユリウス日 
// 【戻値】
//   グレゴリオ暦を格納した配列(0 : 年, 1 : 月, 2 : 日, 3 : 時, 4 : 分, 5 : 秒) 
//////////////////////////////////////////////////
FUNCTION JDToYMD(JD)
	DIM x0 = INT(JD + 68570)
	DIM x1 = INT(x0 / 36524.25)
	DIM x2 = x0 - INT(36524.25 * x1 + 0.75)
	DIM x3 = INT((x2 + 1) / 365.2425)
	DIM x4 = x2 - INT(365.25 * x3) + 31
	DIM x5 = INT(INT(x4) / 30.59)
	DIM x6 = INT(INT(x5) / 11)
	
	DIM t2 = x4 - INT(30.59 * x5)
	DIM t1 = x5 - 12 * x6 + 2
	DIM t0 = 100 * (x1 - 49) + x3 + x6
	
	IFB t1 = 2 AND t2 > 28 THEN
		SELECT TRUE
			CASE t0 MOD 100 = 0 AND t0 MOD 400 = 0
				t2 = 29
			CASE t0 MOD 4 = 0
				t2 = 29
			DEFAULT
				t2 = 28
		SELEND
	ENDIF
	
	DIM tm = 86400 * (JD - INT(JD))
	DIM t3 = INT(tm / 3600)
	DIM t4 = INT((tm - 3600 * t3) / 60)
	DIM t5 = INT(tm - 3600 * t3 - 60 * t4)
	DIM t[] = t0, t1, t2, t3, t4, t5
	
	RESULT = SLICE(t)
FEND

//////////////////////////////////////////////////
// 【引数】
//   JC : ユリウス世紀 
// 【戻値】
//   月黄経(λmoon) 
//////////////////////////////////////////////////
FUNCTION longitudeMoon(JC)
	DIM A[63] = 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0004, 0.0004, 0.0005, 0.0005, 0.0005, 0.0006, 0.0006, 0.0007, 0.0007, 0.0007, 0.0007, 0.0008, 0.0009, 0.0011, 0.0012, 0.0016, 0.0018, 0.0021, 0.0021, 0.0021, 0.0022, 0.0023, 0.0024, 0.0026, 0.0027, 0.0028, 0.0037, 0.0038, 0.004, 0.004, 0.004, 0.005, 0.0052, 0.0068, 0.0079, 0.0085, 0.01, 0.0107, 0.011, 0.0125, 0.0154, 0.0304, 0.0347, 0.0409, 0.0458, 0.0533, 0.0571, 0.0588, 0.1144, 0.1851, 0.2136, 0.6583, 1.274, 6.2888, 481267.8809 * JC, 218.3162
	DIM k[63] = 2322131, 4067, 549197, 1808933, 349472, 381404, 958465, 12006, 39871, 509131, 1745069, 1908795, 2258267, 111869, 27864, 485333, 405201, 790672, 1403732, 858602, 1920802, 1267871, 1856938, 401329, 341337, 71998, 990397, 818536, 922466, 99863, 1379739, 918399, 1934, 541062, 1781068, 133, 1844932, 1331734, 481266, 31932, 926533, 449334, 826671, 1431597, 1303870, 489205, 1443603, 75870, 513197.9, 445267.1, 441199.8, 854535.2, 1367733.1, 377336.3, 63863.5, 966404, 35999, 954397.7, 890534.2, 413335.3, 477198.9, 0, 0
	DIM θ0[63] = 191, 70, 220, 58, 337, 354, 340, 187, 223, 242, 24, 90, 156, 38, 127, 186, 50, 114, 98, 129, 186, 249, 152, 274, 16, 85, 357, 151, 163, 122, 17, 182, 145, 259, 21, 29, 56, 283, 205, 107, 323, 188, 111, 315, 246, 142, 52, 41, 222.5, 27.9, 47.4, 148.2, 280.7, 13.2, 124.2, 276.5, 87.53, 179.93, 145.7, 10.74, 44.963, 0, 0
	DIM λmoon[63]
	FOR n = 0 TO 60
		DIM ang = normalizeAngle(k[n] * JC + θ0[n])
		λmoon[n] = A[n] * COS(degToRad(ang))
	NEXT
	λmoon[61] = normalizeAngle(A[61])
	λmoon[62] = normalizeAngle(A[62])
	RESULT = normalizeAngle(CALCARRAY(λmoon, CALC_ADD))
FEND

//////////////////////////////////////////////////
// 【引数】
//   JC : ユリウス世紀 
// 【戻値】
//   太陽黄経(λsun) 
//////////////////////////////////////////////////
FUNCTION longitudeSun(JC)
	DIM A[18] = 0.0004, 0.0004, 0.0005, 0.0005, 0.0006, 0.0007, 0.0007, 0.0007, 0.0013, 0.0015, 0.0018, 0.0018, 0.0020, 0.0200, -0.0048*JC, 1.9147, 36000.7695*JC, 280.4659
	DIM k[18] = 31557.0, 29930.0, 2281.0, 155.0, 33718.0, 9038.0, 3035.0, 65929.0, 22519.0, 45038.0, 445267.0, 19.0, 32964.0, 71998.1, 35999.05, 35999.05, 0, 0
	DIM θ0[18] = 161.0, 48.0, 221.0, 118.0, 316.0, 64.0, 110.0, 45.0, 352.0, 254.0, 208.0, 159.0, 158.0, 265.1, 267.52, 267.52, 0, 0
	DIM λsun[18]
	FOR n = 0 TO 15
		DIM ang = normalizeAngle(k[n] * JC + θ0[n])
		λsun[n] = A[n] * COS(degToRad(ang))
	NEXT
	λsun[16] = normalizeAngle(A[16])
	λsun[17] = normalizeAngle(A[17])
	RESULT = normalizeAngle(CALCARRAY(λsun, CALC_ADD))
FEND

//////////////////////////////////////////////////
// 【引数】
//   deg : 度数法 
// 【戻値】
//   度単位の角度を0~360度に正規化 
//////////////////////////////////////////////////
FUNCTION normalizeAngle(deg)
	SELECT TRUE
		CASE deg >= 360
			deg = deg - INT(deg / 360) * 360
		CASE deg < 0
			deg = deg + INT(ABS(deg / 360) + 1) * 360
	SELEND
	RESULT = deg
FEND

//////////////////////////////////////////////////
// 【引数】
//   year : 年 
//   month : 月 
//   day : 日 
//   hour : 時 
//   minute : 分 
//   second : 秒 
// 【戻値】
//   ユリウス日 
//////////////////////////////////////////////////
FUNCTION YMDToJD(year, month, day, hour = 0, minute = 0, second = 0)
	IFB month < 3 THEN
		year = year - 1
		month = month + 12
	ENDIF
	DIM JD = INT(year * 365.25)
	JD = JD + INT(year / 400)
	JD = JD - INT(year / 100)
	JD = JD + INT((month - 2) * 30.59)
	JD = JD + 1721088
	JD = JD + day
	DIM t = second / 3600
	t = t + minute / 60
	t = t + hour
	t = t / 24
	JD = JD + t
	RESULT = JD
FEND
  1. GETTIME
  2. YMDToJD
  3. longitudeSun
  4. longitudeMoon
  5. ABS
  6. JDToYMD
  7. RESIZE
  8. JOIN
結果
2020/5/7