tm_model2d.m

来自「Back Projection imaging through wall」· M 代码 · 共 306 行 · 第 1/2 页

% (interior non-PML nodes have sigx=sigz=0, Kx=Kz=1)
sigx = sigxmax.*xdel.^m;
sigz = sigzmax.*zdel.^m;
Kx = 1 + (Kxmax-1)*xdel.^m;
Kz = 1 + (Kzmax-1)*zdel.^m;

% determine FDTD update coefficients
Ca  = (1-dt*sig./(2*ep))./(1+dt*sig./(2*ep));
Cbx  = (dt./ep)./((1+dt*sig./(2*ep))*24*dx.*Kx);
Cbz  = (dt./ep)./((1+dt*sig./(2*ep))*24*dz.*Kz);
Cc = (dt./ep)./(1+dt*sig./(2*ep));
Dbx = (dt./(mu.*Kx*24*dx));
Dbz = (dt./(mu.*Kz*24*dz));
Dc = dt./mu;
Bx = exp(-(sigx./Kx + alpha)*(dt/ep0));
Bz = exp(-(sigz./Kz + alpha)*(dt/ep0));
Ax = (sigx./(sigx.*Kx + Kx.^2*alpha + 1e-20).*(Bx-1))./(24*dx);
Az = (sigz./(sigz.*Kz + Kz.^2*alpha + 1e-20).*(Bz-1))./(24*dz);

% clear unnecessary PML variables as they take up lots of memory
clear sigmax xdel zdel Kx Kz sigx sigz


% ------------------------------------------------------------------------------------------------
% RUN THE FDTD SIMULATION
% ------------------------------------------------------------------------------------------------

disp('Beginning FDTD simulation...')

% initialize gather matrix where data will be stored
gather = zeros(fix((numit-1)/outstep)+1,nrec,nsrc);

% loop over number of sources
for s=1:nsrc

    % zero all field matrices
    Ey = zeros(nx-1,nz-1);          % Ey component of electric field
    Hx = zeros(nx-1,nz);            % Hx component of magnetic field
    Hz = zeros(nx,nz-1);            % Hz component of magnetic field
    Eydiffx = zeros(nx,nz-1);       % difference for dEy/dx
    Eydiffz = zeros(nx-1,nz);       % difference for dEy/dz
    Hxdiffz = zeros(nx-1,nz-1);     % difference for dHx/dz
    Hzdiffx = zeros(nx-1,nz-1);     % difference for dHz/dx
    PEyx = zeros(nx-1,nz-1);        % psi_Eyx (for PML)
    PEyz = zeros(nx-1,nz-1);        % psi_Eyz (for PML)
    PHx = zeros(nx-1,nz);           % psi_Hx (for PML)
    PHz = zeros(nx,nz-1);           % psi_Hz (for PML)
    
    % time stepping loop
    for it=1:numit
        
        % update Hx component...
        
        % determine indices for entire, PML, and interior regions in Hx and property grids
        i = 2:nx-2;  j = 3:nz-2;                % indices for all components in Hx matrix to update
        k = 2*i;  l = 2*j-1;                    % corresponding indices in property grids
        kp = k((k<=kpmlLin | k>=kpmlRin));      % corresponding property indices in PML region
        lp = l((l<=lpmlTin | l>=lpmlBin));
        ki = k((k>kpmlLin & k<kpmlRin));        % corresponding property indices in interior (non-PML) region
        li = l((l>lpmlTin & l<lpmlBin));    
        ip = kp./2;  jp = (lp+1)./2;            % Hx indices in PML region
        ii = ki./2;  ji = (li+1)./2;            % Hx indices in interior (non-PML) region
        
        % update to be applied to the whole Hx grid
        Eydiffz(i,j) = -Ey(i,j+1) + 27*Ey(i,j) - 27*Ey(i,j-1) + Ey(i,j-2);
	    Hx(i,j) = Hx(i,j) - Dbz(k,l).*Eydiffz(i,j);
        
        % update to be applied only to the PML region
        PHx(ip,j) = Bz(kp,l).*PHx(ip,j) + Az(kp,l).*Eydiffz(ip,j);
        PHx(ii,jp) = Bz(ki,lp).*PHx(ii,jp) + Az(ki,lp).*Eydiffz(ii,jp);
        Hx(ip,j) = Hx(ip,j) - Dc(kp,l).*PHx(ip,j);
        Hx(ii,jp) = Hx(ii,jp) - Dc(ki,lp).*PHx(ii,jp);
        
    
        % update Hz component...

        % determine indices for entire, PML, and interior regions in Hz and property grids
        i = 3:nx-2;  j = 2:nz-2;                % indices for all components in Hz matrix to update
        k = 2*i-1;  l = 2*j;                    % corresponding indices in property grids
        kp = k((k<=kpmlLin | k>=kpmlRin));      % corresponding property indices in PML region
        lp = l((l<=lpmlTin | l>=lpmlBin));
        ki = k((k>kpmlLin & k<kpmlRin));        % corresponding property indices in interior (non-PML) region
        li = l((l>lpmlTin & l<lpmlBin));
        ip = (kp+1)./2;  jp = lp./2;            % Hz indices in PML region
        ii = (ki+1)./2;  ji = li./2;            % Hz indices in interior (non-PML) region

        % update to be applied to the whole Hz grid
        Eydiffx(i,j) = -Ey(i+1,j) + 27*Ey(i,j) - 27*Ey(i-1,j) + Ey(i-2,j);
	    Hz(i,j) = Hz(i,j) + Dbx(k,l).*Eydiffx(i,j);
        
        % update to be applied only to the PML region
        PHz(ip,j) = Bx(kp,l).*PHz(ip,j) + Ax(kp,l).*Eydiffx(ip,j);
        PHz(ii,jp) = Bx(ki,lp).*PHz(ii,jp) + Ax(ki,lp).*Eydiffx(ii,jp);
        Hz(ip,j) = Hz(ip,j) + Dc(kp,l).*PHz(ip,j);
        Hz(ii,jp) = Hz(ii,jp) + Dc(ki,lp).*PHz(ii,jp);
        
        
        % update Ey component...
        
        % determine indices for entire, PML, and interior regions in Ey and property grids
        i = 2:nx-2;  j = 2:nz-2;                % indices for all components in Ey matrix to update
        k = 2*i;  l = 2*j;                      % corresponding indices in property grids
        kp = k((k<=kpmlLin | k>=kpmlRin));      % corresponding property indices in PML region
        lp = l((l<=lpmlTin | l>=lpmlBin));
        ki = k((k>kpmlLin & k<kpmlRin));        % corresponding property indices in interior (non-PML) region
        li = l((l>lpmlTin & l<lpmlBin));
        ip = kp./2;  jp = lp./2;                % Ey indices in PML region
        ii = ki./2;  ji = li./2;                % Ey indices in interior (non-PML) region

        % update to be applied to the whole Ey grid
        Hxdiffz(i,j) = -Hx(i,j+2) + 27*Hx(i,j+1) - 27*Hx(i,j) + Hx(i,j-1);
        Hzdiffx(i,j) = -Hz(i+2,j) + 27*Hz(i+1,j) - 27*Hz(i,j) + Hz(i-1,j);
        Ey(i,j) = Ca(k,l).*Ey(i,j) + Cbx(k,l).*Hzdiffx(i,j) - Cbz(k,l).*Hxdiffz(i,j);
        
        % update to be applied only to the PML region
        PEyx(ip,j) = Bx(kp,l).*PEyx(ip,j) + Ax(kp,l).*Hzdiffx(ip,j);
        PEyx(ii,jp) = Bx(ki,lp).*PEyx(ii,jp) + Ax(ki,lp).*Hzdiffx(ii,jp);
        PEyz(ip,j) = Bz(kp,l).*PEyz(ip,j) + Az(kp,l).*Hxdiffz(ip,j);
        PEyz(ii,jp) = Bz(ki,lp).*PEyz(ii,jp) + Az(ki,lp).*Hxdiffz(ii,jp);
        Ey(ip,j) = Ey(ip,j) + Cc(kp,l).*(PEyx(ip,j) - PEyz(ip,j));
        Ey(ii,jp) = Ey(ii,jp) + Cc(ki,lp).*(PEyx(ii,jp) - PEyz(ii,jp));
        
        % add source pulse to Ey at source location
        % (emulates infinitesimal Ey directed line source with current = srcpulse)
        i = srci(s); j = srcj(s);
        Ey(i,j) = Ey(i,j) + srcpulse(it);
    
        
        % plot the Ey wavefield if necessary
        if plotopt(1)==0;   
            if mod(it-1,plotopt(2))==0
                disp(['Source ',num2str(s),'/',num2str(nsrc),', Iteration ',num2str(it),'/',num2str(numit),...
                        ':  t = ',num2str(t(it)*1e9),' ns'])
                figure(1); imagesc(xEy,zEy,Ey'); axis image
                title(['Source ',num2str(s),'/',num2str(nsrc),', Iteration ',num2str(it),'/',num2str(numit),...
                        ':  Ey wavefield at t = ',num2str(t(it)*1e9),' ns']);
                xlabel('Position (m)');  ylabel('Depth (m)');
                caxis([-plotopt(3) plotopt(3)]);
                pause(0.01);
            end
        end
    
        % record the results in gather matrix if necessary
        if mod(it-1,outstep)==0
            tout((it-1)/outstep+1) = t(it);
            for r=1:nrec
                gather((it-1)/outstep+1,r,s) = Ey(reci(r),recj(r));
            end
        end
        
    end
end