Extract phase information from FFT

Hello, I'm rather a newbie to frequency analysis and in order to get a better understanding of it, I have written myself a small example script that first constructs a 5 Hz sinusoidal signal (subplot 1) with a phase shift of pi and then calculates an FFT of this signal:

clear all
close all
clc

SR = 100;       % sampling rate
SI = 1/SR;      % sampling interval
t = 0:SI:1-SI;  % time (s)

f = 5;          % ordinary frequency (number of oscillations, for each second; Hz)
omega = 2*pi*f; % angular frequency (rate of change of the function argument in units of radians / s)
phi = pi;       % phase shift (in radians), where in its cycle the oscillation is at t = 0
A = 1;          % amplitude
s = A * sin(omega*t + phi);
% time domain
figure;
subplot(1,3,1);
plot(t,s);
ylabel('amplitude');
legend(['f=' num2str(f) ' (Hz) ;   omega=' num2str(omega) ' (r/s);  phi=' num2str(phi) ' (r);   A=' num2str(A)] );
title('time domain'); 
% frequency domain
S = fft(s); 
S = S(1:floor(length(S)/2)); % get rid of right-sided spectrum       
S_mag = abs(S);              % get rid of imaginary part
freqs = linspace(0, SR/2, length(S)); % get frequency scaling
subplot(1,3,2);
plot(freqs,S_mag);
title('frequency domain');
% phase information
S_phase = angle(S); % get phase information          
subplot(1,3,3);
stem(freqs,S_phase);

The frequency spectrum (subplot 2) that I get, looks reasonable to me. However, the phase information that I get (subplot 3), does not look like I expected it: For, what I expected was a phase shift of pi for the frequency bins around 5Hz?!

Can anybody explain me what I am doing wrong here, and how I can correctly extract phase information from my FFT output? I would be very thankful for any help.