Babylon.js で天文学に基づいた太陽の位置を表現しようとしてダメだった
こんにちは。やまゆです。
太陽の位置は、天文学的にかなり小さい誤差で算出が可能らしいです。
観測地点の緯度・経度と、日時を入力とすると、 XYZ や azimuth(方位角, 北とか) と inclination(傾斜) で出力することが出来ます。
ちなみに、 mourner/suncalc: A tiny JavaScript library for calculating sun/moon positions and phases. というライブラリもありましたが、興味本位で数式を探してみます。
Position of the Sun - Wikipedia
位置の計算について Wikipedia にそれっぽい数式が載っていたので、これを再現してみようと思います。
地球から太陽を見た時は「傾斜(Declination)」を加味する必要があり、地球は地軸が 23.44 度傾いているので、それを考慮した計算にする必要があるようです。
ばーっと wikipedia を DeepL にかけつつ、計算式の部分を見てみます。
まず最初に、地球には四季が訪れるわけですが、それは太陽からの距離が公転によって若干変わるからだと思います。それを ecliptic longitude と呼ぶみたいです。その軌道の遠心率（正円じゃなく楕円になっている率）は誤差 1 度以下なのでほぼないとみなして、傾斜をこの式で算出できるようです。
function calculateDeclination(N: number): number {
return Math.asin(
Math.sin(-23.44)
* Math.cos(
360 / 365.24
* (N + 10)
+ 360 / Math.PI
* 0.0167
* Math.sin(
360 / 365.24
* (N - 2)
)
)
);
}
ここで N という数値は、 UTC で 1 月 1 日を 0 として、そこから経過した日数(小数点含む)を入れることが可能です。 1 月 1 日 12 時は 0.5 になるでしょう。
function calculateN(date: Date): number {
const d2 = new Date(date);
// 1 月 1 00:00:00 に設定する
d2.setMonth(0, 1);
d2.setHours(0, 0, 0, 0);
// ミリ秒時間が手に入るので「日」に単位変換
const newYear = d2.getTime() / 1000 / 60 / 60 / 24;
const today = date.getTime() / 1000 / 60 / 60 / 24;
const elapsedDate = today - newYear;
return elapsedDate; // Jan. 1 00:00:00 means 0.0
}
また、こちらのページの The formula based on the subsolar point and the atan2 function の部分で、太陽の傾斜と日時、緯度経度から太陽の方角を算出する数式が載っています。
まず、 UTC で 00:00:00 からの経過時間を取得します。
function calculateT(date: Date): number {
const d2 = new Date(date);
d2.setHours(0, 0, 0, 0);
// ミリ秒時間が手に入るので「時間」に単位変換
const startOfDay = d2.getTime() / 1000 / 60 / 60;
const now = date.getTime() / 1000 / 60 / 60;
const elapsedHours = now - startOfDay;
return elapsedHours;
}
また、均時差と呼ばれる時間を計算します。これも 英語版 wikipedia の数式を利用します。
function calculateE(date: Date): number {
const y = date.getUTCFullYear();
const d = Math.floor(calculateN(date));
const D = 6.24004077 + 0.01720197 * (365.25 * (y - 2000) + d);
const delta = -7.659 * Math.sin(D) + 9.863 * Math.sin(2 * D + 3.5932);
return delta;
}
function calculateSunPosition(date: Date, latitude: number, longitude: number): { x: number, y: number, z: number} {
const declinationOfTheSun = calculateDeclination(calculateN(date)); // degree
const TGMT = calculateT(date); // UTC decimal hours since 00:00:00
const Emin = calculateE(date);
const latitudeOfTheSubsolarPoint = declinationOfTheSun;
const longitudeOfTheSubsolarPoint = -15 * (TGMT - 12 + Emin / 60);
const x = Math.cos(latitudeOfTheSubsolarPoint)
* Math.sin(longitudeOfTheSubsolarPoint - longitude);
const y = Math.cos(latitude)
* Math.sin(latitudeOfTheSubsolarPoint)
- Math.sin(latitude)
* Math.cos(latitudeOfTheSubsolarPoint)
* Math.cos(longitudeOfTheSubsolarPoint - longitude);
const z = Math.sin(latitude)
* Math.sin(latitudeOfTheSubsolarPoint)
+ Math.cos(latitude)
* Math.cos(latitudeOfTheSubsolarPoint)
* Math.cos(longitudeOfTheSubsolarPoint - longitude);
return { x, y, z };
}
よし、多分合ってる、と思いつつ Playground に放り込んでみましたが...
時間をお昼にしたのに真っ暗。
恐らく単位の変換が足りてないんだと思います。誰か間違っている所教えて（断念）
Playground はこちら↓↓
Discussion