tmm.m

用传输矩阵的方法来计算光子晶体带隙

function [Qext Qsca Qabs] = TMM(varmLayers,varmFiles,r,m,MaxOrder,WavlenRange)      
%% [Qext Qsca Qabs] = TMM(varmLayers varmFiles r m ...      % Layers
%                         MaxOrder WavlenRange )            % Control
% copyright yhw2005

N = height_y(r);

%% Read refraction index from the input files.
% varm -- Save all the variant refraction indis.
varm = [];
for j = 1:height_y(varmLayers)
    filename = cell2mat(varmFiles(j));
    spe = load(filename);
    varm = [varm; (interp1(spe(:,1)', spe(:,2)', WavlenRange) +...
        interp1(spe(:,1)', spe(:,3)', WavlenRange)*i)];
end
clear spe

%% Main part of this programme,
% where the core matrix M and the scattering coefficients S and T will be
% calculated.

% Initialnize the normalized scattering cross section and extinction cross
% section, Qsca=Csca/(pi*a^2), and Qext=Cext/(pi*a^2).
Qsca=[];
Qext=[];

% Calculate every wavelength.
for wavlen = WavlenRange,
    % Initialization
    Csca_wav = 0;
    Cext_wav = 0;
    for j = 1:height_y(varmLayers)
        m(N+1-varmLayers(j)) = varm(j, find(WavlenRange==wavlen));
    end
    
    % Calculate every order.
    for n = 1:MaxOrder,
        % Scattering Layer
        k = m(1)*2*pi/wavlen;
        rho = k*r(2);
        
        coeff = sqrt(pi/2/rho);
        letBj = besselj(n+1/2,rho);
        letBy = bessely(n+1/2,rho);
        
        psiD = (rho*besselj(n-1+1/2,rho) - n*letBj)*coeff;
        xiD = (rho*bessely(n-1+1/2,rho) - n*letBy)*coeff;
        bj = letBj*coeff;
        by = letBy*coeff;
        M_S = [ bj           by;
                k/rho*psiD   k/rho*xiD ];
        M_T = [ psiD/rho    xiD/rho;
                -k*bj       -k*by    ];
        
        % Generate the transfer matrix from the scattering layer to the core layer.
        for j = 2:N-1,
            % Initialization
            k = m(j)*2*pi/wavlen;
            
            % Transfer matrix for the jth layer at the outer boundary
            rho = k*r(j);
            
            coeff = sqrt(pi/2/rho);
            letBj = besselj(n+1/2,rho);
            letBy = bessely(n+1/2,rho);
        
            psiD = (rho*besselj(n-1+1/2,rho) - n*letBj)*coeff;
            xiD = (rho*bessely(n-1+1/2,rho) - n*letBy)*coeff;
            bj = letBj*coeff;
            by = letBy*coeff;
            M_So = [ bj                 by;
                     k/rho*psiD         k/rho*xiD ];
            M_To = [ psiD/rho       xiD/rho;
                     -k*bj          -k*by ];
                 
            M_S = (M_So(1,1)*M_So(2,2)-M_So(1,2)*M_So(2,1))*...
                [ M_So(2,2)   -M_So(1,2);
                  -M_So(2,1)  M_So(1,1)   ]*M_S;
            M_T = (M_To(1,1)*M_To(2,2)-M_To(1,2)*M_To(2,1))*...
                [ M_To(2,2)   -M_To(1,2);
                  -M_To(2,1)  M_To(1,1)   ]*M_T;
            % Transfer matrix for the jth layer at the inner boundary
            rho = k*r(j+1);
            
            coeff = sqrt(pi/2/rho);
            letBj = besselj(n+1/2,rho);
            letBy = bessely(n+1/2,rho);
            
            psiD = (rho*besselj(n-1+1/2,rho) - n*letBj)*coeff;
            xiD = (rho*bessely(n-1+1/2,rho) - n*letBy)*coeff;
            bj = letBj*coeff;
            by = letBy*coeff;
            M_Si = [ bj             by;
                     k/rho*psiD     k/rho*xiD ];
            M_Ti = [ psiD/rho   xiD/rho;
                     -k*bj      -k*by   ];
            
            M_S = M_Si*M_S;
            M_T = M_Ti*M_T;
            
        end  % layer
        
        % Generate S and T at wavelength=wavlen and order=n.
        k = m(end)*2*pi/wavlen;
        rho = k*r(end);
        
        coeff = sqrt(pi/2/rho);
        letBj = besselj(n+1/2,rho);
        letBy = bessely(n+1/2,rho);
        
        psiD = (rho*besselj(n-1+1/2,rho) - n*letBj)*coeff;
        xiD = (rho*bessely(n-1+1/2,rho) - n*letBy)*coeff;
        bj = letBj*coeff;
        S = (bj*M_S(2,1) - k/rho*psiD*M_S(1,1)) / ...
            (k/rho*psiD*(M_S(1,1)+i*M_S(1,2)) - bj*(M_S(2,1)+i*M_S(2,2)));
        T = (i/rho*psiD*M_T(2,1) + i*k*bj*M_T(1,1)) / ...
            (k*bj*(M_T(1,1)+i*M_T(1,2)) + 1/rho*psiD*(M_T(2,1)+i*M_T(2,2)));
        
        Csca_wav = Csca_wav + (2*n+1)*(S*conj(S)+T*conj(T));
        Cext_wav = Cext_wav + (2*n+1)*real(-S-i*T);
        
    end  % order
    
    % Save the Qsca and Qext at wavelength=wavlen.
    Qsca = [Qsca Csca_wav*2*pi/real(m(1)*2*pi/wavlen)^2/r(end)^2/pi];
    Qext = [Qext Cext_wav*2*pi/real(m(1)*2*pi/wavlen)^2/r(end)^2/pi];
    
end  % wavelength

Qabs = Qext - Qsca;