function [I] = raytrace_with_planes(spheres, planes, lights, cop, dov, xres, yres)
%Does planes in addition to the things that raytrace does. (Not well documented).
%
%Returns an image of the spheres with the light sources
%from the given center of projection and direction of view.
%The spheres list is n by 14 (center, radius, material
%properties).
%the true color image.
I = zeros(xres,yres,3);
%First, set up the Pwfi matrix.
focal_length = .050;
[Pwfi, foc_v] = Pwfi(cop,dov,focal_length,xres,yres);
%Now the flow of the code will be as such:
%for each ray
%   for each sphere
%       detect intersections
%       if smaller t less than current-small t
%           update closest sphere and current-small t
%   if there was a sphere
%       for each light
%           compute specular and diffuse terms to calculate 
%           lighting at that intersection point to light pixel.
%
%And hopefully that should do it.
q = size(spheres);
n_spheres = q(1);
q = size(lights);
n_lights = q(1);
q = size(planes);
n_planes = q(1);
%A ray is given by cop + t*Pwfi.[u v 1]' where u and v are steps
%in xres and yres horizontal and vertical pixels.
e = cop; %just the eye.  Makes for less typing below.
for i = 1:xres
    for j = 1:yres
        min_t = 1e+90; 
        closest_s = 0; closest_p = 0;  s_or_p = 0;
        %remember t=1 is exactly the image plane and so should 
        %be the minimum allowable intersection. s_or_p keeps track
        %of what type of object was closest, -1 for sphere, 1 for plane.
        v = Pwfi*[i j 1]';
        for s = 1:n_spheres
            c = spheres(s,1:3);  %current center.
            r = spheres(s,4);    %radius
            %now cop + t*v is the ray and we need to solve it 
            %simultaneously with the sphere.  The following turns
            %out to be the solution to the quadratic equation.
            t = intersect(e,v,c,r);  %intersects the sphere with the ray, returning t.
            if abs(imag(t(1))) < 1e-20
                if (real(t(1)) < min_t) & (real(t(1)) > 1)
                    min_t = real(t(1));
                    closest_s = s;
                    s_or_p = -1;
                end
            end
            if abs(imag(t(2))) < 1e-20
                if (real(t(2)) < min_t) & (real(t(2)) > 1)
                    min_t = real(t(2));
                    closest_s = s;
                    s_or_p = -1;
                end
            end
        end
        %Now check for closer planes.
        for p = 1:n_planes
            p_on_plane = planes(p,1:3);
            p_normal = planes(p,4:6);
            t = p_intersect(e,v',p_on_plane,p_normal);
            if t > 1 & t < min_t
                min_t = t;
                closest_p = p;
                s_or_p = 1;
            end
        end
        %Now, if there was a closest sphere or plane, let's do some light computation.
        if s_or_p == -1 | s_or_p == 1
        %if closest_s > 0
            %Pick up the material properties.
            sum_l = [0 0 0]; %for the light properties.
            if s_or_p == -1
                am = spheres(closest_s, 5:7);
                dif = spheres(closest_s, 8:10);
                spe = spheres(closest_s, 11:13);
                ex = spheres(closest_s, 14);
            else 
                if s_or_p == 1
                    am = planes(closest_p,7:9);
                    dif = planes(closest_p,10:12);
                    spe = planes(closest_p,13:15);
                    ex = planes(closest_p,16);
                else
                    ['error1'];
                end
            end
            %And the normal of the sphere at the point of intersection.
            %Remember am, dif, and spe are 1x3 vectors, for r,g,b.
            p = cop + min_t*v';  %the point that the ray intersects the sphere or plane.
                                 %I'm reusing variables for different purposes: bad.
            if s_or_p == -1
                n = p - spheres(closest_s,1:3);
            else 
                if s_or_p == 1
                    n = planes(closest_p,4:6);
                end
            end    
            n = n/norm(n);      %the unit normal on the sphere.
            view = cop - p;
            view = view/norm(view); %the unit view direction vector.
            w = xres - j + 1;           %this is for the wierd matlab output.
            %I(w,i,1:3) = am;
            for k = 1:n_lights
                temp2 = zeros(1,1,3);
                l = lights(k,1:3);      %current light.
                sum_l = sum_l + lights(k,4:6);   %the light's color.
                l_dif = lights(k,7:9);  %no idea what a diffuse filter means.
                l_spe = lights(k,10:12);%seems like a softening of the light. 
                l = l - p;
                l = l/norm(l);  %unit light vector from the point.
                %Now we have to do shadow calculation, which is a while loop
                %through all of the other spheres, doing lighting only when 
                %we can.
                s = 1; shad = 0;
                while s <= n_spheres & shad == 0
                    if s ~= closest_s
                        t = intersect(p,l,spheres(s,1:3),spheres(s,4));
                        if (abs(imag(t(1))) < 1e-20 & (real(t(1)) > 0)) | (abs(imag(t(2))) < 1e-20 & (real(t(2)) > 0))
                            shad = 1;
                        end
                    end
                    s = s + 1;
                end
                q = 1;  %go through the planes for shadows.
                if closest_p ~= 1
                    x = [0];
                end
                while q <= n_planes & shad == 0
                    if q ~= closest_p
                        t = p_intersect(p,l,planes(q,1:3),planes(q,4:6));
                        if t > 0 & t < norm(lights(k,1:3) - p)
                            shad = 1;
                        end
                    end
                    q = q + 1;
                end
                if shad == 0
                    do = dot(l,n);
                    if do < 0
                        do = 0;
                    else
                        ref = 2*do*n - l;
                        do2 = dot(view,ref);
                        if do2 < 0
                            do2 = 0;
                        end
                        temp2 = do*dif.*l_dif + (do2^ex)*spe.*l_spe;
                    end
                    %size((I(i,j,1:3)))
                    %I(i,j,1:3) = (I(i,j,1:3)) + do*dif + ((dot(view,ref))^ex)*spe;
                    I(w,i,1) = (I(w,i,1)) + temp2(1);
                    I(w,i,2) = (I(w,i,2)) + temp2(2);
                    I(w,i,3) = (I(w,i,3)) + temp2(3);
                end
            end
            %avg_l = avg_l/n_lights;
            if sum_l(1) > 1
                sum_l(1) = 1;
            end
            if sum_l(2) > 1
                sum_l(2) = 1;
            end
            if sum_l(3) > 1
                sum_l(3) = 1;
            end
            I(w,i,1) = (I(w,i,1)) + am(1)*sum_l(1);
            I(w,i,2) = (I(w,i,2)) + am(2)*sum_l(2);
            I(w,i,3) = (I(w,i,3)) + am(3)*sum_l(3);
            if I(w,i,1) > 1
                I(w,i,1) = 1;
            end
            if I(w,i,2) > 1
                I(w,i,2) = 1;
            end
            if I(w,i,3) > 1
                I(w,i,3) = 1;
            end
        end
    end
end