Octaveで天文計算する際に便利な関数

Octaveを使って天文計算をする際に、便利な関数を紹介します。天文計算と言うと、sinやcos、tanなどの数学関数を利用することがあります。Octaveでは、必要な数学関数や円周率などが用意されています。GNU Fortranでは円周率は、自分で定数として記述していましたが、Octaveでは必要ありませんし、ラジアン値に変換するための定数を記述することもありません。すでに用意されているからです。ここでは、そういった便利な関数を紹介していきましょう。

ラジアン・角度の相互変換

２つの関数が用意されています。

deg2rad … 角度をラジアンに換算する
rad2deg … ラジアンを角度に変換する

deg2rad(316.166396) → 5.518144594364230
rad2deg(5.518144594364230) → 316.1663960000000

sin・cos・tanなどの三角関数の引数を角度で指定

三角関数では、引数はラジアンを指定しなければなりません。そのため、あらかじめ角度をラジアンに換算する必要がありました。Octaveでは、三角関数名にdという修飾子をつけることで、角度で指定することができます。

sind(93.133688) → 0.998504701197285
sin(deg2rad(93.133688)) → 0.998504701197285

sind以外にもcosd tand asind acosd atandなどがありますが、実際に使用するなら、sindやcosdくらいでしょう。電卓と同じ感覚でできるのはいいですね。

自作の関数を作る

関数がなければ作っておくと、あとからコードを記述するのに楽になります。いくつか作ったものを紹介します。

hms2deg

時分秒を角度に換算します。

% hms2deg
function deg = hms2deg(h, m, s)
    deg = 15.0 * (h + m / 60.0 + s / 3600.0);
end

dms2deg

度分秒を角度に換算します。度分秒の場合、マイナスもありますから、sign関数で正負符号を判断しています。計算する際は、abs関数で絶対値にしておかないと、正しく計算されません。

% dms2deg
function deg = dms2deg(d, m, s)
    deg = sign(d) * (abs(d) + m / 60.0 + s / 3600.0);
end

deg2hms

角度から時分秒に換算します。
呼び出し： [h m s] = deg2hms(deg)
戻り値： []で囲まれた変数に値がセットされます。

%deg2hms
function [h m s] = deg2hms(deg)
    deg = deg / 15.0;
    h = floor(deg);
    s = 60.0 * (deg - floor(deg));
    m = floor(s);
    s = 60.0 * (s - m);
end

deg2dms

角度を度分秒に換算します。
呼び出し： [ｄd m s] = deg2dms(deg)
戻り値： []で囲まれた変数に値がセットされます。

%deg2dms
function [d m s] = deg2dms(deg)
    d = floor(abs(deg));
    d = sign(deg) * d;
    s = 60.0 * (abs(deg) - floor(abs(deg)));
    m = floor(s);
    s = 60.0 * (s - m);
end

hms2day

時刻の時分秒を日の少数点に換算します。

% hms2day
function t = hms2day(h, m, s)
    t = h / 24.0 + m / 1440.0 + s / 86400.0;
end

day2hms

日の小数点部分を時刻の時分秒に換算します。秒はround関数で小数点部分を四捨五入しています。
呼び出し： [h m s] = day2hms(deg)
戻り値： []で囲まれた変数に値がセットされます。

%day2hms
function [h m s] = day2hms(t)
    t = 24.0 * (t - floor(t));
    h = floor(t);
    s = 60.0 * (t - floor(t));
    m = floor(s);
    s = round(60.0 * (s - m));
end

引数に対する正負符号については、必要以外は行っていませんので、ご注意ください。

また、時分秒それぞれが範囲を超えた場合、たとえば、秒が60になったといったときは、秒を０にして分を＋１するなどの処理が必要になります。

57