te_model2d.m

Back Projection imaging through wall

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 sigxmax sigzmax 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
    Hy = zeros(nx-1,nz-1);          % Hy component of magnetic field
    Ex = zeros(nx-1,nz);            % Ex component of electric field
    Ez = zeros(nx,nz-1);            % Ez component of electric field
    Hydiffx = zeros(nx,nz-1);       % difference for dHy/dx
    Hydiffz = zeros(nx-1,nz);       % difference for dHy/dz
    Exdiffz = zeros(nx-1,nz-1);     % difference for dEx/dz
    Ezdiffx = zeros(nx-1,nz-1);     % difference for dEz/dx
    PHyx = zeros(nx-1,nz-1);        % psi_Hyx (for PML)
    PHyz = zeros(nx-1,nz-1);        % psi_Hyz (for PML)
    PEx = zeros(nx-1,nz);           % psi_Ex (for PML)
    PEz = zeros(nx,nz-1);           % psi_Ez (for PML)
    
    % time stepping loop
    for it=1:numit
        
        % update Hy component...
        
        % determine indices for entire, PML, and interior regions in Hy and property grids
        i = 2:nx-2;  j = 2:nz-2;                % indices for all components in Hy 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;                % Hy indices in PML region
        ii = ki./2;  ji = li./2;                % Hy indices in interior (non-PML) region

        % update to be applied to the whole Hy grid
        Exdiffz(i,j) = -Ex(i,j+2) + 27*Ex(i,j+1) - 27*Ex(i,j) + Ex(i,j-1);
        Ezdiffx(i,j) = -Ez(i+2,j) + 27*Ez(i+1,j) - 27*Ez(i,j) + Ez(i-1,j);
        Hy(i,j) = Hy(i,j) + Dbz(k,l).*Exdiffz(i,j) - Dbx(k,l).*Ezdiffx(i,j);
        
        % update to be applied only to the PML region
        PHyx(ip,j) = Bx(kp,l).*PHyx(ip,j) + Ax(kp,l).*Ezdiffx(ip,j);
        PHyx(ii,jp) = Bx(ki,lp).*PHyx(ii,jp) + Ax(ki,lp).*Ezdiffx(ii,jp);
        PHyz(ip,j) = Bz(kp,l).*PHyz(ip,j) + Az(kp,l).*Exdiffz(ip,j);
        PHyz(ii,jp) = Bz(ki,lp).*PHyz(ii,jp) + Az(ki,lp).*Exdiffz(ii,jp);
        Hy(ip,j) = Hy(ip,j) + Dc(kp,l).*(PHyz(ip,j) - PHyx(ip,j));
        Hy(ii,jp) = Hy(ii,jp) + Dc(ki,lp).*(PHyz(ii,jp) - PHyx(ii,jp));
        
        
        % update Ex component...
        
        % determine indices for entire, PML, and interior regions in Ex and property grids
        i = 2:nx-2;  j = 3:nz-2;                % indices for all components in Ex 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;            % Ex indices in PML region
        ii = ki./2;  ji = (li+1)./2;            % Ex indices in interior (non-PML) region
        
        % update to be applied to the whole Ex grid
        Hydiffz(i,j) = -Hy(i,j+1) + 27*Hy(i,j) - 27*Hy(i,j-1) + Hy(i,j-2);
	    Ex(i,j) = Ca(k,l).*Ex(i,j) + Cbz(k,l).*Hydiffz(i,j);
        
        % update to be applied only to the PML region
        PEx(ip,j) = Bz(kp,l).*PEx(ip,j) + Az(kp,l).*Hydiffz(ip,j);
        PEx(ii,jp) = Bz(ki,lp).*PEx(ii,jp) + Az(ki,lp).*Hydiffz(ii,jp);
        Ex(ip,j) = Ex(ip,j) + Cc(kp,l).*PEx(ip,j);
        Ex(ii,jp) = Ex(ii,jp) + Cc(ki,lp).*PEx(ii,jp);
        
    
        % update Ez component...

        % determine indices for entire, PML, and interior regions in Ez and property grids
        i = 3:nx-2;  j = 2:nz-2;                % indices for all components in Ez 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;            % Ez indices in PML region
        ii = (ki+1)./2;  ji = li./2;            % Ez indices in interior (non-PML) region

        % update to be applied to the whole Ez grid
        Hydiffx(i,j) = -Hy(i+1,j) + 27*Hy(i,j) - 27*Hy(i-1,j) + Hy(i-2,j);
	    Ez(i,j) = Ca(k,l).*Ez(i,j) - Cbx(k,l).*Hydiffx(i,j);
        
        % update to be applied only to the PML region
        PEz(ip,j) = Bx(kp,l).*PEz(ip,j) + Ax(kp,l).*Hydiffx(ip,j);
        PEz(ii,jp) = Bx(ki,lp).*PEz(ii,jp) + Ax(ki,lp).*Hydiffx(ii,jp);
        Ez(ip,j) = Ez(ip,j) - Cc(kp,l).*PEz(ip,j);
        Ez(ii,jp) = Ez(ii,jp) + Cc(ki,lp).*PEz(ii,jp);
        
        

        % add source pulse to Ex or Ez at source location
        % (emulates infinitesimal Ex or Ez directed line source with current = srcpulse)
        i = srci(s); j = srcj(s);
        if srcloc(s,3)==1
            Ex(i,j) = Ex(i,j) + srcpulse(it);
        elseif srcloc(s,3)==2
            Ez(i,j) = Ez(i,j) + srcpulse(it);
        end
    
        
        % plot the Ex or Ez wavefield if necessary
        if plotopt(1)==1;   
            if mod(it-1,plotopt(3))==0
                disp(['Source ',num2str(s),'/',num2str(nsrc),', Iteration ',num2str(it),'/',num2str(numit),...
                        ':  t = ',num2str(t(it)*1e9),' ns'])
                if plotopt(2)==1
                    figure(1); imagesc(xEx,zEx,Ex'); axis image
                    title(['Source ',num2str(s),'/',num2str(nsrc),', Iteration ',num2str(it),'/',num2str(numit),...
                            ':  Ex wavefield at t = ',num2str(t(it)*1e9),' ns']);
                    xlabel('Position (m)');  ylabel('Depth (m)');
                    caxis([-plotopt(4) plotopt(4)]);
                    pause(0.01);
                elseif plotopt(2)==2
                    figure(1); imagesc(xEz,zEz,Ez'); axis image
                    title(['Source ',num2str(s),'/',num2str(nsrc),', Iteration ',num2str(it),'/',num2str(numit),...
                            ':  Ez wavefield at t = ',num2str(t(it)*1e9),' ns']);
                    xlabel('Position (m)');  ylabel('Depth (m)');
                    caxis([-plotopt(4) plotopt(4)]);
                    pause(0.01);
                end
            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
                if recloc(r,3)==1
                    gather((it-1)/outstep+1,r,s) = Ex(reci(r),recj(r));
                elseif recloc(r,3)==2
                    gather((it-1)/outstep+1,r,s) = Ez(reci(r),recj(r));
                end
            end
        end
        
    end
end