dizhenmonijiluz.m

地球物理学地震子波以及合成记录模拟原代码.

clear;
disp('started!')
path='e:\毕业设计数据\earth';
dir1='模型3';
mkdir(path,dir1);
path=strcat(path,'\',dir1,'\');%整体目录,当模型修改时只需修改这个值即可实现所生成的数据在不同的模型目录下,一个模型下所产生的所有数据都放在此目录下
t=0:0.001:0.02;%s采样间隔为0.1毫秒。
nt=length(t);
fc=35;%子波峰值频率
r(1:nt,1)=0;
r(1,1)=1;
R=(1-2*pi.*pi.*fc.*fc.*t.*t).*exp(-pi.*pi.*fc.*fc.*t.*t);
figure(1)
plot(t,R);
hold on
h1=200;
h2=400;
v1=2000;
v2=3000;
v3=4000;
X=-500:10:500;
dt=0.001;
T=0:0.001:0.6;
m=length(T);
n=length(X);
t1=abs(sqrt(0.25*X.^2+40000))/1000;
t2=abs(sqrt((2*(0.1+4/30))^2+(X.^2)/((0.1*4*10^6+4*9*10^6/30)/(0.1+4/30))));
r1=(v2-v1)./(v2+v1);
r2=((v3-v2)./(v3+v2))*(1-r1.^2);
    tic
    disp('receive');
for(i=1:n)
    if (t1(i)/dt-fix(t1(i)/dt))>=0.5
       t1(i)=fix(t1(i)/dt)+1;
   else t1(i)=fix(t1(i)/dt);
   end 
end
for(i=1:n)
    if (t2(i)/dt-fix(t2(i)/dt))>=0.5
       t2(i)=fix(t2(i)/dt)+1;
   else t2(i)=fix(t2(i)/dt);
   end 
end
R1=zeros(m,n);
for(i=1:n)
    R1(fix(t1(i)),i)=r1;
end
R2=zeros(m,n);
for(i=1:n)
    R2(fix(t2(i)),i)=r2;
end
Re=R1+R2;
srz(1:m+nt-1,1:n)=0;
sr(1:m,1:n)=0;
for(i=1:n)
    srz(1:m+nt-1,i)=conv(R,Re(:,i));
    sr(1:m,i)=srz(1:m,i);
end
for(i=1:n)
    maxsr(i)=max(sr(1:m,i));
end
maxSR=max(maxsr(1:n));
figure(2)
for(i=1:n)
     plot(sr(1:m,i)./maxSR+i,1:m);
     hold on
end 
hold off
ax1 = gca;
set(ax1,'XAxisLocation','top','YAxisLocation','left' ,'XColor','r','YColor','r')
axis(ax1,[1 n 1 m])
axis ij%坐标原点在左上角
path_it=strcat(path,'explore.dat');
fid_it=fopen(path_it,'w');%it代表inverse time,也就是反推出的时间域的值,一个时间段内的值(可以看成各点的接收记录)
fprintf(fid_it,' %.6f',sr);
fclose(fid_it);
clear path_it;
toc
    
disp('end');
return;