function [P, foc_v] = Pifw(cop, focal_length, asize, bsize)
	%This function produces the projection matrix
	%(image from world), a 3x3.  For other programs
	%then, points on the screen will be P.(x-cop)
	%where x is a vector in the world.  This 
	%function is provided a center of projection
	%and a focal length.  35mm film is simulated.
	%The user also indicates the desired resolution,
	%where a is the horizontal and b the vertical.
cop_mag = norm(cop);
cop_unit = cop/cop_mag;
ipp = (cop_mag - focal_length)*cop_unit;
foc_v = ipp-cop;
foc_vu = foc_v/(norm(foc_v));
	%so far then, ipp is the point on the image
	%plane closest to the center of projection,
	%and foc_v is a vector from cop to this point.
z_axis = [0 0 5];
b = z_axis-((dot(z_axis,foc_vu))*foc_vu);
b = 24*(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 = 36*(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 mm. So if multiply
	%a and b by the mm/pixel along their axis, from
	%tail to tip, the resulting projections will be
	%in mm.  Then, finally, P.(x - cop) will give me
	%pixel coordinates (after normalizing).  And I
	%guess that I want pixel 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';
P = inv(Pwfi);