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倍の時間を使っても済まない気がします。

ピンチェンジ割込みを使ってロータリーエンコーダーを読む (Arduino)

 少し前の記事でウォッチドッグタイマー割込みを使ってロータリーエンコーダーの読取りをやってみました。これは、ウォッチドッグタイマーの使い方を試す事例として悪くないと思います。でも、割り込み周期を15msまでしか短く出来ないので、ロータリーエンコーダーをすばやく廻した時に取りこぼしが出てしまうという問題がありました。

 タイマー割込みを使ってロータリーエンコーダーの読取りを行えばこんなことにはならないのですが、タイマーの割り込みが余っていない場合はそうもいきません。

 ということで困っていたのですが、ピンチェンジ割り込み (Pin Change Interrupt) を使えば何とかなりそうなことに気付きました。

 ピンチェンジ割り込みはCPUの機能として存在していますが、Arduinoの言語仕様には出てこないので自分で制御レジスタを設定してやる必要があります。なお、これはArduinoの言語仕様に出てくる外部割込み(INT0, INT1)とは別の割込みになります。

 CPU(ATmega328P)のデーターシートを読むと、ピンチェンジ割込みの詳しい使い方が書かれていますが、物理ピンや割込み信号名の関係がややこしいので一覧表に整理してみました。

▼ピンチェンジ割込み関係ピン、信号対応表
ピンチェンジ割込み関係信号対応表
 今回読もうとしているロータリーエンコーダーは、この表の黄色で示すD9とD10に接続されています。その割込み信号名は PCINT1と PCINT2 で、PCMSK0 レジスタで指定するということが判ります。

 ピンチェンジ割り込みでは、まとまったグループ単位でどこかのピンの状態に変化があった時に割り込みが掛かります。実際にどのピンが変化したのかを特定するためには、プログラムで調べる必要があります。でもロータリーエンコーダーの読取りを行う場合、「二つのピンのどちらかに状態の変化があった」ということが判れば大丈夫なので処理はかなり簡単になります。

 ということで、実際に動作確認するためのプログラムを書いてみました。

▼ピンチェンジ割り込みでロータリーエンコーダーを読むスケッチ
/* ピンチェンジ割込みを使ってロータリーエンコーダーを読む
* A相 = D9, PB1, PCINT1, PCIE0
* B相 = D10, PB2, PCINT2, PCIE0
* 2016/4/8 ラジオペンチ http://radiopench.blog96.fc2.com/
*/
#define ECA 9 // エンコーダーA相 = D9
#define ECB 10 //       B相 = D10

volatile int X, nIRQ;
volatile unsigned long lastMicros;
int data;

void setup() {
pinMode(ECA, INPUT_PULLUP);
pinMode(ECB, INPUT_PULLUP);
Serial.begin(9600);
Serial.println("start");

PCMSK0 |= ((1 << PCINT1) | (1 << PCINT2)); // D9,D10ピンからのピンチェンジ割込みを使う
PCICR |= (1 << PCIE0); // PCIE0グループからの割込み許可
}

void loop() {
if (X != 0) { // エンコーダーの値が変化していたら
data += X;
X = 0;
Serial.print(data); Serial.print(", "); Serial.print(nIRQ); // 動作確認表示
Serial.print(", "); Serial.println(lastMicros);
}
}

ISR(PCINT0_vect) { // PCIINT0グループからの割込み処理
static byte bp = 0; // ビットパターン記録バッファ
static int n = 0;

n++; // 割込み発生回数カウント
bp = bp << 1; // ピンの状態変化をbpに右詰めで記録 0b00ABABAB
if (digitalRead(ECA) == HIGH) bp |= 0x01;
bp = bp << 1;
if (digitalRead(ECB) == HIGH) bp |= 0x01;

bp &= 0x0F; // 下位4ビット残して上位を消す 0b0000ABAB
if (bp == 0b0111) { // このパターンに
X++; // 一致していればインクリメント
nIRQ = n; // 割込み回数を記録し
n = 0; // カウンタをリセット
}

if (bp == 0b1011) {
X--; // 一致していればデクリメント
nIRQ = n;
n = 0;
}

if (n == 0) { // 値の変化があって
if ((micros() - lastMicros) < 50000) { // 50ms以内の更新だったら
X *= 10; // 早送り(10倍速)
}
lastMicros = micros();
}
}
 19、20行目がピンチェンジ割込みを起動するためのレジスタ設定です。ここで設定しておけば、割り込みが掛かるたびに割込み処理ルーチンである32行目の ISR(PCINT0_vect) が実行されます。

 割込み処理ルーチンではパターン一致方式でロータリーエンコーダーの読取りを行っていますが、ちょっとひねりを加えて、素早く廻した時は早送りでカウントが進むようなルーチンを追加してあります。

 実行するとシリアルモニタに状態を表示します。

▼出力例
実行結果
 エンコーダーのカウント値、割込み発生回数、経過時間(μs)の順に表示します。10カウントの時点ですばやく廻したのでカウントが10ステップ間隔で増えています。

 割込み回数は理論的には4になるはずですが、チャッタリングの影響で多少変化しています。なお、スイッチのチャッタリング消しをプログラムで行いたかったのですが、うまい手を思い付かなかったので、回路にチャッタ消しのコンデンサを入れています(下図のC1,C2)。コンデンサの容量は割込み回数を見ながら調整すればいいのですが、今回はそこまではやっていません。

▼回路図
回路図

 ということで、ピンチェンジ割り込みを使ったロータリーエンコーダー読み取りはうまくいったようです。このテクニックはロータリーエンコーダーに限らずArduinoでより複雑な物を作りたい場合に役立つと思います。

 この記事の作成には以下のサイトを参考にさせていただきました。ありがとうございます。
・FIRMLOGICSさん、Arduino のスリープと、Pin Change 割込
・wsnakさん、LEDシーリングライト用赤外線リモコンの製作(9)アイリスオーヤマ用 スリープ対応
・橋本商会さん、ATmega168でピン変化割り込み

ウォッチドッグタイマー割込みを使ってロータリーエンコーダーを読む

 先にお断りしておきますが、この記事の内容ではArduinoを使ってロータリーエンコーダーを完璧に読むことは出来ません。この記事は、ウォッチドッグタイマー割り込みの使い方を解説するために書いています。

 Arduino UNO(ATmega328P)には timer0,timer1,timer2 の3つのタイマーがあります。このうちtimer0はシステムが使っているので、勝手に使うことはできません。いや、使ってもいいのですが、それをやるといろいろ面倒なことが起きるはずです。
 ということで、普通はTimer1とTimer2の割込みを使うことになります。普通のアプリならタイマー割り込みが二つあれば大丈夫でしょう。でももうすこし複雑なことをやらせたい場合、もう一つタイマー割り込みがあると助かります。

 そんなことを考えていたら、居酒屋ガレージ店主さんからコメントで、ATtiny2313を使ったパルスジェネレーターを紹介していただきました。この作品はアセンブラを駆使してCPUの資源をしゃぶりつくして作られているのですが、特に参考になったのは、ウォッチドッグタイマー割込みを使ってデジスイッチの読み取りをやっていることです。

 そうか、その手があったか。時間精度のいらない低頻度の割込みならウォッチドッグタイマーが使えます。このテクニックは組み込み系をやられている人には常識なのかも知れませんが、私にとっては目からウロコでした。

 ということで、早速Arduino UNOでやってみました。

▼回路図
Arduinoにロータリーエンコ-ダー
 Arduino UNOにロータリーエンコーダーを接続します。

▼スケッチ
/* ウォッチドッグタイマー割り込みで、ロータリーエンコーダ
* を読取る実験。(15ms間隔サンプリングなので早回しはダメ)
* 2016/03/31 ラジオペンチ http://radiopench.blog96.fc2.com/
*/

//#include <avr/wdt.h> // ウォッチドッグタイマーを使用

#define ECA 9 // エンコーダーA相 = D9
#define ECB 10 //       B相 = D10

volatile int X;
int data;

void setup() {
pinMode(ECA, INPUT_PULLUP); // ロータリーエンコーダーA
pinMode(ECB, INPUT_PULLUP); // ロータリーエンコーダーB

Serial.begin(9600);
Serial.println("start");
WDTimerStart(0); // 15ms間隔でウォッチドッグタイマー割込み開始
}

void loop() { // メインループ
if (X != 0) { // エンコーダーが動いていたら
data = data + X;
X = 0;
Serial.println(data);
}
}

void WDTimerStart(unsigned int ii) { // ウォッチドッグタイマーをセット。
// 引数はWDTCSRにセットするWDP0-WDP3の値。設定値と動作時間は概略下記
// 0=16ms, 1=32ms, 2=64ms, 3=128ms, 4=250ms, 5=500ms
// 6=1sec, 7=2sec, 8=4sec, 9=8sec
byte bb;
if (ii > 9 ) ii = 9; // 変な値を排除
bb = ii & 7; // 下位3ビットをbbに
if (ii > 7) { // 7以上(8,9)なら
bb |= (1 << 5); // bbの5ビット目(WDP3)を1にする
}
bb |= ( 1 << WDCE );

MCUSR &= ~(1 << WDRF); // MCU Status RegのWatchdog Reset Flag ->0
// start timed sequence
WDTCSR |= (1 << WDCE) | (1 << WDE); // ウォッチドッグ変更許可(WDCEは4サイクルで自動リセット)
// set new watchdog timeout value
WDTCSR = bb; // 制御レジスタを設定
WDTCSR |= _BV(WDIE);
}

ISR(WDT_vect) { // WDT割込み発生時にロータリーエンコーダーを読む
static byte bp = 0; // ビットパターン記録バッファ
bp = bp << 1; // ピンの状態変化をbpに右詰めで記録 0b00ABABAB
if (digitalRead(ECA) == HIGH) bp |= 0x01;
bp = bp << 1;
if (digitalRead(ECB) == HIGH) bp |= 0x01;

bp &= 0x0F; // 下位4ビット残して上位を消す 0b0000ABAB
if (bp == 0b0111) X ++; // 一致していればインクリメント
if (bp == 0b1011) X --; // 一致していればデクリメント
}
 ロータリーエンコーダーを廻すと現在の値をシリアルモニターに表示するプログラムになっています。ロータリーエンコーダーの読取りは、51行目以降の割込み処理ルーチンで行っていますが、アルゴリズムは、以前の記事に書いたビットパターン一致の監視法です。CPUのレジスタを直接いじっているので、6行目の #include は使わないのでコメントアウトしています。

 20行目の WDTimerStart(0); でウォッチドッグタイマー割り込みを起動しています。引数によって割込み間隔が変わり、0で呼ぶと最高速である約15ms間隔で割り込みが発生します。ちなみにこの関数は以前記事にしたdelayWDT関数のプログラムで使った、delayWDT_setup(unsigned int ii) 関数の使い回しです。この割り込みを停止させるために、WDTimerStop()という関数も作った方がいいのですが、とりあえず止めないでも困らないのでまだ書いていません。

 これを動かした結果ですが、割り込みの速度が15ms間隔と遅いので、ロータリーエンコーダーをゆっくり廻さないと取りこぼしが発生します。まあこうなるのは想定の範囲内で、ゆっくりと動かすと正常に値が変化するので割り込みはうまくいっているようです。ちなみに、ロータリーエンコーダーを取りこぼし無く読むには、割込み間隔を1msくらいにする必要があります。

 そういう問題はありますが、ともかくウォッチドックタイマー割込みを使うプログラムはこれでいけそうです。人間が操作するスイッチの読み取りや、バッテリー電圧の監視くらいなら、ウォッチドックタイマの割込みで充分対応出来ると思います。

 あと、同じような機能はネットを探せばライブラリが出てくるので、それを使う手もあります。ただここからは私の意見ですが、Arduinoのユーザーライブラリには古いままメンテナンスされてないものがあって、最新版のIDEではコンパイラが通らないものが結構あります。ということで、よく確認しないでライブラリを使うと、かえってややこしいことになると思います。

◆まとめ
 ウォッチドッグタイマーは、プログラムの暴走など予期しない動きから確実に脱出するために用意されている仕掛けだと思います。でも、さほど信頼性が必要ない場合はこれを割込み源として使うのもアリだと思います。
カレンダー
05 | 2017/06 | 07
- - - - 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 -
プロフィール

ラジオペンチ

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

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