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분 정도는 어긋남이 있지만 괜찮을 것 같습니다.