google-site-verification: google3bd66dd162ef54c7.html

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

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

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

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

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

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

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

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

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

 そんなことで、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倍の時間を使っても済まない気がします。
関連記事

コメントの投稿

管理者にだけ表示を許可する

凄いです。

こんばんは。あっという間に作ってしまったんですね。
すごく難しくて内容は全然理解できませんが、と言うか理解しようとも思うことが出来ない、高度すぎて正直拒絶反応です。自分でもスケッチ焼いて一年の日の出日の入り時刻が表示されました。時刻が時間単位なので悩み、これをどうやってrtcからのデータで変換させるか、格闘中です。 今作っているタイマー装置が超進化したバージョンアップになりそうです。 本当に有難うございました。ありがたく使わさせて頂きます。

re:凄いです。

hiroさんおはようございます。

これ作る動機付けしてくれたのはhiroさんなので感謝してます。

これをどう使うかはゆっくり考えたいと思ってます。

計算式

力作を使わさせて頂いているのですが、いろいろ調べて作った経過日数の計算式に自信が持てません。
ラジオペンチさんのrtcのソースから変数を持ってきて作りました。
int yy, mm, dd, ww;
yy = (dtString[2] & 0xf) * 10 + (dtString[3] & 0xf) + 2000;
mm = (dtString[5] & 0xf) * 10 + (dtString[6] & 0xf);
dd = (dtString[8] & 0xf) * 10 + (dtString[9] & 0xf);

int yyy= yy-1;
int mmm =1+12;

if (mm < 3) {
yy= yy-1 ;
mm= mm+ 12;
}

long ima = (long(365.25 * yy )) + (long(yy / 400)) - (long(yy / 100)) + (long(30.59 * (mm-2))) + dd - 678912;
long kako = (long(365.25 * yyy )) + (long(yyy / 400)) - (long(yyy / 100)) + (long(30.59 * (mmm-2))) + 0 - 678912;
int keikabi = (ima -kako);

今のユリウス通日?を求めて、それを元旦から引く。
こんな感じで良いのでしょうか?ちょっと自信が無くて書き込まさせて頂きました。
もう無理だと諦めていた、日の出日の入りタイマーがラジオペンチさんのプログラムで実現出来る事となり、本当に助かりました。今は試験運用ですが、今後が楽しみです。

re:計算式

hiroさん、今日は。

早速ですが、任意の年月日の間の日数を求める場合は、コメントに書かれたようなやり方でリウス日で計算するのだと思います。

でも今回の関数では年の情報は不要で、1月1日=0とした日付の連番の値だけあれば良いので計算は大幅に簡単になります。
ちなみに、SunRiseTime(Longitude, Latitude, n); の n の値で日付の情報を送りますが、値の範囲は0~364です。

日付の連番を求める方法は、if文を使うとかいろいろありますが、以下のような関数を使うのが一番簡単だと思います。
-----------------------------------
int dayOfYear(int mm, int dd) { // 月日から年間の経過日数を求める
int daySum[] = {0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334}; //前月末の通算日
return daySum[mm - 1] + dd - 1; // 経過日数を返す(1月1日=0)
}
-----------------------------------

dayOfYear の関数を呼べば引数で指定した月日に対する1月1日からの日数を返します。例えば、d = dayOfYear(2, 1); と呼び出すと、d には31が入ります。

なお、この関数は閏年は考慮していないので2月29日は3月1日と同じ値になります。ちなみにSunRiseTime、SunSetTime も閏年を考慮していません。というかエラーチェックはほとんど行っていないのでご注意ください。

ありがとうございます。

了解しました。そうなんですね。難しく考えなくても月の日数は決まっているので配列で簡単に出来てしまうのですね。
閏年も誤差の範囲で特に問題ありませんし、この方法を使わさせて頂きます。

ラジオペンチさんのブログには、本当にお世話になり感謝しています。
arduinoのスケッチもだいぶ分かるようになり、自分のしたい事のヒント・テクニックが、「気温・気圧ロガー 」のスケッチに入っていてとても勉強になりました。
あのスケッチから、いろいろな物を作らさせていただきまた。と言っても私のは稚拙なものばかりですけど。
カレンダー
07 | 2017/08 | 09
- - 1 2 3 4 5
6 7 8 9 10 11 12
13 14 15 16 17 18 19
20 21 22 23 24 25 26
27 28 29 30 31 - -
プロフィール

ラジオペンチ

Author:ラジオペンチ
電子工作を中心としたブログです。たまに近所(東京都稲城市)の話題など。60過ぎて視力や器用さの衰えを感じつつ日々挑戦!
コメントを入れる時にメールアドレスの記入は不要です。なお、非公開コメントは受け付けていません。

記事が気に入ったらクリックを!
最新記事
カテゴリ
最新コメント
リンク
FC2カウンター
検索フォーム
月別アーカイブ
RSSリンクの表示
QRコード
QRコード