MATLABで日の出日の入り

照度計のグラフに日の出と南中と日の入りの線を描いているのだが、当初ハードコードしていたものをプログラムで算出すことにしてみた。

SunCalcというJavascriptのコードを移植してみた。

dayMs = 1000 * 60 * 60 * 24;
J1970 = 2440588;
J2000 = 2451545;

rad  = pi / 180;

toJulian = @(date) date / dayMs - 0.5 + J1970;
fromJulian = @(j) (j + 0.5 - J1970) * dayMs;
toDays = @(date) toJulian(date) - J2000;

e = rad * 23.4397;
declination = @(l, b) asin(sin(b) * cos(e) + cos(b) * sin(e) * sin(l));

J0 = 0.0009;

julianCycle = @(d, lw) round(d - J0 - lw / (2 * pi));
approxTransit = @(Ht, lw, n) J0 + (Ht + lw) / (2 * pi) + n;
hourAngle = @(h, phi, d) acos((sin(h) - sin(phi) * sin(d)) / (cos(phi) * cos(d)));
observerAngle = @(height) -2.076 * sqrt(height) / 60;

solarMeanAnomaly = @(d) rad * (357.5291 + 0.98560028 * d);
solarTransitJ = @(ds, M, L) J2000 + ds + 0.0053 * sin(M) - 0.0069 * sin(2 * L);
eclipticLongitude = @(M) M + (rad * (1.9148 * sin(M) + 0.02 * sin(2 * M) + 0.0003 * sin(3 * M))) + (rad * 102.9372) + pi;


getSetJ = @(h, lw, phi, dec, n, M, L) solarTransitJ(approxTransit(hourAngle(h, phi, dec), lw, n), M, L);

lat = 35.6581;
lng = 139.7414;
height = 0;
now = datetime('now', 'TimeZone', 'UTC');
date = posixtime(now) * 1000;

lw = rad * -lng;
phi = rad * lat;
dh = observerAngle(height);
%d = toDays(date);
d = juliandate(now) - J2000;
n = julianCycle(d, lw);
ds = approxTransit(0, lw, n);
M = solarMeanAnomaly(ds);
L = eclipticLongitude(M);
dec = declination(L, 0);

Jnoon = solarTransitJ(ds, M, L);

solarNoon = fromJulian(Jnoon);
nadir = fromJulian(Jnoon - 0.5);

nan = datetime(solarNoon/1000,'ConvertFrom','posixtime','TimeZone','Asia/Tokyo');

time0 = -0.833;

h0 = (time0 + dh) * rad;

Jset = getSetJ(h0, lw, phi, dec, n, M, L);
Jrise = Jnoon - (Jset - Jnoon);

srise = fromJulian(Jrise);
sset = fromJulian(Jset);

sr = datetime(srise/1000,'ConvertFrom','posixtime', 'TimeZone','Asia/Tokyo');
ss = datetime(sset/1000,'ConvertFrom','posixtime', 'TimeZone','Asia/Tokyo');

disp(sr);
disp(nan);
disp(ss);

functionはThingSpeakのMATLABでは使えないようだったので、無名関数で実装しました。

元のコードは薄明の計算もありますが、このコードは日の出日の入りだけを計算しています。

ユリウス日の計算はMATLABの関数(juliandate)を使ってみました。

緯度経度がThingSpeakのデータに入っている場合はそれを使うのが良いです。

実行してみるとこうなります。

   18-May-2020 04:35:18

   18-May-2020 11:38:52

   18-May-2020 18:42:26

国立天文台のデータと1分程度はずれがあるが大丈夫そうです。