evalfield.m

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

function mF_all = EvalField(vSol, pos, oGd, k, beta, sComponent, mX, mY, isLeaky)
global EPSILON0 MU0 C DELTA ETA
bCancel = 0;
w = k*C;
Z = j*w*MU0;
vISources = [0 cumsum(2*[pos.sources.n])];
sdArray = oGd.sdArray;
mF_all = zeros(size(mX));
pd = 0; %percent done;
hWaitbar = waitbar(pd,['Plotting ' sComponent ' field...'], 'CreateCancelBtn', @cancel_callback);
for iSd = 1:oGd.nSd
    if pos.sources(iSd).n == 0
        continue
    end
    vSd = pos.sources(iSd).iSd;
    for indVsd = 1:length(vSd)
        % iSd_sources points to the sd from which the source for
        % this sd should be taken
        iSd_sources = vSd(indVsd);
        sd = sdArray(iSd);
        Er = sd.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);
        vX = pos.sources(iSd_sources).vX;
        vY = pos.sources(iSd_sources).vY;
        first_magnetic = vISources(iSd_sources);
        first_electric = first_magnetic+pos.sources(iSd_sources).n;
        switch lower(sComponent)
            case 'ez'
                bessel_order = 0;
                first = first_electric;
                scaling = 1;
                bTrans = 0;
            case 'hz'
                bessel_order = 0;
                first = first_magnetic;
                bTrans = 0;
                scaling = ETA^2/ETA_sd2;
            otherwise
                bessel_order = 1;
                bTrans = 1;
                scaling = j*beta/kc;
        end
        mF = zeros(size(mX));
        mB = in_subdomain(mX, mY, sd);
        mX_sd = mX(mB);
        mY_sd = mY(mB);
        if isempty(mX_sd) | isempty(mY_sd)
            continue
        end
        vGreen = pos.sources(iSd_sources).vGreen;
        ang = 2*pi/length(vGreen);
        bFini = any(imag(vGreen));
        for ii = 1:length(vX)
            if bCancel
                mF_all = NaN;
                return
            end
            pd = pd + 1;
            waitbar(pd/vISources(end)*2, hWaitbar);
            for iReflect = 1:length(vGreen)
                mR = sqrt((mX_sd-vX(ii)).^2 + (mY_sd-vY(ii)).^2);
                mFii = 1/(Y*4*j)*kc^2*besselh(bessel_order, 2, kc*mR)*scaling;
                if bTrans
                    switch lower(sComponent)
                        case 'ex'
                            cos_e = (mX_sd-vX(ii))./mR;
                            cos_m = (mY_sd-vY(ii))./mR*Y/(j*beta)*ETA;
                        case 'ey'
                            cos_e = (mY_sd-vY(ii))./mR;
                            cos_m = -(mX_sd-vX(ii))./mR*Y/(j*beta)*ETA;
                        case 'hx'
                            cos_e = -(mY_sd-vY(ii))./mR*Y*ETA/(j*beta);
                            cos_m = (mX_sd-vX(ii))./mR*ETA^2/ETA_sd2;
                        case 'hy'
                            cos_e = (mX_sd-vX(ii))./mR*Y*ETA/(j*beta);
                            cos_m = (mY_sd-vY(ii))./mR*ETA^2/ETA_sd2;
                    end
                    cos_e = cos_e*vGreen(iReflect);
                    cos_m = cos_m*vGreen(iReflect)*(2*bFini-1)^(iReflect-1);
                    mF(mB) = mF(mB) + mFii.*(cos_e*vSol(first_electric+ii)+cos_m*vSol(first_magnetic+ii));
                else
                    bElectric = first == first_electric;
                    mF(mB) = mF(mB) + mFii*vSol(first+ii)*vGreen(iReflect)*(2*(bElectric | bFini)-1)^(iReflect-1);
                end
                if iReflect ~= length(vGreen)
                    [vX(ii), vY(ii)] = reflect(vX(ii), vY(ii), ang, iReflect, oGd.offsetX, oGd.offsetY);
                end
            end
        end
        mF_all(mB) = mF_all(mB) + mF(mB);
    end
end
close(hWaitbar)
if lower(sComponent(1)) == 'h'
    mF_all = mF_all/ETA;
end

    function cancel_callback(varargin)
        bCancel = 1;
        delete(gcf)
    end
end