MATLAB phase of 2D rectangular pulse's Fourier transform
The issue is that the phase value found is between $-\pi$ and $\pi$ (or $0$ and $2\pi$) but it needs to be "unwrapped" to be continuous.

In 1D, the unwrap function will help. Your mileage may vary when applying it to a 2D Fourier transform.

Also, there is a bug in your code that forces the R variable to be real valued:

R=re+imm;
phase=angle(R);


is not going to get you the right phase at all.