%ICONAL resout eq iconale par différences finies dans milieu à gradient de vitesse
% sans interface

clear all

% grille(x,z)
nx = 50; dx = 2;  xl = (nx-1)*dx; xx = 0:dx:xl;
nz = 50; dz = dx; zl = (nz-1)*dz; zz = 0:dz:zl;
% t(x,z)
txz = NaN*ones(nz,nx);


% V0, gradient, pendage
V0 = 500; Vxz = V0*ones(nz,nx);
aa = 10; al = 30*pi/180; ax = aa*sin(al); az = aa*cos(al); for ix = 1:nx, Vxz(:,ix) = V0+ax*(ix-1)*dx+az*zz'; end


% DIFFERENCES FINIES
t1 = dx./Vxz(:,:); t2 = 2*t1.^2;

% onde sphérique dans carré initial de taille ic
ic = 5;
for ia = 1:ic
for ix = 1:ia, txz(ia,ix) = sqrt(((ia-1)^2+(ix-1)^2)*dx^2)/V0; txz(ix,ia) = txz(ia,ix); end
end
txz0 = txz;

% différences finies sur ligne et colonnes ia
for ia = ic:nx-1
% colonne ix = ia
ix = ia;
% temps minimum dans colonne en iz0
iz0 = find(txz(1:ia,ix) == min(txz(1:ia,ix)));
% avance en x+dx, z0
txz(iz0,ix+1) = txz(iz0,ix)+t1(iz0,ix);
% avance en x+dx  pour z>z0
for iz = iz0:ia 
tmp = t2(iz+1,ix+1)-(txz(iz,ix+1)-txz(iz+1,ix)).^2;
if tmp>=0, txz(iz+1,ix+1) = txz(iz,ix)+sqrt(tmp); end
end
% avance en x+dx  pour z<z0
for iz = iz0:-1:2
tmp = t2(iz-1,ix+1)-(txz(iz,ix+1)-txz(iz-1,ix)).^2;
if tmp>=0, txz(iz-1,ix+1) = txz(iz,ix)+sqrt(tmp); end
end

% ligne iz = ia
iz = ia;
% temps minimum dans ligne en ix0
ix0 = find(txz(iz,1:ia) == min(txz(iz,1:ia)));
% descend en z+dz, x0
txz(iz+1,ix0) = txz(iz,ix0)+t1(iz,ix0);
% descend en z+dz  pour x>x0
for ix = ix0:ia 
tmp = t2(iz+1,ix+1)-(txz(iz,ix+1)-txz(iz+1,ix)).^2;
if tmp>=0, txz(iz+1,ix+1) = txz(iz,ix)+sqrt(tmp); end
end
% descend en z+dz  pour x<x0
for ix = ix0:-1:2
tmp = t2(iz+1,ix-1)-(txz(iz+1,ix)-txz(iz,ix-1)).^2;
if tmp>=0, txz(iz+1,ix-1) = txz(iz,ix)+sqrt(tmp); end
end
if ia == fix(nx/2), txz1 = txz; end
end



figure(1), clf,
tl = max(max(txz));  

subplot(2,2,1), box on,  hold on
pcolor(xx,zz,Vxz)
if (aa ~= 0), title(['V = ' num2str(V0) ' +  '  num2str(ax,2)  ' X + '  num2str(az,2)  ' Z M/S'])
else title(['V = ' num2str(V0) ' M/S'])
end
xlabel('X (M)'),ylabel('Z (M)'),set(gca,'Ydir','reverse','DataAspectRatio',[1 1 1]), axis([0 xl 0 zl]), colorbar    % shading interp

subplot(2,2,2),box on,  hold on
txz0(nx,nz) = tl; pcolor(xx,zz,txz0)
title('T(X,Z) SPHERIQUE INITIAL')
xlabel('X (M)'),ylabel('Z (M)'),set(gca,'Ydir','reverse','DataAspectRatio',[1 1 1]), axis([0 xl 0 zl]), colorbar,   % shading interp

subplot(2,2,4),box on,  hold on
txz1(nx,nz) = tl; pcolor(xx,zz,txz1)
title(['T(X,Z) APRES ' num2str(fix(nx/2)) ' AVANCES '])
xlabel('X (M)'),ylabel('Z (M)'),set(gca,'Ydir','reverse','DataAspectRatio',[1 1 1]), axis([0 xl 0 zl]), colorbar,   % shading interp

subplot(2,2,3),box on,  hold on
pcolor(xx,zz,txz)
% fronts d'onde à ti
if (aa ~= 0 && al == 0)
for ti = tl/4:tl/4:tl
at = aa*ti; R = V0*sinh(at)/aa; xt = xx(xx<=R); zt = sqrt(R^2-xt.^2); Z = V0*(cosh(at)-1)/aa;
plot([xt R],[Z+zt Z],'w',[xt R],[Z-zt Z],'w','Linewidth',2)
end
end
title('T(X,Z)')
xlabel('X (M)'),ylabel('Z (M)'),set(gca,'Ydir','reverse','DataAspectRatio',[1 1 1]), axis([0 xl 0 zl]), colorbar,  shading interp



% COMPARAISON TEMPS ARRIVEE NUMERIQUE ET CALCULE
if al == 0
figure(2), clf
% T(ix,z)
ix = [1 fix(nx/4)+1 fix(nx/2)+1 fix(3*nx/4)+1 nx]; xi = (ix-1)*dx;

subplot(1,2,1),box on, grid on, hold on
plot(zz,txz(:,ix(1)),'b',zz,txz(:,ix(2)),'c',zz,txz(:,ix(3)),'g',zz,txz(:,ix(4)),'y',zz,txz(:,ix(5)),'r','LineWidth',2)
% temps direct à V0
if aa == 0, for ii = 1:5, plot(zz,sqrt(xi(ii)^2+zz.^2)/V0,'ko'), end
% temps direct à V0+az
else if al == 0 
for ii = 5:-1:1, 
% paramètre pour rayon entre (0,0) et (x,z) 
ax2 = (aa*xi(ii)).^2; b2 = V0^2+Vxz(:,ix(ii)).^2; c2 = (V0*Vxz(:,ix(ii))).^2; pz = sqrt(4*ax2./((ax2+b2).^2-4*c2));
th0 = asin(V0*pz);  thz = asin(Vxz(:,ix(ii)).*pz); tz = -1/aa*log(tan(th0/2)./tan(thz/2));
% X au point bas des rayons
Xz = Vxz(:,ix(ii))/aa.*sqrt(1-(V0./Vxz(:,ix(ii))).^2); 
plot(zz(xi(ii)<=Xz),tz(xi(ii)<=Xz),'ko',zz(xi(ii)>Xz),-2/aa*log(tan(th0(xi(ii)>Xz)/2))-tz(xi(ii)>Xz),'ks',zz,1/aa*log(1+aa/V0*zz),'ko'), end, end
end
if aa ~= 0, 
legend(['NUM. X = ' num2str(xi(1))],['NUM. X = ' num2str(xi(2))],['NUM. X = ' num2str(xi(3))],['NUM. X = ' num2str(xi(4))],['NUM. X = ' num2str(xi(5)) ],'CAL. RAYON DESCENDANT','CAL. RAYON MONTANT',3)
title(['COMPARAISON T(X, Z=cte) NUMERIQUE ET PAR CALCUL'])
else
legend(['NUM. X = ' num2str(xi(1))],['NUM. X = ' num2str(xi(2))],['NUM. X = ' num2str(xi(3))],['NUM. X = ' num2str(xi(4))],['NUM. X = ' num2str(xi(5)) ],3)
title(['COMPARAISON T(X, Z=cte) NUMERIQUE ET PAR CALCUL'])
end
xlabel('Z (M)'),ylabel('T (S)'), set(gca,'Color',[1. 0.8 0.5]), axis([0 zl 0 tl]), view([90 90])

%T(x,iz)
iz = [1 fix(nz/4)+1 fix(nz/2)+1 fix(3*nz/4)+1 nz]; zi = (iz-1)*dz; ax2 = (aa*xx).^2; 

subplot(1,2,2),box on, grid on, hold on
plot(xx,txz(iz(1),:),'b',xx,txz(iz(2),:),'c',xx,txz(iz(3),:),'g',xx,txz(iz(4),:),'y',xx,txz(iz(5),:),'r','LineWidth',2)
% temps direct à V0 
if aa == 0, for ii = 1:5, plot(xx,sqrt(zi(ii)^2+xx.^2)/V0,'ko'), end
% temps direct à V0+az
else if al == 0
for ii = 1:5, 
% paramètre pour rayon entre (0,0) et (x,z) 
b2 = V0^2+Vxz(iz(ii)).^2; c2 = (V0*Vxz(iz(ii))).^2; pz = sqrt(4*ax2./((ax2+b2).^2-4*c2));
th0 = asin(V0*pz); thz = asin(Vxz(iz(ii))*pz); tz = -1/aa*log(tan(th0/2)./tan(thz/2));
% X au point bas des rayons
Xz = Vxz(iz(ii))/aa*sqrt(1-(V0/Vxz(iz(ii)))^2);
plot(xx(xx<=Xz),tz(xx<=Xz),'ko',xx(xx>Xz), -2/aa*log(tan(th0(xx > Xz)/2))-tz(xx > Xz),'ks',0,1/aa*log(1+aa/V0*zi),'ko'), end, end
end
if aa ~= 0, 
legend(['NUM. Z = ' num2str(zi(1))],['NUM. Z = ' num2str(zi(2))],['NUM. Z = ' num2str(zi(3))],['NUM. Z = ' num2str(zi(4))],['NUM. Z = ' num2str(zi(5)) ],'CAL. RAYON DESCENDANT','CAL. RAYON MONTANT',1)
title(['COMPARAISON T(X, Z=cte) NUMERIQUE ET PAR CALCUL'])
else
legend(['NUM. Z = ' num2str(zi(1))],['NUM. Z = ' num2str(zi(2))],['NUM. Z = ' num2str(zi(3))],['NUM. Z = ' num2str(zi(4))],['NUM. Z = ' num2str(zi(5)) ],1), 
title(['COMPARAISON T(X, Z=cte) NUMERIQUE ET PAR CALCUL'])
end
xlabel('X (M)'),ylabel('T (S)'),set(gca,'Ydir','reverse','Color',[1. 0.8 0.5]), axis([0 xl 0 tl]),
end