class TracerEngine():
    """
    Tracer Engine implements that actual ray tracing. It keeps track of the number
    of objects, and determines which rays intersected which object.
    """
    def __init__(self, parent_assembly):
        """
        Arguments:
        parent_assembly - the highest level assembly
        
        Attributes:
        _asm - the Assembly instance containing the model to trace through.
        tree - a list  used for track parent rays and child rays. Each element
            of the list is a ray bundle created after one iteration of the 
            tracer. Each bundle contains an array listing the parent ray in the
            previous bundle (see ray_bundle.py). When a ray branches, the child
            rays point back to the same index representing that same parent ray.
            Otherwise, the index of each ray points to the ray in the previous
            branch.
        """
        self._asm = parent_assembly
        
    def intersect_ray(self, bundle, surfaces, objects, surf_ownership, \
        ray_ownership, surf_relevancy):
        """
        Finds the first surface intersected by each ray.
        
        Arguments:
        bundle - the RayBundle instance holding incoming rays.
        ownership - an array with the owning object instance for each ray in the
            bundle, or -1 for no ownership.
        
        Returns:
        stack - an s by r boolean array for s surfaces and r rays, stating
            for each surface i=1..s if it is intersected by ray j=1..r
        owned_rays - same size as stack, stating whether ray j was tested at all
            by surface i
        """
        ret_shape = (len(surfaces), bundle.get_num_rays())
        stack = N.zeros(ret_shape)
        owned_rays = N.empty(ret_shape, dtype=N.bool)
        
        # Bounce rays off each object
        for surf_num in xrange(len(surfaces)):
            owned_rays[surf_num] = ((ray_ownership == -1) | \
                (ray_ownership == surf_ownership[surf_num])) & surf_relevancy[surf_num]
            if not owned_rays[surf_num].any():
                continue
            if (~owned_rays[surf_num]).any():
                in_rays = bundle.inherit(owned_rays[surf_num])
            else:
                in_rays = bundle
            stack[surf_num, owned_rays[surf_num]] = \
                surfaces[surf_num].register_incoming(in_rays)
        
        # Raise an error if any of the parameters are negative
        if (stack < -1e-16).any():
            raise ValueError("Parameters must all be positive")
        
        # If parameter == 0, ray does not actually hit object, but originates from there; 
        # so it should be ignored in considering intersections 
        if (stack <= 1e-10).any():
            zeros = N.where(stack <= 1e-6)
            stack[zeros] = N.inf

        # Find the smallest parameter for each ray, and use that as the final one,
        # returns the indices.  If an entire column of the stack is N.inf (the ray misses
        # any surfaces), then take that column to be all False
        stack = ((stack == stack.min(axis=0)) & ~N.isinf(stack))
        
        return stack, owned_rays

    def ray_tracer(self, bundle, reps, min_energy, tree=True):
        """
        Creates a ray bundle or uses a reflected ray bundle, and intersects it
        with all objects, uses intersect_ray(). Based on the intersections,
        generates an outgoing ray in accordance with way the incoming ray
        reflects or refracts off any surfaces.
        
        Arguments:
        bundle - the initial incoming bundle
        reps - stop iteration after this many ray bundles were generated (i.e. 
            after the original rays intersected some surface this many times).
        min_energy - the minimum energy the rays have to have continue tracking
            them; rays with a lower energy are discarded. A float.
        tree - if True, register each bundle in self.tree, otherwise only
            register the last bundle.
        
        Returns: 
        A tuple containing an array of vertices and an array of the the direcitons
        of the last outgoing raybundle (note that the vertices of the new bundle are the 
        intersection points of the previous incoming bundle)
        
        NB: the order of the rays within the arrays may change, but they are tracked
        by the ray tree
        """
        self.tree = RayTree()
        bund = bundle
        if tree is True:
            self.tree.append(bund)
        
        # A list of surfaces and their matching objects:
        surfaces = self._asm.get_surfaces()
        objects = self._asm.get_objects()
        num_surfs = len(surfaces)
        
        surfs_per_obj = [len(obj.get_surfaces()) for obj in objects]
        surfs_until_obj = N.hstack((N.r_[0], N.add.accumulate(surfs_per_obj)))
        surf_ownership = N.repeat(N.arange(len(objects)), surfs_per_obj)
        ray_ownership = -1*N.ones(bund.get_num_rays())
        surfs_relevancy = N.ones((num_surfs, bund.get_num_rays()), dtype=N.bool)
        
        for i in xrange(reps):
            front_surf, owned_rays = self.intersect_ray(bund, surfaces, objects, \
                surf_ownership, ray_ownership, surfs_relevancy)
            outg = []
            record = []
            out_ray_own = []
            new_surfs_relevancy = []
            weak_ray_pos = []
            
            for surf_idx in xrange(num_surfs):
                inters = front_surf[surf_idx, owned_rays[surf_idx]]
                if not any(inters): 
                    surfaces[surf_idx].done()
                    continue
                surfaces[surf_idx].select_rays(N.nonzero(inters)[0])
                new_outg = surfaces[surf_idx].get_outgoing()
                new_record = new_outg
                
                # Fix parent indexing to refer to the full original bundle:
                parents = N.nonzero(owned_rays[surf_idx])[0][new_outg.get_parents()]
                new_outg.set_parents(parents)
        
                # Delete rays with negligible energies
                delete = new_outg.get_energy() <= min_energy
                weak_ray_pos.append(delete)
                if delete.any():
                    new_outg = new_outg.delete_rays(N.nonzero(delete)[0])
                surfaces[surf_idx].done()
                
                # Aggregate outgoing bundles from all the objects
                outg.append(new_outg)
                record.append(new_record)
                
                # Add new ray-ownership information to the total list:
                obj_idx = surf_ownership[surf_idx]
                surf_rel_idx = surf_idx - surfs_until_obj[obj_idx]
                object_owns_outg = objects[obj_idx].own_rays(new_outg, surf_rel_idx)
                out_ray_own.append(N.where(object_owns_outg, obj_idx, -1))
                
                # Add new surface-relevancy information, saying which surfaces
                # of the full list of surfaces must be checked next. This is
                # somewhat memory-intensize and requires optimization.
                surf_relev = N.ones((num_surfs, new_outg.get_num_rays()), dtype=N.bool)
                surf_relev[surf_ownership == obj_idx] = \
                    objects[obj_idx].surfaces_for_next_iteration(new_outg, surf_rel_idx)
                new_surfs_relevancy.append(surf_relev)
            
            bund = concatenate_rays(outg)
            if tree:
                # stores parent branch for purposes of ray tracking
                record = concatenate_rays(record)
                
                if record.get_num_rays() != 0:
                    weak_ray_pos = N.hstack(weak_ray_pos)
                    record = bund + record.inherit(N.nonzero(weak_ray_pos)[0])
                    self.tree.append(record)
                
            if bund.get_num_rays() == 0:
                # All rays escaping
                break
            
            ray_ownership = N.hstack(out_ray_own)
            surfs_relevancy = N.hstack(new_surfs_relevancy)
            
        if not tree:
            # Save only the last bundle. Don't bother moving weak rays to end.
            record = concatenate_rays(record)
            self.tree.append(record)
             
        return bund.get_vertices(), bund.get_directions() def general_axis_rotation(axis,  angle):
    """Generates a rotation matrix around <axis> by <angle>, using the right-hand
    rule.
    Arguments: 
        axis - a 3-component 1D array representing a unit vector
        angle - rotation counterclockwise in radians around the axis when the axis 
            points to the viewer.
    Returns: A 3x3 array representing the matrix of rotation.
    Reference: [1] p.47
    """
    s = sin(angle); c = cos(angle); v = 1 - c
    add = N.array([[0,          -axis[2], axis[1]],  
                            [axis[2],  0,          -axis[0]], 
                            [-axis[1], axis[0],  0        ] ])
    return N.multiply.outer(axis,  axis)*v + N.eye(3)*c + add*s def generate_transform(axis, angle, translation):
    """Generates a transformation matrix                                                      
    Arguments: axis - a 3-component 1D array giving the unit vector to rotate about          
    angle - angle of rotation counter clockwise in radians about the given axis in the 
    parent frame                         
    translation - a 3-component column vector giving the translation along the coordinates
    of the parent object/assembly
    """  
    rot = general_axis_rotation(axis, angle)
    return N.vstack((N.hstack((rot, translation)), N.r_[[0,0,0,1]])) class HemisphereGM(SphericalGM):
    """
    Trims the SphericalGM by only selecting intersection points in the lower
    hemisphere (z < 0).
    """
    def _select_coords(self, coords, prm):
        """
        Select from dual intersections by vetting out rays in the upper
        hemisphere, or if both are below use the default behaviour of choosing
        the first hit.
        """
        select = QuadricGM._select_coords(self, coords, prm) # defaults
        
        coords = N.concatenate((coords, N.ones((2,1,coords.shape[2]))), axis=1)
        local = N.sum(N.linalg.inv(self._working_frame)[None,:,:,None] * \
            coords[:,None,:,:], axis=2)
        bottom_hem = (local[:,2,:] <= 0) & (prm > 0)
        
        select[~N.logical_or(*bottom_hem)] = N.nan
        one_hit = N.logical_xor(*bottom_hem)
        select[one_hit] = N.nonzero(bottom_hem[:,one_hit])[0]
        
        return select
    
    def mesh(self, resolution):
        """
        Represent the surface as a mesh in local coordinates. Uses spherical-
        coordinates bins, i.e. the points are equally distributed by two angles,
        not by x,y.
        
        Arguments:
        resolution - in points per unit length (so the number of points 
            returned is O(A*resolution**3) for area A)
        
        Returns:
        x, y, z - each a 2D array holding in its (i,j) cell the x, y, and z
            coordinate (respectively) of point (i,j) in the mesh.
        """
        rng = N.pi + 1./(self._rad*resolution)
        theta, phi = N.ogrid[
            N.pi/2:rng:1./(self._rad*resolution),
            0:2*rng:1./(self._rad*resolution)]

        x = self._rad * N.sin(theta) * N.cos(phi)
        y = self._rad * N.sin(theta) * N.sin(phi)
        z = N.tile(self._rad * N.cos(theta), (1, x.shape[1]))
        
        return x, y, z class CutSphereGM(SphericalGM):
    """
    Trims the SphericalGM by only selecting intersection points inside a
    bounding volume. This GM still can't supply a mesh of itself.
    """
    def __init__(self, radius=1., bounding_volume=None):
        """
        Arguments:
        radius - of the sphere to cut.
        bunding_volume - an instance BoundaryShape subclass (see 
            ``boundary_shape.py``)
        """
        SphericalGM.__init__(self, radius)
        self._bound = bounding_volume
    
    def find_intersections(self, frame, ray_bundle):
        if self._bound is not None:
            self._bound.transform_frame(frame)
        return SphericalGM.find_intersections(self, frame, ray_bundle)
    
    def _select_coords(self, coords, prm):
        """
        Select from dual intersections by checking which of the  impact points
        is contained in a volume defined at object creation time.
        """
        select = SphericalGM._select_coords(self, coords, prm)
        if self._bound is None:
            return select
        
        #in_bd = self._bound.in_bounds(coords) & (prm > 0)
        in_bd = N.array([self._bound.in_bounds(coords[...,c]) for c in xrange(prm.shape[1])]).T
        in_bd &= prm > 0
        select[~N.logical_or(*in_bd)] = N.nan
        one_hit = N.logical_xor(*in_bd)
        select[one_hit] = N.nonzero(in_bd[:,one_hit])[0]
        
        return select class BoundarySphere(BoundaryShape):
    def __init__(self, location=None,  radius=1.):
        """
        Arguments:
        radius - a float
        location - a 3-component column vector representing the location of the
            sphere's center, passed to the base class.
        
        Attributes:
        _radius - radius of the bounding sphere
        """
        BoundaryShape.__init__(self, location, None)
        self._temp_loc = self._loc
        self._radius = radius
            
    def in_bounds(self, bund_vertices):
        """
        Returns a boolean array for whether or not a ray intersection was within the 
        bounding sphere
        Arguments: bund_vertices - an array of the vertices (see base class)
        """
        return self._radius**2 >= ((bund_vertices - self._temp_loc[:3])**2).sum(axis=1)

    def transform_frame(self, transform):
        """ 
        Transforms the center of the boundary shape into the global coordinates; this occurs
        when the assembly or object containing the surface is transformed
        """
        self._temp_loc = N.dot(transform, N.append(self._loc, N.c_[[1]]))
    
    def bounding_rect_for_plane(self, transform):
        """
        Find a rectangle on the xy plane of a given frame, which contains the
        intersection of the boundary shape and the plane.

        Arguments:
        transform - a 4x4 array, the homog. transf. matrix from the global
            coordinates to the frame whose xy plane intersects the boundary
            shape.

        Returns:
        xmin, xmax, ynin, ymax - of the rect, in the xy plane of the frame,
        """
        cent_proj = N.dot(N.linalg.inv(transform), N.append(self._temp_loc, 1))
        Reff = N.sqrt(self._radius**2 - (self._temp_loc[2] - cent_proj[2])**2)
        return cent_proj[0] - Reff, cent_proj[0] + Reff, \
            cent_proj[1] - Reff, cent_proj[1] + Reff class AssembledObject(Assembly):
    """ Defines an assembly of surfaces as an object. The object has its own set of 
    coordinates such that each surface composing the object can be described in terms of
    the object's coordinate system, and thus the user can rotate or translate the entire
    object together as one piece.
    The object also tracks refractive indices as a ray bundle leaves or enters a new
    material.
    """
    def __init__(self, surfs=None, bounds=None, transform=None):
        """
        Attributes:
        surfaces - a list of Surface objects
        boundaries - a list of Boundary objects that the surfaces are limited 
            by.
        transform - a 4x4 array representing the homogenous transformation 
            matrix of this object relative to the coordinate system of its 
            container
        """
        # Use the supplied values or some defaults:
        if surfs is None:
            self.surfaces = []
        else:
            self.surfaces = surfs
        
        if bounds is None:
            self.boundaries = []
        else:
            self.boundaries = bounds
        
        if transform is None:
            transform = N.eye(4)
        self.set_transform(transform)
    
    def get_surfaces(self):
        return self.surfaces

    def add_surface(self, surface):
        """Adds a surface to the object
        Arguments:  surface - a surface object
        """
        self.surfaces.append(surface)
        self.transform_children()

    def add_boundary(self, boundary):
        """Adds a boundary to the object. Surfaces not enclosed by the boundary
        sphere will not count as hit.
        Arguments: boundary - a spherical boundary objects
        """
        self.boundaries.append(boundary)
        self.transform_children()

    def get_boundaries(self):
        return self.boundaries
    
    def transform_children(self, assembly_transform=N.eye(4)):
        """Transforms an object if the assembly is transformed""" 
        const_t = self.get_transform()
        for member in self.surfaces + self.boundaries:
            member.transform_frame(N.dot(assembly_transform, const_t))
    
    def own_rays(self, rays, surface_id):
        """
        Decide which of the rays continue to propagate inside the object, so
        only the object's surfaces need be checked for intersection.
        This default implementation owns nothing.
        
        Arguments:
        rays - the RayBundle to check. 
        surface_id - the index of the surface which generated this bundle.
        
        Returns:
        a boolean array of length rays.get_num_rays() with False if not owned,
            True if owned.
        """
        return N.zeros(rays.get_num_rays(), dtype=N.bool)
    
    def surfaces_for_next_iteration(self, rays, surface_id):
        """
        Informs the ray tracer that some of the surfaces can be skipped in the
        next ireration for some of the rays.
        This default implementation marks all surfaces as relevant to all rays.
        
        Arguments:
        rays - the RayBundle to check. 
        surface_id - the index of the surface which generated this bundle.
        
        Returns:
        an array of size s by r for s surfaces in this object and r rays,
            stating whether ray i=1..r should be intersected with surface j=1..s
            in the next iteration.
        """
        return N.ones((len(self.surfaces), rays.get_num_rays()), dtype=N.bool) class Assembly(HasFrame):
    """
    Defines an assembly of objects or sub-assemblies.
    
    Attributes:
    _objects - a list of the objects the assembly contains
    _assemblies - a list of the sub assemblies the assembly contains
    """
    def __init__(self, objects=None, subassemblies=None, location=None, rotation=None):
        """
        Arguments:
        objects (optional) - a list of AssembledObject instances that are part
            of this assembly.
        subassemblies (optional) - a list of Assembly instances to be
            transformed together with this assembly.
        location, rotation - passed on to HasFrame.
        """
        if objects is None:
            objects = []
        self._objects = objects
        
        if subassemblies is None:
            subassemblies = []
        self._assemblies = subassemblies
        
        HasFrame.__init__(self, location, rotation)

    def get_local_objects(self):
        """
        Get the list of objects belonging directly to this assembly, without
        querying the child assemblies.
        """
        return self._objects

    def get_assemblies(self):
        return self._assemblies
    
    def get_objects(self):
        """
        Generates a list of AssembledObject instances belonging to this assembly
        or its subassemblies.
        """
        return reduce(operator.add, 
            [asm.get_objects() for asm in self._assemblies] + [self._objects])

    def get_surfaces(self):  
        """
        Generates a list of surface objects out of all the surfaces in the
        objects and subassemblies belonging to this assembly.
        
        The surfaces are guarantied to be in the order that each object returns
        them, and the objects are guarantied to be ordered the same as in 
        self.get_objects()
        """
        return reduce(operator.add, 
            [obj.get_surfaces() for obj in self.get_objects()])

    def add_object(self, object, transform=None):
        """
        Adds an object to the assembly.
        
        Arguments: 
        objects - the AssembledObject to add
        transform - the transformation matrix (as an array object) that describes 
            the object in the coordinate system of the Assembly
        """
        if transform == None:
            transform = N.eye(4)
        self._objects.append(object)
        object.set_transform(transform)
        self.transform_children()

    def add_assembly(self, assembly, transform=None):
        """Adds an assembly to the current assembly.
        
        Arguments:
        assembly - the assembly object to add
        transform - the transformation matrix (as an array object) that describes the 
            new assembly in the coordinate system of the current assembly
        """
        if transform == None:
            transform = N.eye(4)
        self._assemblies.append(assembly)
        assembly.set_transform(transform)
        self.transform_children()

    def set_rotation(self, rotation):
        """
        A recursive version of the parent's set_rotation. Changes the rotation
        part of the assembly's transform, and updates the assembly's children's
        transform accordingly.
        
        Arguments:
        rotation - a 3x3 rotation matrix.
        """
        HasFrame.set_rotation(self, rotation)
        self.transform_children()
    
    def set_location(self, location):
        """
        A recursive version of the parent's set_rotation. Changes the location
        part of the assembly's transform, and updates the assembly's children's
        transform accordingly.
        
        Arguments:
        location - a 3-component location vector.
        """
        HasFrame.set_location(self, location)
        self.transform_children()
    
    def set_transform(self, transform):
        HasFrame.set_transform(self, transform)
        self.transform_children()

    def transform_children(self, assembly_transform=N.eye(4)):
        """
        Transforms the entire assembly
        
        Arguments:
        assembly_transform - the transformation into the parent assembly containing the 
            current assembly
        """
        const_t = self.get_transform()
        for obj in self._assemblies + self._objects:
            obj.transform_children(N.dot(assembly_transform, const_t)) class Assembly(HasFrame):
    def __init__(self, objects=None, subassemblies=None, location=None, rotation=None):
    def get_local_objects(self):
    def get_assemblies(self):
    def get_objects(self):
    def get_surfaces(self):  
    def add_object(self, object, transform=None):
    def add_assembly(self, assembly, transform=None):
    def set_rotation(self, rotation):
    def set_location(self, location):
    def set_transform(self, transform):
    def transform_children(self, assembly_transform=N.eye(4)): class Paraboloid(QuadricGM):
    """Implements the geometry of a circular paraboloid surface"""
    def __init__(self, a=1., b=1.):
        """               
        Arguments: 
        a, b - describe the paraboloid as z = (x/a)**2 + (y/b)**2
            (sorry, legacy)
        
        Private attributes:                                                                  
        a, b - describe the paraboloid as z = a*x**2 + b*y**2
        """ 
        QuadricGM.__init__(self)
        self.a = 1./(a**2)
        self.b = 1./(b**2)

    def _normals(self, hits, directs):
        """
        Finds the normal to the parabola in a bunch of intersection points, by
        taking the derivative and rotating it. Used internally by quadric.
        
        Arguments:
        hits - the coordinates of intersections, as an n by 3 array.
        directs - directions of the corresponding rays, n by 3 array.
        """
        hit = N.dot(N.linalg.inv(self._working_frame), 
            N.vstack((hits.T, N.ones(hits.shape[0]))))
        dir_loc = N.dot(self._working_frame[:3,:3].T, directs.T)
        partial_x = 2*hit[0]*self.a
        partial_y = 2*hit[1]*self.b
        
        local_normal = N.vstack((-partial_x, -partial_y, N.ones_like(partial_x)))
        local_unit = local_normal/N.sqrt(N.sum(local_normal**2, axis=0))
        down = N.sum(dir_loc * local_unit, axis=0) > 0
        local_unit[:,down] *= -1
        normals = N.dot(self._working_frame[:3,:3], local_unit)
        
        return normals  
    
    def get_ABC(self, ray_bundle):
        """
        Determines the variables forming the relevant quadric equation, [1]
        """
        # Transform the the direction and position of the rays temporarily into the
        # frame of the paraboloid for calculations
        d = N.dot(self._working_frame[:3,:3].T, ray_bundle.get_directions())
        v = N.dot(N.linalg.inv(self._working_frame), 
            N.vstack((ray_bundle.get_vertices(), N.ones(d.shape[1]))))[:3]
        
        A = self.a*d[0]**2 + self.b*d[1]**2
        B = 2*self.a*d[0]*v[0] + 2*self.b*d[1]*v[1] - d[2] 
        C = self.a*v[0]**2 + self.b*v[1]**2 - v[2]
        
        return A, B, C class Surface(HasFrame):
    """
    Defines the base of surfaces that interact with rays.
    """
    def __init__(self, geometry, optics, location=None, rotation=None):
        """
        Arguments:
        geometry - a GeometryManager object responsible for finding ray 
            intersections with the surface.
        optics - a callable that gets the geometry manageri, bundle and
            selector, and returns the outgoing ray bundle generated by the
            geometry and bundle.
        location, rotation - passed directly to the HasFrame constructor.
        """
        HasFrame.__init__(self, location, rotation)
        self._geom = geometry
        self._opt = optics
        
    def get_optics_manager(self):
        """
        Returns the optics-manager callable. May be useful for introspection.
        Note that it is a read-only attribute.
        """
        return self._opt
    
    def get_geometry_manager(self):
        """
        Returns the geometry-manager instance. May be useful for introspection.
        Note that it is a read-only attribute.
        """
        return self._geom
    
    def register_incoming(self, ray_bundle):
        """
        Records the incoming ray bundle, and uses the geometry manager to
        return the parametric positions of intersection with the surface along
        the ray.
        
        Arguments:
        ray_bundle - a RayBundle object with at-least its vertices and
            directions specified.
        
        Returns
        A 1D array with the parametric position of intersection along each of
            the rays. Rays that missed the surface return +infinity.
        """
        self._current_bundle = ray_bundle
        return self._geom.find_intersections(self._temp_frame, ray_bundle)
    
    def select_rays(self, idxs):
        """
        Informs the geometry manager that only the specified rays are to be
        used henceforth.
        
        Arguments:
        idxs - an array with indices into the last registered ray bundle,
            marking rays that will be used.
        """
        self._selected = idxs
        self._geom.select_rays(idxs)
    
    def get_outgoing(self):
        """
        Generates a new ray bundle, which is the reflections/refractions of the
        user-selected rays out of the incoming ray-bundle that was previously
        registered.
        
        Returns: 
        a RayBundle object with the new bundle, with vertices on the surface
            and directions according to optics laws.
        """
        return self._opt(self._geom, self._current_bundle, self._selected)
    
    def done(self):
        """
        When this is called, the surface will no longer be queried on the
        results of the latest trace iteration, so it can discard internal
        data to relieve memory pressure.
        """
        if hasattr(self, '_current_bundle'):
            del self._current_bundle
        self._geom.done()
        
    def global_to_local(self, points):
        """
        Transform a set of points in the global coordinates back into the frame
        used during tracing.
        
        Arguments:
        points - a 3 x n array for n 3D points
        
        returns:
        local - a 3 x n array with the respective points in local coordinates.
        """
        return N.dot(N.linalg.inv(self._temp_frame), 
            N.vstack((points, N.ones(points.shape[1]))))
    
    def mesh(self, resolution):
        """
        Represent the surface as a mesh in global coordinates.
        
        Arguments:
        resolution - in points per unit length (so the number of points 
            returned is O(A*resolution**2) for area A)
        
        Returns:
        x, y, z - each a 2D array holding in its (i,j) cell the x, y, and z
            coordinate (respectively) of point (i,j) in the mesh.
        """
        # The geometry manager has the local-coordinates mesh.
        x, y, z = self._geom.mesh(resolution)
        local = N.array((x, y, z, N.ones_like(x)))
        glob = N.tensordot(self._temp_frame, local, axes=([1], [0]))
        
        return glob[:3] class FlatGeometryManager(GeometryManager):
class FiniteFlatGM(FlatGeometryManager):
class RectPlateGM(FiniteFlatGM):
class RoundPlateGM(FiniteFlatGM):
    def find_intersections(self, frame, ray_bundle):
    def select_rays(self, idxs):
    def get_normals(self):
    def get_intersection_points_global(self):
    def done(self):
    def __init__(self):
    def find_intersections(self, frame, ray_bundle):
    def select_rays(self, idxs):
    def __init__(self, width, height):
    def find_intersections(self, frame, ray_bundle):
    def mesh(self, resolution):
    def __init__(self, R):
    def find_intersections(self, frame, ray_bundle):
    def mesh(self, resolution): class Reflective(object):
class AbsorptionAccountant(object):
class ReflectiveReceiver(AbsorptionAccountant):
class AbsorberReflector(Reflective):
class RefractiveHomogenous(object):
class LambertianReflector(object):
class LambertianReceiver(AbsorptionAccountant):
    def __init__(self, absorptivity):
    def __call__(self, geometry, rays, selector):
    def __init__(self, real_optics, absorptivity=1.):
    def reset(self):
    def __call__(self, geometry, rays, selector):
    def get_all_hits(self):
    def __init__(self, absorptivity=1.):
    def __call__(self, geometry, rays, selector):
    def __init__(self, n1, n2):
    def toggle_ref_idx(self, current):
    def __call__(self, geometry, rays, selector):
    def __init__(self, absorptivity):
    def __call__(self, geometry, rays, selector):
    def __init__(self, absorptivity=1.):
        R = N.ones(len(selector))
import unittest
import numpy as N
import math
from tracer.tracer_engine import TracerEngine
from tracer.ray_bundle import RayBundle
from tracer.spatial_geometry import general_axis_rotation, generate_transform
from tracer.sphere_surface import HemisphereGM, CutSphereGM
from tracer.boundary_shape import BoundarySphere
from tracer.object import AssembledObject
from tracer.assembly import Assembly
from tracer import assembly
from tracer.paraboloid import Paraboloid
from tracer.surface import Surface
from tracer import flat_surface
from tracer import optics_callables
# -*- coding: utf-8 -*-



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)

class TestObjectBuilding2(unittest.TestCase):
    """Tests an object composed of two flat surfaces"""
    flat_only = True
    def setUp(self):

        self.assembly = Assembly()