% Fd1d2.m Dielectric slabclear allNN = 200;kc = NN/2;k0 = 100;epsz = 8.85419*10^-12;ddx = .01;dt = ddx/(2*3*10^8);t0 = 50;spread = 15;ex = zeros(1,NN);hy = zeros(1,NN);ca = zeros(1,NN);cb = zeros(1,NN);cond = zeros(1,NN);barr = zeros(1,NN);elm1 = 0.;elm2 = 0.;ehm1 = 0.;ehm2 = 0.;% Initializefor k=1:NNca(k) = 1.;cb(k) = 0.5;endepsilon = 5;sigma = 0.01;eaf  = dt*sigma/(2*epsz*epsilon);for k=80:120cond(k) = sigma;ca(k) = (1 -eaf)/(1+eaf);cb(k) = .5/(epsilon*(1+eaf));barr(k) = .5;endfreq = zeros(2,1);real_pt = zeros(2,NN);imag_pt = zeros(2,NN);amp = zeros(2,NN);freq(1) = 50*1e6;freq(2) = 500*1e6;T = 0;n_step = 1;while n_step > 0n_step = input('How many time steps -->');for n=1:n_stepT = T + 1;ex(1) = elm2;elm2 = elm1elm1 = ex(2)for k=2:NNex(k) = ca(k)*ex(k) + cb(k)*( hy(k-1) - hy(k) );endex(NN) = ehm2;ehm2 = ehm1;ehm1 = ex(NN-1)for k=2:NNfor n=1:2real_pt(n,k) = real_pt(n,k) + cos(2*pi*T*dt*freq(n))*ex(k);imag_pt(n,k) = imag_pt(n,k) + sin(2*pi*T*dt*freq(n))*ex(k);endendpulse = exp(-.5*((t0-T)/spread)^2 );%pulse = sin(2*pi*T/50);pulseif T < 120ex(2) = pulse;endfor k=1:NN-1hy(k) =  hy(k) + .5*( ex(k) - ex(k+1) );endendsubplot(3,1,1)plot(ex)hold onplot(barr,'m--')hold offaxis( [ 1 NN -.5 1.1 ])TT = text(150, .7 , ['T = ',num2str(T)]);set(TT,'fontsize',12)TT = text(87, .3 , [' eps = ',num2str(epsilon)]);set(TT,'fontsize',12)TT = text(87, -.3 , ['sig = ',num2str(sigma)]);set(TT,'fontsize',12)set(gca,'fontsize',12)ylabel('E field')for k=1:NNfor n=1:2amp(n,k) = cond(k)*(real_pt(n,k)^2 + imag_pt(n,k)^2);endendsubplot(6,1,3)plot(amp(1,:))set(gca,'fontsize',12)TT = text(130, 3, ['50 MHz']);set(TT,'fontsize',12)axis( [ 1 NN 0 7 ])set(gca,'XTick',[])ylabel('SAR')subplot(6,1,4)plot(amp(2,:))set(gca,'fontsize',12)TT = text(130, 3, ['500 MHz']);axis( [ 1 NN 0 7 ])set(TT,'fontsize',12)ylabel('SAR')saveas(gcf,'pulse.bmp')end