myfdtd01cylpw.m

Dennis M. Sullivan书中程序3.3.c

function myFDTD = myFDTD01(time)
max_time = 1000;
max_space_x = 60;
max_space_y = 80;
pulse_start_x = 55;
pulse_start_y = 155;
npml = 5;

eps = 1;
pi = 3.14159;
dx = 0.01;
dt = dx/6e8;
epsz = 8.8e-12;
sigma = 0.0;
t0 = 20;
spread = 8;
R = 5;
ic = 30;
jc = 30;

ia = 10;
ja = 10;
ib = max_space_x-ia+1;
jb = max_space_y-ja+1;

ez_inc_low_m1 = 0;
ez_inc_low_m2 = 0;
ez_inc_high_m1 = 0;
ez_inc_high_m2 = 0;

Dz = zeros(max_space_x,max_space_y);
Ez = Dz;
Hx = Ez;
Hy = Ez;
ez_inc = zeros(max_space_y);
hx_inc = ez_inc;
ihx = Dz;
ihy = Dz;
I = Dz;
gb = Dz;
A = ones(max_space_x,max_space_y);
B = Dz;
ga = A;
gi2 = ones(max_space_x);
gi3 = gi2;
fi1 = zeros(max_space_x);
fi2 = gi2;
fi3 = gi2;

gj2 = ones(max_space_y);
gj3 = gj2;
fj1 = zeros(max_space_y);
fj2 = gj2;
fj3 = gj2;

for j = 1:max_space_y
    for i = 1:max_space_x
        xdist = (ic-i);
        ydist = (jc-j);
        dist = (xdist^2 + ydist^2)^0.5;
        if (dist <= R)
            ga(i,j) = 1/(eps + (sigma*dt)/epsz);
            gb(i,j) = (sigma*dt)/epsz;
        end
    end
end

for i = 1:npml;
    xn = 0.33*([(npml - i+1)/npml]^3);
    gi2(i) = 1/(1+xn);
    gi2(max_space_x - i) = 1/(1+xn);
    gi3(i) = (1 - xn)/(1 + xn);
    gi3(max_space_x - i) = (1 - xn)/(1 + xn);
    
    xn = 0.25*([(npml - i + 0.5)/npml]^3);
    fi1(i) = xn;
    fi1(max_space_x - i - 1) = xn;
    fi2(i) = 1/(1 + xn);
    fi2(max_space_x - i - 1) = 1/(1 + xn);
    fi3(i) = (1 - xn)/(1 + xn);
    fi3(max_space_x - i - 1) = (1 - xn)/(1 + xn);
end

for j = 1:npml;
    xn = 0.33*([(npml - j+1)/npml]^3);
    gj2(j) = 1/(1+xn);
    gj2(max_space_y - j ) = 1/(1+xn);
    gj3(j) = (1 - xn)/(1 + xn);
    gj3(max_space_y - j ) = (1 - xn)/(1 + xn);
    
    xn = 0.25*([(npml - j + 0.5)/npml]^3);
    fj1(j) = xn;
    fj1(max_space_y - j - 1) = xn;
    fj2(j) = 1/(1 + xn);
    fj2(max_space_y - j - 1) = 1/(1 + xn);
    fj3(j) = (1 - xn)/(1 + xn);
    fj3(max_space_y - j - 1) = (1 - xn)/(1 + xn);
end

for n = 1:max_time
    
    for j=2:max_space_y
        ez_inc(j)=ez_inc(j)+0.5*(hx_inc(j-1)-hx_inc(j));
    end
    
    ez_inc(1) = ez_inc_low_m2;
    ez_inc_low_m2 = ez_inc_low_m1;
    ez_inc_low_m1 = ez_inc(2);
    
    ez_inc(max_space_y) = ez_inc_high_m2;
    ez_inc_high_m2 = ez_inc_high_m1;
    ez_inc_high_m1 = ez_inc(max_space_y - 1);
           
    for j = 2:max_space_y
        for i = 2:max_space_x
            Dz(i,j) = gi3(i)*gj3(j)*Dz(i,j) + gi2(i)*gj2(j)*0.5*( Hy(i,j) - Hy(i - 1,j) - Hx(i,j) + Hx(i,j-1) );
        end
    end
    
%     pulse = exp ( -.5*( (t0 - n)/spread )^2);
    pulse = sin(2*pi*1000*1e6*dt*n);
    ez_inc(3)=pulse;
           
    for i = ia:ib
        Dz(i,ja) = Dz(i,ja)+0.5*hx_inc(ja-1);
        Dz(i,jb) = Dz(i,jb)-0.5*hx_inc(jb);
    end
    
    for j = 2:max_space_y
        for i = 2:max_space_x
            Ez(i,j) = ga(i,j)*(Dz(i,j) - I(i,j));
            I(i,j) = I(i,j) + gb(i,j)*Ez(i,j);
        end
    end
    
    for j = 1:max_space_y-1
        Ez(1,j) = 0;
        Ez(max_space_x,j) = 0;
    end
    
    for i = 1:max_space_x-1
        Ez(i,1) = 0;
        Ez(i,max_space_y) = 0;
    end
    
    for j=1:max_space_y-1;
        hx_inc(j)=hx_inc(j)+0.5*(ez_inc(j)-ez_inc(j+1));
    end
    
    for j = 1:max_space_y - 1
        for i = 1:max_space_x
            curl_e = Ez(i,j) - Ez(i,j+1);
            ihx(i,j) = ihx(i,j) + fi1(i)*curl_e;
            Hx(i,j) = fj3(j)*Hx(i,j) + fj2(j)*0.5*( curl_e + ihx(i,j) );
        end
    end
    
    for i=ia:ib
        Hx(i,ja-1)=Hx(i,ja-1)+0.5*ez_inc(ja);
        Hx(i,jb)=Hx(i,jb)-0.5*ez_inc(jb);
    end
    
    for j = 1:max_space_y
        for i = 1:max_space_x - 1
            curl_e = Ez(i +1,j) - Ez(i,j);
            ihy(i,j) = ihy(i,j) + fj1(j)*curl_e;
            Hy(i,j) = fi3(i)*Hy(i,j) + fi2(i)*0.5*( curl_e + ihy(i,j) );
        end
    end   
    
    for j=ja:jb
        Hy(ia-1,j)=Hy(ia-1,j)-0.5*ez_inc(j);
        Hy(ib,j)=Hy(ib,j)+0.5*ez_inc(j);
    end
            
    PP=Ez(:,20);
    figure (1)
    mesh(Ez)
    axis([1 max_space_y 1 max_space_x -1.2 1.2])
%     axis([1 max_space_x -1 1 -1 1])
%  plot(PP)
    drawnow
%     pause(0.05)
end