electric field distribution of 1D photonic crystal !!!!

%% 
%% Method:   Transfer Matrix Method (TMM) 
%% 
%% Incidence angle is 0 (not distinguished TE(S) and TM(P))  
%% 
%% every layer correspondence to one (2x2) matrix (AA,BB,CC....) 
%% the total matrix can be written as (AA*BB...)^n      
%% n express the period number of 1-dimensional photonic crystal 
%% 
       
clf; x
clear; 
 
fid1=fopen('normfield.inp','r'); 
wave=fscanf(fid1,'%f',1); 
nn=fscanf(fid1,'%f',1); 
nep=fscanf(fid1,'%f',1); 
 
try 
ne(1:nep)=0; 
nd(1:nep)=fscanf(fid1,'%f',nep); 
ne(1:nep)=fscanf(fid1,'%f',nep); 
status=fclose(fid1); 
catch 
   disp('The number of layer is not fit to the detail result.'); 
   status=fclose(fid1); 
   return 
end 
 
for nj=1:nep 
   kk(nj)=2.0*pi*ne(nj)*wave; 
end 
 
MM=[1,0;0,1]; 
 
for n=1:nep-1 
 
   cc=0.5*(1+kk(n)./kk(n+1)); 
   dd=0.5*(1-kk(n)./kk(n+1)); 
	BB=[0,0;0,0];     
	BB(1,1)=cc.*exp(i*kk(n).*nd(n));     
	BB(1,2)=dd.*exp(-i*kk(n).*nd(n)); 
	BB(2,1)=dd.*exp(i*kk(n).*nd(n)); 
	BB(2,2)=cc.*exp(-i*kk(n).*nd(n));     
 
	MM=BB*MM; 
 
	TT(:,:,n)=MM; 
 
end 
 
a0=1; 
b0=-MM(2,1)/MM(2,2) 
 
x0=0.0; 
m3=1; 
 
for m1=1:nep-2 
   x0=x0+nd(m1); 
   for m2=1:nn 
      nni=real(m2-1); 
      xi=nd(m1+1)*nni/nn; 
      x=x0+xi; 
      eij=(TT(1,1,m1)*a0+TT(1,2,m1)*b0)*exp(i*kk(m1+1).*xi)+... 
         (TT(2,1,m1)*a0+TT(2,2,m1)*b0)*exp(-i*kk(m1+1).*xi); 
    
   eabs=abs(eij); 
   ereal=real(eij); 
   eimag=imag(eij); 
    
RRR(m3)=eabs; 
TTT(m3)=ereal; 
AAA(m3)=eimag; 
lam1(m3)=x; 
 
m3=m3+1; 
end 
 
yy=[lam1;RRR;TTT;AAA]; 
fid1=fopen('normfield.dat','w'); 
fprintf(fid1,'%16.8f%16.8f%16.8f%16.8f\n',yy); 
status=fclose(fid1); 
 
end 
subplot(3,1,1); 
plot(lam1,RRR) 
subplot(3,1,2); 
plot(lam1,TTT) 
subplot(3,1,3); 
plot(lam1,AAA)