fortran
matlab
程序的编写。program seismic implicit none real::lamda1,mu1,r1,dt,dx,dz,vp1,vs1,freq,lamda2,mu2,r2,vp2,vs2 real::c11,c12,c21,c22,c31,c32,d1,d2,d3,pi real,dimension(101)::x,z real,dimension(101,101)::uz1,uz2,uz3,u3,ux1,ux2,ux3 integer::i,j,k,n,c character*10::text,nfile1,nfile2,tp write(*,*)'input ytpe(mid-reflect,len-direct)' read(*,'(a)') tp if(tp.eq.'mid') then print *, 'dx,dz' read(*,*) dx,dz else if(tp.eq.'len') then print *, 'dx,dz' read(*,*) dx,dz end if vp1=1500.0;vs1=1100.0 dt=0.00025 freq=10.0 r1=300.0 lamda1=(vp1**2-2*vs1**2)*r1 mu1=vs1**2*r1 pi=atan(1.0)*4.00 c11=0.5*(lamda1+2.0*mu1)/r1 c21=0.5*(lamda1+mu1)/r1 c31=mu1/(2.0*r1) d1=(dt/dx)**2 d2=dt**2/(4.0*dx*dz) d3=(dt/dz)**2 n=500 ux1=0.0;uz1=0.0 ux2=0.0;uz2=0.0 ux3=0.0;uz3=0.0 do i=2,100 do j=2,100 ux2(i,j)=ux1(i,j)+c11*d1*(ux1(i+1,j)+ux1(i-1,j)-2*ux1(i,j)) & +c21*d2*((uz1(i+1,j+1)-uz1(i-1,j+1))-(uz1(i+1,j-1)-uz1(i-1,j-1))) & +c31*d3*(ux1(i,j+1)+ux1(i,j-1)-2*ux1(i,j)) uz2(i,j)=uz1(i,j)+c11*d1*(uz1(i,j+1)+uz1(i,j-1)-2*uz1(i,j)) & +c21*d2*((ux1(i+1,j+1)-ux1(i-1,j+1))-(ux1(i+1,j-1)-ux1(i-1,j-1))) & +c31*d3*(uz1(i+1,j)+uz1(i-1,j)-2*uz1(i,j)) end do end do ux2(51,51)=0.0 uz2(51,51)=0.0 ux2(50,50)=-3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0) uz2(50,50)=-3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0) ux2(52,50)=3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0) uz2(52,50)=-3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0) ux2(50,52)=-3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0) uz2(50,52)=3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0) ux2(52,52)=3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0) uz2(52,52)=3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0) ux2(50,51)=-3.0*sqrt(2.0)*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0) ux2(52,51)=3.0*sqrt(2.0)*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0) uz2(51,50)=-3.0*sqrt(2.0)*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0) uz2(51,52)=3.0*sqrt(2.0)*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0) do k=2,n do i=2,100 do j=2,100 ux3(i,j)=-ux1(i,j)+2*ux2(i,j)+2*c11*d1*(ux2(i+1,j)+ux2(i-1,j)-2*ux2(i,j)) & +2*c21*d2*((uz2(i+1,j+1)-uz2(i-1,j+1))-(uz2(i+1,j-1)-uz2(i-1,j-1)))& +2*c31*d3*(ux2(i,j+1)+ux2(i,j-1)-2*ux2(i,j)) uz3(i,j)=-uz1(i,j)+2*uz2(i,j)+2*c11*d1*(uz2(i,j+1)+uz2(i,j-1)-2*uz2(i,j)) & +2*c21*d2*((ux2(i+1,j+1)-ux2(i-1,j+1))-(ux2(i+1,j-1)-ux2(i-1,j-1)))& +2*c31*d3*(uz2(i+1,j)+uz2(i-1,j)-2*uz2(i,j)) u3(i,j)=sqrt(ux3(i,j)**2+uz3(i,j)**2) end do end do ux3(51,51)=0.0 uz3(51,51)=0.0 ux3(50,50)=-3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0) uz3(50,50)=-3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0) ux3(52,50)=3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0) uz3(52,50)=-3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0) ux3(50,52)=-3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0) uz3(50,52)=3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0) ux3(52,52)=3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0) uz3(52,52)=3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0) ux3(50,51)=-3.0*sqrt(2.0)*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0) ux3(52,51)=3.0*sqrt(2.0)*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0) uz3(51,50)=-3.0*sqrt(2.0)*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0) uz3(51,52)=3.0*sqrt(2.0)*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0) write(text,'(i3)') k c=mod(k,10) if(c .eq. 0) then nfile1=text(1:3)//'ux.dat' open(88,file=nfile1) do i=1,101 do j=1,101 x(i)=(i-1)*dx z(j)=(j-1)*dzwrite(88,'(1x,i3,f10.4,3x,f10.4,3x,f16.4,3x,f16.4,3x,f16.4)')k,x(i),z(j),ux3(i,j),& uz3(i,j),u3(i,j) end do end do close(88) end if do i=2,100 do j=2,100 ux1(i,j)=ux2(i,j) ux2(i,j)=ux3(i,j) uz1(i,j)=uz2(i,j) uz2(i,j)=uz3(i,j) end do end do end do stop end program seismic
实例采用网格数为100*100，x方向和z方向的间隔dx=dz=1，时间间隔dt=0.00025s， 模拟地质体的纵波速度vp=1500.0m/s，横波速度为vs=1100.0m/s，震源的频率为freq=10.0 Hz，地质体的密度为r=300.0.震源放在模型区域的中心位置，其振动为爆炸震源。通过程序计算结果如下：当T=10*dt时间后，由于处于爆炸的初期，波刚刚向前传播了10*dt*vp的距离。如下图一所示。
当T=200*dt时，波往外传播200*dt*vp=500*0.00025*1500=75m，则波发生反射，如下图三所
有兴趣的可以私话我，一起探讨。