開発環境 メモ帳
実行環境 Internet Explorer 11


eotsim.html
<!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>
 
最終更新:2014年09月06日 22:23