完整程序：

clear all
dx=50;x=-500:dx:500;%炮检距
h=100;V1=1500;
theta=25*pi/180;
V2=V1/sin(theta);
t1=sqrt(x.*x+4*h*h)/V1;%反射波时距曲线
t2=abs(x)./V1;%直达波时距曲线
%折射波时距曲线
xm=2*h*tan(theta);%求盲区
k=1;
for i=1:length(x)
    if x(i)<-xm
        t3_1(k)=2*h*cos(theta)/V1+abs(x(i))/V2;
        x_1(k)=x(i);
        k=k+1;
    end
end
k=1;
for i=1:length(x)
    if x(i)>xm
        t3_2(k)=2*h*cos(theta)/V1+abs(x(i))/V2;
        x_2(k)=x(i);
        k=k+1;
    end
end
plot(x,t1,'b')
hold on
plot(x,t2,'k-.');
plot(x_1,t3_1,'r:');
plot(x_2,t3_2,'r:');
set(gca,'ydir','reverse')
xlabel('Shot-receiver distance(m)');
ylabel('Two-way traveltime(s)');
legend('反射波','直达波','折射波')