Is there a better way than dupliverts to make these scripted monkeys always face "out"?
Instead of dupliverts, we can create it all in one mesh using bmesh.ops.create_monkey().

import bpy
import bmesh
import mathutils

loco_f, loco_v = (-5, 0, 5), (+5, 0, 5)

bpy.ops.object.select_all(action='SELECT')
bpy.ops.object.delete(use_global=False)

bpy.ops.mesh.primitive_ico_sphere_add(location=loco_v, size=4, subdivisions=2)
shape_mesh_v = bpy.context.active_object

bpy.ops.object.mode_set(mode='EDIT')
bpy.ops.mesh.select_all(action='SELECT')
bpy.ops.mesh.delete(type='EDGE_FACE')
bpy.ops.object.mode_set(mode='OBJECT')

bm = bmesh.new()
bm.from_mesh(shape_mesh_v.data)

src = [v for v in bm.verts]

for v in src:
    loc = mathutils.Matrix.Translation(v.co)
    rot = v.normal.to_track_quat('-Y','Z')
    mat = mathutils.Matrix.Identity(3)
    mat.rotate(rot)
    mat = loc * mat.to_4x4()
    bmesh.ops.create_monkey(bm, matrix=mat)

bm.to_mesh(shape_mesh_v.data)
bm.free()


So the key point you are looking for is v.normal.to_track_quat('-Y','Z') to get a quaternion rotation matrix from the vertex normal.

For another way to approach this, you could try animation nodes.



Here we make a list of the location and normal of each vert in an icosphere, we use that list to replicate the suzanne mesh.



EDIT:
After some experimenting with the aim of creating your final goal, I found creating one orbital plane at a time as a unique object seems to be the fastest way, seems to find a balance between object and vertex count. Once they are all created they can all be merged into one mesh for simplicity. While creating each plane, an armature bone can be added and aligned with each orbit allowing easy animation. This is the script I came up with, a simple proof of concept that needs better positioning and timing to be complete.