天文計算に最適?なプログラミング言語 GNU Octave

GNU Octaveは科学技術計算用インタープリタプログラミング言語です。Octaveは基本的にPythonのようにShellをもち、コマンドラインから実行します。また、エディタでプログラムファイルを作成し、それをコマンドラインから実行させることもできます。GUIデスクトップもできるようですが、あまりにも情報が少なく、そこまではわかりません。

Octaveのインストール

Windowsであれば、Octaveは公式のサイトから、インストーラーをダウンロードして、インストールすることができます。ですが、Linuxだと、ディストリビューションがサポートしている、ソフトウェアの管理やソフトウェアの追加と削除を利用してインストールすることにします。

ここでは、Linux Mint 21をのソフトウェアの管理からインストールします。ソフトウェアの管理には、2つのバージョンが用意されています。ひとつは通常版とFlatPak版です。どちらでも構いませんが、ここでは通常版をインストールすることにします。ただし、最新版は7.2ですが、登録されている通常版は6.4といささか古かったので、PPAを利用して最新版をインストールすることにします。

sudo add-apt-repository ppa:ubuntuhandbook1/octave
sudo apt update
sudo apt install octave

なお、警告が表示される場合は、下記の記事を参考にしてください。

Linux Mint 21でPPA登録時に警告が出る問題の修正

インストールが終わると、メニューの「教育・教養」に登録されています。

実行

実行すると、最初にいくつかのダイアログが表示されます。[次]をクリックしていきます。

Octave Welcomeメッセージ
Octave起動画面

コマンドウィンドウからプログラムを入力してもよいのですが、1行ずつ入力しては大変なので、下部の「エディタ」タブをクリックして、プログラムを作成していきます。例として、「マイコン宇宙講座」のサブルーチンMJD、JULIAN、JDATEを作成して、同じ結果になるか確かめてみることにします。

サブルーチン library.m

GNU Fortranで書いたコードをOctaveで書き直したものです。1行目の1;は必ず入力してください。Octaveでは、外部ファイルに記述した関数は、1つのファイルにつき1関数ということになっています。複数の関数を呼び出すには、一度、スクリプトとして実行してから、呼び出すことになります。

また、行末の;(セミコロン)を付け忘れると、実行時に変数の内容が表示されてしまいます。表示させたくないときはセミコロンを忘れないようにしてください。

1;

format long
% ModifiedJulian
function [julian, year, month, day, hour, minute, second] = modifiedJulian(dy, dt)
    year = sign(dy) * floor(abs(dy) / 10000.0);
    month = mod(floor(abs(dy) / 100.0), 100);
    day = mod(abs(dy), 100.0);
    hour = floor(dt / 10000.0);
    minute = mod(floor(dt / 100.0), 100);
    second = mod(dt, 100.0);

    julian = julianDay(year, month, day);
    julian = julian + (hour / 24.0 + minute / 1440.0 + second / 86400.0);
    julian = julian - 0.375;
end

% Julian Day
function julian = julianDay(year, month, day)
    branch = year + (month - 1.0) / 12.0 + day / 365.25;

    if (floor(month) < 3)
        month = month + 12.0;
        year = year - 1.0;
    endif

    if (branch >= 1582.78)
        julian = floor(year * 365.25) + floor(year / 400.0) - floor(year / 100.0) + floor(30.59 * (month - 2.0)) + day + 1721088.5;
    elseif (branch >= 0.0)
        julian = floor(year * 365.25) + floor(30.59 * (month - 2.0)) + day + 1721086.5;
    elseif (year < 0.0)
        julian = sign(year) * floor(abs(year) * 365.25) + floor(30.59 * (month - 2.0)) + day + 1721085.5;
    end
end

% JulianToDate
function [year, month, day] = julianToDate(julian)
    T = [31 28 31 30 31 30 31 31 30 31 30 31];

    jj = julian - 2400000.5;
    year = floor(0.0027379093 * (julian - 2400000.5) + 1858.877);
    month = 1.0;
    day = 0.0;

    julian = julianDay(year, month, day);

    r2 = jj - (julian - 2400000.5);
    if (mod(year, 4.0) == 0 && mod(year, 100.0) ~= 0 || mod(year, 400.0) == 0)
        T(1, 2) = 29
    endif

    r1 = 0.0;
    m = 1;
    while (m < 13)
        if (floor(r2) - r1 - T(1, m) <= 0)
            break;
        endif
        r1 = r1 + T(1, m);
        m = m + 1;
    endwhile

    month = m;
    day = r2 - r1;
    T(1, 2) = 28;

    if (month == 13)
        year = year + 1.0;
        month = month - 12.0;
    endif
end

メインルーチン main.m

% マイコン宇宙講座
% メインJST <-> MJD変換プログラム
% リスト 2-3

printf("\n");
dtime = input("DATE AND TIME (JST)? ");
dy = dtime(1);
dt = dtime(2);

[julian, year, month, day, hour, minute, second] = modifiedJulian(dy, dt);
dtime = hour / 24.0 + minute / 1440.0 + second / 86400.0;

printf("\n");
printf("   INPUT DATE = %4d %2d %7.5f UT\n", year, month, day + dtime - 0.375);
printf("         MJD  = %13.5f UT\n", julian);
printf("\n");

[year, month, day] = julianToDate(julian);
printf("  RETURN DATE = %4d %2d %7.5f UT\n", year, month, day);

jilian = julianDay(year, month, day);
printf("         MJD  = %13.5f UT\n", julian);
printf("\n");

実行

実行は左側のパネルから保存したフォルダをダブルクリックしてから、libray、mainの順に実行します。

library
main

結果

日付と時刻を入力する際は、Octaveの仕様上、[]で必ず囲んでください。ひとであれば、[]で囲む必要はありませんが、複数の場合は[]で囲む必要があります。

DATE AND TIME (JST)? [20221001 120000]

   INPUT DATE = 2022 10 1.12500 UT
         MJD  = 2459853.62500 UT

  RETURN DATE = 2022 10 1.12500 UT
         MJD  = 2459853.62500 UT