clear all; format compact; tic
%function ex3_8c_makeficient
% APPLICATION OF THE FINITE DIFFERENCE METHOD
% This program involves the penetration of a lossless dielectric SPHERE
% by a plane wave. The program provides in the maximum absolute value of
% Ey and Ez during the final half-wave of time-stepping
% Assumption:
% +y-directed incident wave with components Ez and Hx.
% I,J,K,NN correspond to X,Y,Z, and Time.
% IMAX,JMAX,KMAX are the maximum values of x,y,z
% NNMAX is the total number of timesteps.
% NHW represents one half-wave cycle.
% MED is the number of different uniform media sections.
% JS is the j-position of the plane wave front.
% THIS PROGRAM WAS DEVELOPED BY V. BEMMEL [43]
% AND LATER IMPROVED BY D. TERRY
IMAX=19; JMAX=39; KMAX=19;
NMAX=2; NNMAX=500; NHW=40; MED=2; JS=3;
DELTA=3E-3; CL=3.0E8; F=2.5E9;
% Define scatterer dimensions
OI=19.5; OJ=20.0; OK=19.0; RADIUS=15.0;
ER=[1.0,4.0]; % CONSTITUTIVE PARAMETERS
SIG=[0.1,0.0];
% Statement function to compute position w.r.t. center of the sphere
E0=(1E-9)/(36*pi);
U0=(1E-7)*4*pi;
DT=DELTA/(2*CL);
R=DT/E0;
RA=(DT^2)/(U0*E0*(DELTA^2));
RB=DT/(U0*DELTA);
TPIFDT = 2.0*pi*F*DT;
%********************************************************
% EP # 1 - COMPUTE MEDIA PARAMETERS
%********************************************************
CA = 1-R*SIG./ER;
CB = RA./ER;
CBMRB = CB/RB;
% (i) CALCULATE THE REAL/ACTUAL GRID POINTS
% Initialize the media arrays.Index (M) determines which
% medium each point is actually located in and is used to
% index into arrays which determine the constitutive
% parameters of the medium.There are separate M determining
% arrays for EX, EY, and EZ. These arrays correlate the
% integer values of I,J,K to the actual position within
% the lattice. Computing these values now and storing them in these
% arrays as opposed to computing them each time they are
% needed saves a large amount of computation time.
x = 0:(IMAX+1); y = 0:(JMAX+1); z = 0:(KMAX+1);
[Mx,My,Mz]=ndgrid(x,y,z);
IXMED = (sqrt((Mx-OI+.5).^2+(My-OJ).^2+(Mz-OK).^2)<=RADIUS)+1;
IYMED = (sqrt((Mx-OI).^2+(My-OJ+.5).^2+(Mz-OK).^2)<=RADIUS)+1;
IZMED = (sqrt((Mx-OI).^2+(My-OJ).^2+(Mz-OK+.5).^2)<=RADIUS)+1;
% *************************************************
% STEP # 2 - INITIALIZE FIELD COMPONENTS
% *************************************************
% components for output
EY1 = zeros(1,JMAX+2);
EZ1 = zeros(1,JMAX+2);
EX=zeros(IMAX+2,JMAX+2,KMAX+2,NMAX+1);
EY=zeros(IMAX+2,JMAX+2,KMAX+2,NMAX+1);
EZ=zeros(IMAX+2,JMAX+2,KMAX+2,NMAX+1);
HX=zeros(IMAX+2,JMAX+2,KMAX+2,NMAX+1);
HY=zeros(IMAX+2,JMAX+2,KMAX+2,NMAX+1);
HZ=zeros(IMAX+2,JMAX+2,KMAX+2,NMAX+1);
% ********************************************************
% STEP # 3 - USE FD/TD ALGORITHM TO GENERATE
% FIELD COMPONENTS
% ********************************************************
% SINCE ONLY FIELD COMPONENTS AT CURRENT TIME (t) AND PREVIOUS
% TWO TIME STEPS ( t-1 AND t-2) ARE REQUIRED FOR COMPUTATION,
% WE SAVE MEMORY SPACE BY USING THE FOLLOWING INDICES
% NCUR is index in for time t
% NPR1 is index in for t-1
% NPR2 is index in for t-2
% NOTES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% *ind03c.m I incremented the time so it goes 1 2 3, instead of 0 1 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NCUR = 3;
NPR1 = 2;
NPR2 = 1;
for NN = 1: NNMAX % TIME LOOP
if mod(NN,10)==0
disp(['NN = ',num2str(NN)]) % DISPLAY PROGRESS
end
% Next time step - move indices up a notch.
NPR2 = NPR1;
NPR1 = NCUR;
NCUR = mod( NCUR, 3)+1;
for K=0:KMAX % Z LOOP
for J=0:JMAX % Y LOOP
for I=0:IMAX % X LOOP
% (ii)-APPLY SOFT LATTICE TRUNCATION CONDITIONS
%---x=delta/2
if (I==0)
if ((K~=KMAX)&&(K~=0))
HY(0+1,J+1,K+1,NCUR) = (HY(1+1,J+1,K-1+1,NPR2) + HY(1+1,J+1,K+1,NPR2)+ HY(1+1,J+1,K+1+1,NPR2))/3;
HZ(0+1,J+1,K+1,NCUR) = (HZ(1+1,J+1,K-1+1,NPR2) + HZ(1+1,J+1,K+1,NPR2)+ HZ(1+1,J+1,K+1+1,NPR2))/3;
else
if (K==KMAX)
HY(0+1,J+1,KMAX+1,NCUR) = (HY(1+1,J+1,KMAX-1+1,NPR2)+ HY(1+1,J+1,KMAX+1,NPR2))/2;
HZ(0+1,J+1,K+1,NCUR)=( HZ(1+1,J+1,K-1+1,NPR2)+ HZ(1+1,J+1,K+1,NPR2) )/2;
else
HY(0+1,J+1,K+1,NCUR) = ( HY(1+1,J+1,K+1,NPR2)+ HY(1+1,J+1,K+1+1,NPR2))/2;
HZ(0+1,J+1,0+1,NPR2)=(HZ(1+1,J+1,0+1,NPR2)+ HZ(1+1,J+1,1+1,NPR2))/2;
end
end
end
% ---y=0
if (J==0)
EX(I+1,0+1,K+1,NCUR)=EX(I+1,1+1,K+1,NPR2);
EZ(I+1,0+1,K+1,NCUR)=EZ(I+1,1+1,K+1,NPR2);
else
%---y=ymax
if (J==JMAX)
EX(I+1,JMAX+1,K+1,NCUR)=EX(I+1,JMAX-1+1,K+1,NPR2);
EZ(I+1,JMAX+1,K+1,NCUR)=EZ(I+1,JMAX-1+1,K+1,NPR2);
end
end
%---z=0
if(K==0)
if ((I~=0)&&(I~=IMAX))
EX(I+1,J+1,0+1,NCUR) = (EX(I-1+1,J+1,1+1,NPR2) + EX(I+1,J+1,1+1,NPR2)+EX(I+1+1,J+1,1+1,NPR2))/3;
EY(I+1,J+1,0+1,NCUR) = (EY(I-1+1,J+1,1+1,NPR2) + EY(I+1,J+1,1+1,NPR2)+EY(I+1+1,J+1,1+1,NPR2))/3;
else
if (I==0)
EX(0+1,J+1,0+1,NCUR)=(EX(0+1,J+1,1+1,NPR2)+EX(1+1,J+1,1+1,NPR2))/2;
EY(I+1,J+1,0+1,NCUR)=(EY(I+1,J+1,1+1,NPR2)+EY(I+1+1,J+1,1+1,NPR2))/2;
else
EX(I+1,J+1,0+1,NCUR)=(EX(I-1+1,J+1,1+1,NPR2)+EX(I+1,J+1,1+1,NPR2))/2;
EY(I+1,J+1,0+1,NCUR)=(EY(I-1+1,J+1,1+1,NPR2)+EY(I+1,J+1,1+1,NPR2))/2;
end
end
end
% (iii) APPLY FD/TD ALGORITHM
%-----a. HX generation:
HX(I+1,J+1,K+1,NCUR)=HX(I+1,J+1,K+1,NPR1)+RB*(EY(I+1,J+1,K+1+1,NPR1)- EY(I+1,J+1,K+1,NPR1)+EZ(I+1,J+1,K+1,NPR1)-EZ(I+1,J+1+1,K+1,NPR1));
%-----b. HY generation:
HY(I+1,J+1,K+1,NCUR)=HY(I+1,J+1,K+1,NPR1)+RB*(EZ(I+1+1,J+1,K+1,NPR1)- EZ(I+1,J+1,K+1,NPR1)+EX(I+1,J+1,K+1,NPR1)-EX(I+1,J+1,K+1+1,NPR1));
%-----c. HZ generation:
HZ(I+1,J+1,K+1,NCUR)=HZ(I+1,J+1,K+1,NPR1)+RB*(EX(I+1,J+1+1,K+1,NPR1)- EX(I+1,J+1,K+1,NPR1)+EY(I+1,J+1,K+1,NPR1)-EY(I+1+1,J+1,K+1,NPR1));
%---k=kmax ! SYMMETRY
if (K==KMAX)
HX(I+1,J+1,KMAX+1,NCUR)=HX(I+1,J+1,KMAX-1+1,NCUR);
HY(I+1,J+1,KMAX+1,NCUR)=HY(I+1,J+1,KMAX-1+1,NCUR);
end
% -----d. EX generation:
if ((J~=0)&&(J~=JMAX)&&(K~=0))
M = IXMED( I+1, J+1, K+1 );
EX(I+1,J+1,K+1,NCUR) = CA(M)*EX(I+1,J+1,K+1,NPR1) + CBMRB(M)*(HZ(I+1,J+1,K+1,NCUR)-HZ(I+1,J-1+1,K+1,NCUR)+HY(I+1,J+1,K-1+1,NCUR)-HY(I+1,J+1,K+1,NCUR));
end
%-----e. EY generation:
if(K~=0)
M = IYMED( I+1, J+1, K+1 );
if I~=0
EY(I+1,J+1,K+1,NCUR)=CA(M)*EY(I+1,J+1,K+1,NPR1) + CBMRB(M)*(HX(I+1,J+1,K+1,NCUR)-HX(I+1,J+1,K-1+1,NCUR)+HZ(I-1+1,J+1,K+1,NCUR)-HZ(I+1,J+1,K+1,NCUR));
else
EY(I+1,J+1,K+1,NCUR)=CA(M)*EY(I+1,J+1,K+1,NPR1) + CBMRB(M)*(HX(I+1,J+1,K+1,NCUR)-HX(I+1,J+1,K-1+1,NCUR)+ 0 -HZ(I+1,J+1,K+1,NCUR));
end
end
%-----f. EZ generation:
if ((J~=0)&&(J~=JMAX))
M = IZMED( I+1, J+1, K+1 );
% sig(ext)=0 for Ez only from Taflove[14]
if(M==1)
CAM=1;
else
CAM=CA(M);
end
if I~=0
EZ(I+1,J+1,K+1,NCUR)=CAM*EZ(I+1,J+1,K+1,NPR1)+CBMRB(M)*(HY(I+1,J+1,K+1,NCUR)-HY(I-1+1,J+1,K+1,NCUR)+HX(I+1,J-1+1,K+1,NCUR)-HX(I+1,J+1,K+1,NCUR));
else
EZ(I+1,J+1,K+1,NCUR)=CAM*EZ(I+1,J+1,K+1,NPR1)+CBMRB(M)*(HY(I+1,J+1,K+1,NCUR)- 0 +HX(I+1,J-1+1,K+1,NCUR)-HX(I+1,J+1,K+1,NCUR));
end
% (iv) APPLY THE PLANE-WAVE SOURCE
if (J==JS)
EZ(I+1,JS+1,K+1,NCUR) = EZ(I+1,JS+1,K+1,NCUR)+sin( TPIFDT*NN );
end
end
%---i=imax+1/2 ! SYMMETRY
if(I==IMAX)
EY(IMAX+1+1,J+1,K+1,NCUR)=EY(IMAX+1,J+1,K+1,NCUR);
EZ(IMAX+1+1,J+1,K+1,NCUR)=EZ(IMAX+1,J+1,K+1,NCUR);
end
%---k=kmax
if(K==KMAX)
EX(I+1,J+1,KMAX+1+1,NCUR)=EX(I+1,J+1,KMAX-1+1,NCUR);
EY(I+1,J+1,KMAX+1+1,NCUR)=EY(I+1,J+1,KMAX-1+1,NCUR);
end
end % X LOOP
% ********************************************************
% STEP # 4 - RETAIN THE MAXIMUM ABSOLUTE VALUES DURING
% THE LAST HALF-WAVE
% ********************************************************
if ( (K==KMAX)&&(NN>(NNMAX-NHW)) )
TEMP = abs( EY(IMAX+1,J+1,KMAX-1+1,NCUR) );
if (TEMP > EY1(J+1) )
EY1(J+1) = TEMP;
end
TEMP = abs( EZ(IMAX+1,J+1,KMAX+1,NCUR) );
if (TEMP > EZ1(J+1) )
EZ1(J+1) = TEMP;
end
end
end % Y LOOP
end % Z LOOP
end % TIME LOOP
% *******************************************************
toc
figure(3),plot(6:34,EY1(6:34),'.-')
ylabel('Computed |E_y|/|E_i_n_c|')
xlabel('j')
grid on
figure(4),plot(5:34,EZ1(5:34),'.-')
ylabel('Computed |E_z|/|E_i_n_c|')
xlabel('j')
grid on
کد متلب روش های عددی، مثال های کتاب سادیکو، فصل 3، مثال 3-10
,% ,the ,of ,j ,and , ,of the ,is the ,x y , end
درباره این سایت