<!doctype html>
<head>
<title>eotsim</title>
<script>
const PI2 = 2 * Math.PI;
const prps = PI2 / 86400; // planet rotation per second
const epsilon = 1.0e-14;
function Calc()
{
var obliquity = parseFloat(document.getElementById("obliquity").value);
var eccentricity = parseFloat(document.getElementById("eccentricity").value);
DrawObliquity(obliquity);
DrawEccentricity(eccentricity);
DrawGraph(obliquity, eccentricity);
}
function ToRadians(a)
{
return a * Math.PI / 180;
}
function ModAngle(angle)
{
while (angle <= -Math.PI) angle += PI2;
while (Math.PI < angle) angle -= PI2;
return angle;
}
function GetSvg(id)
{
var svg = document.getElementById(id);
while (svg.firstChild) {
svg.removeChild(svg.firstChild);
}
return svg;
}
function Append(svg, type, attr)
{
var e = document.createElementNS("http://www.w3.org/2000/svg", type);
for (var n in attr) {
e.setAttribute(n, attr[n]);
}
svg.appendChild(e);
}
function DrawObliquity(obliquity)
{
var svg = GetSvg("svg_o");
var rotate = "rotate(" + obliquity + ",150,150)";
var attr;
// 惑星
attr = {cx:150, cy:150, r:80, fill:"#88f"};
Append(svg, "circle", attr);
// 地軸
attr = {x1:150, y1:30, x2:150, y2:270, stroke:"#000", transform:rotate};
Append(svg, "line", attr);
// 赤道
attr = {x1:70, y1:150, x2:230, y2:150, stroke:"#f00", transform:rotate};
Append(svg, "line", attr);
// 黄道
attr = {x1:70, y1:150, x2:230, y2:150, stroke:"#ff0"};
Append(svg, "line", attr);
}
function DrawEccentricity(e)
{
var svg = GetSvg("svg_e");
var attr;
var a = 100;
var b = Math.round(Math.sqrt(1 - e * e) * a);
// 楕円軌道
attr = {cx:150, cy:150, rx:a, ry:b, stroke:"#000", fill:"#fff"};
Append(svg, "ellipse", attr);
// 補助線
attr = {x1:150, y1:30, x2:150, y2:270, stroke:"#000", "stroke-dasharray":"5,5"};
Append(svg, "line", attr);
attr = {x1:30, y1:150, x2:270, y2:150, stroke:"#000", "stroke-dasharray":"5,5"};
Append(svg, "line", attr);
// 太陽
var ae = Math.round(a * e);
attr = {cx:150 + ae, cy:150, r:5, fill:"#f84"};
Append(svg, "circle", attr);
}
function DrawGraph(obliquity, e)
{
var svg = GetSvg("svg_g");
var attr;
var data = [];
var maxsec = 1;
for (var day = 0; day < 360; day++) {
// 平均近点角(概算)近日点から次の近日点までの角度
var Ma = PI2 * (day / 360);
var Ta = KeplersEquation(e, Ma); // 真近点角
// 黄径(概算)春分点から次の春分点までの角度
var eclipticLon = Ta;
var ellipseEffect = Math.round(ModAngle(Ma - Ta) / prps);
var obliquityEffect = CalcObliquityEffect(obliquity, eclipticLon);
var equationOfTime = ellipseEffect + obliquityEffect; // 均時差
maxsec = Math.max(maxsec, Math.abs(equationOfTime));
var rec = {};
rec.ellipseEffect = ellipseEffect;
rec.obliquityEffect = obliquityEffect;
rec.equationOfTime = equationOfTime;
data.push(rec);
}
var hour = Math.ceil(maxsec / 3600);
for (var h = -hour; h <= hour; h++) {
var y = 150 - Math.round(150 * h / hour);
attr = {x1:0, y1:y, x2:720, y2:y, stroke:"#888"};
Append(svg, "line", attr);
}
var points_e = [];
var points_o = [];
var points_t = [];
var a = 150 / (3600 * hour);
for (var i = 0; i < data.length; i++) {
console.log(i + ":" + data[i].ellipseEffect + "," + data[i].obliquityEffect);
var x = i * 2;
points_e.push(x + "," + (150 - Math.round(a * data[i].ellipseEffect)));
points_o.push(x + "," + (150 - Math.round(a * data[i].obliquityEffect)));
points_t.push(x + "," + (150 - Math.round(a * data[i].equationOfTime)));
}
attr = {points:points_e.join(" "), stroke:"#00f", fill:"none"};
Append(svg, "polyline", attr);
attr = {points:points_o.join(" "), stroke:"#f00", fill:"none"};
Append(svg, "polyline", attr);
attr = {points:points_t.join(" "), stroke:"#000", fill:"none"};
Append(svg, "polyline", attr);
}
// 傾斜効果の計算
function CalcObliquityEffect(obliquity, eclipticLon)
{
// 地球を基準とした太陽の公転
var x = Math.cos(eclipticLon);
var r = Math.sin(eclipticLon);
var y = Math.cos(ToRadians(obliquity)) * r;
var celestialEquator = Math.atan2(y, x); // 天の赤道上の角度
return Math.round(eclipticLon / prps - celestialEquator / prps);
}
// 漸化式によりケプラー方程式を解く
// M 平均近点角(mean anomaly)
// E 離心近点角(Eccentric anomaly)
// T 真近点角(true anomaly)
function KeplersEquation(e, M)
{
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;
}
var K = Math.sqrt((1 + e) / (1 - e));
var T = Math.atan(K * Math.tan(E / 2)) * 2;
return T;
}
</script>
</head>
<body>
<table>
<tr><td>黄道傾斜角(度)</td><td><input id="obliquity" value="23.44"></td></tr>
<tr><td>離心率</td><td><input id="eccentricity" value="0.01671"></td></tr>
</table>
<br>
<button onclick="Calc()">Calc</button><br>
<svg id="svg_o" width=300 height=300></svg>
<svg id="svg_e" width=300 height=300></svg><br>
<svg id="svg_g" width=720 height=300></svg><br>
</body>