Since this answer was originally posted, the Blender software package has undergone significant revision causing the Python scripts to break. I have modified the scripts so that they should work from v2.79b to v2.93LTS.

Updated Blender scripts

Note: I added an update below to import position and orientation transformations and view the 3D simulation results in Mathematica.

I have used the free program Blender v2.79b to simulate the handling of 100s of complex shapes through a geometrically complex industrial machine with many moving parts including vibrating elements. So, it should be able to handle a "coin flip". I believe Blender still uses the Bullet Physics Engine as its solver. I should warn you that collision simulation can get difficult and that there are a lot of tricks of the trade you must learn to be accurate and fast in a general case.

Blender has a python interface and it can be run as a background task (Bullet also has a python interface, but I am not familiar with its operation). Since Mathematica can create text files with StringTemplate and can execute system commands, we should be able to create a python script to drive a Blender simulation.

Python Script Generation

Blender has a fairly well-documented API and there many resources one can find online to generate python script.

import bpy
from math import pi

# Read current Blender version
version = bpy.app.version

for o in bpy.data.objects:
    if version < (2, 80, 0):
        if o.type == 'MESH' or o.type == 'EMPTY':
            o.select = True
        else:
            o.select = False
    else:
        o.select_set(o.type == 'EMPTY' or o.type == 'MESH')

# Delete all objects in the scene
bpy.ops.object.delete()

# Add the floor
if version < (2, 80, 0):
    bpy.ops.mesh.primitive_cube_add(radius=5, location=(0, 0, 0))
else:
    bpy.ops.mesh.primitive_cube_add(size=10, location=(0, 0, 0))
bpy.ops.transform.resize(value=(1, 1, 0.1))
bpy.ops.rigidbody.objects_add(type='PASSIVE')
boxObj = bpy.context.active_object
boxObj.rigid_body.collision_shape = "BOX"
boxObj.name = "Ground"

# Add the Coin
bpy.ops.mesh.primitive_cylinder_add(radius=1, depth=0.1, location=(0, 0, 3))
bpy.ops.rigidbody.objects_add(type='ACTIVE')
boxObj = bpy.context.active_object
boxObj.rigid_body.collision_shape = "CYLINDER"
bpy.context.object.rigid_body.friction = 0.25
bpy.context.object.rigid_body.restitution = 0.75
boxObj.name = "Coin"
# Set reference to the coin
coin = bpy.data.objects["Coin"]

# Set a reference to the scene
sce = bpy.context.scene
# Set first frame
sce.frame_set(1)
# Set Keyframes
coin.keyframe_insert(data_path="location")
coin.keyframe_insert(data_path="rotation_euler")
bpy.context.object.rigid_body.kinematic = True
bpy.context.object.keyframe_insert('rigid_body.kinematic')

# Advance two frames and add translational and rotational motion
sce.frame_set(4)
# Translate up a little
coin.location.z = 3.45
# Rotate coin predominantly around the x-axis
coin.rotation_euler.x = 1
coin.rotation_euler.y = 0.1
coin.rotation_euler.z = 0.1
# Set Keyframes
coin.keyframe_insert(data_path="location")
coin.keyframe_insert(data_path="rotation_euler")
bpy.context.object.rigid_body.kinematic = False
bpy.context.object.keyframe_insert('rigid_body.kinematic')

# Set frame to the end
sce.frame_set(250)

# Bake rigid body simulation
override = {'scene': bpy.context.scene,
            'point_cache': bpy.context.scene.rigidbody_world.point_cache}
# bake to current frame
bpy.ops.ptcache.bake(override, bake=False)

# Get transformations
tr = coin.matrix_world.translation
eu = coin.matrix_world.to_euler()
print("          X                  Y                  Z                  RX                  RY                  RZ")
print(tr.x, tr.y, tr.z, eu.x, eu.y, eu.z)

if eu.x > pi / 2.:
    print("Coin is heads")
else:
    print("Coin is tails")

The script above does is the following:

1. Selects and deletes all extraneous objects in the scene.
2. Adds a passive cube (note that cubes are better than planes for collisions) for the ground.
3. Adds an active cylindrical coin.
   • Sets friction to 0.25
   • Sets restitution to 0.75
4. Adds keyframe animate to set initial velocity and rotation.
5. Prints out and tests if coin has rotated more than $frac{pi}{2}$ radians to determine heads/tails.

Driving from Mathematica

We can create a parametric model in Mathematica by replacing hard coded parameters with template variables using `` delimiters as in the createCoinFlip function.

createCoinFlip[z_, rx_, ry_, rz_, friction_, restitution_] := 
 StringTemplate["import bpy
from math import pi

# Read current Blender version
version = bpy.app.version

for o in bpy.data.objects:
    if version < (2, 80, 0):
        if o.type == 'MESH' or o.type == 'EMPTY':
            o.select = True
        else:
            o.select = False
    else:
        o.select_set(o.type == 'EMPTY' or o.type == 'MESH')

# Delete all objects in the scene
bpy.ops.object.delete()

# Add the floor
if version < (2, 80, 0):
    bpy.ops.mesh.primitive_cube_add(radius=5, location=(0, 0, 0))
else:
    bpy.ops.mesh.primitive_cube_add(size=10, location=(0, 0, 0))
bpy.ops.transform.resize(value=(1, 1, 0.1))
bpy.ops.rigidbody.objects_add(type='PASSIVE')
boxObj = bpy.context.active_object
boxObj.rigid_body.collision_shape = "BOX"
boxObj.name = "Ground"

# Add the Coin
bpy.ops.mesh.primitive_cylinder_add(radius=1, depth=0.1, 
location=(0, 0, 3))
bpy.ops.rigidbody.objects_add(type='ACTIVE')
boxObj = bpy.context.active_object
boxObj.rigid_body.collision_shape = "CYLINDER"
bpy.context.object.rigid_body.friction = `friction`
bpy.context.object.rigid_body.restitution = `restitution`
boxObj.name = "Coin"
# Set reference to the coin
coin = bpy.data.objects["Coin"]

# Set a reference to the scene
sce = bpy.context.scene
# Set first frame
sce.frame_set(1)
# Set Keyframes
coin.keyframe_insert(data_path="location")
coin.keyframe_insert(data_path="rotation_euler")
bpy.context.object.rigid_body.kinematic = True
bpy.context.object.keyframe_insert('rigid_body.kinematic')

# Advance two frames and add translational and rotational motion
sce.frame_set(4)
# Translate up a little
coin.location.z = `z`
# Rotate coin predominantly around the x-axis
coin.rotation_euler.x = `rx`
coin.rotation_euler.y = `ry`
coin.rotation_euler.z = `rz`
# Set Keyframes
coin.keyframe_insert(data_path="location")
coin.keyframe_insert(data_path="rotation_euler")
bpy.context.object.rigid_body.kinematic = False
bpy.context.object.keyframe_insert('rigid_body.kinematic')

# Set frame to the end
sce.frame_set(250)

# Bake rigid body simulation
override = {'scene': bpy.context.scene,
            'point_cache': 
bpy.context.scene.rigidbody_world.point_cache}
# bake to current frame
bpy.ops.ptcache.bake(override, bake=False)

# Get transformations
tr = coin.matrix_world.translation
eu = coin.matrix_world.to_euler()
print("
          X                  Y                  Z                  RX                  R
Y                  RZ")
print(tr.x, tr.y, tr.z, eu.x, eu.y, eu.z)

if eu.x > pi / 2.:
    print("Coin flip result is heads")
else:
    print("Coin flip result is tails")
"][<|"z" -> z, "rx" -> rx, "ry" -> ry, "rz" -> rz, 
"friction" -> friction, "restitution" -> restitution|>]

Blender will send a lot of information to standard out. We can parse this output with Find to extract a line of interest. Putting it all together, the following will create a python script, run Blender in the background, and parse the output.

fileName = "coinflip.py";
file = OpenWrite[fileName];
WriteString[file, createCoinFlip[3.95, 1, 0.1, 0.1, 0.25, 0.75]];
Close[file];
outputfile = CreateFile[];
Run["blender --background --python coinflip.py >>" <> outputfile];
stext = OpenRead[outputfile];
Find[stext, "Coin"]
Close[stext];
DeleteFile[outputfile]
(* Coin is tails *)

You can visualize the results of the simulation by removing the "--background" and repeating the step above.

fileName = "coinflip.py";
file = OpenWrite[fileName];
WriteString[file, createCoinFlip[3.45, 1, 0.1, 0.1, 0.25, 0.75]];
Close[file];
outputfile = CreateFile[];
Run["blender --python coinflip.py >>" <> outputfile];
stext = OpenRead[outputfile];
Find[stext, "Coin"]
Close[stext];
DeleteFile[outputfile]

If you left click anywhere on the screen and hit the play button, then you should see the following:

Rendering in Blender

You can take advantage of Blender's photorealistic rendering capability to create a nice animation if desired.

Update: Additional Post Processing in Mathematica

Blender is geared more towards the artist whereas Mathematica is geared more towards the physicist. We can find synergy when we combine the strengths of both tools.

What follows is a simple example of how to conduct additional post-processing on a Blender simulation in Mathematica.

First, let's modify the python generation script to give the positions and orientations of the coin at each frame (we will insert a string "PosRot" to identify the appropriate lines).

createCoinFlipTransform[z_, rx_, ry_, rz_, friction_, restitution_] :=
  StringTemplate["import bpy
from math import pi

# Read current Blender version
version = bpy.app.version

for o in bpy.data.objects:
    if version < (2, 80, 0):
        if o.type == 'MESH' or o.type == 'EMPTY':
            o.select = True
        else:
            o.select = False
    else:
        o.select_set(o.type == 'EMPTY' or o.type == 'MESH')

# Delete all objects in the scene
bpy.ops.object.delete()

# Add the floor
if version < (2, 80, 0):
    bpy.ops.mesh.primitive_cube_add(radius=5, location=(0, 0, 0))
else:
    bpy.ops.mesh.primitive_cube_add(size=10, location=(0, 0, 0))
bpy.ops.transform.resize(value=(1, 1, 0.1))
bpy.ops.rigidbody.objects_add(type='PASSIVE')
boxObj = bpy.context.active_object
boxObj.rigid_body.collision_shape = "BOX"
boxObj.name = "Ground"

# Add the Coin
bpy.ops.mesh.primitive_cylinder_add(radius=1, depth=0.1, 
location=(0, 0, 3))
bpy.ops.rigidbody.objects_add(type='ACTIVE')
cylObj = bpy.context.active_object
cylObj.rigid_body.collision_shape = "CYLINDER"
bpy.context.object.rigid_body.friction = `friction`
bpy.context.object.rigid_body.restitution = `restitution`
cylObj.name = "Coin"
# Set reference to the coin
coin = bpy.data.objects["Coin"]

# Set a reference to the scene
sce = bpy.context.scene
# Set first frame
sce.frame_set(1)
# Set Keyframes
coin.keyframe_insert(data_path="location")
coin.keyframe_insert(data_path="rotation_euler")
bpy.context.object.rigid_body.kinematic = True
bpy.context.object.keyframe_insert('rigid_body.kinematic')

# Advance two frames and add translational and rotational motion
sce.frame_set(4)
# Translate up a little
coin.location.z = `z`
# Rotate coin predominantly around the x-axis
coin.rotation_euler.x = `rx`
coin.rotation_euler.y = `ry`
coin.rotation_euler.z = `rz`
# Set Keyframes
coin.keyframe_insert(data_path="location")
coin.keyframe_insert(data_path="rotation_euler")
bpy.context.object.rigid_body.kinematic = False
bpy.context.object.keyframe_insert('rigid_body.kinematic')

# Set frame to the end
sce.frame_set(250)

# Bake rigid body simulation
override = {'scene': bpy.context.scene,
            'point_cache': 
bpy.context.scene.rigidbody_world.point_cache}
# bake to current frame
bpy.ops.ptcache.bake(override, bake=False)

# Get transformations
tr = coin.matrix_world.translation
eu = coin.matrix_world.to_euler()

for i in range(250):
    sce.frame_set(i)
    tr = coin.matrix_world.translation
    eu = coin.matrix_world.to_euler()
    print("PosRot",tr.x, tr.y, tr.z, eu.x , eu.y , eu.z )
"][<|"z" -> z, "rx" -> rx, "ry" -> ry, "rz" -> rz, 
"friction" -> friction, "restitution" -> restitution|>]

We can extract the positions and orientations of the simulation with the following code.

fileName = "coinflip.py";
file = OpenWrite[fileName];
WriteString[file, createCoinFlipTransform[4, -Pi 0.75, 0.1, 0.1, 0.25, 0.75]];
Close[file];
outputfile = CreateFile[];
Run["blender --background --python coinflip.py >>" <> outputfile];
stext = OpenRead[outputfile];
data = ToExpression@StringSplit[#] & /@ FindList[stext, "PosRot"];
{tx, ty, tz, rx, ry, rz} = Transpose@data[[All, {2, 3, 4, 5, 6, 7}]];
Close[stext];
DeleteFile[outputfile]

We can define a cuboid and cylinder that have the same dimensions as the Blender simulation and we can create a transformation function with the following code.

box = {Cuboid[{-5, -5, -0.5}, {5, 5, 0.5}]};
cyl = {Cylinder[{{0, 0, -0.05}, {0, 0, 0.05}}, 1], 
   AbsolutePointSize[10], 
   Opacity[1], {Black, Point[{0, 0, 0}]}, {Red, 
    Point[{1, 0, 0}]}, {Green, Point[{0, 1, 0}]}, {Blue, 
    Point[{0, 0, 1}]}};
m = IdentityMatrix[4];
m[[1 ;; 3, 1 ;; 3]] = EulerMatrix[{a, b, c}, {1, 2, 3}];
m[[1 ;; 3, -1]] = {x, y, z};
transform[a_, b_, c_, x_, y_, z_] = TransformationFunction[m];

Now, we can combine plots of position and orientation (or other quantities like angular momentum) into a Manipulate[] function.

Manipulate[
 Column[{Row[{ListPlot[{tx[[1 ;; i]], ty[[1 ;; i]], tz[[1 ;; i]]}, 
      Filling -> Axis, ImageSize -> {200, 200}, PlotRange -> All, 
      PlotLegends -> {"tx", "ty", "tz"}],
     ListPlot[{rx[[1 ;; i]], ry[[1 ;; i]], rz[[1 ;; i]]}, 
      Filling -> Axis, ImageSize -> {200, 200}, PlotRange -> All, 
      PlotLegends -> {"rx", "ry", "rz"}]}], 
   Graphics3D[{{Opacity[0.75], Red, box}, 
     GeometricTransformation[{Opacity[.85], Yellow, cyl}, 
      transform[rx[[i]], ry[[i]], rz[[i]], tx[[i]], ty[[i]], 
       tz[[i]]]]}, SphericalRegion -> True,  Boxed -> False, 
    ImageSize -> {400, 400}]}], {i, 1, 250, 1}]