e_sources.m

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

% calculate TM fields. 
% each line in the matrix represents a matching-point (MP).
function [mHet, mEet, mEez] = e_sources(testing, sources, Y, kc, beta)
global EPSILON0 MU0 C

vSx = testing.vX;
vSy = testing.vY;
vNx = testing.vNx;
vNy = testing.vNy;
vEx = sources.vX;
vEy = sources.vY;

[mSx, mEx] = ndgrid(vSx, vEx);
[mSy, mEy] = ndgrid(vSy, vEy);
[mNx, mEx] = ndgrid(vNx, vEx);
[mNy, mEy] = ndgrid(vNy, vEy);

mR = sqrt((mSx-mEx).^2 + (mSy-mEy).^2);

%if abs(kc)>TOL
    % electric field radial component
%   mEr = 1/(4*Y)*beta*kc*besselh(1,2,kc*mR);
    mEr = 1/(4*Y)*beta*kc*hankel2(1,kc*mR);
    % electric field z component
%    mEez = 1/(Y*4*j)*kc^2*besselh(0,2,kc*mR);
    mEez = 1/(Y*4*j)*kc^2*hankel2(0,kc*mR);
% else
%      % electric field radial component
%      mEr = j/(4*Y)*beta/pi*2./mR;
%      % electric field z component
%      mEez = 1/(Y*4*j)*j*2/pi*log(1.781/2*mR)*kc^2;
% end
% magnetic field azimuthal component
mHphi = Y/(j*beta)*mEr;
% calculate tangential components
mDx = mSx-mEx;
mDy = mSy-mEy;
mAlpha = 1./mR.*(mDx.*mNx + mDy.*mNy);
mBeta = 1./mR.*(-mDx.*mNy + mDy.*mNx);
mEet = mEr.*mBeta;
mHet = mHphi.*mAlpha;