Why is the phase of the FFT not exactly 0 degrees for a cosine and 90 degrees for sine wave?

Question

Let's say I have two signals. The first is a cosine wave and the second is a sine wave. Each oscillates at 0.01 Hz. The sample rate is 1 Hz and the length of time series is 1000 seconds. Each has an amplitude of 1.

My understanding is that an FFT of these two signals should recover a spike in the amplitude spectrum of 0.5 at 0.01 Hz (and similarly at the corresponding negative frequency). This all makes sense and I can get this to work as expected.

But the phase of the FFT is a bit perplexing. My expectation is that the cosine wave will give a phase of 0° at 0.01 Hz, and the sine wave will give a phase of 90° at 0.01 Hz. However, the result I get gives a phase of 3.6° for cosine, and 86.4° for sine. (The negative frequencies are complex conjugates).

Why can the FFT recover the precise amplitudes, but can't do so with the phases? Is there some reason for this? Is it just some sort of numerical or indexing issue or is there some deeper reason? Is my Fourier frequency list incorrect and I'm not actually sampling the spectrum at exactly 0.01 Hz?

MATLAB code to replicate is below.

Any help is appreciated.

f = 0.01; %signal frequency
fs = 1; %sample rate
dt = 1./fs; %length of sample
t = (dt:dt:1000)'; %time vector

df = fs/length(t); %frequency spacing
fAxis = (0:df:(fs-df)) - (fs-mod(length(t),2)*df)/2; %frequency axis with negative freqs

b1 = cos(2*pi*f*t); %first time series signal
b2 = sin(2*pi*f*t); %second time series signal

%FFT each signal and scale
B1 = fftshift((fft(b1)./length(fAxis)));
B2 = fftshift((fft(b2)./length(fAxis)));

%Find index of 0.01 frequency
indf = find(abs(fAxis-0.01)<10^-9);

%Magnitudes look okay. Each returns 0.5
abs(B1(indf)) 
abs(B2(indf))

%But phases???
angle(B1(indf))*180/pi
angle(B2(indf))*180/pi

Answer 1

Your time vector is slightly off, it should start at t=0, and end at t=999 (this latter modification keeps the length of the FFT at 1000 to keep a resolution $d_f = 1/1000$, allowing for the frequency component at $f = 0.01\texttt{Hz}$ to fall in one bin):

t = (0:dt:1000-1/fs)'; %time vector

Result:

ans =

   8.8369e-14

ans =

  -90.0000