边界的声波逆时偏移

初步的代码，只保存边界的声波逆时偏移 RTM

%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