calc_field.m

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

function mF = calc_field(vSol, pos, oGd, k0, neff, sComponent, strClass, x, y, varargin)
% field = CALC_FIELD(solution, positions, oGd, k0, Neff, sComponent,...
%                 sClass, x, y).
% Calculates a field component (or longitudinal component of Poynting's
% vector). Argument solution is the source amplitudes, positions is the
% positions of sources and testing points, oGd, the geometry descriptor, k0
% the free-space wave number, and Neff the (real) effective index. The component
% plotted is defined by sComponent, which is a two letter string with the
% following allowed values: 'Ex', 'Ey', 'Ez', 'Hx', 'Hy', 'Hz', and 'Sz'.
% Argument sClass is the class string (see CUT_POS). Set to 'none' if no
% symmetry is required. Arguments x and y are vectors to the x and y
% coordinates where the field is evaluated.
% If Neff is complex, the function must be called with one more argument:
% field = CALC_FIELD(solution, positions, oGd, k0, Neff, sComponent,...
%                 sClass, x, y, 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.
beta = neff*k0;
if ~isreal(beta)
    isLeaky = varargin{1};
else
    isLeaky = 0;
end

if ~strcmp(strClass, 'none')
    xmin = x(1);
    xmax = x(end);
    ymin = y(1);
    ymax = y(end);
    if strClass(1) == 'f'
        r(1) = sqrt(x(1)^2+y(1)^2);
        r(2) = sqrt(x(1)^2+y(end)^2);
        r(3) = sqrt(x(end)^2+y(1)^2);
        r(4) = sqrt(x(end)^2+y(end)^2);
        r = max(r);
    else
        r = max([abs(x), abs(y)]);
    end
    points = round(length(x)*2*r/(x(end)-x(1)));
    x = linspace(-r, r, points);
    y = x;
end

[mX, mY] = meshgrid(x, y);
if sComponent == 'Sz'
    mHx = calc_field(vSol, pos, oGd, k0, neff, 'Hx', strClass, x, y, isLeaky);
    if isnan(mHx)
        mF = NaN;
        return
    end
    mHy = calc_field(vSol, pos, oGd, k0, neff, 'Hy', strClass, x, y, isLeaky);
    if isnan(mHy)
        mF = NaN;
        return
    end
    mEx = calc_field(vSol, pos, oGd, k0, neff, 'Ex', strClass, x, y, isLeaky);
    if isnan(mEx)
        mF = NaN;
        return
    end
    mEy = calc_field(vSol, pos, oGd, k0, neff, 'Ey', strClass, x, y, isLeaky);
    if isnan(mEy)
        mF = NaN;
        return
    end
    mF = 1/2*(mEx.*conj(mHy)-mEy.*conj(mHx));
else
    if ~strcmp(strClass, 'none')
        [mX, mY] = cut_quad(mX, mY);
        if (strClass(1) == 'f') && (lower(sComponent(2)) == 'x' | lower(sComponent(2)) == 'y')
            sCompx = 'Ex';
            sCompy = 'Ey';
            if lower(sComponent(1)) == 'h'
                sCompx = 'Hx';
                sCompy = 'Hy';
            end
            mFx = EvalField(vSol, pos, oGd, k0, beta, sCompx, mX, mY, isLeaky);
            if isnan(mFx)
                mF = NaN;
                return
            end
            mFy = EvalField(vSol, pos, oGd, k0, beta, sCompy, mX, mY, isLeaky);
            if isnan(mFy)
                mF = NaN;
                return
            end
            mF = expand_quad_trans(mFx, mFy, sComponent, strClass);
        else
            mF = EvalField(vSol, pos, oGd, k0, beta, sComponent, mX, mY, isLeaky);
            mF = expand_quad(mF, sComponent, strClass);
        end
    else
        mF = EvalField(vSol, pos, oGd, k0, beta, sComponent, mX, mY, isLeaky);
    end
end
if isnan(mF)
    mF = NaN;
    return
end
if ~strcmp(strClass, 'none')
    iXmin = floor((xmin/r/2+1/2)*points+1);
    iXmax = ceil((xmax/r/2+1/2)*points+1);
    iYmin = floor((ymin/r/2+1/2)*points+1);
    iYmax = ceil((ymax/r/2+1/2)*points+1);
    mF = mF(iYmin:iYmax, iXmin:iXmax);
end