Thursday 7 March 2013

Raspberry Pi - Minecraft - Planetary Gravity Sim

First off I cant take the credit for this, I came across this great post and video titled "A naive simulator of gravity, written in Python" by a guy called Thanassis Tsiodras which describes and provides the code for a planetary gravity simulator, and I 're-used' the code, changing a few bits to migrate this to use the minecraft api.

I don't have a complete grasp of all the code (maths has never been my strong point), but I managed to reverse engineer Thanassis's code to understand what I needed to change to 'downsize' the universe so it could be rendered within the Minecraft world.

The original code throws 30 random planets around a large sun and lets gravity do its job, usual a good percentage of the planets get drawn into the sun or collide with each other, those that don't usually end up in elliptical orbits around the sun.

This was a bit hardcore for the Pi to cope with rendering in Minecraft, so I simplified the code to initially create 3 planets that would sit in stable orbits; I then introduced a rouge planet into the system which throws everything out and the whole system is destroyed in a few cycles when the planets collide and eventually all fall into the sun.


Its really easy to change the planets in the simulation, just head to the line in the code which is headed with a comment #CREATE THE PLANETS HERE, and add as many as you want, giving them a position (x,y) and a velocity (x,y).


Download and run
You can download the code direct from git-hub, so run minecraft, open/create a world and follow the instructions:

sudo apt-get install git-core
cd ~
git clone https://github.com/martinohanlon/minecraft-planets.git
cd minecraft-planets
python minecraft-planets.py


The code
If you want to have a go yourself:

Create a directory for the program

mkdir ~/minecraft-planets

Copy the python api class library from minecraft to the programs directory

cp -r ~/mcpi/api/python/mcpi ~/minecraft-planets/minecraft

Create minecraft-planets.py python program

nano ~/minecraft-planets/minecraft-planets.py

or open Idle and save minecraft-planets.py to the minecraft-planets directory

Code
https://github.com/martinohanlon/minecraft-planets.git

#!/usr/bin/env python
"""
An updated version of Thanassis Tsiodras's gravity simulator which uses
minecraft as the method of rendering the planets!

Martin O'Hanlon
www.stuffaboutcode.com

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

An improved version of my Python-based gravity simulator, using Runge-Kutta
4th order solution of the differential equations - coded during Xmas 2012.
Happy holidays, everyone!

I've always been fascinated by space - ever since I read 'The Family of
the Sun', when I was young. And I always wanted to simulate what I've read
about Newton's gravity law, and see what happens in... a universe of my own
making :-)

So: The following code 'sprays' some 'planets' randomly, around a sun,
inside a 900x600 window (the values are below, change them at will).
Afterwards, it applies a very simple set of laws:

- Gravity, inversely proportional to the square of the distance, and linearly
  proportional to the product of the two masses
- Elastic collissions of two objects if they are close enough to touch: a
  merged object is then created, that maintains the momentum (mass*velocity)
  and the mass of the two merged ones.
- This updated version of the code is using the RK4 solution of the velocity/
  acceleration differential equation, and is in fact based on the excellent
  blog of Glenn Fiedler (http://gafferongames.com)

Use the numeric keypad's +/- to zoom in/out, and press SPACE to toggle
showing/hiding the orbits trace.

Blog post at:

    http://users.softlab.ntua.gr/~ttsiod/gravity.html
    http://ttsiodras.github.com/gravity.html

Thanassis Tsiodras
ttsiodras@gmail.com
"""

import sys
import math
import random
from collections import defaultdict

#MaOH
#import the minecraft.py module from the minecraft directory
import minecraft.minecraft as minecraft
#import minecraft block module
import minecraft.block as block
#import time, so delays can be used
import time

# The window size
# MaOH - smaller window or system, as in minecraft the blocks are much bigger!
WIDTH, HEIGHT = 50, 50
WIDTHD2, HEIGHTD2 = WIDTH/2., HEIGHT/2.

# The density of the planets - used to calculate their mass
# from their volume (i.e. via their radius)
DENSITY = 0.001

# The gravity coefficient - it's my universe, I can pick whatever I want :-)
GRAVITYSTRENGTH = 1.e4

# The global list of planets
g_listOfPlanets = []

# MaOH - Z POSITION IN MINECRAFT WORLD
ZPOSITION = 15

class State:
    """Class representing position and velocity."""
    def __init__(self, x, y, vx, vy):
        self._x, self._y, self._vx, self._vy = x, y, vx, vy

    def __repr__(self):
        return 'x:{x} y:{y} vx:{vx} vy:{vy}'.format(
            x=self._x, y=self._y, vx=self._vx, vy=self._vy)


class Derivative:
    """Class representing velocity and acceleration."""
    def __init__(self, dx, dy, dvx, dvy):
        self._dx, self._dy, self._dvx, self._dvy = dx, dy, dvx, dvy

    def __repr__(self):
        return 'dx:{dx} dy:{dy} dvx:{dvx} dvy:{dvy}'.format(
            dx=self._dx, dy=self._dy, dvx=self._dvx, dvy=self._dvy)


class Planet:
    """Class representing a planet. The "_st" member is an instance of "State",
    carrying the planet's position and velocity - while the "_m" and "_r"
    members represents the planet's mass and radius."""
    #MaOH - changed the parameters of the planet and so I can pass in an initial state
    def __init__(self, blockType, initialState=None):
        if initialState != None:
            self._st = initialState
        else:
            # otherwise pick a random position and velocity
            self._st = State(
               float(random.randint(0, WIDTH)),
               float(random.randint(0, HEIGHT)),
               float(random.randint(0, 40)/100.)-0.2,
               float(random.randint(0, 40)/100.)-0.2)
            
        self._r = 0.55
        
        self.setMassFromRadius()
        self._merged = False
        #MaOH - create a last vec3 which is representative of where the planet is now
        self._lastVec3 = minecraft.Vec3(roundNo(self._st._x), roundNo(self._st._y), ZPOSITION)
        self._blockType = blockType

    def __repr__(self):
        return repr(self._st)

    def acceleration(self, state, unused_t):
        """Calculate acceleration caused by other planets on this one."""
        ax = 0.0
        ay = 0.0
        for p in g_listOfPlanets:
            if p is self or p._merged:
                continue  # ignore ourselves and merged planets
            dx = p._st._x - state._x
            dy = p._st._y - state._y
            dsq = dx*dx + dy*dy  # distance squared
            dr = math.sqrt(dsq)  # distance
            force = GRAVITYSTRENGTH*self._m*p._m/dsq if dsq>1e-10 else 0.
            # Accumulate acceleration...
            ax += force*dx/dr
            ay += force*dy/dr
        return (ax, ay)

    def initialDerivative(self, state, t):
        """Part of Runge-Kutta method."""
        ax, ay = self.acceleration(state, t)
        return Derivative(state._vx, state._vy, ax, ay)

    def nextDerivative(self, initialState, derivative, t, dt):
        """Part of Runge-Kutta method."""
        state = State(0., 0., 0., 0.)
        state._x = initialState._x + derivative._dx*dt
        state._y = initialState._y + derivative._dy*dt
        state._vx = initialState._vx + derivative._dvx*dt
        state._vy = initialState._vy + derivative._dvy*dt
        ax, ay = self.acceleration(state, t+dt)
        return Derivative(state._vx, state._vy, ax, ay)

    def updatePlanet(self, t, dt):
        """Runge-Kutta 4th order solution to update planet's pos/vel."""
        a = self.initialDerivative(self._st, t)
        b = self.nextDerivative(self._st, a, t, dt*0.5)
        c = self.nextDerivative(self._st, b, t, dt*0.5)
        d = self.nextDerivative(self._st, c, t, dt)
        dxdt = 1.0/6.0 * (a._dx + 2.0*(b._dx + c._dx) + d._dx)
        dydt = 1.0/6.0 * (a._dy + 2.0*(b._dy + c._dy) + d._dy)
        dvxdt = 1.0/6.0 * (a._dvx + 2.0*(b._dvx + c._dvx) + d._dvx)
        dvydt = 1.0/6.0 * (a._dvy + 2.0*(b._dvy + c._dvy) + d._dvy)
        self._st._x += dxdt*dt
        self._st._y += dydt*dt
        self._st._vx += dvxdt*dt
        self._st._vy += dvydt*dt

    def setMassFromRadius(self):
        """From _r, set _m: The volume is (4/3)*Pi*(r^3)..."""
        self._m = DENSITY*4.*math.pi*(self._r**3.)/3.

    def setRadiusFromMass(self):
        """Reversing the setMassFromRadius formula, to calculate radius from
        mass (used after merging of two planets - mass is added, and new
        radius is calculated from this)"""
        self._r = (3.*self._m/(DENSITY*4.*math.pi))**(0.3333)

    #MaOH - functions to draw planet in minecraft
    def drawPlanet(self, mc):
        newVec3 = minecraft.Vec3(roundNo(self._st._x), roundNo(self._st._y), ZPOSITION)
        #has the planet moved
        if (newVec3.x != self._lastVec3.x) or (newVec3.y != self._lastVec3.y) or (newVec3.z != self._lastVec3.z):
            # clear last block
            self.drawPlanetInMC(mc, self._lastVec3, self._r, block.AIR)
            # draw new planet
            self.drawPlanetInMC(mc, newVec3, self._r, self._blockType)
            self._lastVec3 = newVec3

    def clearPlanet(self, mc):
        self.drawPlanetInMC(mc, self._lastVec3, self._r, block.AIR)
        
    def drawPlanetInMC(self, mc, vec3, radius, blockType):
        # if the diameter is greater than 1, create a sphere, otherwise its just a block!
        if radius * 2.0 > 1.5:
            #round up radius
            radius = int(radius + 0.5)
            for x in range(radius*-1,radius):
                for y in range(radius*-1, radius):
                    for z in range(radius*-1,radius):
                        if x**2 + y**2 + z**2 < radius**2:
                            mc.setBlock(vec3.x + x, vec3.z + z, vec3.y + y, blockType)
        else:
            mc.setBlock(vec3.x, vec3.z, vec3.y, blockType)



def roundNo(number):
    return int(number + 0.5)

def main():

    #Connect to minecraft by creating the minecraft object
    # - minecraft needs to be running and in a game
    mc = minecraft.Minecraft.create()

    mc.postToChat("Hi, Minecraft Planetary Gravity Simulator")
    mc.postToChat("www.stuffaboutcode.com")
    
    #clear an area
    mc.setBlocks(0 - WIDTH, ZPOSITION - 5, 0 - HEIGHT, WIDTH, ZPOSITION + 5, HEIGHT, block.AIR)
    time.sleep(5)
    
    global g_listOfPlanets, PLANETS

    # And God said: Let there be lights in the firmament of the heavens...
    g_listOfPlanets = []
    #MaOH - changed to create a fixed collection of planets, rather than random
    #CREATE THE PLANETS HERE Planet(typeOfBlock, State(x position, y position, x velocity, y velocity))
    # nice planets in sensible orbits
    g_listOfPlanets.append(Planet(block.GOLD_BLOCK, State(-10, 0, 0, 0.2)))
    g_listOfPlanets.append(Planet(block.DIAMOND_BLOCK, State(10, 0, 0, -0.2)))
    g_listOfPlanets.append(Planet(block.SNOW_BLOCK, State(20, 0, 0, 0.15)))
    # rouge planet 
    #g_listOfPlanets.append(Planet(block.GOLD_ORE, State(-12, -12, 0, 0.15)))

    def planetsTouch(p1, p2):
        dx = p1._st._x - p2._st._x
        dy = p1._st._y - p2._st._y
        dsq = dx*dx + dy*dy
        dr = math.sqrt(dsq)
        return dr<=(p1._r + p2._r)

    #MaOH - the sun is now the centre of the minecraft universe
    sun = Planet(block.GLOWING_OBSIDIAN, State(0, 0, 0, 0))
    
    sun._m *= 100
    sun.setRadiusFromMass()
    g_listOfPlanets.append(sun)

    #MaOH - draw the sun
    sun.drawPlanetInMC(mc, sun._lastVec3, sun._r, sun._blockType)

    #merge planets that are inside the sun
    for p in g_listOfPlanets:
        if p is sun:
            continue
        if planetsTouch(p, sun):
            p._merged = True  # ignore planets inside the sun

    # t and dt are unused in this simulation, but are in general, 
    # parameters of engine (acceleration may depend on them)
    t, dt = 0., 1.

    while True:
        t += dt
        for p in g_listOfPlanets:
            if not p._merged:  # for planets that have not been merged, draw a
                # circle based on their radius, but take zoom factor into account
                #MaOH - draw planet
                p.drawPlanet(mc)

        # Update all planets' positions and speeds (should normally double
        # buffer the list of planet data, but turns out this is good enough :-)
        for p in g_listOfPlanets:
            if p._merged or p is sun:
                continue
            # Calculate the contributions of all the others to its acceleration
            # (via the gravity force) and update its position and velocity
            p.updatePlanet(t, dt)

        # See if we should merge the ones that are close enough to touch,
        # using elastic collisions (conservation of total momentum)
        for p1 in g_listOfPlanets:
            if p1._merged:
                continue
            for p2 in g_listOfPlanets:
                if p1 is p2 or p2._merged:
                    continue
                if planetsTouch(p1, p2):
                    if p1._m < p2._m:
                        p1, p2 = p2, p1  # p1 is the biggest one (mass-wise)
                    p2._merged = True
                    #MaOH - if the planet is merged, get rid of it
                    p2.clearPlanet(mc)
                    if p1 is sun:
                        continue  # No-one can move the sun :-)
                    newvx = (p1._st._vx*p1._m+p2._st._vx*p2._m)/(p1._m+p2._m)
                    newvy = (p1._st._vy*p1._m+p2._st._vy*p2._m)/(p1._m+p2._m)
                    p1._m += p2._m  # maintain the mass (just add them)
                    p1.setRadiusFromMass()  # new mass --> new radius
                    p1._st._vx, p1._st._vy = newvx, newvy

if __name__ == "__main__":
    try:
        main()
    except KeyboardInterrupt:
        print "stopped"

16 comments:

  1. VERY interesting. (Though I assume you meant Rogue planet, not Rouge.)

    ReplyDelete
    Replies
    1. O dear, my poor spelling, good job I have skills in other areas

      Delete
  2. I want to get this up so i can see how the script works, but im so confused

    ReplyDelete
    Replies
    1. Ive updated the post with instructions on how you can download the code and run it, hopefully that will make it easier to get up and running.

      What you are trying to do? Whats confusing you? Let me know and I'll try and help.

      Delete
  3. Hello,

    I know this may not be the proper place to put a request for help but I have tried to search for this and have come up with nothing.

    I have tried the usual scripts for clearing the world (fill with AIR -127x -64y -127z by 127x 64y 127z, floor at -64y), and this works except for one little itty bitty nagging problem. The world shifts, 0x0y0z is no longer the center. In one try the X min/max is -240 and 16, Z is -113 and 141, y remains -64 and unknown. I keep shifting the min/max for the world and it moves again. I have tried with freshly created worlds and the problem persists.

    Can anyone tell me how to create or where to download a centered world?
    Is there some file I can edit to fix this?

    ReplyDelete
    Replies
    1. I think what you are saying is that when try and clear a whole world using setBlocks, the centre of the world isnt always the same.

      The only suggestion I can make is to go 'over the edges' i.e. clear a bigger area than you need. you can set a block outside the boundry of the world, it just wont show up.

      Delete
    2. This comment has been removed by the author.

      Delete
  4. (python script names are from the web, barely modified, mainly co-ordinates, commented out server references since only access local Minecraft, added line with minecraft.Minecraft.create())
    The Minecraft Pi I am using is the stock version that came with my RPi B+ Raspian.

    I started fresh.
    Deleted the world
    created fresh world
    Using playerposition.py script (I know the co-ordinates are also displayed at the top, but I wanted comparison):
    I confirmed the x/z size of the world (-127, -127 by 127, 127)
    y min at -64 (lowest floor), y death at -128.
    Ran my combo clearall/playerposition.py script:
    it cleared the world from -254, -254, -254 by 254, 254, 254
    (8 times the size of a normal world)
    created 3 axis walls:
    (0, -64, -127, 0, 64, 127)
    (-127, 0, -127, 127, 0, 127) (non-death floor)
    (-127, -64, 0, 127, 64, 0)
    spawned (funny thing, I was inside the wall, whacka, whacka, whacka)

    I then went to confirm the size of the world again. While the overall size was (-127, -64, -127, 127, ?, 127) (? seems to be positive infinity), the z coordinate for the axis walls stopped at 126 (x/z, y=0, and y/z, x=0). You can walk over the edge and die, or fly around without being stopped, this behavior is seen in every world I have tried. This says to me the z axis shifted by -1 (to -128) this time, which can't be traveled to, and you can't set any blocks with a z=127.
    This shifting gets worse the more times you clear the world, and not always just the z axis (can be x or both), but doesn't affect y. One x was min/max -240/-16, z was -113/141 as stated previously.

    ReplyDelete
    Replies
    1. I've written a new script to delete/build the world (from freshly generated) using the old reliable method: one block at a time, from the top to the bottom (y: 64 to -64). It is really really slow (changing 8388225 blocks individually, and I haven't overclocked my CPU).

      I'll let you know how it goes.

      Delete
  5. Well, the initial test took 8+ hours to run.

    ReplyDelete
  6. Looks like any form of programmatic solution is out. Every routine, regardless of how simplistic, when ran on the whole 'world' results in an axis shift.

    Used co-ordinates: (x and z) -127 to 127, (y) -63 to 64 (-64 wasn't used because it would have resulted in no 'floor'. Co-ordinates confirmed before altering world using program(s), and after program(s) ran. Minimum shift is 1 on the x/z axes, and up to a maximum of 120. Y axis appears to be unaltered.

    ReplyDelete
  7. Hello again,

    Well I finally found something that works. I used a variation of another's code which I will post below. Instructions are commented inside. As long as the initial world's corner points have variations on -127 and 127, I have gotten this to work every time. It only takes about 5 minutes to run and all the boundries are marked (solid bottom, and lines for the corners and top edges).

    Code will be in following post due to size.

    ReplyDelete
  8. #!/usr/bin/env python

    #original idea and instructions on ghex editing file to rename 'world'
    #http://winkleink.blogspot.co.uk/2013/02/raspberry-pi-minecraft-and-python.html
    #expanded for maximum area and dimensional stability-
    #using the original code resulted in x and z axes randomly shifting on my RPi.
    #When finished, you will have a mostly empty world: a floor (the lowest bedrock level),
    #4 corner markers, 4 upper boundary markers, and a (0,0,0) Diamond block marker.
    #Or you can change column and center markers as desired.


    #Instructions/Notes
    #1) Clean out the world folder/directory
    # home/pi/.minecraft/games/com.mojang/minecraftWorlds
    # delete or move world variants from folder, if any
    #2) Start Minecraft PI (server) Program
    # create a new world
    #3) walk/fly to 4 corners to confirm co-ordinates
    # x and z should be various combinations of -127 and 127
    # (-127,y,-127), (-127,y,127), (127,y,-127), (127,y,127)
    # If not, restart from step 1.
    # This seems to lock the corners.
    # stay near one of these points. DO NOT GO BACK TO THE CENTER.
    #4) Run this script in Idle (not Idle 3) or a command line may also work
    # This takes several minutes to run, be patient.
    # Some of the world may clear as you look around and it is re-rendered.
    # if you were walking, you will drop to bottom bedrock 'floor' at level -64
    # (your Y co-ordinatewill be -63 since you are standing on the floor)
    # check all 4 corners again, a bedrock column should be at each
    #5) Exit world/game/server
    #6) Restart Minecraft PI (server) Program
    # reload world
    #7) Recheck the 4 columns to make sure they are still there
    #8) If all is well, exit the world, make a backup, rename if you can.
    # Renaming the folder will allow creating multiple worlds (changes line above 'Creative')
    # Use the instructions on the above webpage to edit 'level.dat'/change world name
    # (changes 1st line, 5 characters only)
    # You can change worlds by swiping then to the left (or right, if possible) - as needed -
    # on the World Select screen.
    # If not, restart from step 1, this may take several tries.
    #
    #Interesting tidbits:
    #If you fly to the top, you can land outside the world and walk around,
    #but you can't build anything.
    #Change the corner and boundary block type to a flowing type of block and see what happens!!


    #import the minecraft.py module from the minecraft directory
    import mcpi.minecraft as minecraft
    #import minecraft block module
    import mcpi.block as block


    if __name__ == "__main__":
    mc = minecraft.Minecraft.create()
    #clear world area
    # I know the co-ordinates are overkill but this seems the only
    # way to clear Everything without the x/z co-ordinates shifting
    mc.setBlocks(-254, -63, -254, 254, 63, 254, block.AIR.id)

    #mark x,y,z center (0,0,0)
    mc.postToChat("Putting a diamond block at 0,0,0")
    mc.setBlock(0,0,0,block.DIAMOND_BLOCK.id)

    #apparently the setblocks routine goes from positive to negative
    #originally I used -127, but the column was always 1 off
    #on the x, z, or both axes; changing to -128 fixed this.
    blockid=block.LAPIS_LAZULI_BLOCK.id
    #blockid=block.WATER.id
    mc.postToChat("Marking 4 Corners")
    mc.setBlocks(-128,-63,-128,-128, 63,-128,blockid)
    mc.setBlocks(-128,-63, 127,-128, 63, 127,blockid)
    mc.setBlocks( 127,-63, 127, 127, 63, 127,blockid)
    mc.setBlocks( 127,-63,-128, 127, 63,-128,blockid)

    #upper boundry markers
    mc.postToChat("Marking 4 Upper Boundary Markers")
    mc.setBlocks(-128, 63,-128, 127, 63,-128,blockid)
    mc.setBlocks( 127, 63,-128, 127, 63, 127,blockid)
    mc.setBlocks( 127, 63, 127,-128, 63, 127,blockid)
    mc.setBlocks(-128, 63, 127,-128, 63,-128,blockid)

    ReplyDelete
  9. This comment has been removed by the author.

    ReplyDelete
  10. Note: It looks like the indention are gone from above. The all lines after If are supposed to be indented.

    After it completes, if the program is still running, you can terminate the program.

    ReplyDelete

Note: only a member of this blog may post a comment.