fdtd(pml).m

a PML FDTD. computational domain

close all; 
clear all;
%arrays for pressure and velocity
spatialWidth=100; 
temporalWidth=2;
p=zeros(spatialWidth,temporalWidth);
v=zeros(spatialWidth,temporalWidth);



length=2;%waveguide length
dx=length/spatialWidth; %spatial resolution

%Temporal resolution (courant limit)
c=340;
dt=dx/c;

%density and modulous
rho=1.21;
k=rho*c^2;

duration=5; %Simulation Run for 5secs
iterations=duration/dt;

excitationPoint=spatialWidth/2; % ectiation point half way along the waveguide

xPML=10;

for n=2:iterations
    t=n*dt;
    for i=2: spatialWidth-1 
        
        if i>(spatialWidth-xPML)
            xi=xPML-(spatialWidth-i);
            a0=log(10)/(K*dt);
            a1=a0*(xi/xPML)^2;
            a2=a0*((xi-1/2)/xPML)^2;
            p(i,1)=exp(-(a1*K)*dt)*p(i,2)-(1-exp(-(a1*K)*dt))/(a1*K)*K*1/dx*(v(i+1,2)-v(i,2));
            v(i,1)=exp(-(a2*K)*dt)*v(i,2)-(1-exp(-(a2*K)*dt))/(a2*K)*(1/rho)*1/dx*(p(i,1)-p(i-1,1));
        else
            p(i,1)=p(i,2)-K*dt/dx*(v(i+1,2)-v(i,2));
            if i==excitationPoint
                %p(i,1)=cos(2*pi*500*t);
                sigma=0.0005;
                t0=3*sigma;
                fs=exp( - ((t-t0)/0.0005)^2  );
                p(i,1)=p(i,1)+fs;

            end
            v(i,1)=v(i,2)-(1/rho)*dt/dx*(p(i,1)-p(i-1,1));
        end
        p(i,2)=p(i,1);
        v(i,2)=v(i,1);
    end
    plot(p(:,2)) 
    axis([0 spatialWidth -2 2]); 
    frame = getframe();
end