Project Lovecraft Problem 3

[danbaron]

Problem 3

Numerically estimate (or calculate) the volume of the part of Sphere 0 which is not
intersected by any of Spheres 1-6.

0. x^2 + y^2 + z^2 = 64

1. (x - 8*(2^0.5))^2 + y^2 + z^2 = 64

2. (x + 8*(2^0.5))^2 + y^2 + z^2 = 64

3. x^2 + (y - 8*(2^0.5))^2 + z^2 = 64

4. x^2 + (y + 8*(2^0.5))^2 + z^2 = 64

5. x^2 + y^2 + (z - 8*(2^0.5))^2 = 64

6. x^2 + y^2 + (z + 8*(2^0.5))^2 = 64

DTB.

[danbaron]

Here is the solution for this problem.

Dan

------------------------------------------------------------------------------------

REVISION: 2011-12-11

Now, I see that I should have explained what "a" is in the calculus formula (below).


V = 4*pi*R^3/3 - 12*pi*(2*R^3/3 - R^2*a + a^3/3).

'For, R = 8, and, a = 4*Sqr(2),

'V = 650.4773974..
In a cross section showing the xy plane, Spheres 0-4 will be displayed, and will appear as circles with radius 8.

Say, a straight line is drawn through the two points of intersection of Circles 0 and 1.

Say, a line segment is drawn perpendicular from that line to the center of Circle 0.

The length of that segment, is, "a".

You can do the same thing with Circles, 0 and 2, 0 and 3, and, 0 and 4.

And, you can also do the same thing in the xz and yz planes, with two other sets of perimeter circles.

------------------------------------------------------------------------------------


'-----------------------------------------------------------------------------------
Uses "console"
Uses "math"
Uses "file"
'-----------------------------------------------------------------------------------

'Will use the Monte Carlo Method.

'Is named after the casinos of Monte Carlo.

'The method is usually easy to program for a particular problem. On the other
'hand, its accuracy is proportional to the number of trials that are simulated.
'In other words, the accuracy is proportional to how long the program runs.

'The accuracy of this method also greatly depends on the ability of the
'pseudo-random number generator to produce a uniform distribution, and the length
'of the cycle before the distribution begins to repeat.

'Will generate random points in Sphere 0, and count the number of points which
'are not also contained in the other spheres. The problem is symmetric about
'three planes, x = 0, y = 0, z = 0. Therefore, only need to generate points in
'the first octant.

'pv = partial volume of Sphere 0
'tv = total volume of Sphere 0
'np = number of points in Sphere 0, which are not in any of the other spheres
'tp = total number of points

'pv = tv * np / tp

'-----------------------------------------------------------------------------------

'RESULTS

'I ran the program twice, once using thinBasic's pseudo-random number generator
'(RNG), and once using an RNG I found in a book. Both times, I generated 10
'million points. Here are the results.

'thinBasic RNG

'# of points volume
'01000000 651.063192
'02000000 650.383335
'03000000 650.807978
'04000000 650.832641
'05000000 650.553621
'06000000 650.283608
'07000000 650.204409
'08000000 650.367518
'09000000 650.599946
'10000000 650.608953

'book RNG

'# of points volume
'01000000 650.850871
'02000000 650.992419
'03000000 650.614243
'04000000 650.992955
'05000000 650.506867
'06000000 650.661426
'07000000 650.633341
'08000000 650.688145
'09000000 650.572065
'10000000 650.478772

'-----------------------------------------------------------------------------------

'I also solved the problem using calculus. (I guess in most cases, the Monte
'Carlo Method is used for problems which cannot be solved using calculus.)

'V = 4*pi*R^3/3 - 12*pi*(2*R^3/3 - R^2*a + a^3/3).

'For, R = 8, and, a = 4*Sqr(2),

'V = 650.4773974..

'(The volume of the entire sphere, is just, 4*pi*R^3/3, or,

'V = 2144.6605848..)

'-----------------------------------------------------------------------------------

'One thing to note, is that for a three dimensional Monte Carlo Method problem,
'like this one, the range of the random points generated, must be a box. I used,

'x = %radius * Rnd
'y = %radius * Rnd
'z = %radius * Rnd.

'This forms a cube with sides 8. For each point generated, I determined if the
'point was contained in Sphere 0. If so, then I continued using it in the
'calculations. If not, then I discarded that point, and generated another one.

'If, instead, you use spherical coordinates, and generate points as shown below,
'then, the calculations will be incorrect. It is not apparent, but, x, y, and z,
'in points generated as shown below, are not independent. I think the reason is
'that, each point has been "forced" to reside within Sphere 0.

'r = Rnd * %radius
'phi = Rnd * Pi / 2
'theta = Rnd * Pi / 2

'x = r * Sin(phi) * Cos(theta)
'y = r * Sin(phi) * Sin(theta)
'z = r * Cos(phi)

'-----------------------------------------------------------------------------------

%radius = 8
%tottrials = 10 ^ 07
%updateinterval = 10 ^ 06

'sphere volume = 4/3*pi*r^3

%spherevol = 4 / 3 * Pi * %radius ^ 3 As Ext
%poffset = %radius * Sqr(2) As Ext
%noffset = - %poffset As Ext

'Global filestring As String = "Lovecraft3thin.txt"
Global filestring As String = "Lovecraft3C.txt"

Global outfile As DWord
Global origin(6, 3) As Ext

'-----------------------------------------------------------------------------------

Function tbmain()

Ext phi, theta
Ext x, y, z
DWord i

DWord totpoints = 0
DWord numgoodpoints = 0
Quad idum0 = -1234

origin(1, 1) = %noffset
origin(1, 2) = 0
origin(1, 3) = 0

origin(2, 1) = %poffset
origin(2, 2) = 0
origin(2, 3) = 0

origin(3, 1) = 0
origin(3, 2) = %noffset
origin(3, 3) = 0

origin(4, 1) = 0
origin(4, 2) = %poffset
origin(4, 3) = 0

origin(5, 1) = 0
origin(5, 2) = 0
origin(5, 3) = %noffset

origin(6, 1) = 0
origin(6, 2) = 0
origin(6, 3) = %poffset

outfile = FILE_Open(filestring, "OUTPUT")
'Randomize(789)

Do

'x = %radius * Rnd
'y = %radius * Rnd
'z = %radius * Rnd

x = %radius * ran2(idum0)
y = %radius * ran2(idum0)
z = %radius * ran2(idum0)

If pointinsphere(x, y, z) Then
totpoints += 1
Else
Iterate Do
End If

numgoodpoints += nointersection(x, y, z)

If Mod(totpoints, %updateinterval) = 0 Then writeoutput(numgoodpoints, totpoints)

If totpoints = %tottrials Then Exit Do

Loop

FILE_Close(outfile)

WaitKey

End Function

'-----------------------------------------------------------------------------------

Function pointinsphere(px As Ext, py As Ext, pz As Ext) As Integer
Ext Static r

r = Sqr(px ^ 2 + py ^ 2 + pz ^ 2)
If r > %radius Then Return FALSE
Return TRUE

End Function

'-----------------------------------------------------------------------------------

Function nointersection(px As Ext, py As Ext, pz As Ext) As Integer
Ext Static d
Integer i

For i = 2 To 6 Step 2
d = Sqr((px - origin(i, 1)) ^ 2 + (py - origin(i, 2)) ^ 2 + (pz - origin(i, 3)) ^ 2)
If d < %radius Then Return FALSE
Next
Return TRUE

End Function

'-----------------------------------------------------------------------------------

Function writeoutput(ngp As DWord, j As DWord)
Ext partvol = %spherevol * ngp / j
String numfor = Format$(j, "00000000")
String volfor = Format$(partvol, "000.000000")
String allstring = numfor & " " & volfor

FILE_LinePrint(outfile, allstring)
Console_WriteLine(allstring)

End Function

'-----------------------------------------------------------------------------------

Sub ran2(ByRef idum As Quad) As Ext
'pseudo random EXT generator
'from "Numerical Recipes in C", 2nd edition, p.282.
'0 < RAN2() < 1

'Set idum negative to seed function.
'Then, do not alter its value between function calls.

Static maxj As Quad = -10
Static minj As Quad = 100
Static im1 As Quad = 2147483563
Static im2 As Quad = 2147483399
Static am As Ext = 1.0 / im1
Static imm1 As Quad = im1 - 1
Static ia1 As Quad = 40014
Static ia2 As Quad = 40692
Static iq1 As Quad = 53668
Static iq2 As Quad = 52774
Static ir1 As Quad = 12211
Static ir2 As Quad = 3791
Static ntab As Quad = 32
Static ndiv As Quad = 1 + imm1 / ntab
Static ndive As Ext = ndiv
Static eps As Ext = 10^(-18)
Static rnmx As Ext = 1.0 - eps
Static idum2 As Quad = 123456789
Static iy As Quad = 0
Static iv(ntab) As Quad
Static temp As Ext
Static scale As Ext = 32.0 / 33.0
Local je As Ext
Local j, k As Quad

If idum <= 0 Then
idum = Abs(idum)
If idum < 1 Then idum = 1
idum2 = idum
For j = ntab + 8 To 1 Step -1
k = idum / iq1
idum = ia1 * (idum - k * iq1) - k * ir1
If idum < 0 Then idum += im1
If j <= ntab Then iv(j) = idum
Next j
iy = iv(1)
End If

k = idum / iq1
idum = ia1 * (idum - k * iq1) - k * ir1
If idum < 0 Then idum += im1
k = idum2 / iq2
idum2 = ia2 * (idum2 - k * iq2) - k * ir2
If idum2 < 0 Then idum2 += im2
je = 1.0 + iy / ndive
j = scale * je
iy = iv(j) - idum2
iv(j) = idum
If iy < 1 Then iy += imm1
temp = am * iy
If temp > rnmx Then
Return rnmx
Else
Return temp
End If

End Function

'-----------------------------------------------------------------------------------

[ErosOlmi]

Numerical Recipes in C
The Art of Scientific Computing
William H. Press - Brian P. Flannery
Saul A. Teukolsky - William T. Vetterling
Cambridge University Press
http://www.nr.com/nrccover.gif
Dan,

thanks a lot for your script very well documented.
I have the book you mentioned but I have first edition
Worth to add into thinBasic arsenal your translated pseudo random generator. I think in the book it is the function called "Portable Random Number Generator"

Thanks again.
Eros

[danbaron]

In the second edition (hardcover), under the section, "Portable Random Number Generators" (page 278 ), there are three versions, "ran1()", "ran2()", and "ran3()". ran2() appears on page 282. On page 281, the authors state that they will pay a reader $1000 if he can find a flaw in ran2(). Beneath the function heading, it states that the period for ran2() is greater than 10^18. I translated it as best I could. I had to put an extra variable in the function, "Static scale As Ext = 32.0 / 33.0". Without that factor, "scale", no matter what I did, I would get an array index violation for "iv(ntab)". The indices should go from 0 to 31, for the function as it appears in the book. For thinBasic, the indices become from 1 to 32. No matter what I tried, without "scale", I would get an index violation, saying I illegally accessed either iv(0) or iv(33). "scale" fixed it, but I don't really understand what was happening.

You don't need to thank me. You put in the section for math. I thank you for doing so.

Dan

[ErosOlmi]

Ok, added Ran2 function into Math module.
Will be present in next release.
With a compiled Ran2 function execution speed is improved quite a lot.

I was thinking about to create a general form of

Function pointinsphere(px As Ext, py As Ext, pz As Ext) As Integer
Function nointersection(px As Ext, py As Ext, pz As Ext) As Integer