天体の位置を計算してみた

天体の位置を計算するための方法が、海上保安庁から提供されています。

[天文・暦情報](https://www1.kaiho.mlit.go.jp/KOHO/)

計算してみましょう！！

今回は面倒なので太陽だけですが、太陽、金星、火星、木星、土星、月、といった、代表的な天体が含まれています。

天体の計算が、なぜ国立天文台ではなく、海上保安庁なのか?
GPSなんてものがなかった時代、広い海の上で自分の居場所を知るには、天測航法に頼る必要がありました。そのため、各国の海軍やそれに相当する機関が天文観測に注力していました。
測位は現在でこそGPS等の衛星ベースの技術にその座を奪われ、その運用も米では宇宙軍が、日本では内閣府が行っていますが、かつての航法システム、例えばオメガやロランといったシステムは、各国の海軍や沿岸警備隊、日本では海上保安庁といった、海に関わる機関が運用していました。

## ソースコード

```C#
class CelestialPositionCalculation
{
    public enum Type
    {
        SunRa,
        SunDec,
        SunAu
    }

    public struct Data
    {
        public int Year { get; set; }
        public int a { get; set; }
        public int b { get; set; }
        public double[] Value { get; set; }

        public Data(int Year, int a, int b, double[] Value)
        {
            this.Year = Year;
            this.a = a;
            this.b = b;
            this.Value = Value;
        }
    }

    private static readonly ReadOnlyDictionary<Type, Data[]> dataDic = generateValuesetDic();
    private static readonly ReadOnlyDictionary<int, int> deltaTDic = generateDeltaTdic();

    public static double Calc(DateTimeOffset Date, Type title)
    {
        double t = CalcT(Date);

        Data page;
        if (!searchPage(title, Date, t, out page))
        {
            return (double.NaN);
        }

        int a = page.a;
        int b = page.b;
        double theta = Math.Acos((2.0 * t - (a + b)) / (b - a));
        double f = page.Value.Select((tmp, i) => tmp * Math.Cos(i * theta)).Sum();

        return (f);
    }

    public static double CalcT(DateTimeOffset Date)
    {
        DateTime UTC = Date.UtcDateTime;

        int P = UTC.Month - 1;
        int Q = (int)Math.Floor((UTC.Month + 7) / 10.0);
        int Y = (int)Math.Floor(UTC.Year / 4.0 - Math.Floor(UTC.Year / 4.0) + 0.77);
        int S = (int)Math.Floor(P * 0.55 - 0.33);
        int T = 30 * P + Q * (S - Y) + P * (1 - Q) + UTC.Day;
        double F = UTC.Hour / 24.0 + UTC.Minute / 1440.0 + UTC.Second / 86400.0;
        double t = T + F + deltaTDic[UTC.Year] / 86400.0;

        return (t);
    }

    private static bool searchPage(Type title, DateTimeOffset Date, double t, out Data page)
    {
        if (dataDic.ContainsKey(title))
        {
            var tmp = dataDic[title].Where(a => a.Year == Date.UtcDateTime.Year && a.a <= t && t <= a.b);

            if (tmp.Any())
            {
                page = tmp.First();
                return (true);
            }
        }

        page = new Data();
        return (false);
    }

    private static ReadOnlyDictionary<Type, Data[]> generateValuesetDic()
    {
        return (
            new ReadOnlyDictionary<Type, Data[]>(new Dictionary<Type, Data[]>()
            {
                {
                    Type.SunRa, new Data[] {
                        new Data(2020, 0, 122, new[] {
                            22.71187, 3.929287, -0.104071, 0.035631, 0.006885,
                            -0.002385, 0.000138, 0.00005, -0.000056, 0.000014,
                            -0.000054, -0.000018, 0.000062, 0.000008, -0.000027,
                            -0.000002, 0.000003, 0, }),
                        new Data(2020, 121, 245, new[] {
                            6.6489, 4.134433, -0.04333, -0.04038, 0.005603,
                            0.003175, -0.000428, -0.000194, 0.000005, 0.000005,
                            -0.000046, 0.000041, 0.000056, -0.000028, -0.000029,
                            0.00001, 0.000007, -0.000004, }),
                        new Data(2020, 244, 367, new[] {
                            14.575618, 4.056804, 0.152148, 0.012255, -0.01366,
                            -0.002129, 0.000302, 0.000116, 0.000024, -0.000018,
                            0, 0.000069, 0.000002, -0.000045, -0.000002,
                            0.000014, -0.000001, -0.000005, }),
                    } },
                {
                    Type.SunDec, new Data[] {
                        new Data(2020, 0, 122, new [] {
                            -5.75369, 20.14227, 1.7549, -1.00916, 0.01301,
                            0.00848, -0.00361, 0.00073, 0.00003, -0.00011,
                            -0.00025, -0.00014, 0.00027, 0.00016, -0.00012,
                            -0.00008, 0.00001, 0.00002, }),
                        new Data(2020, 121, 245, new []
                        {
                            17.14314, -3.5435, -5.77996, 0.22272, 0.1592,
                            -0.01133, -0.00458, 0.00086, 0.00031, 0.00019,
                            -0.00009, -0.00001, -0.00003, -0.00011, 0.00004,
                            0.00007, -0.00001, -0.00004, }),
                        new Data(2020, 244, 367, new []
                        {
                            -10.75539, -16.72221, 3.56387, 0.97158, -0.02871,
                            -0.02364, -0.00538, 0.00036, 0.00024, 0.00008,
                            0.00011, -0.00034, -0.00003, 0.00028, -0.00001,
                            -0.00012, 0.00001, 0.00004, }),
                    } },
                {
                    Type.SunAu, new Data[] {
                        new Data(2020, 0, 122, new []
                        {
                            0.993193, 0.012806, 0.002311, -0.000662, -0.000048,
                            0.000004, 0.000002, -0.000014, -0.000002, -0.000004,
                            -0.000004, 0.000017, 0.000004, -0.000011, -0.000002,
                            0.000004, 0.000001, -0.000001, }),
                        new Data(2020, 121, 245, new []
                        {
                            1.012388, 0.000983, -0.004213, -0.000046, 0.000088,
                            -0.000009, -0.000003, -0.000012, 0.000006, -0.000001,
                            0.000007, 0.000015, -0.000009, -0.000011, 0.000004,
                            0.000004, -0.000002, -0.000001, }),
                        new Data(2020, 244, 367, new []
                        {
                            0.994543, -0.013788, 0.00185, 0.000719, -0.000052,
                            -0.000012, -0.000006, 0, 0.000009, 0,
                            0.000014, 0, -0.000017, 0, 0.000008,
                            0, -0.000002, 0, }),
                    } },
            }));
    }

    private static ReadOnlyDictionary<int, int> generateDeltaTdic()
    {
        return (new ReadOnlyDictionary<int, int>(new Dictionary<int, int>() {
                { 2020, 70 }
            }));
    }
}
```

定数をソースコードに埋め込んでいるため、そのあたりが複雑ですが、計算自体は結構簡単です。

## 使用例

<details><summary>Chartでグラフ化のサンプル</summary><div>

```C#
Font font = new Font("", 20, FontStyle.Bold, GraphicsUnit.Pixel);

DateTimeOffset start = new DateTimeOffset(new DateTime(2020, 1, 1), new TimeSpan());
DateTimeOffset end = new DateTimeOffset(new DateTime(2021, 1, 1), new TimeSpan());

Chart chart = new Chart();
chart.Size = new Size(1920, 1080);

ChartArea area = new ChartArea();
chart.ChartAreas.Add(area);

Legend legend = new Legend();
chart.Legends.Add(legend);
legend.Docking = Docking.Top;

Series marker = new Series();
chart.Series.Add(marker);
marker.IsVisibleInLegend = false;
marker.ChartType = SeriesChartType.Point;
marker.MarkerSize = 20;
marker.MarkerStyle = MarkerStyle.Circle;
marker.Font = font;
marker.LabelBackColor = Color.LightGray;

Series ra = new Series("R.A. [hour]");
chart.Series.Add(ra);
ra.ChartType = SeriesChartType.Line;
ra.BorderWidth = 3;

Series dec = new Series("Dec. [deg]");
chart.Series.Add(dec);
dec.ChartType = SeriesChartType.Line;
dec.BorderWidth = 3;

Series au = new Series("Dist. [AU]");
chart.Series.Add(au);
au.YAxisType = AxisType.Secondary;
au.ChartType = SeriesChartType.Line;
au.BorderWidth = 3;

area.AxisX.Minimum = start.UtcDateTime.ToOADate();
area.AxisX.Maximum = end.UtcDateTime.ToOADate();
area.AxisX.Title = "UTC";
area.AxisX.TitleFont = font;
area.AxisX.LabelStyle.Format = "yyyy/MM";
area.AxisX.IntervalType = DateTimeIntervalType.Months;
area.AxisX.Interval = 1;
area.AxisX.MinorGrid.Enabled = true;
area.AxisX.MinorGrid.IntervalType = DateTimeIntervalType.Weeks;
area.AxisX.MinorGrid.Interval = 1;
area.AxisX.MinorGrid.LineColor = Color.LightGray;

area.AxisY.Title = "[hour], [deg]";
area.AxisY.TitleFont = font;
area.AxisY.Minimum = -25;
area.AxisY.Maximum = +25;
area.AxisY.Interval = 5;
area.AxisY.MinorGrid.Enabled = true;
area.AxisY.MinorGrid.Interval = 1;
area.AxisY.MinorGrid.LineColor = Color.LightGray;

area.AxisY2.Title = "[AU]";
area.AxisY2.TitleFont = font;
area.AxisY2.Minimum = 0.98;
area.AxisY2.Maximum = 1.03;
area.AxisY2.Interval = 0.005;

for (DateTimeOffset date = start; date < end; date += new TimeSpan(1, 0, 0))
{
    double oad = date.UtcDateTime.ToOADate();

    ra.Points.AddXY(oad, CelestialPositionCalculation.Calc(date, CelestialPositionCalculation.Type.SunRa) % 24);
    dec.Points.AddXY(oad, CelestialPositionCalculation.Calc(date, CelestialPositionCalculation.Type.SunDec));
    au.Points.AddXY(oad, CelestialPositionCalculation.Calc(date, CelestialPositionCalculation.Type.SunAu));
}

{
    DateTimeOffset date;

    date = new DateTimeOffset(new DateTime(2020, 3, 20, 12, 50, 0));
    marker.Points.Add(new DataPoint(date.UtcDateTime.ToOADate(),
        CelestialPositionCalculation.Calc(date, CelestialPositionCalculation.Type.SunDec))
    { Label = "春分の日", });

    date = new DateTimeOffset(new DateTime(2020, 6, 21, 6, 44, 0));
    marker.Points.Add(new DataPoint(date.UtcDateTime.ToOADate(),
        CelestialPositionCalculation.Calc(date, CelestialPositionCalculation.Type.SunDec))
    { Label = "夏至", });

    date = new DateTimeOffset(new DateTime(2020, 9, 22, 22, 31, 0));
    marker.Points.Add(new DataPoint(date.UtcDateTime.ToOADate(),
        CelestialPositionCalculation.Calc(date, CelestialPositionCalculation.Type.SunDec))
    { Label = "秋分の日", });

    date = new DateTimeOffset(new DateTime(2020, 12, 21, 19, 2, 0));
    marker.Points.Add(new DataPoint(date.UtcDateTime.ToOADate(),
        CelestialPositionCalculation.Calc(date, CelestialPositionCalculation.Type.SunDec))
    { Label = "冬至", });
}

Bitmap bmp = new Bitmap(chart.Width, chart.Height);
chart.DrawToBitmap(bmp, new Rectangle() { Size = bmp.Size });
bmp.Save("log.png");
```
</div></details>

![log.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/264452/fb0bac84-76d3-f74e-91c6-155c8cfc9d87.png)

太陽の赤緯がゼロを通過してプラスに変わったところが春分、そのときの太陽の方向が赤経ゼロで、その日が春分の日です。同様に、赤緯がマイナスに変わる日が秋分の日です。
太陽の赤緯が最大になる日が夏至、最小になる日が冬至です。

また、太陽と地球の距離は、1月上旬を最小とし、7月上旬を最大としています。
太陽から離れれば離れるほど、エネルギーは小さくなりますから、エネルギー的に言えば、実は夏より冬のほうが、太陽から受け取るエネルギーは大きいのです(エネルギーの比でおよそ7%程度)。そのため、人工衛星の太陽電池の能力は、夏至(ほぼ太陽から一番遠い時期)の発電量で表されていたりします。

## その他

　ちょうど1年前にも同じようなエントリ書いてるんだな…… 書いてから気がついた。。。