-
Notifications
You must be signed in to change notification settings - Fork 2
/
qreki.rb
666 lines (591 loc) · 28.7 KB
/
qreki.rb
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
require "qreki/version"
require "date"
# =========================================================================
# Based on 「旧暦計算サンプルプログラム」
# Copyright (C) 1993,1994 by H.Takano
# http://www.vector.co.jp/soft/dos/personal/se016093.html
# =========================================================================
module Qreki
class Srk
attr_accessor :year, :month, :day
end
class Qrk
attr_accessor :year, :uruu, :month, :day, :rokuyou, :sekki
end
$k = Math::PI / 180.0 # Pi
$tz = - (9.0 / 24.0) # タイムゾーンオフセット
$rm_sun0 = 0 # 太陽黄経
#=========================================================================
# 旧暦算出
#
# 引数 .... 日付
# 戻り値
# q.year : 旧暦年
# q.uruu : 平月／閏月 flag .... 平月:false 閏月:true
# q.month : 旧暦月
# q.day : 旧暦日
# q.rokuyou : 六曜
# q.sekki : 二十四節気
#=========================================================================
def self.calc(year, month, day)
calc_from_date(Date.new(year, month, day))
end
def self.calc_from_date(tm)
array = calc_kyureki(tm.year, tm.month, tm.day)
qreki = Qrk.new
qreki.year = array[0]
qreki.uruu = array[1]
qreki.month = array[2]
qreki.day = array[3]
qreki.rokuyou = rokuyou(tm.year, tm.month, tm.day)
qreki.sekki = sekki(tm.year, tm.month, tm.day)
qreki
end
# =========================================================================
# 新暦に対応する、旧暦を求める。
#
# 呼び出し時にセットする変数：
# year : 計算する年（JSTユリウス日）
# mon : 計算する月（JSTユリウス日）
# day : 計算する日（JSTユリウス日）
# 戻り値：
# this.year : 旧暦年
# this.uruu : 平月／閏月 flag .... 平月:false 閏月:true
# this.month : 旧暦月
# this.day : 旧暦日
# =========================================================================
def self.calc_kyureki(year, mon, day)
chu = Array.new(5, [])
m = Array.new(5).map{Array.new(3,0)}
saku = Array.new
kyureki = Array.new
tm = ymdt2jd(year, mon, day, 0, 0, 0)
# -----------------------------------------------------------------------
# 計算対象の直前にあたる二分二至の時刻を求める
# -----------------------------------------------------------------------
chu[0] = calc_chu(tm, 90)
# -----------------------------------------------------------------------
# 上で求めた二分二至の時の太陽黄経をもとに朔日行列の先頭に月名をセット
# -----------------------------------------------------------------------
m[0][0] = ( $rm_sun0 / 30.0 ).to_i + 2
# -----------------------------------------------------------------------
# 中気の時刻を計算（３回計算する）
# chu[i]:中気の時刻
# -----------------------------------------------------------------------
(1..3).each do |i|
chu[i] = calc_chu(chu[i - 1] + 32, 30)
end
# -----------------------------------------------------------------------
# 計算対象の直前にあたる二分二至の直前の朔の時刻を求める
# -----------------------------------------------------------------------
saku[0] = calc_saku(chu[0])
# -----------------------------------------------------------------------
# 朔の時刻を求める
# -----------------------------------------------------------------------
(1..4).each do |i|
saku[i] = calc_saku(saku[i - 1] + 30)
# 前と同じ時刻を計算した場合（両者の差が26日以内）には、初期値を
# +33日にして再実行させる。
if (saku[i - 1].to_i.abs - saku[i].to_i) <= 26
saku[i] = calc_saku(saku[i - 1] + 35)
end
end
# -----------------------------------------------------------------------
# saku[1]が二分二至の時刻以前になってしまった場合には、朔をさかのぼり過ぎ
# たと考えて、朔の時刻を繰り下げて修正する。
# その際、計算もれ（saku[4]）になっている部分を補うため、朔の時刻を計算
# する。（近日点通過の近辺で朔があると起こる事があるようだ...？）
# -----------------------------------------------------------------------
if saku[1].to_i <= chu[0].to_i
(0..4).each do |i|
saku[i] = saku[i + 1]
end
saku[4] = calc_saku(saku[3] + 35)
# -----------------------------------------------------------------------
# saku[0]が二分二至の時刻以後になってしまった場合には、朔をさかのぼり足
# りないと見て、朔の時刻を繰り上げて修正する。
# その際、計算もれ（saku[0]）になっている部分を補うため、朔の時刻を計算
# する。（春分点の近辺で朔があると起こる事があるようだ...？）
# -----------------------------------------------------------------------
elsif saku[0].to_i > chu[0].to_i
[4,3,2,1].each do |i|
saku[i] = saku[i - 1]
end
saku[0] = calc_saku(saku[0] - 27)
end
# -----------------------------------------------------------------------
# 閏月検索Ｆｌａｇセット
# （節月で４ヶ月の間に朔が５回あると、閏月がある可能性がある。）
# lap=false:平月 lap=true:閏月
# -----------------------------------------------------------------------
lap = (saku[4].to_i <= chu[3].to_i)
# -----------------------------------------------------------------------
# 朔日行列の作成
# m[i][0] ... 月名（1:正月 2:２月 3:３月 ....）
# m[i][1] .... 閏フラグ（false:平月 true:閏月）
# m[i][2] ...... 朔日のjd
# -----------------------------------------------------------------------
m[0][1] = false
m[0][2] = saku[0].to_i
(1..4).each do |i|
if lap && i > 1
if (chu[i - 1] <= saku[i - 1].to_i) || (chu[i - 1] >= saku[i].to_i)
m[i-1][0] = m[i-2][0]
m[i-1][1] = true
m[i-1][2] = saku[i - 1].to_i
lap = false
end
end
m[i][0] = m[i-1][0] + 1
m[i][0] -= 12 if m[i][0] > 12
m[i][2] = saku[i].to_i
m[i][1] = false
end
# -----------------------------------------------------------------------
# 朔日行列から旧暦を求める。
# -----------------------------------------------------------------------
state = 0
i_tmp = 0
(1..4).each do |i|
i_tmp = i
if tm.to_i < m[i][2].to_i
state = 1
break
elsif tm.to_i == m[i][2].to_i
state = 2
break
end
end
i_tmp -= 1 if state == 0 || state == 1
i = i_tmp
kyureki[1] = m[i][1]
kyureki[2] = m[i][0]
kyureki[3] = tm.to_i - m[i][2].to_i + 1
# -----------------------------------------------------------------------
# 旧暦年の計算
# （旧暦月が10以上でかつ新暦月より大きい場合には、
# まだ年を越していないはず...）
# -----------------------------------------------------------------------
a = jd2ymdt(tm)
kyureki[0] = a[0]
if kyureki[2] > 9 && kyureki[2] > a[1]
kyureki[0] -= 1
end
[kyureki[0], kyureki[1], kyureki[2], kyureki[3]]
end
# =========================================================================
# 直前の二分二至／中気の時刻を求める
#
# パラメータ
# tm ............ 計算基準となる時刻（JSTユリウス日）
# longitude ..... 求める対象（90:二分二至,30:中気））
# 戻り値
# 求めた時刻（JSTユリウス日）を返す。
# グローバル変数 $rm_sun0 に、その時の太陽黄経をセットする。
#
# ※ 引数、戻り値ともユリウス日で表し、時分秒は日の小数で表す。
# 力学時とユリウス日との補正時刻=0.0secと仮定
# =========================================================================
def self.calc_chu(tm, longitude)
# -----------------------------------------------------------------------
# 時刻引数を小数部と整数部とに分解する（精度を上げるため）
# -----------------------------------------------------------------------
tm1 = tm.to_i
tm2 = tm - tm1 + $tz # JST -> UTC
# -----------------------------------------------------------------------
# 直前の二分二至の黄経 λsun0 を求める
# -----------------------------------------------------------------------
t = (tm2 + 0.5) / 36525.0 + (tm1 - 2451545.0) / 36525.0
rm_sun = longitude_sun(t)
$rm_sun0 = longitude * (rm_sun / longitude.to_f).to_i
# -----------------------------------------------------------------------
# 繰り返し計算によって直前の二分二至の時刻を計算する
# （誤差が±1.0 sec以内になったら打ち切る。）
# -----------------------------------------------------------------------
delta_t1 ||= 0.0
delta_t2 ||= 1.0
while (delta_t1 + delta_t2).abs > (1.0 / 86400.0)
# -------------------------------------------------------------------
# λsun(t) を計算
# t = (tm + .5 - 2451545) / 36525;
# -------------------------------------------------------------------
t = (tm2 + 0.5) / 36525.0 + (tm1 - 2451545.0) / 36525.0
rm_sun = longitude_sun(t)
# -------------------------------------------------------------------
# 黄経差 Δλ＝λsun －λsun0
# -------------------------------------------------------------------
delta_rm = rm_sun - $rm_sun0
# -------------------------------------------------------------------
# Δλの引き込み範囲（±180°）を逸脱した場合には、補正を行う
# -------------------------------------------------------------------
if delta_rm > 180.0
delta_rm -= 360.0
elsif delta_rm < -180.0
delta_rm += 360.0
end
# -------------------------------------------------------------------
# 時刻引数の補正値 Δt
# delta_t = delta_rm * 365.2 / 360;
# -------------------------------------------------------------------
delta_t1 = (delta_rm * 365.2 / 360.0).to_i
delta_t2 = (delta_rm * 365.2 / 360.0) - delta_t1
# -------------------------------------------------------------------
# 時刻引数の補正
# tm -= delta_t;
# -------------------------------------------------------------------
tm1 = tm1 - delta_t1
tm2 = tm2 - delta_t2
if tm2 < 0
tm1 -= 1.0
tm2 += 1.0
end
end
# -----------------------------------------------------------------------
# 戻り値の作成
# 時刻引数を合成し、戻り値（JSTユリウス日）とする
# -----------------------------------------------------------------------
tm2 + tm1 - $tz
end
# =========================================================================
# 直前の朔の時刻を求める
#
# 呼び出し時にセットする変数
# tm ........ 計算基準の時刻（JSTユリウス日）
# 戻り値
# 朔の時刻
#
# ※ 引数、戻り値ともJSTユリウス日で表し、時分秒は日の小数で表す。
# 力学時とユリウス日との補正時刻=0.0secと仮定
# =========================================================================
def self.calc_saku(tm)
# -----------------------------------------------------------------------
# ループカウンタのセット
# -----------------------------------------------------------------------
lc = 1
# -----------------------------------------------------------------------
# 時刻引数を小数部と整数部とに分解する（精度を上げるため）
# -----------------------------------------------------------------------
tm1 = tm.to_i
tm2 = tm - tm1 + $tz # JST -> UTC
# -----------------------------------------------------------------------
# 繰り返し計算によって朔の時刻を計算する
# （誤差が±1.0 sec以内になったら打ち切る。）
# -----------------------------------------------------------------------
delta_t1 ||= 0.0
delta_t2 ||= 1.0
while ( delta_t1 + delta_t2 ).abs > ( 1.0 / 86400.0 )
lc += 1
# -------------------------------------------------------------------
# 太陽の黄経λsun(t) ,月の黄経λmoon(t) を計算
# t = (tm + .5 - 2451545) / 36525;
# -------------------------------------------------------------------
t = (tm2 + 0.5) / 36525.0 + (tm1 - 2451545.0) / 36525.0
rm_sun = longitude_sun(t)
rm_moon = longitude_moon(t)
# -------------------------------------------------------------------
# 月と太陽の黄経差Δλ
# Δλ＝λmoon－λsun
# -------------------------------------------------------------------
delta_rm = rm_moon - rm_sun
# -------------------------------------------------------------------
# ループの１回目（lc=1）で delta_rm < 0 の場合には引き込み範囲に
# 入るように補正する
# -------------------------------------------------------------------
if lc==1 && delta_rm < 0
delta_rm = normalization_angle(delta_rm)
# -------------------------------------------------------------------
# 春分の近くで朔がある場合（0 ≦λsun≦ 20）で、
# 月の黄経λmoon≧300 の場合には、
# Δλ＝ 360 － Δλ と計算して補正する
# -------------------------------------------------------------------
elsif rm_sun >= 0 && rm_sun <= 20 && rm_moon >= 300
delta_rm = normalization_angle(delta_rm)
delta_rm = 360 - delta_rm
# -------------------------------------------------------------------
# Δλの引き込み範囲（±40°）を逸脱した場合には、補正を行う
# -------------------------------------------------------------------
elsif delta_rm.abs > 40
delta_rm = normalization_angle(delta_rm)
end
# -------------------------------------------------------------------
# 時刻引数の補正値 Δt
# delta_t = delta_rm * 29.530589 / 360;
# -------------------------------------------------------------------
delta_t1 = ( delta_rm * 29.530589 / 360.0 ).to_i
delta_t2 = ( delta_rm * 29.530589 / 360.0 ) - delta_t1
# -------------------------------------------------------------------
# 時刻引数の補正
# tm -= delta_t;
# -------------------------------------------------------------------
tm1 = tm1 - delta_t1
tm2 = tm2 - delta_t2
if( tm2 < 0 )
tm1 -= 1.0
tm2 += 1.0
end
# -------------------------------------------------------------------
# ループ回数が15回になったら、初期値 tm を tm-26 とする。
# -------------------------------------------------------------------
if( lc == 15 && (delta_t1 + delta_t2).abs > (1 / 86400.0) )
tm1 = (tm - 26).to_i
tm2 = 0
# -------------------------------------------------------------------
# 初期値を補正したにも関わらず、振動を続ける場合には初期値を答えとし
# て返して強制的にループを抜け出して異常終了させる。
# -------------------------------------------------------------------
elsif( lc > 30 && (delta_t1 + delta_t2).abs > (1 / 86400.0) )
tm1 = tm
tm2 = 0
break
end
end
# -----------------------------------------------------------------------
# 戻り値の作成
# 時刻引数を合成し、戻り値（ユリウス日）とする
# -----------------------------------------------------------------------
tm2 + tm1 - $tz
end
# =========================================================================
# 角度の正規化を行う。すなわち引数の範囲を ０≦θ＜３６０ にする。
# =========================================================================
def self.normalization_angle(angle)
if angle >= 0.0
return angle - 360.0 * ( angle / 360.0 ).to_i
else
return 360.0 + angle - 360.0 * ( angle / 360.0 ).to_i
end
end
# =========================================================================
# 太陽の黄経 λsun(t) を計算する（t は力学時）
# =========================================================================
def self.longitude_sun(t)
# -----------------------------------------------------------------------
# 摂動項の計算
# -----------------------------------------------------------------------
th = 0.0004 * Math.cos( $k * normalization_angle( 31557.0 * t + 161.0 ) )
th += 0.0004 * Math.cos( $k * normalization_angle( 29930.0 * t + 48.0 ) )
th += 0.0005 * Math.cos( $k * normalization_angle( 2281.0 * t + 221.0 ) )
th += 0.0005 * Math.cos( $k * normalization_angle( 155.0 * t + 118.0 ) )
th += 0.0006 * Math.cos( $k * normalization_angle( 33718.0 * t + 316.0 ) )
th += 0.0007 * Math.cos( $k * normalization_angle( 9038.0 * t + 64.0 ) )
th += 0.0007 * Math.cos( $k * normalization_angle( 3035.0 * t + 110.0 ) )
th += 0.0007 * Math.cos( $k * normalization_angle( 65929.0 * t + 45.0 ) )
th += 0.0013 * Math.cos( $k * normalization_angle( 22519.0 * t + 352.0 ) )
th += 0.0015 * Math.cos( $k * normalization_angle( 45038.0 * t + 254.0 ) )
th += 0.0018 * Math.cos( $k * normalization_angle( 445267.0 * t + 208.0 ) )
th += 0.0018 * Math.cos( $k * normalization_angle( 19.0 * t + 159.0 ) )
th += 0.0020 * Math.cos( $k * normalization_angle( 32964.0 * t + 158.0 ) )
th += 0.0200 * Math.cos( $k * normalization_angle( 71998.1 * t + 265.1 ) )
th -= 0.0048 * Math.cos( $k * normalization_angle( 35999.05 * t + 267.52 ) ) * t
th += 1.9147 * Math.cos( $k * normalization_angle( 35999.05 * t + 267.52 ) )
# -----------------------------------------------------------------------
# 比例項の計算
# -----------------------------------------------------------------------
ang = normalization_angle( 36000.7695 * t )
ang = normalization_angle( ang + 280.4659 )
th = normalization_angle( th + ang )
th
end
# =========================================================================
# 月の黄経 λmoon(t) を計算する（t は力学時）
# =========================================================================
def self.longitude_moon(t)
# -----------------------------------------------------------------------
# 摂動項の計算
# -----------------------------------------------------------------------
th = 0.0003 * Math.cos( $k * normalization_angle( 2322131.0 * t + 191.0 ) )
th += 0.0003 * Math.cos( $k * normalization_angle( 4067.0 * t + 70.0 ) )
th += 0.0003 * Math.cos( $k * normalization_angle( 549197.0 * t + 220.0 ) )
th += 0.0003 * Math.cos( $k * normalization_angle( 1808933.0 * t + 58.0 ) )
th += 0.0003 * Math.cos( $k * normalization_angle( 349472.0 * t + 337.0 ) )
th += 0.0003 * Math.cos( $k * normalization_angle( 381404.0 * t + 354.0 ) )
th += 0.0003 * Math.cos( $k * normalization_angle( 958465.0 * t + 340.0 ) )
th += 0.0004 * Math.cos( $k * normalization_angle( 12006.0 * t + 187.0 ) )
th += 0.0004 * Math.cos( $k * normalization_angle( 39871.0 * t + 223.0 ) )
th += 0.0005 * Math.cos( $k * normalization_angle( 509131.0 * t + 242.0 ) )
th += 0.0005 * Math.cos( $k * normalization_angle( 1745069.0 * t + 24.0 ) )
th += 0.0005 * Math.cos( $k * normalization_angle( 1908795.0 * t + 90.0 ) )
th += 0.0006 * Math.cos( $k * normalization_angle( 2258267.0 * t + 156.0 ) )
th += 0.0006 * Math.cos( $k * normalization_angle( 111869.0 * t + 38.0 ) )
th += 0.0007 * Math.cos( $k * normalization_angle( 27864.0 * t + 127.0 ) )
th += 0.0007 * Math.cos( $k * normalization_angle( 485333.0 * t + 186.0 ) )
th += 0.0007 * Math.cos( $k * normalization_angle( 405201.0 * t + 50.0 ) )
th += 0.0007 * Math.cos( $k * normalization_angle( 790672.0 * t + 114.0 ) )
th += 0.0008 * Math.cos( $k * normalization_angle( 1403732.0 * t + 98.0 ) )
th += 0.0009 * Math.cos( $k * normalization_angle( 858602.0 * t + 129.0 ) )
th += 0.0011 * Math.cos( $k * normalization_angle( 1920802.0 * t + 186.0 ) )
th += 0.0012 * Math.cos( $k * normalization_angle( 1267871.0 * t + 249.0 ) )
th += 0.0016 * Math.cos( $k * normalization_angle( 1856938.0 * t + 152.0 ) )
th += 0.0018 * Math.cos( $k * normalization_angle( 401329.0 * t + 274.0 ) )
th += 0.0021 * Math.cos( $k * normalization_angle( 341337.0 * t + 16.0 ) )
th += 0.0021 * Math.cos( $k * normalization_angle( 71998.0 * t + 85.0 ) )
th += 0.0021 * Math.cos( $k * normalization_angle( 990397.0 * t + 357.0 ) )
th += 0.0022 * Math.cos( $k * normalization_angle( 818536.0 * t + 151.0 ) )
th += 0.0023 * Math.cos( $k * normalization_angle( 922466.0 * t + 163.0 ) )
th += 0.0024 * Math.cos( $k * normalization_angle( 99863.0 * t + 122.0 ) )
th += 0.0026 * Math.cos( $k * normalization_angle( 1379739.0 * t + 17.0 ) )
th += 0.0027 * Math.cos( $k * normalization_angle( 918399.0 * t + 182.0 ) )
th += 0.0028 * Math.cos( $k * normalization_angle( 1934.0 * t + 145.0 ) )
th += 0.0037 * Math.cos( $k * normalization_angle( 541062.0 * t + 259.0 ) )
th += 0.0038 * Math.cos( $k * normalization_angle( 1781068.0 * t + 21.0 ) )
th += 0.0040 * Math.cos( $k * normalization_angle( 133.0 * t + 29.0 ) )
th += 0.0040 * Math.cos( $k * normalization_angle( 1844932.0 * t + 56.0 ) )
th += 0.0040 * Math.cos( $k * normalization_angle( 1331734.0 * t + 283.0 ) )
th += 0.0050 * Math.cos( $k * normalization_angle( 481266.0 * t + 205.0 ) )
th += 0.0052 * Math.cos( $k * normalization_angle( 31932.0 * t + 107.0 ) )
th += 0.0068 * Math.cos( $k * normalization_angle( 926533.0 * t + 323.0 ) )
th += 0.0079 * Math.cos( $k * normalization_angle( 449334.0 * t + 188.0 ) )
th += 0.0085 * Math.cos( $k * normalization_angle( 826671.0 * t + 111.0 ) )
th += 0.0100 * Math.cos( $k * normalization_angle( 1431597.0 * t + 315.0 ) )
th += 0.0107 * Math.cos( $k * normalization_angle( 1303870.0 * t + 246.0 ) )
th += 0.0110 * Math.cos( $k * normalization_angle( 489205.0 * t + 142.0 ) )
th += 0.0125 * Math.cos( $k * normalization_angle( 1443603.0 * t + 52.0 ) )
th += 0.0154 * Math.cos( $k * normalization_angle( 75870.0 * t + 41.0 ) )
th += 0.0304 * Math.cos( $k * normalization_angle( 513197.9 * t + 222.5 ) )
th += 0.0347 * Math.cos( $k * normalization_angle( 445267.1 * t + 27.9 ) )
th += 0.0409 * Math.cos( $k * normalization_angle( 441199.8 * t + 47.4 ) )
th += 0.0458 * Math.cos( $k * normalization_angle( 854535.2 * t + 148.2 ) )
th += 0.0533 * Math.cos( $k * normalization_angle( 1367733.1 * t + 280.7 ) )
th += 0.0571 * Math.cos( $k * normalization_angle( 377336.3 * t + 13.2 ) )
th += 0.0588 * Math.cos( $k * normalization_angle( 63863.5 * t + 124.2 ) )
th += 0.1144 * Math.cos( $k * normalization_angle( 966404.0 * t + 276.5 ) )
th += 0.1851 * Math.cos( $k * normalization_angle( 35999.05 * t + 87.53 ) )
th += 0.2136 * Math.cos( $k * normalization_angle( 954397.74 * t + 179.93 ) )
th += 0.6583 * Math.cos( $k * normalization_angle( 890534.22 * t + 145.7 ) )
th += 1.2740 * Math.cos( $k * normalization_angle( 413335.35 * t + 10.74 ) )
th += 6.2888 * Math.cos( $k * normalization_angle( 477198.868 * t + 44.963) )
#-----------------------------------------------------------------------
# 比例項の計算
#-----------------------------------------------------------------------
ang = normalization_angle( 481267.8809 * t )
ang = normalization_angle( ang + 218.3162 )
th = normalization_angle( th + ang )
th
end
#=========================================================================
# 年月日、時分秒（世界時）からユリウス日（JD）を計算する
#
# ※ この関数では、グレゴリオ暦法による年月日から求めるものである。
# （ユリウス暦法による年月日から求める場合には使用できない。）
#=========================================================================
def self.ymdt2jd(year, mon, day, hour, min, sec)
if mon < 3
year -= 1
mon += 12
end
jd = ( 365.25 * year ).to_i
jd += ( year / 400.0 ).to_i
jd -= ( year / 100.0 ).to_i
jd += ( 30.59 * ( mon - 2 ) ).to_i
jd += 1721088
jd += day
t = sec / 3600.0
t += min / 60.0
t += hour
t = t / 24.0
jd += t
jd
end
#=========================================================================
# ユリウス日（JD）から年月日、時分秒（世界時）を計算する
#
# 戻り値の配列TIME[]の内訳
# TIME[0] ... 年 TIME[1] ... 月 TIME[2] ... 日
# TIME[3] ... 時 TIME[4] ... 分 TIME[5] ... 秒
#
# ※ この関数で求めた年月日は、グレゴリオ暦法によって表されている。
#
#=========================================================================
def self.jd2ymdt(jd)
time = []
x0 = ( jd + 68570.0 ).to_i
x1 = ( x0 / 36524.25 ).to_i
x2 = x0 - ( 36524.25 * x1 + 0.75 ).to_i
x3 = ( ( x2+1 ) / 365.2425 ).to_i
x4 = x2 - ( 365.25 * x3 ).to_i + 31.0
x5 = ( x4.to_i / 30.59 ).to_i
x6 = ( x5.to_i / 11.0 ).to_i
time[2] = x4 - ( 30.59 * x5 ).to_i
time[1] = x5 - 12 * x6 + 2
time[0] = 100 * ( x1 - 49 ) + x3 + x6
if time[1] == 2 && time[2] > 28
if time[0] % 100 == 0 && time[0] % 400 == 0
time[2] = 29
elsif time[0] % 4 ==0
time[2]=29
else
time[2]=28
end
end
tm = 86400.0 * ( jd - jd.to_i )
time[3] = ( tm / 3600.0 ).to_i
time[4] = ( (tm - 3600.0 * time[3] ) / 60.0 ).to_i
time[5] = ( tm - 3600.0 * time[3] - 60 * time[4] ).to_i
time
end
#=========================================================================
# 六曜算出関数
#
# 引数 .... 計算対象となる年月日 $year $mon $day
# 戻り値 .... 六曜 (大安 赤口 先勝 友引 先負 仏滅)
#=========================================================================
def self.rokuyou(year, mon, day)
rokuyou = %w(大安 赤口 先勝 友引 先負 仏滅)
q_yaer, uruu, q_mon, q_day = calc_kyureki(year, mon, day)
rokuyou[ (q_mon + q_day) % 6 ]
end
#=========================================================================
# 今日が２４節気かどうか調べる
#
# 引数 .... 計算対象となる年月日 $year $mon $day
# 戻り値 .... ２４節気の名称
#=========================================================================
def self.sekki(year, mon, day)
#-----------------------------------------------------------------------
# ２４節気の定義
#-----------------------------------------------------------------------
sekki24 = %w(春分 清明 穀雨 立夏 小満 芒種 夏至 小暑 大暑 立秋 処暑 白露
秋分 寒露 霜降 立冬 小雪 大雪 冬至 小寒 大寒 立春 雨水 啓蟄)
tm = ymdt2jd(year, mon, day, 0, 0, 0)
#-----------------------------------------------------------------------
# 時刻引数を分解する
#-----------------------------------------------------------------------
tm1 = tm.to_i
tm2 = tm - tm1
tm2 -= 9.0 / 24.0
t = (tm2 + 0.5) / 36525.0 + (tm1 - 2451545.0) / 36525.0
#今日の太陽の黄経
rm_sun_today = longitude_sun(t)
tm += 1
tm1 = tm.to_i
tm2 = tm - tm1
tm2 -= 9.0 / 24.0
t = (tm2 + 0.5) / 36525.0 + (tm1 - 2451545.0) / 36525.0
# 明日の太陽の黄経
rm_sun_tommorow = longitude_sun(t)
#
rm_sun_today0 = 15.0 * (rm_sun_today / 15.0).to_i
rm_sun_tommorow0 = 15.0 * (rm_sun_tommorow / 15.0).to_i
if rm_sun_today0 != rm_sun_tommorow0
return sekki24[rm_sun_tommorow0 / 15]
else
return ''
end
end
def self.shunbun(year)
day = (20.8431 + 0.242194 * ( year - 1980 ) - ( year - 1980 ) / 4).to_i
sreki = Srk.new
sreki.year = year
sreki.month = 3
sreki.day = day
sreki
end
def self.shuubun(year)
day = (23.2488 + 0.242194 * ( year - 1980 ) - ( year - 1980 ) / 4).to_i
sreki = Srk.new
sreki.year = year
sreki.month = 9
sreki.day = day
sreki
end
end