function [Pwfi, foc_v] = Pwfi(cop, dov, focal_length, asize, bsize)
%This is a modified version of the code for Program 1.
%
%This function produces the projection matrix
%(image from world), a 3x3.  For other programs
%then, points in space will be cop + t*Pwfi.x
%where x is a vector in pixel coordinates.  This 
%function is provided a center of projection
%and a focal length and direction of view.  35mm 
%film is simulated. The user also indicates the 
%desired resolution, where a is the horizontal 
%and b the vertical. units are meters.
foc_v = focal_length*(dov/norm(dov));
foc_vu = foc_v/(norm(foc_v));
%foc_v is a vector from the center of 
%projection to the he point on the image
%plane closest to the cop. foc_vu is unit length.
z_axis = [0 0 5];
b = z_axis-((dot(z_axis,foc_vu))*foc_vu);
b = .036*(b/(norm(b)));
%the b image plane vector is the 24mm projection
%of the z-axis onto the image plane. the
%horizontal screen axis a is the cross of foc_v
%and b, of magnitude 36mm.
a = cross(foc_vu,b);
a = .036*(a/(norm(a)));
o = foc_v - .5*b - .5*a;
%Now, we have o, b, and a, all in mm.  But the
%the screen coords are given in pixels, with
%the result (in the world) in m. So if multiply
%a and b by the mm/pixel along their axis, from
%tail to tip, the resulting projections will be
%in m.  Then, finally, t*P.x (given pixel coords)
%will give me meter coordinates.
%And I guess that I want meter coordinates, so...the
%below makes my images bxa.
a = (1/asize)*a;
b = (1/bsize)*b;
Pwfi = zeros(3);
Pwfi(1:3,1) = a';
Pwfi(1:3,2) = b';
Pwfi(1:3,3) = o';