smt_solve.m

来自「利用电磁场的源激发方法来计算光子晶体波导例如光子晶体光纤」· M 代码 · 共 52 行

function [vErr, varargout] = smt_solve(Neff, oGd, pos, k0, varargin)
% error = SMT_SOLVE(Neff, oGd, positions, k0)
% Solves the geometry defined in oGd.
% Sources and testing points are distributed according to positions.
% Free-space wave number: k0.
% Effective index: Neff, may be a vector
% Returns the error in continuity conditions.
% [error, solution] = SMT_SOLVE(Neff, oGd, pos, k0)
% Returns the vector of source amplitudes as well.
% If Neff is a vector, solution is a matrix. Each column of the matrix
% corresponds to one entry of the Neff vector.
% If Neff is complex, the function must be called with one more argument:
% [error, solution] = SMT_SOLVE(Neff, oGd, pos, k0, isLeaky).
% isLeaky determines the behavior at infinity. A value of 1 means that
% phase propagation is outwards (to infinity). A value of 0 means that
% phase propagation is inwards (from infinity), as in plasmons.
persistent guess
nTries = 4;
vBeta = Neff*k0;
vErr = zeros(1,length(vBeta));
mSol = zeros(2*sum([pos.sources.n]), length(vBeta));
if ~isreal(Neff)
    isLeaky = varargin{1};
else
    isLeaky = 0; % if Neff is real isLeaky is irrelevant
end
for ind = 1:length(vBeta)
    beta = vBeta(ind);
    [mZ, mBNegate] = imp_matrix(oGd, pos, k0, beta, isLeaky);
    if (length(guess) ~= size(mZ, 2))
        guess = [];
    end
    err = -1;
    iTry = 1;
    while (err == -1) & (iTry < nTries)
        if iTry > 1
            guess = [];
            disp('Arnoldi convergence failure')
        end
        [err, guess] = non_trivial(mZ, mBNegate, guess);
        iTry = iTry + 1;
    end
    vErr(ind) = bc_error(mZ, mBNegate, guess);

    if nargout == 2
        mSol(:,ind) = guess(:);
    end
end
if nargout == 2
    varargout{1} = mSol;
end