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分程度はずれがあるが大丈夫そうです。
Author And Source
この問題について(MATLABで日の出日の入り), 我々は、より多くの情報をここで見つけました https://qiita.com/yamori813/items/99ce6a8f9b8a0e9be9e8著者帰属:元の著者の情報は、元のURLに含まれています。著作権は原作者に属する。
Content is automatically searched and collected through network algorithms . If there is a violation . Please contact us . We will adjust (correct author information ,or delete content ) as soon as possible .