rotationMatrixFromGravity

CMRotationMatrix rotationMatrixFromGravity(float x, float y, float z)
{	
	// The Z axis of our rotated frame is opposite gravity
	vec3f_t zAxis = vec3f_normalize(vec3f_init(-x, -y, -z));
	
	// The Y axis of our rotated frame is an arbitrary vector perpendicular to gravity
	// Note that this convention will have problems as zAxis.x approaches +/-1 since the magnitude of
	// [0, zAxis.z, -zAxis.y] will approach 0
	vec3f_t yAxis = vec3f_normalize(vec3f_init(0, zAxis.z, -zAxis.y));
	
	// The X axis is just the cross product of Y and Z
	vec3f_t xAxis = vec3f_crossProduct(yAxis, zAxis);
	
	// each array is a row
	CMRotationMatrix mat = {
		[ xAxis.x, xAxis.y, xAxis.z ],
		[ yAxis.x, yAxis.y, yAxis.z ],
		[ zAxis.x, zAxis.y, zAxis.z ] };

	return mat;
}

Here's some octave code that I used to prove to myself that it works:

#1. Use an acceleration reading to compute the rotation matrix per the gist
#    The value of g was obtained from my accelerometer at rest.

octave:33> g = [1.010040589617603, 0.0065919980468153935, 0.09399700918607135]
g =

   1.0100406   0.0065920   0.0939970

octave:34> rawZ = -g
rawZ =

  -1.0100406  -0.0065920  -0.0939970

octave:35> zAxis = rawZ/norm(rawZ)
zAxis =

  -0.9956766  -0.0064983  -0.0926603

octave:36> rawY = [0, zAxis(3), -zAxis(2)]
rawY =

   0.000000  -0.092660   0.006498

octave:37> yAxis = rawY / norm(rawY)
yAxis =

   0.00000  -0.99755   0.06996

octave:38> xAxis = cross(yAxis, zAxis)
xAxis =

   0.092888  -0.069656  -0.993237

octave:39> R = [xAxis; yAxis; zAxis]
R =

   0.09289   0.00000  -0.99568
  -0.06966  -0.99755  -0.00650
  -0.99324   0.06996  -0.09266

#2. Compute the product of the rotation matrix with the original acceleration reading, and you should get around [0, 0, 1].  
#    Or in english: If your accelerometer reading is is the same as what you've said means "earth reference" then 
#    The rotated value of that reading should look like a object laying flat and parallel to the ground.

octave:40> R*g'  # g' because we need a column vector here.
ans =

   2.2987e-04
  -7.7542e-02
  -1.0115e+00

Conclusion: It works!

function R = rotMatrix(g)
  rawZ  = -g;
  zAxis = rawZ/norm(rawZ);
  rawY  = [0, zAxis(3), -zAxis(2)];
  yAxis = rawY / norm(rawY);
  xAxis = cross(yAxis, zAxis);
  R     = [xAxis; yAxis; zAxis];
endfunction

g = [1.0100406, 0.0065920, 0.0939970];
g4 = [0.70700, 0.70700, 0.001];

Applying a rotation. Useful if the sensor rotates in place without translating. Requires gyros and accelerometers:

octave> g1 = [ 1, 0, 0 ]
octave> r  = rotMatrix(g1)   # calculate the initial rotation matrix, converting device => N,E,D coordinates

octave> function r = d2r(deg)
  r = deg * pi / 180;
endfunction

octave> function R = degRotZMatrix(deg)
  R = [ cos(d2r(deg)), -sin(d2r(deg)), 0
        sin(d2r(deg)),  cos(d2r(deg)), 0
        0            ,  0            , 1 ];
endfunction

#In the next data-frame the device has rotated.  The device => N,E,D mapping needs updated
octave> g2 = [ 0.707, 0.707, 0 ]  # accelerometer readings after a 45deg clockwise rotation about Z as determined by gyros.
                                  # This should map to [ 0, 0, -1 ]

octave> newRotationMatrix = r * degRotZMatrix(-45)  # Clockwise is negative
octave> newRotationMatrix * g2'                     # This should eq N,E,D [ 0, 0 , -1 ]

ans =

   9.1881e-02
   6.9958e-05
  -9.9562e-01

Yay! it works!