Module: MkSunmoon::Compute

Included in:
Sunmoon
Defined in:
lib/mk_sunmoon/compute.rb

Class Method Summary collapse

Class Method Details

.compute_angle_ecliptic(jy, t, lambda, beta) ⇒ Object

時刻(t)における黄経、黄緯(λ(jy),β(jy))の天体の方位角(ang)計算

@param: lambda (天体の黄経(λ(T)(度))) @param: beta (天体の黄緯(β(T)(度))) @param: jy (ユリウス年) @param: t (時刻(0.xxxx日)) @return: angle (角度(xx.x度))



513
514
515
516
 
# File 'lib/mk_sunmoon/compute.rb', line 513

def compute_angle_ecliptic(jy, t, lambda, beta)
  res = eclip2equat(jy, lambda, beta)
  return compute_angle_equator(jy, t, *res)
end

.compute_angle_equator(jy, t, alpha, delta) ⇒ Object

時刻(t)における赤経、赤緯(α(jy),δ(jy))(度)の天体の方位角(ang)計算

@param: jy (ユリウス年) @param: t (時刻 (0.xxxx日)) @param: alpha (天体の赤経(α(jy)(度))) @param: delta (天体の赤緯(δ(jy)(度))) @return: angle (角度(xx.x度))



527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
 
# File 'lib/mk_sunmoon/compute.rb', line 527

def compute_angle_equator(jy, t, alpha, delta)
  time_sidereal = compute_sidereal_time(jy, t)
  hour_angle = time_sidereal - alpha
  a_0  = -1.0 * Math.cos(Const::PI_180 * delta) \
              * Math.sin(Const::PI_180 * hour_angle)
  a_1  = Math.sin(Const::PI_180 * delta) \
       * Math.cos(Const::PI_180 * @lat)  \
       - Math.cos(Const::PI_180 * delta) \
       * Math.sin(Const::PI_180 * @lat)  \
       * Math.cos(Const::PI_180 * hour_angle)
  angle  = Math.atan(a_0 / a_1) / Const::PI_180
  angle += 360.0 if a_1 > 0.0 && angle < 0.0
  angle += 180.0 if a_1 < 0.0
  return angle
end

.compute_beta_moon(jy) ⇒ Object

月視黄緯の計算

@param: jy (ユリウス年(JST)) @return: lambda



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
 
# File 'lib/mk_sunmoon/compute.rb', line 313

def compute_beta_moon(jy)
  bm  =  0.0005 * Math.sin(Const::PI_180 * norm_angle(307.0   +    19.4    * jy))
  bm +=  0.0026 * Math.sin(Const::PI_180 * norm_angle( 55.0   +    19.34   * jy))
  bm +=  0.0040 * Math.sin(Const::PI_180 * norm_angle(119.5   +     1.33   * jy))
  bm +=  0.0043 * Math.sin(Const::PI_180 * norm_angle(322.1   +    19.36   * jy))
  bm +=  0.0267 * Math.sin(Const::PI_180 * norm_angle(234.95  +    19.341  * jy))
  bt  =  0.0003 * Math.sin(Const::PI_180 * norm_angle(234.0   + 19268.0    * jy))
  bt +=  0.0003 * Math.sin(Const::PI_180 * norm_angle(146.0   +  3353.3    * jy))
  bt +=  0.0003 * Math.sin(Const::PI_180 * norm_angle(107.0   + 18149.4    * jy))
  bt +=  0.0003 * Math.sin(Const::PI_180 * norm_angle(205.0   + 22642.7    * jy))
  bt +=  0.0004 * Math.sin(Const::PI_180 * norm_angle(147.0   + 14097.4    * jy))
  bt +=  0.0004 * Math.sin(Const::PI_180 * norm_angle( 13.0   +  9325.4    * jy))
  bt +=  0.0004 * Math.sin(Const::PI_180 * norm_angle( 81.0   + 10242.6    * jy))
  bt +=  0.0004 * Math.sin(Const::PI_180 * norm_angle(238.0   + 23281.3    * jy))
  bt +=  0.0004 * Math.sin(Const::PI_180 * norm_angle(311.0   +  9483.9    * jy))
  bt +=  0.0005 * Math.sin(Const::PI_180 * norm_angle(239.0   +  4193.4    * jy))
  bt +=  0.0005 * Math.sin(Const::PI_180 * norm_angle(280.0   +  8485.3    * jy))
  bt +=  0.0006 * Math.sin(Const::PI_180 * norm_angle( 52.0   + 13617.3    * jy))
  bt +=  0.0006 * Math.sin(Const::PI_180 * norm_angle(224.0   +  5590.7    * jy))
  bt +=  0.0007 * Math.sin(Const::PI_180 * norm_angle(294.0   + 13098.7    * jy))
  bt +=  0.0008 * Math.sin(Const::PI_180 * norm_angle(326.0   +  9724.1    * jy))
  bt +=  0.0008 * Math.sin(Const::PI_180 * norm_angle( 70.0   + 17870.7    * jy))
  bt +=  0.0010 * Math.sin(Const::PI_180 * norm_angle( 18.0   + 12978.66   * jy))
  bt +=  0.0011 * Math.sin(Const::PI_180 * norm_angle(138.3   + 19147.99   * jy))
  bt +=  0.0012 * Math.sin(Const::PI_180 * norm_angle(148.2   +  4851.36   * jy))
  bt +=  0.0012 * Math.sin(Const::PI_180 * norm_angle( 38.4   +  4812.68   * jy))
  bt +=  0.0013 * Math.sin(Const::PI_180 * norm_angle(155.4   +   379.35   * jy))
  bt +=  0.0013 * Math.sin(Const::PI_180 * norm_angle( 95.8   +  4472.03   * jy))
  bt +=  0.0014 * Math.sin(Const::PI_180 * norm_angle(219.2   +   299.96   * jy))
  bt +=  0.0015 * Math.sin(Const::PI_180 * norm_angle( 45.8   +  9964.00   * jy))
  bt +=  0.0015 * Math.sin(Const::PI_180 * norm_angle(211.1   +  9284.69   * jy))
  bt +=  0.0016 * Math.sin(Const::PI_180 * norm_angle(135.7   +   420.02   * jy))
  bt +=  0.0017 * Math.sin(Const::PI_180 * norm_angle( 99.8   + 14496.06   * jy))
  bt +=  0.0018 * Math.sin(Const::PI_180 * norm_angle(270.8   +  5192.01   * jy))
  bt +=  0.0018 * Math.sin(Const::PI_180 * norm_angle(243.3   +  8206.68   * jy))
  bt +=  0.0019 * Math.sin(Const::PI_180 * norm_angle(230.7   +  9244.02   * jy))
  bt +=  0.0021 * Math.sin(Const::PI_180 * norm_angle(170.1   +  1058.66   * jy))
  bt +=  0.0022 * Math.sin(Const::PI_180 * norm_angle(331.4   + 13377.37   * jy))
  bt +=  0.0025 * Math.sin(Const::PI_180 * norm_angle(196.5   +  8605.38   * jy))
  bt +=  0.0034 * Math.sin(Const::PI_180 * norm_angle(319.9   +  4433.31   * jy))
  bt +=  0.0042 * Math.sin(Const::PI_180 * norm_angle(103.9   + 18509.35   * jy))
  bt +=  0.0043 * Math.sin(Const::PI_180 * norm_angle(307.6   +  5470.66   * jy))
  bt +=  0.0082 * Math.sin(Const::PI_180 * norm_angle(144.9   +  3713.33   * jy))
  bt +=  0.0088 * Math.sin(Const::PI_180 * norm_angle(176.7   +  4711.96   * jy))
  bt +=  0.0093 * Math.sin(Const::PI_180 * norm_angle(277.4   +  8845.31   * jy))
  bt +=  0.0172 * Math.sin(Const::PI_180 * norm_angle(  3.18  + 14375.997  * jy))
  bt +=  0.0326 * Math.sin(Const::PI_180 * norm_angle(328.96  + 13737.362  * jy))
  bt +=  0.0463 * Math.sin(Const::PI_180 * norm_angle(172.55  +   698.667  * jy))
  bt +=  0.0554 * Math.sin(Const::PI_180 * norm_angle(194.01  +  8965.374  * jy))
  bt +=  0.1732 * Math.sin(Const::PI_180 * norm_angle(142.427 +  4073.3220 * jy))
  bt +=  0.2777 * Math.sin(Const::PI_180 * norm_angle(138.311 +    60.0316 * jy))
  bt +=  0.2806 * Math.sin(Const::PI_180 * norm_angle(228.235 +  9604.0088 * jy))
  bt +=  5.1282 * Math.sin(Const::PI_180 * norm_angle( 93.273 +  4832.0202 * jy + bm))
  return bt
end

.compute_diff_moon(jy) ⇒ Object

月の視差の計算

@param: jy (ユリウス年) @return: t (出入時刻(0.xxxx日))



375
376
377
378
379
380
381
382
383
384
385
386
 
# File 'lib/mk_sunmoon/compute.rb', line 375

def compute_diff_moon(jy)
  t  =  0.0003 * Math.sin(Const::PI_180 * norm_angle(227.0  +  4412.0   * jy))
  t +=  0.0004 * Math.sin(Const::PI_180 * norm_angle(194.0  +  3773.4   * jy))
  t +=  0.0005 * Math.sin(Const::PI_180 * norm_angle(329.0  +  8545.4   * jy))
  t +=  0.0009 * Math.sin(Const::PI_180 * norm_angle(100.0  + 13677.3   * jy))
  t +=  0.0028 * Math.sin(Const::PI_180 * norm_angle(  0.0  +  9543.98  * jy))
  t +=  0.0078 * Math.sin(Const::PI_180 * norm_angle(325.7  +  8905.34  * jy))
  t +=  0.0095 * Math.sin(Const::PI_180 * norm_angle(190.7  +  4133.35  * jy))
  t +=  0.0518 * Math.sin(Const::PI_180 * norm_angle(224.98 +  4771.989 * jy))
  t +=  0.9507 * Math.sin(Const::PI_180 * norm_angle(90.0))
  return t
end

.compute_dist_sun(jy) ⇒ Object

太陽の距離の計算

@param: jy (ユリウス年(JST)) @return: distance



218
219
220
221
222
223
224
225
226
227
 
# File 'lib/mk_sunmoon/compute.rb', line 218

def compute_dist_sun(jy)
  r_sun  = 0.000007 * Math.sin(Const::PI_180 * norm_angle(156.0 +  329.6  * jy))
  r_sun += 0.000007 * Math.sin(Const::PI_180 * norm_angle(254.0 +  450.4  * jy))
  r_sun += 0.000013 * Math.sin(Const::PI_180 * norm_angle( 27.8 + 4452.67 * jy))
  r_sun += 0.000030 * Math.sin(Const::PI_180 * norm_angle( 90.0))
  r_sun += 0.000091 * Math.sin(Const::PI_180 * norm_angle(265.1 +  719.98 * jy))
  r_sun += (0.007256 - 0.0000002 * jy) * Math.sin(Const::PI_180 * norm_angle(267.54 + 359.991 * jy))
  r_sun  = 10.0 ** r_sun
  return r_sun
end

.compute_dt(year, month, day) ⇒ Object

ΔT の計算

@param: year @param: month @param: day @return: dt



61
62
63
64
65
 
# File 'lib/mk_sunmoon/compute.rb', line 61

def compute_dt(year, month, day)
  ymd = sprintf("%04d%02d%02d", year, month, day)
  tm = MkTime.new(ymd)
  return tm.dt
end

.compute_height_ecliptic(jy, t, lambda, beta) ⇒ Object

時刻(t)における黄経、黄緯(λ(jy),β(jy))の天体の高度(height)計算

@param: jy (ユリウス年) @param: t (時刻(0.xxxx日)) @param: lambda (天体の黄経(λ(T)(度))) @param: beta (天体の黄緯(β(T)(度))) @return: height (高度(xx.x度))



552
553
554
555
 
# File 'lib/mk_sunmoon/compute.rb', line 552

def compute_height_ecliptic(jy, t, lambda, beta)
  res = eclip2equat(jy, lambda, beta)
  return compute_height_equator(jy, t, *res)
end

.compute_height_equator(jy, t, alpha, delta) ⇒ Object

時刻(t)における赤経、赤緯(α(jy),δ(jy))(度)の天体の高度(height)計算

@param: jy (ユリウス年) @param: t (時刻 (0.xxxx日)) @param: alpha (天体の赤経α(jy)(度)) @param: delta (天体の赤緯δ(jy)(度)) @return: height (高度(xx.x度))



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
 
# File 'lib/mk_sunmoon/compute.rb', line 566

def compute_height_equator(jy, t, alpha, delta)
  time_sidereal = compute_sidereal_time(jy, t)
  sidereal = time_sidereal - alpha
  height  = Math.sin(Const::PI_180 * delta) \
          * Math.sin(Const::PI_180 * @lat)  \
          + Math.cos(Const::PI_180 * delta) \
          * Math.cos(Const::PI_180 * @lat)  \
          * Math.cos(Const::PI_180 * sidereal)
  height  = Math.asin(height) / Const::PI_180

  # 大気差補正
  # [ 以下の内、3-2の計算式を採用 ]
  # # 1. 日月出没計算 by「菊池さん」による計算式
  # #   [ http://kikuchisan.net/ ]
  # h = 0.0167 / Math.tan( Const::PI_180 * ( height + 8.6 / ( height + 4.4 ) ) )

  # # 2. 中川用語集による計算式 ( 5度 - 85度用 )
  # #   [ http://www.es.ris.ac.jp/~nakagawa/term_collection/yogoshu/ll/ni.htm ]
  # h  = 58.1      / Math.tan( height )
  # h -=  0.07     / Math.tan( height ) ** 3
  # h +=  0.000086 / Math.tan( height ) ** 5
  # h *= 1 / 3600.0

  # # 3-1. フランスの天文学者ラドー(R.Radau)の平均大気差と１秒程度の差で大気差を求めることが可能
  # # ( 標準的大気(気温10ﾟC，気圧1013.25hPa)の場合 )
  # # ( 視高度30ﾟ以上 )
  # h  = ( 58.294  / 3600.0 ) * Math.tan( Const::PI_180 * ( 90.0 - height ) )
  # h -= (  0.0668 / 3600.0 ) * Math.tan( Const::PI_180 * ( 90.0 - height ) ) ** 3

  # 3-2. フランスの天文学者ラドー(R.Radau)の平均大気差と１秒程度の差で大気差を求めることが可能
  # ( 標準的大気(気温10ﾟC，気圧1013.25hPa)の場合 )
  # ( 視高度 4ﾟ以上 )
  a = Math.tan(Const::PI_180 * (90.0 - height))
  h  = (58.76   + \
       (-0.406  + \
       (- 0.0192) \
       * a) * a) * a * 1 / 3600.0

  # # 3-3. さらに、上記の大気差(3-1,3-2)を気温、気圧を考慮する
  # # ( しかし、気温・気圧を考慮してもさほど変わりはない )
  # pres = 1013.25 # <= 変更
  # temp = 30.0    # <= 変更
  # h *= pres / 1013.25
  # h *= 283.25 / ( 273.15 + temp )

  return height + h
end

.compute_hour_angle_diff(alpha, delta, time_sidereal, height, div) ⇒ Object

出入点(k)の時角(tk)と天体の時角(t)との差(dt=tk-t)を計算する

@param: alpha(R.A.) (赤経) @param: delta(Decl.) (赤緯) @param: time_sidereal(恒星時Θ(度)) @param: height (出没高度(度)) @param: div (0: 出, 1: 入, 2: 南中) @return: hour angle difference (時角の差　dt)



458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
 
# File 'lib/mk_sunmoon/compute.rb', line 458

def compute_hour_angle_diff(alpha, delta, time_sidereal, height, div)
  if div == 2
    tk = 0
  else
    tk  = Math.sin(Const::PI_180 * height)
    tk -= Math.sin(Const::PI_180 * delta) * Math.sin(Const::PI_180 * @lat)
    tk /= Math.cos(Const::PI_180 * delta) * Math.cos(Const::PI_180 * @lat)
    # 出没点の時角
    return 0.0 if tk.abs > 1
    tk = Math.acos(tk) / Const::PI_180
    # tkは出のときマイナス、入のときプラス
    tk = -tk if div == 0 && tk > 0
    tk = -tk if div == 1 && tk < 0
  end
  # 天体の時角
  t = time_sidereal - alpha
  dt = tk - t
  # dtの絶対値を180°以下に調整
  while dt >  180.0; dt -= 360.0; end
  while dt < -180.0; dt += 360.0; end
  return dt
end

.compute_lambda_moon(jy) ⇒ Object

月視黄経の計算

@param: jy (ユリウス年(JST)) @return: lambda



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
 
# File 'lib/mk_sunmoon/compute.rb', line 235

def compute_lambda_moon(jy)
  am  = 0.0006 * Math.sin(Const::PI_180 * norm_angle( 54.0   +    19.3    * jy))
  am += 0.0006 * Math.sin(Const::PI_180 * norm_angle( 71.0   +     0.2    * jy))
  am += 0.0020 * Math.sin(Const::PI_180 * norm_angle( 55.0   +    19.34   * jy))
  am += 0.0040 * Math.sin(Const::PI_180 * norm_angle(119.5   +     1.33   * jy))
  rm  = 0.0003 * Math.sin(Const::PI_180 * norm_angle(280.0   + 23221.3    * jy))
  rm += 0.0003 * Math.sin(Const::PI_180 * norm_angle(161.0   +    40.7    * jy))
  rm += 0.0003 * Math.sin(Const::PI_180 * norm_angle(311.0   +  5492.0    * jy))
  rm += 0.0003 * Math.sin(Const::PI_180 * norm_angle(147.0   + 18089.3    * jy))
  rm += 0.0003 * Math.sin(Const::PI_180 * norm_angle( 66.0   +  3494.7    * jy))
  rm += 0.0003 * Math.sin(Const::PI_180 * norm_angle( 83.0   +  3814.0    * jy))
  rm += 0.0004 * Math.sin(Const::PI_180 * norm_angle( 20.0   +   720.0    * jy))
  rm += 0.0004 * Math.sin(Const::PI_180 * norm_angle( 71.0   +  9584.7    * jy))
  rm += 0.0004 * Math.sin(Const::PI_180 * norm_angle(278.0   +   120.1    * jy))
  rm += 0.0004 * Math.sin(Const::PI_180 * norm_angle(313.0   +   398.7    * jy))
  rm += 0.0005 * Math.sin(Const::PI_180 * norm_angle(332.0   +  5091.3    * jy))
  rm += 0.0005 * Math.sin(Const::PI_180 * norm_angle(114.0   + 17450.7    * jy))
  rm += 0.0005 * Math.sin(Const::PI_180 * norm_angle(181.0   + 19088.0    * jy))
  rm += 0.0005 * Math.sin(Const::PI_180 * norm_angle(247.0   + 22582.7    * jy))
  rm += 0.0006 * Math.sin(Const::PI_180 * norm_angle(128.0   +  1118.7    * jy))
  rm += 0.0007 * Math.sin(Const::PI_180 * norm_angle(216.0   +   278.6    * jy))
  rm += 0.0007 * Math.sin(Const::PI_180 * norm_angle(275.0   +  4853.3    * jy))
  rm += 0.0007 * Math.sin(Const::PI_180 * norm_angle(140.0   +  4052.0    * jy))
  rm += 0.0008 * Math.sin(Const::PI_180 * norm_angle(204.0   +  7906.7    * jy))
  rm += 0.0008 * Math.sin(Const::PI_180 * norm_angle(188.0   + 14037.3    * jy))
  rm += 0.0009 * Math.sin(Const::PI_180 * norm_angle(218.0   +  8586.0    * jy))
  rm += 0.0011 * Math.sin(Const::PI_180 * norm_angle(276.5   + 19208.02   * jy))
  rm += 0.0012 * Math.sin(Const::PI_180 * norm_angle(339.0   + 12678.71   * jy))
  rm += 0.0016 * Math.sin(Const::PI_180 * norm_angle(242.2   + 18569.38   * jy))
  rm += 0.0018 * Math.sin(Const::PI_180 * norm_angle(  4.1   +  4013.29   * jy))
  rm += 0.0020 * Math.sin(Const::PI_180 * norm_angle( 55.0   +    19.34   * jy))
  rm += 0.0021 * Math.sin(Const::PI_180 * norm_angle(105.6   +  3413.37   * jy))
  rm += 0.0021 * Math.sin(Const::PI_180 * norm_angle(175.1   +   719.98   * jy))
  rm += 0.0021 * Math.sin(Const::PI_180 * norm_angle( 87.5   +  9903.97   * jy))
  rm += 0.0022 * Math.sin(Const::PI_180 * norm_angle(240.6   +  8185.36   * jy))
  rm += 0.0024 * Math.sin(Const::PI_180 * norm_angle(252.8   +  9224.66   * jy))
  rm += 0.0024 * Math.sin(Const::PI_180 * norm_angle(211.9   +   988.63   * jy))
  rm += 0.0026 * Math.sin(Const::PI_180 * norm_angle(107.2   + 13797.39   * jy))
  rm += 0.0027 * Math.sin(Const::PI_180 * norm_angle(272.5   +  9183.99   * jy))
  rm += 0.0037 * Math.sin(Const::PI_180 * norm_angle(349.1   +  5410.62   * jy))
  rm += 0.0039 * Math.sin(Const::PI_180 * norm_angle(111.3   + 17810.68   * jy))
  rm += 0.0040 * Math.sin(Const::PI_180 * norm_angle(119.5   +     1.33   * jy))
  rm += 0.0040 * Math.sin(Const::PI_180 * norm_angle(145.6   + 18449.32   * jy))
  rm += 0.0040 * Math.sin(Const::PI_180 * norm_angle( 13.2   + 13317.34   * jy))
  rm += 0.0048 * Math.sin(Const::PI_180 * norm_angle(235.0   +    19.34   * jy))
  rm += 0.0050 * Math.sin(Const::PI_180 * norm_angle(295.4   +  4812.66   * jy))
  rm += 0.0052 * Math.sin(Const::PI_180 * norm_angle(197.2   +   319.32   * jy))
  rm += 0.0068 * Math.sin(Const::PI_180 * norm_angle( 53.2   +  9265.33   * jy))
  rm += 0.0079 * Math.sin(Const::PI_180 * norm_angle(278.2   +  4493.34   * jy))
  rm += 0.0085 * Math.sin(Const::PI_180 * norm_angle(201.5   +  8266.71   * jy))
  rm += 0.0100 * Math.sin(Const::PI_180 * norm_angle( 44.89  + 14315.966  * jy))
  rm += 0.0107 * Math.sin(Const::PI_180 * norm_angle(336.44  + 13038.696  * jy))
  rm += 0.0110 * Math.sin(Const::PI_180 * norm_angle(231.59  +  4892.052  * jy))
  rm += 0.0125 * Math.sin(Const::PI_180 * norm_angle(141.51  + 14436.029  * jy))
  rm += 0.0153 * Math.sin(Const::PI_180 * norm_angle(130.84  +   758.698  * jy))
  rm += 0.0305 * Math.sin(Const::PI_180 * norm_angle(312.49  +  5131.979  * jy))
  rm += 0.0348 * Math.sin(Const::PI_180 * norm_angle(117.84  +  4452.671  * jy))
  rm += 0.0410 * Math.sin(Const::PI_180 * norm_angle(137.43  +  4411.998  * jy))
  rm += 0.0459 * Math.sin(Const::PI_180 * norm_angle(238.18  +  8545.352  * jy))
  rm += 0.0533 * Math.sin(Const::PI_180 * norm_angle( 10.66  + 13677.331  * jy))
  rm += 0.0572 * Math.sin(Const::PI_180 * norm_angle(103.21  +  3773.363  * jy))
  rm += 0.0588 * Math.sin(Const::PI_180 * norm_angle(214.22  +   638.635  * jy))
  rm += 0.1143 * Math.sin(Const::PI_180 * norm_angle(  6.546 +  9664.0404 * jy))
  rm += 0.1856 * Math.sin(Const::PI_180 * norm_angle(177.525 +   359.9905 * jy))
  rm += 0.2136 * Math.sin(Const::PI_180 * norm_angle(269.926 +  9543.9773 * jy))
  rm += 0.6583 * Math.sin(Const::PI_180 * norm_angle(235.700 +  8905.3422 * jy))
  rm += 1.2740 * Math.sin(Const::PI_180 * norm_angle(100.738 +  4133.3536 * jy))
  rm += 6.2887 * Math.sin(Const::PI_180 * norm_angle(134.961 +  4771.9886 * jy + am))
  rm += norm_angle(218.3161 + 4812.67881 * jy)
  return rm
end

.compute_lambda_sun(jy) ⇒ Object

太陽視黄経の計算

@param: jy (ユリウス年(JST)) @return: lambda



188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
 
# File 'lib/mk_sunmoon/compute.rb', line 188

def compute_lambda_sun(jy)
  rm  = 0.0003 * Math.sin(Const::K * norm_angle(329.7  +   44.43  * jy))
  rm += 0.0003 * Math.sin(Const::K * norm_angle(352.5  + 1079.97  * jy))
  rm += 0.0004 * Math.sin(Const::K * norm_angle( 21.1  +  720.02  * jy))
  rm += 0.0004 * Math.sin(Const::K * norm_angle(157.3  +  299.30  * jy))
  rm += 0.0004 * Math.sin(Const::K * norm_angle(234.9  +  315.56  * jy))
  rm += 0.0005 * Math.sin(Const::K * norm_angle(291.2  +   22.81  * jy))
  rm += 0.0005 * Math.sin(Const::K * norm_angle(207.4  +    1.50  * jy))
  rm += 0.0006 * Math.sin(Const::K * norm_angle( 29.8  +  337.18  * jy))
  rm += 0.0007 * Math.sin(Const::K * norm_angle(206.8  +   30.35  * jy))
  rm += 0.0007 * Math.sin(Const::K * norm_angle(153.3  +   90.38  * jy))
  rm += 0.0008 * Math.sin(Const::K * norm_angle(132.5  +  659.29  * jy))
  rm += 0.0013 * Math.sin(Const::K * norm_angle( 81.4  +  225.18  * jy))
  rm += 0.0015 * Math.sin(Const::K * norm_angle(343.2  +  450.37  * jy))
  rm += 0.0018 * Math.sin(Const::K * norm_angle(251.3  +    0.20  * jy))
  rm += 0.0018 * Math.sin(Const::K * norm_angle(297.8  + 4452.67  * jy))
  rm += 0.0020 * Math.sin(Const::K * norm_angle(247.1  +  329.64  * jy))
  rm += 0.0048 * Math.sin(Const::K * norm_angle(234.95 +   19.341 * jy))
  rm += 0.0200 * Math.sin(Const::K * norm_angle(355.05 +  719.981 * jy))
  rm += (1.9146 - 0.00005 * jy) * Math.sin(Const::K * norm_angle(357.538 + 359.991 * jy))
  rm += norm_angle(280.4603 + 360.00769 * jy)
  return norm_angle(rm)
end

.compute_moon(div) ⇒ Object

月の出・入・南中の計算

@param: div(0: 出, 1: 入, 2: 南中) @return: [time_str, hour_angle]



97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
 
# File 'lib/mk_sunmoon/compute.rb', line 97

def compute_moon(div)
  time_val = compute_time_moon(div)
  if time_val == 0.0
    time_str = "--:--:--"
    angle    = "---"
  else
    time_str = val2hhmmss(time_val * 24.0)
    jy = (@jd + time_val + @dt / 86400.0 - 2451545.0) / 365.25
    lambda = compute_lambda_moon(jy)
    beta   = compute_beta_moon(jy)
    if div == 2
      angle = compute_height_ecliptic(jy, time_val, lambda, beta)
    else
      angle = compute_angle_ecliptic(jy, time_val, lambda, beta)
    end
  end
  return [time_str, angle]
end

.compute_sidereal_time(jy, t) ⇒ Object

恒星時Θ(度)の計算

@param: jy (ユリウス年(JST)) @param: t (時刻(0.xxxx日)) @param: longitude @return: sidereal time



439
440
441
442
443
444
445
446
 
# File 'lib/mk_sunmoon/compute.rb', line 439

def compute_sidereal_time(jy, t)
  val  = 325.4606       + \
        (360.007700536  + \
        (  0.00000003879) \
        * jy) * jy        \
        + 360.0 * t + @lon
  return norm_angle(val)
end

.compute_sun(div) ⇒ Object

日の出・入・南中の計算

@param: div(0: 出, 1: 入, 2: 南中) @return: [time_str, hour_angle]



73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
 
# File 'lib/mk_sunmoon/compute.rb', line 73

def compute_sun(div)
  time_val = compute_time_sun(div)
  if time_val == 0.0
    time_str = "--:--:--"
    angle    = "---"
  else
    time_str = val2hhmmss(time_val * 24.0)
    jy = (@jd + time_val + @dt / 86400.0 - 2451545.0) / 365.25
    lambda = compute_lambda_sun(jy)
    if div == 2
      angle = compute_height_ecliptic(jy, time_val, lambda, 0.0)
    else
      angle = compute_angle_ecliptic(jy, time_val, lambda, 0.0)
    end
  end
  return [time_str, angle]
end

.compute_time_moon(div) ⇒ Object

月の出/月の入/月の南中計算

@param: div (0: 出, 1: 入, 2: 南中) @return: time (0.xxxx日)



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
 
# File 'lib/mk_sunmoon/compute.rb', line 156

def compute_time_moon(div)
  rev = 1    # 補正値初期値
  t   = 0.5  # 逐次計算時刻(日)初期設定

  while rev.abs > Const::CONVERGE
    jy = (@jd + t + @dt / 86400.0 - 2451545.0) / 365.25
    lambda = compute_lambda_moon(jy)
    beta   = compute_beta_moon(jy)
    alpha, delta = eclip2equat(jy, lambda, beta)
    unless div == 2  # 南中のときは計算しない
      diff = compute_diff_moon(jy)
      height = -1 * Const::REFRACTION - @incl + diff
    end
    time_sidereal = compute_sidereal_time(jy, t)
    hour_angle_diff = compute_hour_angle_diff(
      alpha, delta, time_sidereal, height, div
    )
    # 仮定時刻に対する補正値
    rev = hour_angle_diff / 347.8
    t += rev
  end
  # 月の出/月の入りがない場合は 0 とする
  t = 0.0 if t < 0.0 || t >= 1.0
  return t
end

.compute_time_sun(div) ⇒ Object

日の出/日の入/日の南中計算

@param: div (0: 出, 1: 入, 2: 南中) @return: time (0.xxxx日)



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
 
# File 'lib/mk_sunmoon/compute.rb', line 122

def compute_time_sun(div)
  rev = 1    # 補正値初期値
  t   = 0.5  # 逐次計算時刻(日)初期設定

  while rev.abs > Const::CONVERGE
    jy = (@jd + t + @dt / 86400.0 - 2451545.0) / 365.25
    lambda = compute_lambda_sun(jy)
    dist   = compute_dist_sun(jy)
    alpha, delta = eclip2equat(jy, lambda, 0.0)
    unless div == 2  # 南中のときは計算しない
      r = 0.266994 / dist
      diff = 0.0024428 / dist
      height = -1 * r - Const::REFRACTION - @incl + diff
    end
    time_sidereal = compute_sidereal_time(jy, t)
    hour_angle_diff = compute_hour_angle_diff(
      alpha, delta, time_sidereal, height, div
    )
    return 0.0 if hour_angle_diff == 0.0
    # 仮定時刻に対する補正値
    rev = hour_angle_diff / 360.0
    t += rev
  end
  # 日の出・入がない場合は 0 とする
  t = 0.0 if t < 0.0 || t >= 1.0
  return t
end

.eclip2equat(jy, lambda, beta) ⇒ Object

黄道座標 -> 赤道座標変換

@param: jy (ユリウス年(JST)) @param: lambda (黄経) @param: beta (黄緯) @return: [alpha(R.A.), delta(Decl.)]



415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
 
# File 'lib/mk_sunmoon/compute.rb', line 415

def eclip2equat(jy, lambda, beta)
  eps = (23.439291 - 0.000130042 * jy) * Const::PI_180
  lambda *= Const::PI_180
  beta   *= Const::PI_180
  a  =      Math.cos(beta) * Math.cos(lambda)
  b  = -1 * Math.sin(beta) * Math.sin(eps)
  b +=      Math.cos(beta) * Math.sin(lambda) * Math.cos(eps)
  c  =      Math.sin(beta) * Math.cos(eps)
  c +=      Math.cos(beta) * Math.sin(lambda) * Math.sin(eps)
  alpha  = b / a
  alpha  = Math.atan(alpha) / Const::PI_180
  alpha += 180 if a < 0
  delta  = Math.asin(c) / Const::PI_180
  return [alpha, delta]
end

.gc2jd(year, month, day, hour = 0, min = 0, sec = 0) ⇒ Object

Gregorian Calendar -> Julian Day

  • フリーゲルの公式を使用する

    JD

    int( 365.25 × year )

    + int( year / 400 )

    • int( year / 100 )

    + int( 30.59 ( month - 2 ) ) + day + 1721088

    ※上記の int( x ) は厳密には、x を超えない最大の整数

    (

@param: year @param: month @param: day @param: hour @param: minute @param: second @return: jd ( ユリウス日 )



28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
 
# File 'lib/mk_sunmoon/compute.rb', line 28

def gc2jd(year, month, day, hour = 0, min = 0, sec = 0)
  # 1月,2月は前年の13月,14月とする
  if month < 3
    year  -= 1
    month += 12
  end
  # 日付(整数)部分計算
  jd  = (365.25 * year).truncate
  jd += (year / 400.0).truncate
  jd -= (year / 100.0).truncate
  jd += (30.59 * (month - 2)).truncate
  jd += day
  jd += 1721088.125
  # 時間(小数)部分計算
  t  = sec / 3600.0
  t += min / 60.0
  t += hour
  t  = t / 24.0
  return jd + t
end

.norm_angle(angle) ⇒ Object

角度の正規化

@param: angle @return: angle



394
395
396
397
398
399
400
401
402
403
404
405
 
# File 'lib/mk_sunmoon/compute.rb', line 394

def norm_angle(angle)
  if angle < 0
    angle1  = angle * (-1)
    angle2  = (angle1 / 360.0).truncate
    angle1 -= 360 * angle2
    angle1  = 360 - angle1
  else
    angle1  = (angle / 360.0).truncate
    angle1  = angle - 360.0 * angle1
  end
  return angle1
end

.val2hhmmss(val) ⇒ Object

時刻：数値->時刻：時分変換 (xx.xxxx -> hh:mm:ss)

@param: time value (時刻, xx.xxxx日) @return: time string (時刻, hh:mm:ss)



487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
 
# File 'lib/mk_sunmoon/compute.rb', line 487

def val2hhmmss(val)
  val_h = val.truncate            # 整数部(時)
  val_2 = val - val_h             # 小数部
  val_m = (val_2 * 60).truncate   # (分)計算
  val_3 = val_2 - (val_m / 60.0)  # (秒)計算
  val_s = (val_3 * 60 * 60).round
  if val_s == 60
    val_s = 0
    val_m +=1
  end
  if val_m == 60
    val_m = 0
    val_h += 1
  end
  return sprintf("%02d:%02d:%02d", val_h, val_m, val_s)
end