|
uchinyan |
|
#25809への追加：
少し長くなりますが、一応、プログラムを書いておきましょう。 北極点に点 A、経度 0 度の子午線上に点 B、球面の任意の位置に点 C をとって、点 B、C を極座標で表します。なお、球の中心を O とします。 点 B は、経度 0 度上なので、緯度の自由度だけでこれを y とします。点 C は、経度を a、緯度を z にしています。 これら、a, y, z は自由に等確率で与えることができると仮定していますが、さすがに、これは大丈夫でしょう。 次に、球面上の幾何学の知識を使います。私は、今回始めてWebで勉強したので、ひょっとしたら間違っているかもしれません．．．^^; さて、球面三角形ABC を考えます。すると、∠A は a になることが分かります。 ただし、球面三角形は、測地線の短い方を取るとしているので、a がπを超える場合は、∠A = 2π - a とします。 さらに、∠B = b、∠C = c とおきます。 球面三角形の各辺は、その辺となる測地線の中心角で表すので、∠AOB = y、∠AOC = z ですが、さらに、∠BOC = x とします。 球面三角法の公式から、 cos(x) = cos(y)cos(z) + sin(x)sin(y)cos(a) cos(b)sin(x) = cos(y)sin(z) - sin(y)cos(z)cos(a) cos(c)sin(x) = cos(z)sin(y) - sin(z)cos(y)cos(a) がいえ、b, c が求まります。 球面三角形の面積 s は、やはり、球面上の幾何学から、s = a + b + c - π と求まります。 後は、この面積を足しあげて、球の表面積を 1 としているので 4πで割って、取ってきた三角形の個数で平均します。 なお、プログラムの限界による桁あふれなどで計算がエラーになることがあるので、それを防ぐための工夫、DO LOOP の辺り、をしています。 2005/09/16 8:00am 頃追加。 RANDOMIZEの位置を少し変えました。この方が、ランダム性が増すようです。 2005/09/16 11:30am 頃追加。 ほとんど不要ですが、安全性を増すために、さらに工夫しました。 とはいえ、計算結果に改善は見られません。まぁ、乱数ですから、当然ですね。むしろ、揺れが大きくなったかなぁ。 100000 : .124803944399878 200000 : .125444581427515 300000 : .125101556380867 400000 : .125008553196104 500000 : .124829641436551 600000 : .125150645416347 700000 : .125171487607043 800000 : .125133544380509 900000 : .124987211093014 1000000 : .125020609081723 しかしいずれにせよ、1/8 = 0.125 の辺りをフラフラしている感じはあります。 以下、プログラムです。 LET n = 1000000 FOR k = 0 TO n STEP 100000 RANDOMIZE LET sum = 0 FOR i = 1 TO k LET x = 0 LET s = 10 LET t = 10 DO WHILE (ABS(s) > 1) OR (ABS(t) > 1) DO WHILE SIN(x) = 0 LET a = 2 * PI * RND IF a > PI THEN LET a = 2 * PI - a END IF LET y = PI * RND LET z = PI * RND LET x = ACOS(COS(y) * COS(z) + SIN(y) * SIN(z) * COS(a)) LOOP LET s = (COS(y) * SIN(z) - SIN(y) * COS(z) * COS(a))/SIN(x) LET t = (COS(z) * SIN(y) - SIN(z) * COS(y) * COS(a))/SIN(x) LOOP LET b = ACOS(s) LET c = ACOS(t) LET s = a + b + c - PI LET sum = sum + s NEXT i IF k <> 0 THEN LET p = sum/(4 * PI) LET p = p/k PRINT k;":"; p END IF NEXT k END |
|
ネコの住む家
9月16日（金） 11:28:27
MAIL:uchi@sco.bekkoame.ne.jp 25813 |