Whip's Inertia Tensor

/// <summary>
/// This computes an inertial tensor relative to a reference block
/// Help provided by Equinox :)
/// </summary>
/// <param name="reference">Reference block to base tensor off of</param>
/// <returns>Inertia tensor of the grid</returns>
MatrixD CalculateInertiaTensor(IMyShipController reference)
{
    MatrixD inertiaTensor = new MatrixD();
    var shipMass = reference.CalculateShipMass().PhysicalMass;

    var localRotation = new Matrix();
    double gridSize = reference.CubeGrid.GridSize; //get size of blocks on grid
    reference.Orientation.GetMatrix(out localRotation);
    localRotation.Translation = gridSize * reference.Position;
    var localInverse = Matrix.Invert(localRotation);
    var centerOfMass = Vector3D.Transform(reference.CenterOfMass, MatrixD.Invert(reference.WorldMatrix)); //convert COM to block relative
    var gridMinimum = Vector3D.Transform(reference.CubeGrid.Min * gridSize, localInverse); //get minimum grid vector
    var gridMaximum = Vector3D.Transform(reference.CubeGrid.Max * gridSize, localInverse); //get maximum grid vector

    var boundingBox = new BoundingBoxD(gridMinimum, gridMaximum); //Get bounding box of this ship
    var corners = boundingBox.GetCorners(); //get corners of bounding box

    foreach (var thisCorner in corners)
    {
        var diagonal = centerOfMass - thisCorner;
        var offset =  diagonal / 2;
        var x = diagonal.X;
        var y = diagonal.Y;
        var z = diagonal.Z;
        var localMass = shipMass / 8;
        var I_xx = localMass / 12 * (y * y + z * z);
        var I_yy = localMass / 12 * (x * x + z * z);
        var I_zz = localMass / 12 * (y * y + x * x);
        var localTensor = new MatrixD(I_xx, 0, 0,
                                      0, I_yy, 0,
                                      0, 0, I_zz);
        inertiaTensor += localTensor + localMass * (offset.LengthSquared() * MatrixD.Identity - OuterProduct(offset, offset));
    }
    return inertiaTensor;
}