日の出計算 その2 2012/10/24

Javaです。ここんとこずーと日の出の計算をおこなっていました。

計算結果を国立天文台の値とくらべてみます。


 結果は以下です。1分違いの不一致が188個でした。どうでしょうか。使えない値でしょうか。左の値が、国立天文台の時間です。

年月日 時分 比較
2012/1/1 6:50 6:50
2012/1/2 6:51 6:50 ×
2012/1/3 6:51 6:50 ×
2012/1/4 6:51 6:50 ×
2012/1/5 6:51 6:51
2012/1/6 6:51 6:51
2012/1/7 6:51 6:51
2012/1/8 6:51 6:51
2012/1/9 6:51 6:51
2012/1/10 6:51 6:51
2012/1/11 6:51 6:50 ×
2012/1/12 6:51 6:50 ×
2012/1/13 6:51 6:50 ×
2012/1/14 6:50 6:50
2012/1/15 6:50 6:50
2012/1/16 6:50 6:50
2012/1/17 6:50 6:49 ×
2012/1/18 6:49 6:49
2012/1/19 6:49 6:49
2012/1/20 6:49 6:48 ×
2012/1/21 6:48 6:48
2012/1/22 6:48 6:47 ×
2012/1/23 6:47 6:47
2012/1/24 6:47 6:46 ×
2012/1/25 6:46 6:46
2012/1/26 6:46 6:45 ×
2012/1/27 6:45 6:45
2012/1/28 6:45 6:44 ×
2012/1/29 6:44 6:43 ×
2012/1/30 6:43 6:43
2012/1/31 6:42 6:42
2012/2/1 6:42 6:41 ×
2012/2/2 6:41 6:41
2012/2/3 6:40 6:40
2012/2/4 6:39 6:39
2012/2/5 6:39 6:38 ×
2012/2/6 6:38 6:37 ×
2012/2/7 6:37 6:36 ×
2012/2/8 6:36 6:35 ×
2012/2/9 6:35 6:35
2012/2/10 6:34 6:34
2012/2/11 6:33 6:33
2012/2/12 6:32 6:32
2012/2/13 6:31 6:31
2012/2/14 6:30 6:30
2012/2/15 6:29 6:28 ×
2012/2/16 6:28 6:27 ×
2012/2/17 6:27 6:26 ×
2012/2/18 6:26 6:25 ×
2012/2/19 6:24 6:24
2012/2/20 6:23 6:23
2012/2/21 6:22 6:22
2012/2/22 6:21 6:21
2012/2/23 6:20 6:19 ×
2012/2/24 6:19 6:18 ×
2012/2/25 6:17 6:17
2012/2/26 6:16 6:16
2012/2/27 6:15 6:14 ×
2012/2/28 6:14 6:13 ×
2012/2/29 6:12 6:12
2012/3/1 6:11 6:11
2012/3/2 6:10 6:09 ×
2012/3/3 6:08 6:08
2012/3/4 6:07 6:07
2012/3/5 6:06 6:05 ×
2012/3/6 6:04 6:04
2012/3/7 6:03 6:03
2012/3/8 6:02 6:01 ×
2012/3/9 6:00 6:00
2012/3/10 5:59 5:58 ×
2012/3/11 5:57 5:57
2012/3/12 5:56 5:56
2012/3/13 5:55 5:54 ×
2012/3/14 5:53 5:53
2012/3/15 5:52 5:51 ×
2012/3/16 5:50 5:50
2012/3/17 5:49 5:49
2012/3/18 5:48 5:47 ×
2012/3/19 5:46 5:46
2012/3/20 5:45 5:44 ×
2012/3/21 5:43 5:43
2012/3/22 5:42 5:41 ×
2012/3/23 5:40 5:40
2012/3/24 5:39 5:39
2012/3/25 5:38 5:37 ×
2012/3/26 5:36 5:36
2012/3/27 5:35 5:34 ×
2012/3/28 5:33 5:33
2012/3/29 5:32 5:31 ×
2012/3/30 5:30 5:30
2012/3/31 5:29 5:29
2012/4/1 5:28 5:27 ×
2012/4/2 5:26 5:26
2012/4/3 5:25 5:24 ×
2012/4/4 5:23 5:23
2012/4/5 5:22 5:21 ×
2012/4/6 5:21 5:20 ×
2012/4/7 5:19 5:19
2012/4/8 5:18 5:17 ×
2012/4/9 5:16 5:16
2012/4/10 5:15 5:15
2012/4/11 5:14 5:13 ×
2012/4/12 5:12 5:12
2012/4/13 5:11 5:11
2012/4/14 5:10 5:09 ×
2012/4/15 5:08 5:08
2012/4/16 5:07 5:07
2012/4/17 5:06 5:05 ×
2012/4/18 5:05 5:04 ×
2012/4/19 5:03 5:03
2012/4/20 5:02 5:01 ×
2012/4/21 5:01 5:00 ×
2012/4/22 5:00 4:59 ×
2012/4/23 4:58 4:58
2012/4/24 4:57 4:57
2012/4/25 4:56 4:55 ×
2012/4/26 4:55 4:54 ×
2012/4/27 4:54 4:53 ×
2012/4/28 4:52 4:52
2012/4/29 4:51 4:51
2012/4/30 4:50 4:50
2012/5/1 4:49 4:49
2012/5/2 4:48 4:48
2012/5/3 4:47 4:46 ×
2012/5/4 4:46 4:45 ×
2012/5/5 4:45 4:44 ×
2012/5/6 4:44 4:43 ×
2012/5/7 4:43 4:42 ×
2012/5/8 4:42 4:42
2012/5/9 4:41 4:41
2012/5/10 4:40 4:40
2012/5/11 4:39 4:39
2012/5/12 4:38 4:38
2012/5/13 4:38 4:37 ×
2012/5/14 4:37 4:36 ×
2012/5/15 4:36 4:35 ×
2012/5/16 4:35 4:35
2012/5/17 4:35 4:34 ×
2012/5/18 4:34 4:33 ×
2012/5/19 4:33 4:33
2012/5/20 4:32 4:32
2012/5/21 4:32 4:31 ×
2012/5/22 4:31 4:31
2012/5/23 4:31 4:30 ×
2012/5/24 4:30 4:30
2012/5/25 4:30 4:29 ×
2012/5/26 4:29 4:28 ×
2012/5/27 4:29 4:28 ×
2012/5/28 4:28 4:28
2012/5/29 4:28 4:27 ×
2012/5/30 4:27 4:27
2012/5/31 4:27 4:26 ×
2012/6/1 4:27 4:26 ×
2012/6/2 4:26 4:26
2012/6/3 4:26 4:25 ×
2012/6/4 4:26 4:25 ×
2012/6/5 4:25 4:25
2012/6/6 4:25 4:25
2012/6/7 4:25 4:25
2012/6/8 4:25 4:24 ×
2012/6/9 4:25 4:24 ×
2012/6/10 4:25 4:24 ×
2012/6/11 4:25 4:24 ×
2012/6/12 4:25 4:24 ×
2012/6/13 4:25 4:24 ×
2012/6/14 4:25 4:24 ×
2012/6/15 4:25 4:24 ×
2012/6/16 4:25 4:24 ×
2012/6/17 4:25 4:24 ×
2012/6/18 4:25 4:24 ×
2012/6/19 4:25 4:25
2012/6/20 4:25 4:25
2012/6/21 4:25 4:25
2012/6/22 4:26 4:25 ×
2012/6/23 4:26 4:25 ×
2012/6/24 4:26 4:26
2012/6/25 4:27 4:26 ×
2012/6/26 4:27 4:26 ×
2012/6/27 4:27 4:27
2012/6/28 4:28 4:27 ×
2012/6/29 4:28 4:27 ×
2012/6/30 4:28 4:28
2012/7/1 4:29 4:28 ×
2012/7/2 4:29 4:29
2012/7/3 4:30 4:29 ×
2012/7/4 4:30 4:30
2012/7/5 4:31 4:30 ×
2012/7/6 4:31 4:31
2012/7/7 4:32 4:31 ×
2012/7/8 4:32 4:32
2012/7/9 4:33 4:32 ×
2012/7/10 4:33 4:33
2012/7/11 4:34 4:34
2012/7/12 4:35 4:34 ×
2012/7/13 4:35 4:35
2012/7/14 4:36 4:35 ×
2012/7/15 4:37 4:36 ×
2012/7/16 4:37 4:37
2012/7/17 4:38 4:37 ×
2012/7/18 4:39 4:38 ×
2012/7/19 4:39 4:39
2012/7/20 4:40 4:40
2012/7/21 4:41 4:40 ×
2012/7/22 4:41 4:41
2012/7/23 4:42 4:42
2012/7/24 4:43 4:42 ×
2012/7/25 4:44 4:43 ×
2012/7/26 4:44 4:44
2012/7/27 4:45 4:45
2012/7/28 4:46 4:45 ×
2012/7/29 4:47 4:46 ×
2012/7/30 4:47 4:47
2012/7/31 4:48 4:48
2012/8/1 4:49 4:48 ×
2012/8/2 4:50 4:49 ×
2012/8/3 4:50 4:50
2012/8/4 4:51 4:51
2012/8/5 4:52 4:52
2012/8/6 4:53 4:52 ×
2012/8/7 4:54 4:53 ×
2012/8/8 4:54 4:54
2012/8/9 4:55 4:55
2012/8/10 4:56 4:55 ×
2012/8/11 4:57 4:56 ×
2012/8/12 4:57 4:57
2012/8/13 4:58 4:58
2012/8/14 4:59 4:59
2012/8/15 5:00 4:59 ×
2012/8/16 5:01 5:00 ×
2012/8/17 5:01 5:01
2012/8/18 5:02 5:02
2012/8/19 5:03 5:02 ×
2012/8/20 5:04 5:03 ×
2012/8/21 5:05 5:04 ×
2012/8/22 5:05 5:05
2012/8/23 5:06 5:06
2012/8/24 5:07 5:06 ×
2012/8/25 5:08 5:07 ×
2012/8/26 5:08 5:08
2012/8/27 5:09 5:09
2012/8/28 5:10 5:09 ×
2012/8/29 5:11 5:10 ×
2012/8/30 5:11 5:11
2012/8/31 5:12 5:12
2012/9/1 5:13 5:12 ×
2012/9/2 5:14 5:13 ×
2012/9/3 5:14 5:14
2012/9/4 5:15 5:15
2012/9/5 5:16 5:15 ×
2012/9/6 5:17 5:16 ×
2012/9/7 5:17 5:17
2012/9/8 5:18 5:18
2012/9/9 5:19 5:18 ×
2012/9/10 5:20 5:19 ×
2012/9/11 5:20 5:20
2012/9/12 5:21 5:21
2012/9/13 5:22 5:21 ×
2012/9/14 5:23 5:22 ×
2012/9/15 5:23 5:23
2012/9/16 5:24 5:24
2012/9/17 5:25 5:25
2012/9/18 5:26 5:25 ×
2012/9/19 5:27 5:26 ×
2012/9/20 5:27 5:27
2012/9/21 5:28 5:28
2012/9/22 5:29 5:28 ×
2012/9/23 5:30 5:29 ×
2012/9/24 5:30 5:30
2012/9/25 5:31 5:31
2012/9/26 5:32 5:31 ×
2012/9/27 5:33 5:32 ×
2012/9/28 5:33 5:33
2012/9/29 5:34 5:34
2012/9/30 5:35 5:35
2012/10/1 5:36 5:35 ×
2012/10/2 5:37 5:36 ×
2012/10/3 5:37 5:37
2012/10/4 5:38 5:38
2012/10/5 5:39 5:38 ×
2012/10/6 5:40 5:39 ×
2012/10/7 5:41 5:40 ×
2012/10/8 5:41 5:41
2012/10/9 5:42 5:42
2012/10/10 5:43 5:43
2012/10/11 5:44 5:43 ×
2012/10/12 5:45 5:44 ×
2012/10/13 5:46 5:45 ×
2012/10/14 5:46 5:46
2012/10/15 5:47 5:47
2012/10/16 5:48 5:48
2012/10/17 5:49 5:49
2012/10/18 5:50 5:49 ×
2012/10/19 5:51 5:50 ×
2012/10/20 5:52 5:51 ×
2012/10/21 5:53 5:52 ×
2012/10/22 5:54 5:53 ×
2012/10/23 5:54 5:54
2012/10/24 5:55 5:55
2012/10/25 5:56 5:56
2012/10/26 5:57 5:57
2012/10/27 5:58 5:58
2012/10/28 5:59 5:59
2012/10/29 6:00 5:59 ×
2012/10/30 6:01 6:00 ×
2012/10/31 6:02 6:01 ×
2012/11/1 6:03 6:02 ×
2012/11/2 6:04 6:03 ×
2012/11/3 6:05 6:04 ×
2012/11/4 6:06 6:05 ×
2012/11/5 6:07 6:06 ×
2012/11/6 6:08 6:07 ×
2012/11/7 6:09 6:08 ×
2012/11/8 6:10 6:09 ×
2012/11/9 6:11 6:10 ×
2012/11/10 6:12 6:11 ×
2012/11/11 6:13 6:12 ×
2012/11/12 6:14 6:13 ×
2012/11/13 6:15 6:14 ×
2012/11/14 6:16 6:15 ×
2012/11/15 6:17 6:16 ×
2012/11/16 6:18 6:17 ×
2012/11/17 6:19 6:18 ×
2012/11/18 6:20 6:19 ×
2012/11/19 6:21 6:20 ×
2012/11/20 6:22 6:21 ×
2012/11/21 6:23 6:22 ×
2012/11/22 6:24 6:23 ×
2012/11/23 6:25 6:24 ×
2012/11/24 6:25 6:25
2012/11/25 6:26 6:26
2012/11/26 6:27 6:27
2012/11/27 6:28 6:28
2012/11/28 6:29 6:29
2012/11/29 6:30 6:30
2012/11/30 6:31 6:31
2012/12/1 6:32 6:32
2012/12/2 6:33 6:32 ×
2012/12/3 6:34 6:33 ×
2012/12/4 6:35 6:34 ×
2012/12/5 6:36 6:35 ×
2012/12/6 6:36 6:36
2012/12/7 6:37 6:37
2012/12/8 6:38 6:38
2012/12/9 6:39 6:38 ×
2012/12/10 6:40 6:39 ×
2012/12/11 6:40 6:40
2012/12/12 6:41 6:41
2012/12/13 6:42 6:41 ×
2012/12/14 6:42 6:42
2012/12/15 6:43 6:43
2012/12/16 6:44 6:43 ×
2012/12/17 6:44 6:44
2012/12/18 6:45 6:45
2012/12/19 6:46 6:45 ×
2012/12/20 6:46 6:46
2012/12/21 6:47 6:46 ×
2012/12/22 6:47 6:47
2012/12/23 6:48 6:47 ×
2012/12/24 6:48 6:48
2012/12/25 6:49 6:48 ×
2012/12/26 6:49 6:48 ×
2012/12/27 6:49 6:49
2012/12/28 6:50 6:49 ×
2012/12/29 6:50 6:49 ×
2012/12/30 6:50 6:50
2012/12/31 6:50 6:50



不一致:188

ソースコード

import static java.lang.Math.*;

import java.text.SimpleDateFormat;
import java.util.Calendar;
import java.util.Date;

/**
 *
 * @author nakawakashigeto
 *
 */
public class A20121021 {

    public static void main(String[] args) {
        Calendar calendar = Calendar.getInstance();
        calendar.set(2012, 1 - 1, 1, 0, 0, 0);
      
      
        // 検証のために1年分の出力
        while(true){
            int yyyy = calendar.get(Calendar.YEAR);
            int m = calendar.get(Calendar.MONTH) + 1;
            int d = calendar.get(Calendar.DAY_OF_MONTH);
            sunrise("東京", yyyy, m, d, 35.6581 , 139.7414);
            // 一日進める
            calendar.add(Calendar.DAY_OF_YEAR, 1);
            if(2013 <= calendar.get(Calendar.YEAR)) break;
        }
    }
  
    /**
     * 日付表現クラス
     *
     */
    static class YMDH {
        int Y, M, D, H;

        public YMDH(int yyyy, int mm, int dd, int hh) {
            this.Y = yyyy;
            this.M = mm;
            this.D = dd;
            this.H = hh;
        }

        public String toString() {
            StringBuilder builder = new StringBuilder();
            builder.append("y:[").append(Y).append("] ").append("m:[")
                    .append(M).append("] ").append("d:[").append(D)
                    .append("] ").append("h:[").append(H).append("] ");
            return builder.toString();
        }
    }

    /**
     * 角度を表現するクラス
     */
    static class Degree {
        double deg;
        public Degree(double argdeg) {
            this.deg = argdeg;
        }
        public double toR() {
            return toRadians(deg);
        }
        public String toString() {
            StringBuilder builder = new StringBuilder();
            builder.append("deg:[").append(deg).append("] ");
            return builder.toString();
        }
    }

    /**
     * 経緯度を表現するクラス
     *
     */
    static class LongitudeAndLatitude {
        /** 緯度 */
        Degree latitude;
        /** 経度 */
        Degree longitude;

        public LongitudeAndLatitude(double a, double b) {
            this.latitude = new Degree(a);
            this.longitude = new Degree(b);
        }

        public String toString() {
            StringBuilder builder = new StringBuilder();
            builder.append("lat:[").append(latitude.deg).append("] ")
                    .append("lng[").append(longitude.deg).append("]");
            return builder.toString();
        }
    }

    /**
     * 太陽の赤経、赤緯、距離、恒星時をまとめて扱うためのクラス
     *
     */
    static class Sun {
        /** ra:赤経 α right ascension */
        Degree ra;
        /** δ:赤緯 デルタ declination */
        Degree dec;
        /** r:距離 */
        double r;
        /** th:恒星時 θ シータ sidereal time */
        Degree th;

        public String toString() {
            StringBuilder builder = new StringBuilder();
            builder.append("a:[").append(ra).append("] ").append("δ:[")
                    .append(dec).append("] ").append("r:[").append(r)
                    .append("] ").append("th:[").append(th).append("] ");
            return builder.toString();
        }
    }

    /**
     * 日の出計算を行う
     *
     * @param location
     * @param Y
     * @param M
     * @param D
     * @param longitude
     * @param latitude
     */
    static void sunrise(String location, int Y, int M, int D, double longitude,
            double latitude) {

        YMDH ymdh = new YMDH(Y, M, D, 6);
        LongitudeAndLatitude longitudeAndLatitude = new LongitudeAndLatitude(
                longitude, latitude);
        double d = calcSunrise(ymdh, longitudeAndLatitude);

        // 補正なぜこれが必要になるのか?
        // 書籍には少数点の値のみでよいと書いてある
        if(d >= 1)
            d -= floor(d);
        // System.out.println(d);
        d = d * 24.;
        double d2 = d % 1;

        String s = location + "\t" + ymdh.Y + "年" + ymdh.M + "月"
                + ymdh.D + "日" + (int) floor(d) + "時" + (int) floor(d2 * 60)
                + "分 \t";
        String s2 = String.format("%d:%02d", (int) floor(d),  (int) floor(d2 * 60));
        System.out.println(s + s2);
    }

    /**
     * 仮定時刻に対する日の出計算を行い、補正値を返す。
     *
     * @param ymdh
     * @param longitudeAndLatitude
     * @return
     */
    static double calcSunrise(YMDH ymdh, LongitudeAndLatitude longitudeAndLatitude) {

        // (1) 日の出時刻を仮定し、日単位の時刻変数dに換算する
        // d:時刻変数
        double d = ymdh.H / 24.;
      
        // (2) 太陽の赤経、赤緯、恒星時を求める
        Sun sun = a2(ymdh.Y, ymdh.M, ymdh.D, ymdh.H, d,
                longitudeAndLatitude.longitude.deg);

        // (3) 出没高度k
        double k = k(sun.r);

        // (4) 時角tk
        double tk = tk(sun.dec.deg, longitudeAndLatitude.latitude.deg, k);

        tk = toDegrees(acos(tk));

        // (5) 恒星時から太陽の時角tを計算する
        double t = t(sun.ra.deg, sun.th.deg);
        //
        if (0 < tk)
            tk *= -1;
      
        // (6) 仮定時刻に対する補正値Δdの計算をする
        double dd = deltaD(tk, t);
        // System.out.print("ra:" + (int)sun.ra.deg +"\tt:"+ t + "\ttk:" + tk + "\td:" + d + " dd:" + dd + "   ");
        // 補正値を加えて新しく仮定時刻を決定する
        double newd = d + dd;//  補正値

        return newd;
    }

    /**
     *
     * @param Y
     * @param M
     * @param D
     * @param H
     * @param λ2
     * @return
     */
    static Sun a2(int Y, int M, int D, int H, double d, double λ2) {

            // T 時刻変数
            double T = 時刻変数T(Y, M, D);

            double λs = λs(T);

            double q = q(T);

            double r = r(q);

            double e = toRadians(23.439291) - toRadians(0.000130042) * T;

            double tana = tan(λs) * cos(e);
            double sinδ = sin(λs) * sin(e);

            double _λs = toDegrees(λs) % 360.;
            tana = atan(tana);

            if (0 <= _λs && _λs < 90) {
                tana += toRadians(360);
            } else if (90 <= _λs && _λs < 180) {
                tana += toRadians(180);
            } else if (180 <= _λs && _λs < 270) {
                tana += toRadians(180);
            } else if (270 <= _λs && _λs < 360) {
                tana += toRadians(360);
            }

            // System.out.print((int)toDegrees(tana)+"\t");
            sinδ = asin(sinδ);
          
            // 恒星時を求める
            double th = siderealTimeθ(T, d, λ2);

            Sun sun = new Sun();
            {

                sun.ra = new Degree(toDegrees(tana));
                sun.dec = new Degree(toDegrees(sinδ));
                sun.r = r;
                sun.th = new Degree(th);
            }
            // System.out.print(sun+"\t");
            return sun;
    }

    /**
     * 仮定時刻に対する補正値Δdの計算
     *
     * @param r
     * @return
     */
    static double deltaD(double tk, double t) {
        double deltaD = 0;
        deltaD = (tk - t) / 360.;
        return deltaD;
    }

    /**
     * 太陽の出没高度kを計算する
     *
     * @param r
     * @return
     */
    static double k(double r) {
        // 太陽の視半径
        final double S = 0.266994 / r;
        // 大気差
        final double R = 0.585556;
        // 視差
        final double Π = 0.0024428 / r;
        //
        final double E0 = 0;

        double k = -S - E0 - R + Π;
        return k;
    }

    /**
     * 時角を求める φ
     *
     * @return
     */
    static double tk(double delta, double phi, double k) {
        double tk = 0;
        double sink = sin(toRadians(k));
        double cosDelta = cos(toRadians(delta));
        double sinDelta = sin(toRadians(delta));
        double cosPhi = cos(toRadians(phi));
        double sinPhi = sin(toRadians(phi));
        tk = (sink - sinDelta * sinPhi) / (cosDelta * cosPhi);
        return tk;
    }

    /**
     *
     * @param delta
     * @param phi
     * @param k
     * @return
     */
    static double tk(Degree delta, Degree phi, Degree k) {
        double tk = 0;
        double sink = sin(k.toR());
        double cosDelta = cos(delta.toR());
        double sinDelta = sin(delta.toR());
        double cosPhi = cos(phi.toR());
        double sinPhi = sin(phi.toR());
        tk = (sink - sinDelta * sinPhi) / (cosDelta * cosPhi);
        return tk;
    }

    /**
     * 恒星時から太陽の時角tを計算する
     *
     * @param delta
     * @param phi
     * @param k
     * @return
     */
    static double t(double a, double th) {
        double t = 0;
        t = th - a;
        return t;
    }

    /**
     * 距離 r
     * @param q
     * @return
     */
    static double r(double q) {
        return pow(10., q);
    }

    /**
     * 恒星時
     *
     * @param T
     * @param d
     * @param λ
     * @return
     */
    static double siderealTimeθ(double T, double d, double λ) {
        double th = 0;

        th = 325.4606 + 360.007700536 * T + 0.00000003879 * pow(T, 2) + 360.
                * d + λ;

        double th2 = th;

        if (th2 < 0) {
            while (th2 < 0) {
                th2 += 360;
            }
        } else if (360 <= th2) {
            while (360 <= th2) {
                th2 -= 360;
            }
        }
        return th2;
    }

    static double q(double T) {
        double[] L = new double[6];
        L[0] = (0.007256 - 0.0000002 * T)
                * sin(toRadians(267.54) + toRadians(359.991) * T);
        L[1] = (0.000091) * sin(toRadians(265.1) + toRadians(719.98) * T);
        L[2] = (0.000030) * sin(toRadians(90.0));
        L[3] = (0.000013) * sin(toRadians(27.8) + toRadians(4452.67) * T);
        L[4] = (0.000007) * sin(toRadians(254) + toRadians(450.4) * T);
        L[5] = (0.000007) * sin(toRadians(156) + toRadians(329.6) * T);
        double q = 0;
        for (double d : L) {
            q += d;
        }
        return q;
    }

    static double λs(double T) {
        double[] L = new double[20];
        L[0] = toRadians(280.4603) + toRadians(360.00769) * T;
        L[1] = (toRadians(1.9146) - toRadians(0.00005))
                * sin(toRadians(357.538) + toRadians(359.991) * T);
        L[2] = (toRadians(0.0200))
                * sin(toRadians(355.05) + toRadians(719.981) * T);
        L[3] = (toRadians(0.0048))
                * sin(toRadians(234.95) + toRadians(19.341) * T);
        L[4] = (toRadians(0.0020))
                * sin(toRadians(247.1) + toRadians(329.64) * T);
        L[5] = (toRadians(0.0018))
                * sin(toRadians(297.8) + toRadians(4452.67) * T);
        L[6] = (toRadians(0.0018))
                * sin(toRadians(251.3) + toRadians(0.20) * T);
        L[7] = (toRadians(0.0015))
                * sin(toRadians(343.2) + toRadians(450.37) * T);
        L[8] = (toRadians(0.0013))
                * sin(toRadians(81.4) + toRadians(225.18) * T);
        L[9] = (toRadians(0.0008))
                * sin(toRadians(132.5) + toRadians(659.29) * T);
        L[10] = (toRadians(0.0007))
                * sin(toRadians(153.3) + toRadians(90.38) * T);
        L[11] = (toRadians(0.0007))
                * sin(toRadians(206.8) + toRadians(30.35) * T);
        L[12] = (toRadians(0.0006))
                * sin(toRadians(29.8) + toRadians(337.18) * T);
        L[13] = (toRadians(0.0005))
                * sin(toRadians(207.4) + toRadians(1.50) * T);
        L[14] = (toRadians(0.0005))
                * sin(toRadians(291.2) + toRadians(22.81) * T);
        L[15] = (toRadians(0.0004))
                * sin(toRadians(234.9) + toRadians(315.56) * T);
        L[16] = (toRadians(0.0004))
                * sin(toRadians(157.3) + toRadians(299.30) * T);
        L[17] = (toRadians(0.0004))
                * sin(toRadians(21.1) + toRadians(720.02) * T);
        L[18] = (toRadians(0.0003))
                * sin(toRadians(352.5) + toRadians(1079.97) * T);
        L[19] = (toRadians(0.0003))
                * sin(toRadians(329.7) + toRadians(44.43) * T);

        double λs = 0;
        for (double d : L) {
            λs += d;
        }
        return λs;
    }

    static double 時刻変数T(int Y, int M, int D) {
        double deltaT = deltaT(Y, M);
        double T = 0.;
        double K = 経過日数K(Y, M, D, deltaT);
        T = K / 365.25;
        return T;
    }

    static double 経過日数K(int Y, int M, int D, double deltaT) {
        // 1月と2月は前年の13月、14月として考える
        if (M == 1 || M == 2) {
            M += 12;
            Y -= 1;
        }
        // J2000.0からの経過日数を求める
        // (2000 + Y)年M月D日として計算する
        Y = Y - 2000;

        return 経過日数K2(Y, M, D, deltaT);
    }

    static double 経過日数K2(int Y, int M, int D, double deltaT) {

        double K = 0;
        double KK = 0;
        KK = 365 * Y + 30 * M + D - 33.875 + Math.floor(3 * (M + 1) / 5.)
                + Math.floor(Y / 4.);
        K = KK + 6 / 24. + deltaT / 86400;
        return K;
    }

    static double deltaT(int yyyy, int mm) {
        // y = year + (month - 0.5)/12
        double y = yyyy + (mm - 0.5) / 12.;

        if (1961 <= y && y <= 1986) {
            // Between years 1961 and 1986, calculate:
            // ΔT = 45.45 + 1.067*t - t^2/260 - t^3 / 718
            // where: t = y - 1975
            double t = y - 1975;
            double dT = 45.45 + 1.067 * t - Math.pow(t, 2) / 260.
                    - Math.pow(t, 3) / 718.;
            return dT;
        }

        if (1986 <= y && y <= 2005) {
            // Between years 1986 and 2005, calculate:
            // ΔT = 63.86 + 0.3345 * t - 0.060374 * t^2 + 0.0017275 * t^3 +
            // 0.000651814 * t^4
            // + 0.00002373599 * t^5
            // where: t = y - 2000
            double t = y - 2000;
            double dT = 63.86 + 0.3345 * t - 0.060374 * Math.pow(t, 2)
                    + 0.0017275 * Math.pow(t, 3) + 0.000651814 * Math.pow(t, 4)
                    + 0.00002373599 * Math.pow(t, 5);
            return dT;
        }

        if (2005 <= y && y <= 2050) {
            // Between years 2005 and 2050, calculate:
            // ΔT = 62.92 + 0.32217 * t + 0.005589 * t^2
            // where: t = y - 2000
            double t = y - 2000;
            double dT = 62.92 + 0.32217 * t + 0.005589 * Math.pow(t, 2);
            return dT;
        }

        return 0;
    }

}

: