DOS1D.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% this program calculates and plots the total density of photonic modes (DOS) allowed in a 1D photonic
%%% crystal with the unit cell made of two homogeneous dielectric layers;
%%% the calculation is separated in two cases: TE(s) and TM(p)
%%% polarizations; DOS is calculated via numerical integration of symbolic
%%% expressions

%%% Author: Cazimir-Gabriel Bostan, 2008 (Bucharest, Romania)
%%%	http://cgbostan.evonet.ro
%%%	[e-mail]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tic
clear all
warning off

%%% refractive indices of the layers
n1 = 1.4;
n2 = 3.45;

%%% normalized thicknesses of the layers
d1 = n2/(n1+n2);
d2 = 1-d1;

N=1000; % no. of discretization points
u = linspace(0.15, 0.45, N); % vector of normalized frequency
Ds = []; Dp = []; % DOS for "s" and "p" pol.

for j=1:N
    o=u(j);
    
    % for "s" polarization
    Fs = @(x) x.*((n1^2*d1./sqrt(n1^2-x.^2)).*(sin(2*pi*d1*o.*sqrt(n1^2-x.^2)).*cos(2*pi*d2*o.*sqrt(n2^2-x.^2))+0.5.*(sqrt(n1^2-x.^2)./sqrt(n2^2-x.^2)+sqrt(n2^2-x.^2)./sqrt(n1^2-x.^2)).*sin(2*pi*d1*o.*sqrt(n1^2-x.^2)).*cos(2*pi*d2*o.*sqrt(n2^2-x.^2)))+...
        + (n2^2*d2./sqrt(n2^2-x.^2)).*(sin(2*pi*d2*o.*sqrt(n2^2-x.^2)).*cos(2*pi*d1*o.*sqrt(n1^2-x.^2))+0.5.*(sqrt(n2^2-x.^2)./sqrt(n1^2-x.^2)+sqrt(n1^2-x.^2)./sqrt(n2^2-x.^2)).*sin(2*pi*d2*o.*sqrt(n2^2-x.^2)).*cos(2*pi*d1*o.*sqrt(n1^2-x.^2)))-...
        - (1/o)*((n1^2-n2^2)^2/(4*pi)).*x.^2.*sin(2*pi*d1*o.*sqrt(n1^2-x.^2)).*sin(2*pi*d2*o.*sqrt(n2^2-x.^2))./(sqrt(n1^2-x.^2).*sqrt(n2^2-x.^2)).^3)./(1-(cos(2*pi*d1*o.*sqrt(n1^2-x.^2)).*cos(2*pi*d2*o.*sqrt(n2^2-x.^2))-0.5.*(sqrt(n1^2-x.^2)./sqrt(n2^2-x.^2)+...
        + sqrt(n2^2-x.^2)./sqrt(n1^2-x.^2)).*sin(2*pi*d1*o.*sqrt(n1^2-x.^2)).*sin(2*pi*d2*o.*sqrt(n2^2-x.^2))).^2);
    
    Ds(j)=abs(quadl(Fs,0,n1-0.001));
    
    % for "p" polarization
    Fp = @(x) x.*((n1^2*d1./sqrt(n1^2-x.^2)).*(sin(2*pi*d1*o.*sqrt(n1^2-x.^2)).*cos(2*pi*d2*o.*sqrt(n2^2-x.^2))+0.5.*((n2^2/n1^2).*sqrt(n1^2-x.^2)./sqrt(n2^2-x.^2)+(n1^2/n2^2).*sqrt(n2^2-x.^2)./sqrt(n1^2-x.^2)).*sin(2*pi*d1*o.*sqrt(n1^2-x.^2)).*cos(2*pi*d2*o.*sqrt(n2^2-x.^2)))+...
        + (n2^2*d2./sqrt(n2^2-x.^2)).*(sin(2*pi*d2*o.*sqrt(n2^2-x.^2)).*cos(2*pi*d1*o.*sqrt(n1^2-x.^2))+0.5.*((n1^2/n2^2).*sqrt(n2^2-x.^2)./sqrt(n1^2-x.^2)+(n2^2/n1^2).*sqrt(n1^2-x.^2)./sqrt(n2^2-x.^2)).*sin(2*pi*d2*o.*sqrt(n2^2-x.^2)).*cos(2*pi*d1*o.*sqrt(n1^2-x.^2)))-...
        - (1/o)*((n1^2-n2^2)^2*(1/n1^2+1/n2^2-1)/(4*pi)).*x.^2.*sin(2*pi*d1*o.*sqrt(n1^2-x.^2)).*sin(2*pi*d2*o.*sqrt(n2^2-x.^2))./(sqrt(n1^2-x.^2).*sqrt(n2^2-x.^2)).^3)./(1-(cos(2*pi*d1*o.*sqrt(n1^2-x.^2)).*cos(2*pi*d2*o.*sqrt(n2^2-x.^2))-0.5.*((n2^2/n1^2).*sqrt(n1^2-x.^2)./sqrt(n2^2-x.^2)+...
        + (n1^2/n2^2).*sqrt(n2^2-x.^2)./sqrt(n1^2-x.^2)).*sin(2*pi*d1*o.*sqrt(n1^2-x.^2)).*sin(2*pi*d2*o.*sqrt(n2^2-x.^2))).^2);
    
    Dp(j)=abs(quadl(Fp,0,n1-0.0001)); 
    
    display(sprintf('Calculation for o[%d] is finished',j));
    
end

plot(u,Ds); xlabel('Normalized frequency'); ylabel('DOS for s-pol (a.u.)');
figure
plot(u,Dp,'r'); xlabel('Normalized frequency'); ylabel('DOS for p-pol (a.u.)');

toc