Project Lovecraft Problem 3

**danbaron** · 02-03-2010, 07:30

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** · 05-03-2010, 08:54

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** · 05-03-2010, 11:25

[info]

Numerical Recipes in C
The Art of Scientific Computing
William H. Press - Brian P. Flannery
Saul A. Teukolsky - William T. Vetterling
Cambridge University Press

[/info]
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** · 05-03-2010, 22:30

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** · 05-03-2010, 23:15

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
[code=thinbasic]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[/code]