 |
Example - Gravity
* The "\" symbol has been used at the end of lines to indicate that the line was too long for the display and has been
wrapped to the line below in order to fit within the width of the web page. These lines would need to be joined in the normal
program editor.
I. Introduction
This program models the dynamics of several bodies under gravity, such as the planets in our solar system. A finite element array is used to simulate the system.
II. Mathematics and Physics Involved
As Newton discovered, the attractive force exerted by a point mass m2 on another point mass m1 is  , where G is a constant, and r1 and r2 are the appropriate position vectors.
By considering the sum over all infinitesimal particles in two spheres, it can be shown that the force exerted by sphere 2 on sphere 1 is the same (with r1 and r2 representing the positions of the centres of mass) if the two spheres do not coincide, but if they do, the force is equal to a constant multiplied by the distance between the two centres of mass.
The force on sphere 1 due to more than one other body is obtained by summing the force due to each other body. Note that the net gravitational force exerted by a body on itself is zero.
III. How The Program Works
1) Initialisation
The finite element array is initialised after the user sets the number of points and the time period over which to simulate. The initial positions and velocities are determined randomly each time to give a different problem. The formula that is used to calculate the gravitational forces makes use of the Summation() method to sum over all points.
' set the initial position and velocity of each
' point, but make sure that the arithmetic mean of
' each is zero
' initial x position
FOR i%% = 1 TO n%%:init%(i%%) = (200 * RND (1)) - 100:NEXT i%%
meanval% = Summation("init%(i%%)","i%%",1,n%%) / n%%
FOR i%% = 1 TO n%%:f.x.setinitialvalue(160 + init%(i%%) - meanval% \
,0,i%%,i%%):NEXT i%%
' initial y position
FOR i%% = 1 TO n%%:init%(i%%) = (200 * RND (1)) - 100:NEXT i%%
meanval% = Summation("init%(i%%)","i%%",1,n%%) / n%%
FOR i%% = 1 TO n%%:f.y.setinitialvalue(160 + init%(i%%) - meanval% \
,0,i%%,i%%):NEXT i%%
' initial velocity in x direction
FOR i%% = 1 TO n%%:init%(i%%) = (24 * RND (1)) - 12:NEXT i%%
meanval% = Summation("init%(i%%)","i%%",1,n%%) / n%%
FOR i%% = 1 TO n%%:f.x.setinitialvalue(init%(i%%) - meanval%,1,i%%,\
i%%):NEXT i%%
' initial velocity in y direction
FOR i%% = 1 TO n%%:init%(i%%) = (24 * RND (1)) - 12:NEXT i%%
meanval% = Summation("init%(i%%)","i%%",1,n%%) / n%%
FOR i%% = 1 TO n%%:f.y.setinitialvalue(init%(i%%) - meanval%,1,i%%, \
i%%):NEXT i%%
fx$ = "if( ((x()-x(0,i%%))^2)+((y()-y(0,i%%))^2)<=4 ,(x(0,i%%)-x())/\
8,"
fx$ = fx$ + "((x(0,i%%)-x())/((((x(0,i%%)-x())^2)+((y(0,i%%)-y())^2)\
)^1.5)))"
fy$ = "if( ((x()-x(0,i%%))^2)+((y()-y(0,i%%))^2)<=4 ,(y(0,i%%)-y())/\
8,"
fy$ = fy$ + "((y(0,i%%)-y())/((((x(0,i%%)-x())^2)+((y(0,i%%)-y())^2)\
)^1.5)))"
fx2$ = "5000*Superbase.Summation(fx$,i$,1,n%%)"
fy2$ = "5000*Superbase.Summation(fy$,i$,1,n%%)"
i$ = "i%%"
f.x.setformula(fx2$,2)
f.y.setformula(fy2$,2)
f.advanceto(0.0)
Because the calculation speed of the positions in the finite element array varies due to their proximity to each other, a blank file is set up to store the positions of the elements as they are calculated using the CREATE MEMORY command to create a volatile file. This will later be used to play back the simulation at a constant speed.
fname$ = "points"
CREATE MEMORY fname$
ADD "t;num"
ADD "element;nmi"
ADD "xpos;nmi"
ADD "ypos;nmi"
MAKE fname$
FILE fname$
2) Updating of the Positions
Each time the array is advanced, the points are moved according to the values in the finite element array. The positions are then stored in the volatile file.
FOR i%% = 1 TO n%%
x% = f.x(0,i%%):x% = IF (x% < 10,10, IF (x% > 310,310,x%))
y% = f.y(0,i%%):y% = IF (y% < 40,40, IF (y% > 340,340,y%))
x$ = "forms.f.e" + STR$ (i%%,"000") + ".move(x%,y%,x%,y%)":\
EXECUTE x$
IF forms.f.TraceCheck.Value = - 1 THEN
IF lastx001% <> 0 THEN
SET fc = forms.f.Add("","line")
x$ = "fc.move(x%,y%,lastx" + STR$ (i%%,"000") + "%,lasty"\
+ STR$(i%%,"000") + "%)":EXECUTE x$
x$ = "fc.forecolor = " + STR$ (i%%,"000"):EXECUTE x$
fc.tag = "traceline"'so we can delet it later
END IF
x$ = "lastx" + STR$ (i%%,"000") + "% = x%":EXECUTE x$
x$ = "lasty" + STR$ (i%%,"000") + "% = y%":EXECUTE x$
END IF
forms.f.page.SetDirtyRect(x% - 5,Y% - 5,x% + 5,y% + 5)
forms.f.refresh()
BLANK
t.fname$ = f.t
element = i%%
xpos = x%
ypos = y%
STORE
NEXT i%%
3) Replay
When the simulation is finished, the user is offered the chance to watch the replay.
WHILE (replay%% <> 0)
i%% = 0
SELECT FIRST
t% = NOW + 50
WHILE NOT EOF (fname$)
i%% = i%% + 1
x$ = "forms.f.e" + STR$ (element,"000") + ".move(xpos,ypos,\
xpos,ypos)":EXECUTE x$
IF i%% > n%% THEN
forms.f.page.SetDirtyRect(xpos - 5,ypos - 5,xpos + 5,ypos\
+ 5)
ELSE
forms.f.page.SetDirtyRect(9,39,302,302)
END IF
forms.f.refresh()
IF ((i%% MOD n%%) = 0) THEN
SET STATUS "Time: " + STR$ (t.fname$,"999.000")
IF (t% > NOW ) THEN
t% = (t% - NOW ) / 1000
IF t% <= 0 THEN t% = .1
WAIT FOR t%
END IF
t% = NOW + 50
END IF
SELECT NEXT
WEND
REQUEST "Click OK to replay again","",1,replay%%
WEND
The points are displayed in a box as below:
|
