function rct=rc(theta,t)
global tc
% RC(THETA,T) calcule le rayon de la lava: il est constant avant 
% d'atteindre le temps critique; puis devient une sinusoïde corrigée
% par une exponentielle.

% Rayon constant
pct=0.3*legendreQ3(cos((theta+0.003)/2))-0.25*legendreP(1,cos(theta));

% Dépendance au temps
ft=zeros(size(t));
i=find(t<tc);
ft=sqrt(2)*exp(-3*(t-tc)).*cos(6*pi*(t-tc)+pi/4);
ft(i)=t(i)./(t(i)+eps);

% le rayon final
rct=0.8*(1+pct.*ft);
