clc;clear;
%C:\Users\Administrator\Documents\MATLAB\3d_data
dizhi='E:\serialresult\';%%
mmm='_20181001';
Tsteps=60;
ie=60;je=60;ke=60;
is=ie/2;js=je/2;ks=ke/2;
cc=3e8;f0=1.5e9;
deltex=0.01;
deltet=deltex/2/cc;
t0=20;tau=6;
% nmax=floor(/deltet);
rate=1;
len_xishu=Tsteps/rate;
JJ=1;
gax=ones(ie,je,ke);gay=gax;gaz=gax;
dx=zeros(ie,je,ke);dy=dx;dz=dx;
ex=zeros(ie,je,ke);ey=ex;ez=ex;
hx=zeros(ie,je,ke);hy=hx;hz=hx;
gaz(is,js,16:45)=0;gaz(is,js,ks)=1;
data_ex=zeros(len_xishu,ie,je,ke);data_ey=zeros(len_xishu,ie,je,ke);data_ez=zeros(len_xishu,ie,je,ke);
data_hx=zeros(len_xishu,ie,je,ke); data_hy=zeros(len_xishu,ie,je,ke); data_hz=zeros(len_xishu,ie,je,ke);
for t=1:Tsteps
for k=2:ke
for j=2:je
for i=2:ie
dx(i,j,k)=dx(i,j,k)+0.5*(hz(i,j,k)-hz(i,j-1,k)...
-hy(i,j,k)+hy(i,j,k-1));
dy(i,j,k)=dy(i,j,k)+0.5*(hx(i,j,k)-hx(i,j,k-1)...
-hz(i,j,k)+hz(i-1,j,k));
dz(i,j,k)=dz(i,j,k)+0.5*(hy(i,j,k)-hy(i-1,j,k)...
-hx(i,j,k)+hx(i,j-1,k));
end
end
end
pulse=exp(-0.5*((t0-t)/tau)^2);
dz(is,js,ks)=pulse;
% pulse *cos(2*pi*t)
for k=2:ke-1
for j=2:je-1
for i=2:ie-1
ex(i,j,k)=gax(i,j,k)*dx(i,j,k);
ey(i,j,k)=gay(i,j,k)*dy(i,j,k);
ez(i,j,k)=gaz(i,j,k)*dz(i,j,k);
end
end
end
for k=2:ke-1
for j=2:je-1
for i=2:ie
hx(i,j,k)=hx(i,j,k)+0.5*(ey(i,j,k+1)...
-ey(i,j,k)-ez(i,j+1,k)+ez(i,j,k));
end
end
end
for k=2:ke-1
for j=2:je
for i=2:ie-1
hy(i,j,k)=hy(i,j,k)+0.5*(ez(i+1,j,k)...
-ez(i,j,k)-ex(i,j,k+1)+ex(i,j,k));
end
end
end
for k=2:ke
for j=2:je-1
for i=2:ie-1
hz(i,j,k)=hz(i,j,k)+0.5*(ex(i,j+1,k)...
-ex(i,j,k)-ey(i+1,j,k)+ez(i,j,k));
end
end
end
% if t==30
% for i=1:ie
% for j=1:je
% ezkout(i,j)=ez(i,j,ks-5);
% end
% end
% [i,j]=meshgrid(1:1:ie,1:1:ke);
% % subplot(2,2,1);
% mesh(i,j,ezkout)
% surf(i,j,ezkout)
% axis([1,ie,1,je,-0.02,0.03]);
% title('n=30');
%end
% if t==40
% for i=1:ie
% for j=1:je
% ezkout(i,j)=ez(i,j,ks-5);
% end
% end
% [i,j]=meshgrid(1:1:ie,1:1:ke);
% % subplot(2,2,2);
% mesh(i,j,ezkout)
% surf(i,j,ezkout)
% axis([1,ie,1,je,-0.02,0.03]);
% title('n=40');
% end
% if t==50
% for i=1:ie
% for j=1:je
% ezkout(i,j)=ez(i,j,ks-5);
% end
% end
% [i,j]=meshgrid(1:1:ie,1:1:ke);
% subplot(2,2,3);
% mesh(i,j,ezkout)
% surf(i,j,ezkout)
% axis([1,ie,1,je,-0.02,0.03]);
% title('n=50');
% end
% if t==60
% for i=1:ie
% for j=1:je
% ezkout(i,j)=ez(i,j,ks-5);
% end
% end
% [i,j]=meshgrid(1:1:ie,1:1:ke);
% subplot(2,2,4);
% mesh(i,j,ezkout)
% surf(i,j,ezkout)
% axis([1,ie,1,je,-0.02,0.03]);
% title('n=60');
% end
x=1:1:je;y=1:1:ke;
if mod(t,rate)==0
JJ=t/rate;
data_ex(JJ,:,:,:)=ex(1:ie,x,y);data_ey(JJ,:,:,:)=ey(1:ie,x,y);data_ez(JJ,:,:,:)=ez(1:ie,x,y);
data_hx(JJ,:,:,:)=hx(1:ie,x,y);data_hy(JJ,:,:,:)=hy(1:ie,x,y);data_hz(JJ,:,:,:)=hz(1:ie,x,y);
end
save([dizhi,'_','ex',mmm,'.mat'],'data_ex');
save([dizhi,'_','ey',mmm,'.mat'],'data_ey');
save([dizhi,'_','ez',mmm,'.mat'],'data_ez');
save([dizhi,'_','hx',mmm,'.mat'],'data_hx');
save([dizhi,'_','hy',mmm,'.mat'],'data_hy');
save([dizhi,'_','hz',mmm,'.mat'],'data_hz');
end
评论0