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