Arduinoで日の出・日の入り時刻を計算

　先日コメントで Arduino で日の出・日の入り時刻の計算が出来ないかという問い合わせを頂きました。日の出や日の入り時刻を計算で求めることが出来れば、いろいろなことに応用出来るので面白いテーマです。

　ということで、調べてみたのですが、これ結構厄介でした。ネットにいろんな情報があるのですが、この話の教科書的な書物として、以下の本が有名なので図書館で借りてきました。

▼日の出日の入りの計算　（長沢　工著）

　実はこの本にも私が求めているズバリの答えは書いてありませんでした。でも、読んでみると頭の中がかなり整理出来ました。日の出、日の入り時刻を求める方法には大きく以下の三つのアプローチがあるようです。なお、ここで前提としている課題は、日付と経緯度が判っている時に、日の出・日の入り時刻を知る方法です。

1) 基準地点の日の出・日の入りデーターを持っておき、経緯度の差から補正する方法。
2) 地球と太陽の相対位置を軌道計算から求める方法。（ユリウス暦日が出てきます）
3) 太陽位置を表す近似式を使う方法。

1) 項の方法は簡単な計算で出来そうです。元のデーターは適当に間引いて持っておき、使う時に補間すればあまりデーター量が増えなくても済みそうです。ただ、経緯度が大きく違うと計算誤差が大きくなりそうです。

2) 項の方法は精度の高い計算結果が得られますが計算がやっかいです。それに軌道要素は少しずつ変わるので、最新の理科年表の値を参照する必要があり、マイコンのプログラムとして埋め込むのは難しそうです。

3) 項の方法で使う近似式には大きく2種類あるようです。高精度な近似式(17次くらいの式です）ともう少し精度が低い近似式(3次くらいの式）です。高精度な近似式を使えば秒オーダーの精度が出るようですが、常に最新の近似式を使う必要があり、マイコンに埋め込むことは難しいです。ちなみに海上保安庁の海洋情報部から近似式で使う定数が毎年発表されています。

　一方で、精度の低い近似式ではいつも同じ係数で計算するのでマイコンのプログラムに埋め込むことが出来ます。ただ、問題はどのくらいの精度で計算出来るか判らない点です。また、Arduino の浮動小数点の値は 32ビットで計算されるのでＰＣと比べて精度が低く、誤差が出やすくなります。ということで、このあたりは実際にプログラムを作って確かめるしか無さそうです。

　そんなことで、3) 項の精度の低い近似式を使う方法で計算プログラムを作りました。作成に当たっては下記のサイトを参考にさせて頂きました。ありがとうございます。
http://k-ichikawa.blog.enjoy.jp/etc/HP/js/sunRise/srs.html 　日の出・日の入時刻の年間変化（JavaScript版）
http://www.iot-kyoto.com/satoh/2016/01/22/post-99/ 　日の出時間と日没時間の計算方法～プログラム編

/*
  Arduinoで日の出・日の入り時刻を計算
  近似式を使って日の出・日の入り時刻を求める。
  by ラジオペンチ, 2017/5/10, http://radiopench.blog96.fc2.com/
  参考サイト
  http://k-ichikawa.blog.enjoy.jp/etc/HP/js/sunRise/srs.html
  http://www.iot-kyoto.com/satoh/2016/01/22/post-99/
*/

#include <math.h>  // 数学ライブラリ　https://www.arduino.cc/en/Math/H

#define  M_PI   3.14159265358979323846    // πの値
#define DEG(a)  ((a) * 180 / M_PI)        // ラジアンを度に変換するマクロ
#define RAD(a)  ((a) * M_PI /180)         // 度をラジアンに変換するマクロ

float t1;                                 // 日の入り時刻（単位：時）
float t2;                                 // 日の入り時刻（単位：時）
float tm;                                 // 南中時刻（単位：時）
float tL;                                 // 昼の長さ（単位：時）

float Longitude = 139.7414;               // 経度（東経を入力）
float Latitude = 35.6581;                 // 緯度（北緯を入力）

void setup() {
  Serial.begin(115200);
}

void loop() {
  for (int n = 0; n < 365; n++) {              // 日付連番 0から364日の範囲を計算
    t1 = SunRiseTime(Longitude, Latitude, n);  // 経緯度、日付連番から日の出時刻を求める
    t2 = SunSetTime(Longitude, Latitude, n);   // 経緯度、日付連番から日の入り出時刻を求める
    tm = (t1 + t2) / 2;                        // 太陽南中時刻
    tL = (t2 - t1);                            // 昼の時間
    Serial.print(n + 1), Serial.print(", ");
    Serial.print(t1, 4), Serial.print(", ");
    Serial.print(t2, 4), Serial.print(", ");
    Serial.print(tm, 4), Serial.print(", ");
    Serial.print(tL, 4), Serial.println();
  }
  for (;;) {}                                  // ここで停止
}

float SunRiseTime(float x, float y, int n) {   // 日の出時刻を求める関数
  float d, e, t;
  y = RAD(y);                                  // 緯度をラジアンに変換
  d = dCalc(n);                                // 太陽赤緯を求める
  e = eCalc(n);                                // 均時差を求める
  // 太陽の時角幅を求める（視半径、大気差などを補正 (-0.899度)　）
  t = DEG(acos( (sin(RAD(-0.899)) - sin(d) * sin(y)) / (cos(d) * cos(y)) ) );
  return ( -t + 180.0 - x + 135.0) / 15.0 - e; // 日の出時刻を返す
}

float SunSetTime(float x, float y, int n) {    // 日の入り時刻を求める関数
  float d, e, t;
  y = RAD(y);                                  // 緯度をラジアンに変換
  d = dCalc(n);                                // 太陽赤緯を求める
  e = eCalc(n);                                // 均時差を求める
  // 太陽の時角幅を求める（視半径、大気差などを補正 (-0.899度)　）
  t = DEG(acos( (sin(RAD(-0.899)) - sin(d) * sin(y)) / (cos(d) * cos(y)) ) );
  return ( t + 180.0 - x + 135.0) / 15.0 - e;  // 日の入り時刻を返す
}

float dCalc(int n) {                           // 近似式で太陽赤緯を求める
  float d, w;
  w = (n + 0.5) * 2 * M_PI / 365;              // 日付をラジアンに変換
  d = + 0.33281
      - 22.984 * cos(w) - 0.34990 * cos(2 * w) - 0.13980 * cos(3 * w)
      + 3.7872 * sin(w) + 0.03250 * sin(2 * w) + 0.07187 * sin(3 * w);
  return RAD(d);                               // 赤緯を返す（単位はラジアン）
}

float eCalc(int n) {                           // 近似式で均時差を求める
  float e, w;
  w = (n + 0.5) * 2 * M_PI / 365;              // 日付をラジアンに換算
  e = + 0.0072 * cos(w) - 0.0528 * cos(2 * w) - 0.0012 * cos(3 * w)
      - 0.1229 * sin(w) - 0.1565 * sin(2 * w) - 0.0041 * sin(3 * w);
  return e;                                    // 均一時差を返す（単位は時）
}

プログラムの解説
・Arduino IDE は1.8.1を使いました。
・10, 12行目　三角関数の acos を使うために念のために入れましたが、コメントアウトしても大丈夫でした。Arduino IDE に元々入ってるのかもしれません。
・21, 22行目　使う場所の経度、緯度を入れます。なお、デフォルトの値は暦の計算で使われる東京の座標です。
・30, 31行目　SunRiseTime と SunSetTime の関数で日の出と日の入り時刻を求めています。一つの関数でまとめて日の出/日の入り時刻を得るやり方もあり、その方が無駄な計算をしないで済みます。でも判り易さを優先してこのような作りにしてみました。なお、この関数一つの一回の実行時間は約 2.3ms でした(計算する値で多少変化します）。

▼実行結果（先頭部分）

　一年分の日の出時刻、日の入り時刻、太陽南中時刻、昼の長さを出力します。

▼グラフ

　Arduino IDE のシリアルプロッタでグラフにすると判り易いです。なお n をプロットすると見づらくなるので、34行目をコメントアウトしています。

　ということで日の出・日の入り時刻が計算出来るようになりました。気になるのはその精度です。

▼エクセルで計算精度確認


　計算精度確認のために、結果を天文台の暦のページから発表されている日の出、日の入り時刻と比較してみました。

　なお、比較に使ったデーターは、アナレンマを描く基本検討をやった時に入手した2016年のデーターです。この年は閏年なので366日あって比較対象としてはあまり良くないです。でも悪い条件で比較しておけばこれより悪くなることは無いでしょうからこのままいきます。（というのは言い訳で、データーを取り直すのが面倒です）

　長々と書いてきましたが、以下のグラフがこの記事の結論です。グラフの横軸は一年間を表し、縦軸が天文台の暦のページで発表されている日の出・日の入り時刻と、このプログラムの計算結果との差です。

▼日の出時刻の計算誤差


▼日の入り時刻の計算誤差

　ほとんどのデーターは±60秒以内に入っています。日の出・日の入りがこれくらいの精度で計算出来ていれば問題無いと思います。天文台から発表される日の出・日の入り時刻は1分単位なので、プラスマイナス30秒の丸め誤差が入っているのですが、その様子が見えています。別の言い方をすると、プロットされた点は60秒の帯状の範囲に分布していますが、この帯の中心線が正確な日の出・日の入り時刻です。

◆まとめ
　Arduinoを使った日の出・日の入り時刻の計算プログラムが出来ました。実際に何かのアプリで使うためには、日付と時刻の情報が必要になり、RTCなどと組み合わせて使うことになると思います。でもともかく一歩前進です。

　日の出・日の入り時刻を求める計算では、三角関数の計算を沢山やらないといけません。これは8ビットのマイコンには重い仕事なので、実行時間がかなり長くなるのではないかと心配しました。しかし、やってみるとミリセカンドオーダーで計算されたのには驚きました。
　ちなみに昔の8ビットのマイコン、例えば Z80 に BASICインタープリターで同じ計算をやらせたら、どれくらいの実行時間になるか興味深いところです。たぶんこの10倍の時間を使っても済まない気がします。