imp_matrix.m

利用电磁场的源激发方法来计算光子晶体波导例如光子晶体光纤

function [mZ, varargout] = imp_matrix(oGd, pos, k, beta, isLeaky)
global EPSILON0 MU0 C ETA

w = k*C;
sdArray = oGd.sdArray;
mZ = zeros(sum([pos.testing.n])*4, sum([pos.sources.n])*2);
mBNegate = logical(mZ);
% these vectors give the edge indices of each region in the matrix
vISources = [0 cumsum(2*[pos.sources.n])];
vITesting = [0 cumsum([pos.testing.n])];
indEz = size(mZ,1)/4;
indHt = 2*indEz;
indHz = 3*indEz;
Z = j*w*MU0;
for iSd = 1:oGd.nSd
    if pos.sources(iSd).n == 0
        continue
    end
    Er = sdArray(iSd).Er;
    if ischar(Er)
        Er = feval(Er, 2*pi/k);
    end
    Y = j*EPSILON0*Er*w;
    ETA_sd2 = Z/Y;
    isHankel = (iSd == 1);
    kc = calc_kc(k, beta, Er, isHankel, isLeaky);
    % testing points on the 'inside' curve
    curve = sdArray(iSd).insideCurve;
    if ~isempty(curve)
        iCurve = pos.testing(curve.iCurve).iCurve;
        [Fe, Ee, Be, Fm, Em, Bm] = green_fun(pos.testing(iCurve), pos.sources(iSd), Y, kc, beta, oGd);
        mZeros = zeros(size(Be));
        rows = vITesting(iCurve)+1:vITesting(iCurve+1);
        cols = vISources(iSd)+1:vISources(iSd+1);
        mZ(rows, cols) = [-ETA*Fm Ee];
        mZ(rows+indEz, cols) = [mZeros Be];
        mZ(rows+indHt, cols) = [ETA^2/ETA_sd2*Em ETA*Fe];
        mZ(rows+indHz, cols) = [ETA^2/ETA_sd2*Bm mZeros];
    end
    % testing points on the 'outside' curves
    for iOutsideCurve = 1:length(sdArray(iSd).outsideCurveArray)
        curve = sdArray(iSd).outsideCurveArray(iOutsideCurve);
        if ~isempty(curve)
            if pos.testing(curve.iCurve).n == 0
                continue
            end
            iCurve = pos.testing(curve.iCurve).iCurve;
            [Fe, Ee, Be, Fm, Em, Bm] = green_fun(pos.testing(iCurve), pos.sources(iSd), Y, kc, beta, oGd);
            mZeros = zeros(size(Be));
            rows = vITesting(iCurve)+1:vITesting(iCurve+1);
            cols = vISources(iSd)+1:vISources(iSd+1);
            mZ(rows, cols) = -[-ETA*Fm Ee];
            mZ(rows+indEz, cols) = -[mZeros Be];
            mZ(rows+indHt, cols) = -[ETA^2/ETA_sd2*Em ETA*Fe];
            mZ(rows+indHz, cols) = -[ETA^2/ETA_sd2*Bm mZeros];
            mBNegate([rows rows+indEz rows+indHt rows+indHz], cols) = 1;
        end
    end
end

%mZ = sparse(mZ);
if nargout == 2
    varargout(1) = {mBNegate};
end