The problem is that there is not an exact whole number of cycles of the sine wave in the original time series. If you look at the values of x(1) and x(end), you'll find they are the same, because the final point is at time t=1, when a new cycle is just starting.
Unless you have an exactly whole number of cycles of a sine wave present in the initial series, there won't be a sample point in the frequency domain lying at exactly the frequency of the wave. In your case, the sample at X(101) (with its complex conjugate at X(902)) is close to, but not exactly at, the frequency f=100. This slight shift in frequency doesn't affect the amplitude too much, but it has a big effect on the phase.
A minimal change to get the answer you expect is to reduce the length of the input by 1, for example like this:
which will give a whole number of cycles provided f is an integer. (As Walter Roberson points out, linspace may be preferable - but rounding errors aren't in fact the problem in this case.) This change gives the phase as expected, as you can see by printing the value of P(101).
There's a general lesson: if you have real-world data, you are very unlikely to have an exactly whole number of cycles in the sequence, and so you can't get the phase easily from the FFT.
Best Answer