<!doctype html>
<head>
<title>EquationOfTime</title>
<script>
const ns = "http://www.w3.org/2000/svg";
const PI2 = 2 * Math.PI;
const solarYear = 365.24219; // 太陽年(day) 365d5h48m45s
const anomalisticYear = 365.259643; // 近点年(day) 365d6h13m53.1552s
const e = 0.01671022; // 離心率(Orbital eccentricity)
const earthRotationPerSec = PI2 / 86400; // 地球が1秒間に回転する角度
const K = Math.sqrt((1 + e) / (1 - e));
const epsilon = 1.0e-14;
// 太陽の視角(0.533deg)と大気差(35m8s)による昼の長さの補正
const revision = Math.sin(Radians(0.533 / 2 + (35 * 60 + 8) / 3600));
var ST; // 標準時
var equationOfTimeRev; // 南中時補正
var latRad;
var ra;
var rx;
var svg;
var data;
function Calc()
{
var lat = parseFloat(document.getElementById("lat").value);
var lon = parseFloat(document.getElementById("lon").value);
ST = parseFloat(document.getElementById("ST").value);
var year = parseInt(document.getElementById("year").value);
equationOfTimeRev = 240 * lon - 3600 * ST;
latRad = Radians(lat);
ra = Math.sin(latRad) * revision;
rx = Math.cos(latRad) * revision;
data = [];
var date = new Date(year, 0, 1);
for (var i = 0; date.getFullYear() == year; i++) {
CalcDay(date);
date.setDate(date.getDate() + 1);
}
DrawGraph();
DrawTable();
}
function CalcDay(date)
{
// 修正ユリウス日
var y = date.getFullYear();
var m = date.getMonth() + 1;
var d = date.getDate();
if (m < 3) {
y--;
m += 12;
}
var MJD = Math.floor(365.25 * y) + Math.floor(y / 400) - Math.floor(y / 100) +
Math.floor(30.59 * (m - 2)) + d - 678912;
MJD += (12 - ST) / 24; // 標準時における12:00
// 黄道傾斜角
var T = (MJD - 51554.5) / 36525; // 2000/1/1 12:00(UT)からのユリウス世紀(36525日)
var obliquity = (84381.406 - 46.836769 * T - 0.00059 * T * T + 0.001813 * T * T * T) / 3600;
// 平均近点角（概算）近日点から次の近日点までの角度
var Ma = ModAngle(PI2 * ((MJD / anomalisticYear - 0.1242853) % 1.0));
var out = {};
KeplersEquation(Ma, out);
var Ta = out.T; // 真近点角
// 春分点(vernal equinox)の真近点角
var MJDv = (Math.floor(MJD / solarYear - 0.3399541) + 0.3399541) * solarYear;
var Mv = PI2 * ((MJDv / anomalisticYear - 0.1242853) % 1.0);
KeplersEquation(Mv, out);
var Tv = out.T;
// 黄径（概算）春分点から次の春分点までの角度
var eclipticLon = ModAngle(Ta - Tv);
// 楕円効果と傾斜効果
var ellipseEffect = Math.round(Ma / earthRotationPerSec - Ta / earthRotationPerSec);
var obliquityEffect = CalcObliquityEffect(obliquity, eclipticLon);
var equationOfTime = ellipseEffect + obliquityEffect; // 均時差
var transit = 43200 - equationOfTime - equationOfTimeRev; // 南中時
// 天球上の太陽軌道の高さと半径
var solarDecl = Radians(Math.sin(eclipticLon) * obliquity); // 太陽の赤緯
var solarAlt = Math.sin(solarDecl);
var solarRad = Math.cos(solarDecl);
// 天球上の太陽軌道と地平面の交点＝日出・日没
var x = ((solarAlt + ra) * -Math.tan(latRad) - rx) / solarRad;
var halfDaytime = Math.round(43200 * Math.acos(x) / Math.PI);
var rec = {};
rec.date = new Date(date.getTime());
rec.rising = transit - halfDaytime;
rec.transit = transit;
rec.setting = transit + halfDaytime;
data.push(rec);
}
function Radians(x)
{
return x * Math.PI / 180;
}
function ModAngle(angle)
{
while (angle <= -Math.PI) angle += PI2;
while (Math.PI < angle) angle -= PI2;
return angle;
}
// 傾斜効果の計算
function CalcObliquityEffect(obliquity, eclipticLon)
{
// 地球を基準とした太陽の公転
var x = Math.cos(eclipticLon);
var r = Math.sin(eclipticLon);
var y = Math.cos(Radians(obliquity)) * r;
var celestialEquator = Math.atan2(y, x); // 天の赤道上の角度
return Math.round(eclipticLon / earthRotationPerSec - celestialEquator / earthRotationPerSec);
}
// 漸化式によりケプラー方程式を解く
// M 平均近点角(mean anomaly)
// E 離心近点角(Eccentric anomaly)
// T 真近点角(true anomaly)
function KeplersEquation(M, out)
{
var E0 = M; // 初項
var E;
for (var i = 0; ; ) {
E = M + e * Math.sin(E0);
if ((E0 - epsilon < E) && (E < E0 + epsilon)) {
break;
}
if (10 <= ++i) {
// 計算打ち切り
break;
}
E0 = E;
}
out.E = E;
out.T = Math.atan(K * Math.tan(E / 2)) * 2;
}
//------------------------------------------------------------------------------
function DrawGraph()
{
svg = document.getElementById("svg");
while (svg.firstChild) {
svg.removeChild(svg.firstChild);
}
for (var i = 0; i < data.length; i++) {
var d = data[i];
var color = (d.date.getMonth() & 1) ? "#888" : "#ccc";
DrawRect(i, 0, d.rising, color);
DrawRect(i, d.rising, d.transit, "#8fc");
DrawRect(i, d.transit, d.setting, "#fc8");
DrawRect(i, d.setting, 86400, color);
}
for (var i = 1; i < 24; i++) {
var x = i * 30;
var l = document.createElementNS(ns, "line");
l.setAttribute("x1", x);
l.setAttribute("y1", 0);
l.setAttribute("x2", x);
l.setAttribute("y2", 365);
l.setAttribute("stroke", (i % 6) ? "white" : "red");
svg.appendChild(l);
}
}
function DrawRect(y, x1, x2, color)
{
var r = document.createElementNS(ns, "rect");
r.setAttribute("x", Math.floor(x1 / 120));
r.setAttribute("y", y);
r.setAttribute("width", Math.floor((x2 - x1) / 120));
r.setAttribute("height", 1);
r.setAttribute("stroke", color);
svg.appendChild(r);
}
//------------------------------------------------------------------------------
function DrawTable()
{
var div = document.getElementById("div");
while (div.firstChild) {
div.removeChild(div.firstChild);
}
var str = Math.floor(equationOfTimeRev / 60) + "m" + (equationOfTimeRev % 60) + "s";
var elm;
elm = document.createTextNode("南中時補正:" + str);
div.appendChild(elm);
elm = document.createElement("br");
div.appendChild(elm);
var table = document.createElement("table");
table.border = 1;
div.appendChild(table);
var thead = table.createTHead();
var row = thead.insertRow();
row.bgColor = "KHAKI";
AppendCell(row, "日付");
AppendCell(row, "日出");
AppendCell(row, "南中");
AppendCell(row, "日没");
var tbody = table.createTBody();
for (var i = 0; i < data.length; i++) {
var d = data[i];
var row = tbody.insertRow();
row.bgColor = (i & 1) ? "WHEAT" : "WHITE";
AppendCell(row, (d.date.getMonth() + 1) + "/" + d.date.getDate());
AppendCell(row, HourMin(d.rising));
AppendCell(row, HourMin(d.transit));
AppendCell(row, HourMin(d.setting));
}
}
function AppendCell(row, str)
{
var cell = row.insertCell();
var text = document.createTextNode(str);
cell.appendChild(text);
return cell;
}
function HourMin(sec)
{
sec += 30;
var hour = Math.floor(sec / 3600);
var min = Math.floor(sec / 60) % 60;
return hour + ":" + ("0" + min).slice(-2);
}
</script>
</head>
<body>
<table>
<tr><td>緯度</td><td><input id="lat" value="35.45"></td></tr>
<tr><td>経度</td><td><input id="lon" value="139.65"></td></tr>
<tr><td>標準時</td><td><input id="ST" value="9"></td></tr>
<tr><td>年</td><td><input id="year" value="2014"></td></tr>
</table>
<br>
<a href="http://eco.mtk.nao.ac.jp/koyomi/dni/">各地のこよみ(日の出入り、月の出入り、南中時)</a><br>
<br>
<button onclick="Calc()">Calc</button><br>
<br>
<svg id="svg" height=366></svg><br>
<br>
<div id="div"></div>
</body>
