最短路径法弯曲射线追踪的MATLAB代码.zip

clc clear close all % Moser曲射线追踪主程序 % 算法： % 作者：shangxiang % 时间：2021.3.9 % xmin = 0; xmax = 6; ymin = 0; ymax = 7; dx = 1; dy = 1; ndot = 5; % 内插节点数 nex = (xmax - xmin)/dx; ney = (ymax - ymin)/dy; exy = nex*ney; nxydot = (nex+1)*ndot*ney + (ney+1)*ndot*nex; S2d = (1/4000)*ones(ney,nex); % 慢度 S2d(4:4,2:5) = 1/2400; % 画网格图 [xx,yy] = meshgrid(xmin:dx:xmax,ymin:dy:ymax); plot(xx,yy,'k:',xx',yy','k:');hold on axis([xmin - .2*dx, xmax + .2*dx, ymin - .2*dy, ymax + .2*dy]); xlabel('x/m'); ylabel('y/m'); set(gca,'FontSize',14); % 画路径显示图，iplot = 0,不画，= 其他,则画。 iplot = 1; % Floyd最短路径算法计算最短路径 [d,r,pz] = sMoserRoadsub(xmin,xmax,ymin,ymax,dx,dy,S2d,ndot); a = 1; b = nxydot - 35; while 1 if (r(a,b) ~= b) plot(pz(a,1),pz(a,2),'bo'); c = a; a = r(a,b); plot([pz(a,1),pz(c,1)],[pz(a,2),pz(c,2)],'b-','LineWidth',.5); else plot(pz(a,1),pz(a,2),'bo'); plot(pz(b,1),pz(b,2),'bo'); plot([pz(a,1),pz(b,1)],[pz(a,2),pz(b,2)],'b-','LineWidth',.5); break; end end shortest_time = d(a,b); Fresnel_zone = 1/(2*.24*1e3) + shortest_time; clear a b a = 1; b = nxydot - 35; for i = 1:nxydot S_shortest_time = d(a,i); R_shortest_time = d(b,i); SR_time = S_shortest_time + R_shortest_time; if SR_time <= Fresnel_zone % plot(pz(i,1),pz(i,2),'r*'); plot([pz(a,1),pz(i,1)],[pz(a,2),pz(i,2)],'b-','LineWidth',.5); plot([pz(i,1),pz(b,1)],[pz(i,2),pz(b,2)],'b-','LineWidth',.5); end end figure(2) zzS = ones(ney+1,nex+1); zzS(1:ney,1:nex) = S2d; pcolor(xx,yy,zzS); colormap('jet'); h = colorbar; h.Label.String = 'Slowness (s/m)'; caxis([0.0001 0.0030]); set(gca,'FontSize',14);