E0=1
nu=0.042
nus=0.02772
nui=0.01428
v=3.5
C=zeros(13,61)
E= zeros(13,61);
for t= 1:13,
for r=40:100,
alpha=v*t/r;
y=exp(-alpha^2)*alpha^2;
F=trapz(y); % calcul de F: partie integral dans l'expression de C
C(t,r)=E0*(1-(1+nus*v*t)*exp(-nus*v*t))/(4/sqrt(pi)*F);% calcul de C
if (t>r/v)% heaviside vaut 1 si t>r/v et 0 si t<r/v
E(t,r)=E0*exp(-nu*v*t)*(dirac(t-r/v)/(4*pi*v*r^2)+ nus*1/(4*pi*v*r*t)*log((1+r/(v*t)))/(1-r/(v*t))+C*1*(3*nus/(4*pi*v*r*t))^1.5*exp((-3*nus*r^2)/(4*v*t)-nui*v*t));
elseif (t<r/v)
E(t,r)=E0*exp(-nu*v*t)*(dirac(t-r/v)/(4*pi*v*r^2));
end
end
plot(log10(4*pi*r*r*E(t,r)),r);
end
Partager