How do I get the phase at each point from a sweep sine function?

Hello, I am trying to extract the phase from a sweep sine function, in order to find the phase shift between two sweep sine function, at each point. Given the function

S(t) = cos(2*pi*f(t)*t + Phi(t))

where f(t) = a*t + fmin and Phi(t) = pi*t (as a simple example)

I am trying to find Phi(t).

I have been using the following code to test this but so far no results.

%constant definition------------------------------------
Fmin = 1; Fmax = 10;
Tmin = 0; Tmax = 1;
fech = 1000;
alpha = (Fmax - Fmin)/(Tmax - Tmin);

%definition of f and Phi--------------------------------
t = linspace(Tmin,Tmax,(Tmax-Tmin)*fech);
f = alpha*t + Fmin; 
Phi = 1*pi*t;

S1 = cos(2*pi*f.*t+Phi);
S2 = cos(2*pi*f.*t);

%Phase calculation using the fourier transform----------
S1FFT = fftshift(fft(S1));
S2FFT = fftshift(fft(S2));

Phi1 = atan2(imag(S1FFT),real(S1FFT));
Phi2 = atan2(imag(S2FFT),real(S2FFT));

Dphi = Phi2-Phi1;

Is this feasable and if yes as anyone got any ideas on how to do this?

Thank you very much Sincerely yours, Jean