The original question:
Distribute N points according to a homogeneous Poisson point process in a two-dimensional space of unit area.
That doesn't make sense if N is a fixed number. For the definition of the Poisson point process, the N has to be a Poisson random variable with its mean related to the area/size of the simulation region. This is non-negotiable. But if you fix N=n to some natural number (that is, in probability language, you condition on N=n), you then get a binomial point process.
As Nguyen Anh suggested above (and I wrote in another response), you treat this problem in polar coordinates.
For each point, you uniformly choose a random angular coordinate (or component) on the interval 
. For the radial coordinate (or component), you first choose a random number on the unit interval 
, and then you --- and this is very important -- square root that random number and multiply it by the radius of the disk. You have to take the square root to assure uniform points (because the area of a disk/sector is proportional to the radius squared, not the radius).
. For the radial coordinate (or component), you first choose a random number on the unit interval , and then you --- and this is very important -- square root that random number and multiply it by the radius of the disk. You have to take the square root to assure uniform points (because the area of a disk/sector is proportional to the radius squared, not the radius).I recently blogged about simulating these two point processes. The binomial point process on a rectangle is here. The Poisson point process on a rectangle and on a disk.
The code for a binomial process on a disk would be.
%Simulation window parameters
r=1; %radius of disk
xx0=0; yy0=0; %centre of disk
%Simulate binomial point process
pointsNumber=100;
theta=2*pi*(rand(pointsNumber,1)); %angular coordinates
rho=r*sqrt(rand(pointsNumber,1)); %radial coordinates
%Convert from polar to Cartesian coordinates
[xx,yy]=pol2cart(theta,rho); %x/y coordinates of Poisson points
%Shift centre of disk to (xx0,yy0)
xx=xx+xx0;
yy=yy+yy0;
%Plotting
scatter(xx,yy);
xlabel('x');ylabel('y');
axis square;
