天体の位置を計算するための方法が、海上保安庁から提供されています。
計算してみましょう！！
今回は面倒なので太陽だけですが、太陽、金星、火星、木星、土星、月、といった、代表的な天体が含まれています。
天体の計算が、なぜ国立天文台ではなく、海上保安庁なのか?
GPSなんてものがなかった時代、広い海の上で自分の居場所を知るには、天測航法に頼る必要がありました。そのため、各国の海軍やそれに相当する機関が天文観測に注力していました。
測位は現在でこそGPS等の衛星ベースの技術にその座を奪われ、その運用も米では宇宙軍が、日本では内閣府が行っていますが、かつての航法システム、例えばオメガやロランといったシステムは、各国の海軍や沿岸警備隊、日本では海上保安庁といった、海に関わる機関が運用していました。
ソースコード
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 }
}));
}
}
定数をソースコードに埋め込んでいるため、そのあたりが複雑ですが、計算自体は結構簡単です。
使用例
Chartでグラフ化のサンプル
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");
太陽の赤緯がゼロを通過してプラスに変わったところが春分、そのときの太陽の方向が赤経ゼロで、その日が春分の日です。同様に、赤緯がマイナスに変わる日が秋分の日です。
太陽の赤緯が最大になる日が夏至、最小になる日が冬至です。
また、太陽と地球の距離は、1月上旬を最小とし、7月上旬を最大としています。
太陽から離れれば離れるほど、エネルギーは小さくなりますから、エネルギー的に言えば、実は夏より冬のほうが、太陽から受け取るエネルギーは大きいのです(エネルギーの比でおよそ7%程度)。そのため、人工衛星の太陽電池の能力は、夏至(ほぼ太陽から一番遠い時期)の発電量で表されていたりします。
その他
ちょうど1年前にも同じようなエントリ書いてるんだな…… 書いてから気がついた。。。