def regular_square_bundle(num_rays, center, direction, width):
	"""
	Generate a ray bundles whose rays are equally spaced along a square grid,
	and all pointing in the same direction.
	
	Arguments:
	num_rays - number of rays to generate.
	center - a column 3-array with the 3D coordinate of the disk's center
	direction - a 1D 3-array with the unit direction vector for the bundle.
	width - of the square of starting points.
	
	Returns: 
	A RayBundle object with the above charachteristics set.
	"""
	rot = rotation_to_z(direction)
	directions = N.tile(direction[:,None], (1, num_rays))
	range = N.s_[-width:width:float(2*width)/N.sqrt(num_rays)]
	xs, ys = N.mgrid[range, range]
	vertices_local = N.array([xs.flatten(),  ys.flatten(),  N.zeros(len(xs.flatten()))])
	vertices_global = N.dot(rot,  vertices_local)

	rayb = RayBundle()
	rayb.set_vertices(vertices_global + center)
	rayb.set_directions(directions)
	return rayb

class TestTraceProtocol1(unittest.TestCase):
    """ 
    Tests intersect_ray and the bundle driver with a single flat surface, not rotated, with 
    a single interation 
    """
    def setUp(self):
        dir = N.array([[1,1,-1],[-1,1,-1],[-1,-1,-1],[1,-1,-1]]).T/math.sqrt(3)
        position = N.c_[[0,0,1],[1,-1,1],[1,1,1],[-1,1,1]]
        self._bund = RayBundle(position, dir, energy=N.ones(4))
        
        self.assembly = Assembly()
        object = AssembledObject()
        object.add_surface(Surface(FlatGeometryManager(), opt.perfect_mirror))
        self.assembly.add_object(object)
        self.engine = TracerEngine(self.assembly)
        
    def test_intersect_ray1(self):
        surfaces = self.assembly.get_surfaces()
        objects = self.assembly.get_objects()
        surfs_per_obj = [len(obj.get_surfaces()) for obj in objects]
        surf_ownership = N.repeat(N.arange(len(objects)), surfs_per_obj)
        ray_ownership = -1*N.ones(self._bund.get_num_rays())
        surfs_relevancy = N.ones((len(surfaces), self._bund.get_num_rays()), dtype=N.bool)
        
        params = self.engine.intersect_ray(self._bund, surfaces, objects, \
            surf_ownership, ray_ownership, surfs_relevancy)[0]
        self.failUnless(params.all())

    def test_ray_tracer(self):
        """Ray tracer after one iteration returns what the surface would have"""
        params = self.engine.ray_tracer(self._bund,1,.05)[0]
        correct_pts = N.zeros((3,4))
        correct_pts[:2,0] = 1
        
        N.testing.assert_array_almost_equal(params, correct_pts)

def edge_rays_bundle(num_rays,  center,  direction,  radius, ang_range, flux=None, radius_in=0.):

	radius = float(radius)
	radius_in = float(radius_in)
	a = edge_rays_directions(num_rays, ang_range)
		
	# Rotate to a frame in which <direction> is Z:
	perp_rot = rotation_to_z(direction)
	directions = N.sum(perp_rot[...,None] * a[None,...], axis=1)
	# Locations:
	# See [1]
	xi1 = random.uniform(size=num_rays)
	thetas = random.uniform(high=2.*N.pi, size=num_rays)
	rs = N.sqrt(radius_in**2.+xi1*(radius**2.-radius_in**2.))
	xs = rs * N.cos(thetas)
	ys = rs * N.sin(thetas)

	# Rotate locations to the plane defined by <direction>:
	vertices_local = N.vstack((xs, ys, N.zeros(num_rays)))
	vertices_global = N.dot(perp_rot, vertices_local)

	rayb = RayBundle(vertices=vertices_global + center, directions=directions)
	if flux != None:
		rayb.set_energy(N.pi*(radius**2.-radius_in**2.)/num_rays*flux*N.ones(num_rays))
	return rayb

class TestTraceProtocol2(unittest.TestCase):
    """
    Tests intersect_ray with a flat surface rotated around the x axis 45 degrees
    """
    def setUp(self):
        ns = -1/N.sqrt(2)
        dir = N.c_[[0,0,1],[0,0,-1],[0,ns,ns]]
        position = N.c_[[0,0,1],[0,1,2],[0,0,1]]
        self._bund = RayBundle(position, dir, energy=N.ones(3))
        
    def test_intersect_ray2(self):
        rot = general_axis_rotation([1,0,0],N.pi/4)
        surface = Surface(FlatGeometryManager(), opt.perfect_mirror, rotation=rot)
        assembly = Assembly()
        object = AssembledObject()
        object.add_surface(surface)
        assembly.add_object(object)
        
        engine = TracerEngine(assembly)
        surfaces = assembly.get_surfaces()
        objects = assembly.get_objects()
        surfs_per_obj = [len(obj.get_surfaces()) for obj in objects]
        surf_ownership = N.repeat(N.arange(len(objects)), surfs_per_obj)
        ray_ownership = -1*N.ones(self._bund.get_num_rays())
        surfs_relevancy = N.ones((len(surfaces), self._bund.get_num_rays()), dtype=N.bool)
        
        params = engine.intersect_ray(self._bund, surfaces, objects, \
            surf_ownership, ray_ownership, surfs_relevancy)[0]
        correct_params = N.array([[False, True, False]])

        N.testing.assert_array_almost_equal(params, correct_params)

class ExampleScene(TracerScene):
    source_y = t_api.Range(0., 5., 2.)
    source_z = t_api.Range(0., 5., 1.)

    def __init__(self):
        # The energy bundle we'll use for now:
        nrm = 1/(N.sqrt(2))
        direct = N.c_[[0,-nrm, nrm],[0,0,-1]]
        position = N.tile(N.c_[[0, self.source_y, self.source_z]], (1, 2))
        self.bund = RayBundle(vertices=position, directions=direct, energy=N.r_[1, 1])

        # The assembly for ray tracing:
        rot1 = N.dot(G.rotx(N.pi/4)[:3,:3], G.roty(N.pi)[:3,:3])
        surf1 = rect_one_sided_mirror(width=10, height=10)
        surf1.set_rotation(rot1)
        surf2 = rect_one_sided_mirror(width=10, height=10)
        self.assembly = Assembly(objects=[surf1, surf2])

        TracerScene.__init__(self, self.assembly, self.bund)

    @t_api.on_trait_change('_scene.activated')
    def initialize_camere(self):
        self._scene.mlab.view(0, -90)
        self._scene.mlab.roll(0)

    @t_api.on_trait_change('source_y, source_z')
    def bundle_move(self):
        position = N.tile(N.c_[[0, self.source_y, self.source_z]], (1, 2))
        self.bund.set_vertices(position)
        self.plot_ray_trace()

    view = tui.View(
        tui.Item('_scene', editor=SceneEditor(scene_class=MayaviScene),
            height=400, width=300, show_label=False),
        tui.HGroup('-', 'source_y', 'source_z'))

class TestInterface(unittest.TestCase):
    def setUp(self):
        self.num_rays = 10
        dir = N.tile(N.c_[[0, 0, -1]], (1, self.num_rays))
        theta = N.linspace(0, 2*N.pi, self.num_rays, endpoint=False)
        position = N.vstack((N.cos(theta), N.sin(theta), N.ones(self.num_rays)))
        
        self._bund = RayBundle(position, dir)
        
        # The boundary is positioned to create a bottom hemisphere.
        boundary = BoundarySphere(radius=4., location=N.r_[0., 0., -4*N.sqrt(3)/2.])
        self.gm = CutSphereGM(2., boundary)
        self.prm = self.gm.find_intersections(N.eye(4), self._bund)
    
    def test_find_intersections(self):
        """The correct parametric locations are found for cut sphere geometry"""
        self.failUnlessEqual(self.prm.shape, (self.num_rays,))
        N.testing.assert_array_almost_equal(self.prm, 1 + 2*N.sin(N.pi/3))
    
    def test_get_normals(self):
        """Cut sphere surface returns center-pointing normals"""
        self.gm.select_rays(N.arange(self.num_rays))
        n = self.gm.get_normals()
        N.testing.assert_array_almost_equal(n[-1,0], n[-1,1:])
        N.testing.assert_array_almost_equal(self._bund.get_vertices()[:2],
            -n[:2]/N.sqrt((n[:2]**2).sum(axis=0)))
    
    def test_inters_points_global(self):
        """Cut sphere returns correct intersections"""
        self.gm.select_rays(N.arange(self.num_rays))
        pts = self.gm.get_intersection_points_global()
        N.testing.assert_array_equal(pts[:2], self._bund.get_vertices()[:2])
        N.testing.assert_array_almost_equal(pts[2], -2*N.sin(N.pi/3))

class TestParabolicDish(unittest.TestCase):
    def setUp(self):
        pos = N.zeros((3,4))
        pos[0] = N.r_[0, 0.5, 2, -2]
        pos[2] = 2.
        dir = N.tile(N.c_[[0,0,-1]], (1,4))

        self.bund = RayBundle()
        self.bund.set_vertices(pos)
        self.bund.set_directions(dir)

        self.surf = Surface(ParabolicDishGM(2., 1.), opt.perfect_mirror)
    
    def test_selection_at_origin(self):
        """Simple dish rejects missing rays"""
        misses = N.isinf(self.surf.register_incoming(self.bund))
        N.testing.assert_array_equal(misses, N.r_[False, False, True, True])
    
    def test_transformed(self):
        """Translated and rotated dish rejects missing rays"""
        trans = generate_transform(N.r_[1., 0., 0.], N.pi/4., N.c_[[0., 0., 1.]])
        
        self.surf.transform_frame(trans)
        misses = N.isinf(self.surf.register_incoming(self.bund))
        N.testing.assert_array_equal(misses, N.r_[False, False, True, True])
    
    def test_mesh(self):
        """Parabolic dish mesh looks OK"""
        p = ParabolicDishGM(5, 3)
        x, y, z = p.mesh(5)
        
        N.testing.assert_array_almost_equal(z, p.a*(x**2 + y**2))
        self.failIf(N.any(x**2 + y**2 > 6.25))

 def runTest(self):
     pos = N.array([[0, 1.5], [0, -1.5], [1, 0], [-1, 0], [0.1, 0.1], [-0.1, 0.6]])
     bund = RayBundle()
     bund.set_vertices(N.vstack((pos.T, N.ones(pos.shape[0]))))
     bund.set_directions(N.tile(N.c_[[0,0,-1]], (1,6)))
     surf = Surface(HexagonalParabolicDishGM(2., 1.), opt.perfect_mirror)
     
     misses = N.isinf(surf.register_incoming(bund))
     N.testing.assert_array_equal(misses, N.r_[True, True, True, True, False, False])

    def setUp(self):
        dir = N.c_[[0., 0, -1], [0, 1, -1], [0, 11, -2], [0, 1, 0]]
        dir /= N.sqrt(N.sum(dir**2, axis=0))
        position = N.c_[[0., 0, 1], [0, -1, 1], [0, -11, 2], [0, 1, 1]]

        bund = RayBundle()
        bund.set_vertices(position)
        bund.set_directions(dir)
        self.bund = bund
        self.correct = N.r_[1., N.sqrt(2), N.sqrt(11**2 + 4)]

    def test_up_down(self):
        """Rays coming from below are absorbed, from above reflected"""
        going_down = N.c_[[1, 1, -1], [-1, 1, -1], [-1, -1, -1], [1, -1, -1]] / N.sqrt(3)
        going_up = going_down.copy()
        going_up[2] = 1 / N.sqrt(3)
        
        pos_up = N.c_[[0,0,1], [1,-1,1], [1,1,1], [-1,1,1]]
        pos_down = pos_up.copy()
        pos_down[2] = -1

        bund = RayBundle()
        bund.set_directions(N.hstack((going_down, going_up)))
        bund.set_vertices(N.hstack((pos_up, pos_down)))
        bund.set_energy(N.tile(100, 8))
        bund.set_ref_index(N.tile(1, 8))
        
        gm = FlatGeometryManager()
        prm = gm.find_intersections(N.eye(4), bund)
        absref = optics_callables.AbsorberReflector(0.)
        selector = N.arange(8)
        gm.select_rays(selector)
        outg = absref(gm, bund, selector)
        
        e = outg.get_energy()
        N.testing.assert_array_equal(e[:4], 100)
        N.testing.assert_array_equal(e[4:], 0)

class TestHomogenizer(unittest.TestCase):
    def setUp(self):
        """A homogenizer transforms a bundle correctly"""
        hmg = rect_homogenizer(5., 3., 10., 0.9)
        self.engine = TracerEngine(hmg)
        self.bund = RayBundle()
        
        # 4 rays starting somewhat above (+z) the homogenizer
        pos = N.zeros((3,4))
        pos[2] = N.r_[11, 11, 11, 11]
        self.bund.set_vertices(pos)
        
        # One ray going to each wall:
        dir = N.c_[[1, 0, -1], [-1, 0, -1], [0, 1, -1], [0, -1, -1]]/N.sqrt(2)
        self.bund.set_directions(dir)
        
        # Laborious setup details:
        self.bund.set_energy(N.ones(4)*4.)
        self.bund.set_ref_index(N.ones(4))
    
    def test_first_hits(self):
        """Test bundle enters homogenizer correctly"""
        v, d = self.engine.ray_tracer(self.bund, 1, 0.05)
        
        out_dirs = N.c_[[-1, 0, -1], [1, 0, -1], [0, -1, -1], [0, 1, -1]]/N.sqrt(2)
        N.testing.assert_array_almost_equal(d, out_dirs)
        
        out_hits = N.c_[
            [2.5, 0, 8.5], 
            [-2.5, 0, 8.5], 
            [0, 1.5, 9.5], 
            [0, -1.5, 9.5]]
        N.testing.assert_array_almost_equal(v, out_hits)

def solar_disk_bundle(num_rays,  center,  direction,  radius, ang_range, flux=None, radius_in=0., angular_span=[0.,2.*N.pi], procs=1):
    """
    Generates a ray bundle emanating from a disk, with each surface element of 
    the disk having the same ray density. The rays all point at directions uniformly 
    distributed between a given angle range from a given direction.
    Setting of the bundle's energy is left to the caller.
    
    Arguments:
    num_rays - number of rays to generate.
    center - a column 3-array with the 3D coordinate of the disk's center
    direction - a 1D 3-array with the unit average direction vector for the
        bundle.
    radius - of the disk.
    ang_range - in radians, the maximum deviation from <direction>.
    flux - if not None, the ray bundle's energy is set such that each ray has
        an equal amount of energy, and the total energy is flux*pi*radius**2
    radius_in - Inner radius if the disc is pierced
    angular_span - wedge of the disc to consider
    
    Returns: 
    A RayBundle object with the above characteristics set.
    """

	# FIXME why should 'center' be a column vector... that's just annoying.

    radius = float(radius)
    radius_in = float(radius_in)
    a = pillbox_sunshape_directions(num_rays, ang_range)
        
    # Rotate to a frame in which <direction> is Z:
    perp_rot = rotation_to_z(direction)
    directions = N.sum(perp_rot[...,None] * a[None,...], axis=1)
    # Locations:
    # See [1]
    xi1 = random.uniform(size=num_rays)
    thetas = random.uniform(low=angular_span[0], high=angular_span[1], size=num_rays)
    rs = N.sqrt(radius_in**2.+xi1*(radius**2.-radius_in**2.))
    xs = rs * N.cos(thetas)
    ys = rs * N.sin(thetas)

    # Rotate locations to the plane defined by <direction>:
    vertices_local = N.vstack((xs, ys, N.zeros(num_rays)))
    vertices_global = N.dot(perp_rot, vertices_local)

    rayb = RayBundle(vertices=vertices_global + center, directions=directions)
    if flux != None:
        rayb.set_energy(N.pi*(radius**2.-radius_in**2.)/num_rays*flux*N.ones(num_rays))
    else:
        rayb.set_energy(N.ones(num_rays)/num_rays/procs)
    return rayb

class TestTraceProtocol6(unittest.TestCase):
    """
    Tests a spherical surface
    """
    def setUp(self):
        surface1 = Surface(HemisphereGM(2.), opt.perfect_mirror,
            rotation=general_axis_rotation(N.r_[1,0,0], N.pi/2.))
        surface2 = Surface(HemisphereGM(2.), opt.perfect_mirror, 
            location=N.array([0,-2,0]), 
            rotation=general_axis_rotation(N.r_[1,0,0], -N.pi/2.))
        
        self._bund = RayBundle()
        self._bund.set_directions(N.c_[[0,1,0]])
        self._bund.set_vertices(N.c_[[0,-1,0]])
        self._bund.set_energy(N.r_[[1]]) 
        self._bund.set_ref_index(N.r_[[1]])

        assembly = Assembly()
        object1 = AssembledObject()
        object2 = AssembledObject()
        object1.add_surface(surface1)
        object2.add_surface(surface2)
        assembly.add_object(object1)
        assembly.add_object(object2)

        self.engine = TracerEngine(assembly)
        
    def test_ray_tracers1(self):
        params = self.engine.ray_tracer(self._bund, 1, .05)[0]
        correct_params = N.c_[[0,2,0]]

        N.testing.assert_array_almost_equal(params,correct_params)

class TestObjectBuilding1(unittest.TestCase):
    """Tests an object composed of sphere surfaces"""
    def setUp(self):
        self.assembly = Assembly()
        surface1 = Surface(HemisphereGM(3.), optics_callables.perfect_mirror, 
            location=N.array([0,0,-1.]),
            rotation=general_axis_rotation(N.r_[1,0,0], N.pi))
        surface2 = Surface(HemisphereGM(3.), optics_callables.perfect_mirror, 
            location=N.array([0,0,1.]))
        
        self.object = AssembledObject()
        self.object.add_surface(surface1)
        self.object.add_surface(surface2)
        self.assembly.add_object(self.object)

        dir = N.c_[[0,0,1.],[0,0,1.]]
        position = N.c_[[0,0,-3.],[0,0,-1.]]
    
        self._bund = RayBundle(position, dir, energy=N.ones(2))
    
    def test_object(self):
        """Tests that the assembly heirarchy works at a basic level"""
        self.engine = TracerEngine(self.assembly)

        inters = self.engine.ray_tracer(self._bund,1,.05)[0]
        correct_inters = N.c_[[0,0,2],[0,0,-2]]

        N.testing.assert_array_almost_equal(inters, correct_inters)
    
    def test_translation(self):
        """Tests an assembly that has been translated"""
        trans = N.array([[1,0,0,0],[0,1,0,0],[0,0,1,1],[0,0,0,1]])
        self.assembly.transform_children(trans)

        self.engine = TracerEngine(self.assembly)

        params =  self.engine.ray_tracer(self._bund,1,.05)[0]
        correct_params = N.c_[[0,0,3],[0,0,-1]]

        N.testing.assert_array_almost_equal(params, correct_params)

    def test_rotation_and_translation(self):
        """Tests an assembly that has been translated and rotated"""
        self._bund = RayBundle()
        self._bund.set_vertices(N.c_[[0,-5,1],[0,5,1]])
        self._bund.set_directions(N.c_[[0,1,0],[0,1,0]])
        self._bund.set_energy(N.r_[[1,1]])
        self._bund.set_ref_index(N.r_[[1,1]])

        trans = generate_transform(N.r_[[1,0,0]], N.pi/2, N.c_[[0,0,1]])
        self.assembly.transform_children(trans)

        self.engine = TracerEngine(self.assembly)

        params =  self.engine.ray_tracer(self._bund,1,.05)[0]
        correct_params = N.c_[[0,-2,1]]

        N.testing.assert_array_almost_equal(params, correct_params)

 def setUp(self):
     self.num_rays = 10
     dir = N.tile(N.c_[[0, 0, -1]], (1, self.num_rays))
     theta = N.linspace(0, 2*N.pi, self.num_rays, endpoint=False)
     position = N.vstack((N.cos(theta), N.sin(theta), N.ones(self.num_rays)))
     self._bund = RayBundle(position, dir)
     
     self.gm = HemisphereGM(radius=2.)
     self.prm = self.gm.find_intersections(N.eye(4), self._bund)

 def test_paraxial_ray(self):
     """A paraxial ray in reflected correctly"""
     bund = RayBundle()
     bund.set_vertices(N.c_[[0.01, 0., 2.]])
     bund.set_directions(N.c_[[0., 0., -1.]])
     bund.set_energy(N.r_[100.])
     bund.set_ref_index(N.r_[1])
     
     self.engine.ray_tracer(bund, 15, 10.)
     non_degenerate = self.engine.tree[-1].get_energy() > 10
     v = self.engine.tree[-1].get_vertices()[:,non_degenerate]
     d = self.engine.tree[-1].get_directions()[:,non_degenerate]
     # Not high equality demanded, because of spherical aberration.
     N.testing.assert_array_almost_equal(v, N.c_[[-0.01, 0., 1.5]], 2)
     N.testing.assert_array_almost_equal(d, N.c_[[0., 0., 1.]], 2)

 def setUp(self):
     dir = N.array([[1,1,-1],[-1,1,-1],[-1,-1,-1],[1,-1,-1]]).T/math.sqrt(3)
     position = N.c_[[0,0,1],[1,-1,1],[1,1,1],[-1,1,1]]
     self._bund = RayBundle(position, dir, energy=N.ones(4))
     
     self.assembly = Assembly()
     object = AssembledObject()
     object.add_surface(Surface(FlatGeometryManager(), opt.perfect_mirror))
     self.assembly.add_object(object)
     self.engine = TracerEngine(self.assembly)

def triangular_bundle(num_rays, A, AB, AC, direction, ang_range=N.pi/2., flux=None, procs=1):
	"""
	Triangular ray-casting surface anchored on the point A.
	Arguments:
	- num_rays: the number of rays 
	- A: The first summit of the triangle and its anchor point.
	- AB and AC the vertices of the sides of the triangle in its plane of reference.
	- direction: The direction at which the source is pointing
	- ang_range: the angular range of the rays emitted by the source

	Returns: 
	- A ray bundle object for tracing
	"""
	# Triangle ray vertices:
	# Declare random numbers:
	r1 = N.vstack(N.random.uniform(size=num_rays))
	r2 = N.vstack(N.random.uniform(size=num_rays))
	# Define points in a local referential where A is at [0,0] on a z=0 plane.
	sqrtr1 = N.sqrt(r1)
	Plocs = sqrtr1*(1.-r2)*AB+r2*sqrtr1*AC # Triangle point picking

	vertices_local = N.array([Plocs[:,0], Plocs[:,1], N.zeros(num_rays)])

	# Bring everything back to the global referential:
	rot = rotation_to_z(direction)
	vertices_global = N.dot(rot, vertices_local)+N.vstack(A)
	
	# Local referential directions:
	a = pillbox_sunshape_directions(num_rays, ang_range)
	# Rotate to a frame in which <direction> is Z:
	directions = N.sum(rot[...,None] * a[None,...], axis=1)

	rayb = RayBundle()

	rayb.set_vertices(vertices_global)
	rayb.set_directions(directions)

	l1 = N.sqrt(N.sum(AB**2))
	l2 = N.sqrt(N.sum(AC**2))
	l3 = N.sqrt(N.sum((-AB+AC)**2))
	s = (l1+l2+l3)/2.
	area = N.sqrt(s*(s-l1)*(s-l2)*(s-l3))
	if flux != None:
		rayb.set_energy(N.ones(num_rays)*flux*area/float(num_rays))
	else:
		rayb.set_energy(N.ones(num_rays)/float(num_rays)/procs)

	return rayb

def single_ray_source(position, direction, flux=None):
	'''
	Establishes a single ray source originating from a definned point on a defined exact 
	direction for the purpose of testing single ray behviours.

	Arguments:
	position - column 3-array with the ray's starting position.
	direction - a 1D 3-array with the unit average direction vector for the
				bundle.
	flux - if not None, the energy transported by the ray.

	Returns:
	A Raybundle object with the corresponding characteristics.
	'''
	directions = N.tile(direction[:,None],1)
	directions /= N.sqrt(N.sum(directions**2, axis=0))
	singray = RayBundle(vertices = position, directions = directions)
	singray.set_energy(flux*N.ones(1))
	return singray

class TestTraceProtocol3(unittest.TestCase):
    """
    Tests intersect_ray and the bundle driver with two rotated planes, with a single iteration
    """
    def setUp(self):
        self.x = 1/(math.sqrt(2))
        dir = N.c_[[0,self.x,-self.x],[0,1,0]]
        position = N.c_[[0,0,1],[0,0,1]]
        self._bund = RayBundle(position, dir, energy=N.ones(2))

        rot1 = general_axis_rotation([1,0,0],N.pi/4)
        rot2 = general_axis_rotation([1,0,0],N.pi/(-4))
        surf1 = Surface(FlatGeometryManager(), opt.perfect_mirror, rotation=rot1)
        surf2 = Surface(FlatGeometryManager(), opt.perfect_mirror, rotation=rot2)
        
        self.assembly = Assembly()
        object = AssembledObject()
        object.add_surface(surf1)
        object.add_surface(surf2)
        self.assembly.add_object(object)

        self.engine = TracerEngine(self.assembly)
        
    def test_intersect_ray(self):
        surfaces = self.assembly.get_surfaces()
        objects = self.assembly.get_objects()
        surfs_per_obj = [len(obj.get_surfaces()) for obj in objects]
        surf_ownership = N.repeat(N.arange(len(objects)), surfs_per_obj)
        ray_ownership = -1*N.ones(self._bund.get_num_rays())
        surfs_relevancy = N.ones((len(surfaces), self._bund.get_num_rays()), dtype=N.bool)
        
        params = self.engine.intersect_ray(self._bund, surfaces, objects, \
            surf_ownership, ray_ownership, surfs_relevancy)[0]
        correct_params = N.array([[True, True],[False, False]])

        N.testing.assert_array_almost_equal(params,correct_params)    

    def test_ray_tracer1(self):
        params = self.engine.ray_tracer(self._bund, 1,.05)[0]
        correct_params = N.c_[[0,.5,.5],[0,1,1]]

        N.testing.assert_array_almost_equal(params,correct_params)