%RTM
clear;
clc
tic
nt=1201;????%?number?of?time?steps

isnap=10;????%?snapshot?sampling

c=zeros（150150）;
nz=150;
nx=150;
for?i=1:150
????for?j=1:150
????????c（ij）=2000;
????end
end

v=c;

coeff=[?-3.02359?1.75?-0.291667?0.0648148?-0.0132576?0.00212121?-0.000226625?0.0000118929];%tra
p=zeros（[nz?nx]）;??pold=p;?d2px=p;?d2pz=p;pboundarynew=p;

M1=zeros（nznx3）;
s11=zeros（nznx3）;
s22=zeros（nznx3）;

dx=15;??%?calculate?space?increment
x=（0:（nx-1））*dx;
z=（0:（nz-1））*dx;??%?initialize?space?coordinates
dt=0.001;?%?calculate?time?step?from?stability?criterion
h=dx;
tau=dt;
src=zeros（1nt）;

src（1:size（ricker（0.00120）））=ricker（0.00120）;
zs=22;
for?shot=1:1
????
????c=zeros（150150）;
????nz=150;
????nx=150;
????for?i=1:150
????????for?j=1:150
????????????c（ij）=2000;
????????end
????end
????v=c;

????if?shot==1
????????xs=35;
????elseif?shot==2
????????xs=75;
????else
????????xs=105;
????end
????
????seis_record1=zeros（ntnx）;
????
????for?it=1:nt
????????%正演?雙程波的波場
????????for?k=2:nz-2
????????????for?j=2:nx-2
????????????????if?（k<8?||?j<8?||?k>nz?-7?||?j>nx-7）
????????????????????d2px（kj）=-2*p（kj）/dx^2;
????????????????????d2px（kj）=d2px（kj）+（p（k+1j）+p（k-1j））/dx^2;
????????????????????d2pz（kj）=-2*p（kj）/dx^2;
????????????????????d2pz（kj）=d2pz（kj）+（p（kj+1）+p（kj-1））/dx^2;
????????????????????
????????????????else
????????????????????d2px（kj）=coeff（1）*p（kj）/dx^2;
????????????????????for?ii=2:length（coeff）
????????????????????????d2px（kj）=d2px（kj）+coeff（ii）*（p（k-ii+1j）+p（k+ii-1j））/dx^2;
????????????????????end
????????????????????
????????????????????d2pz（kj）=coeff（1）*p（kj）/dx^2;
????????????????????for?ii=2:length（coeff）
????????????????????????d2pz（kj）=d2pz（kj）+coeff（ii）*（p（kj-ii+1）+p（kj+ii-1））/dx^2;
????????????????????end
????????????????end
????????????????
????????????????%正演?雙程波的波場
????????????end
????????end
????????pnew=2*p-pold+c.*c.*（d2px+d2pz）*dt^2;
????????
????????%單程波的波場
????????for?kkk=1:20
????????????%上邊界
????????????k=20-kkk+1;
????????????for?j=22+1-kkk:nx-21+kkk-1?%對于每一行，左右各有兩個角點先不算
????????????????pboundarynew（kj）=（2*h*tau^2*v（kj）*（（pold（kj）?-?pold（k+1j）?+?pnew（k+1j））/（2*h*tau）?-?（pold（kj）?-?2*p（k+1j）?-?2*p（kj）?+?pold（k+1j）?+?pnew（k+1j））/（2*tau^2*v（kj））?+?...
????????????????????（v（kj）*（pold（kj+1）?-?2*pold（kj）?+?pold（kj-1）?-?2*pnew（k+1j）?+?pnew（k+1j+1）?+?pnew（k+1j-1）））/（4*h^2）））/（h?+?tau*v（kj））;
????????????end
????????????
????????????%底邊界
????????????k=nz-20+kkk;
????????????for?j=22+1-kkk:nx-21+kkk-1?%對于每一行，左右各有兩個角點先不算
????????????????pboundarynew（kj）=?（2*h*tau^2*v（kj）*（（pold（kj）?-?pold（k-1j）?+??pnew（k-1j））/（2*h*tau）?-?（pold（kj）?-?2*p（k-1j）-?2*p（kj）?+?pold（k-1j）?...
????????????????????+?pnew（k-1j））/（2*tau^2*v（kj））?+?（v（kj）*（pold（kj+1）?-?2*pold（kj）?+?pold（kj-1）?-?2*pnew（k-1j）?+?pnew（k-1j+1）+?pnew（k-1j-1）））/（4*h^2）））/（h?+?tau*v（kj））;
????????????end
??????????

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件??????27025??2013-04-18?08:50??total8bian.m

-----------?---------??----------?-----??----

????????????????27025????????????????????1