calc_att.m

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

function [Neff_i, mSz] = calc_att(vSol, pos, oGd, k0, neff, x, y)
% [neff_i, Sz] = CALC_ATT(solution, positions, oGd, k0, Neff, x, y)
% Calculates the attenuation of a mode due to radiation by a perturbational
% method.
% 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 effective index. The points on
% which Poynting's vector is evaluated is determined by x and y.
% Their range should include all the boundaries, but should not be too 
% far away from the outermost boundary. 
% The output is the imaginary part of the effective index, and the longitudinal
% component of Poynting's vector. 
% The function may be used iteratively, first with a real effective index,
% and on the next iteration, with the imaginary part found previously.
% Note: at the moment, symmetry is not supported. The effective index may
% be found with the aid of symmetry, but the solution must then be
% recomputed, at the known effecitve index, without symmetry.
[mX, mY] = meshgrid(x, y);
bottom = y(1);
left = x(1);
right = x(end);
top = y(end);

[rows, cols] = size(mX);
beta = neff*k0;
% calculate Sz:

mHx = EvalField(vSol, pos, oGd, k0, beta, 'Hx', mX, mY, 1);
mHy = EvalField(vSol, pos, oGd, k0, beta, 'Hy', mX, mY, 1);
mEx = EvalField(vSol, pos, oGd, k0, beta, 'Ex', mX, mY, 1);
mEy = EvalField(vSol, pos, oGd, k0, beta, 'Ey', mX, mY, 1);
mSz = mEx.*conj(mHy)-mEy.*conj(mHx);

% calculate Sy1:
% (power into the bottom edge)

mX = x;
mY = repmat(bottom, 1, cols);
mHz = EvalField(vSol, pos, oGd, k0, beta, 'Hz', mX, mY, 1);
mHx = EvalField(vSol, pos, oGd, k0, beta, 'Hx', mX, mY, 1);
mEz = EvalField(vSol, pos, oGd, k0, beta, 'Ez', mX, mY, 1);
mEx = EvalField(vSol, pos, oGd, k0, beta, 'Ex', mX, mY, 1);
mSy1 = mEz.*conj(mHx)-mEx.*conj(mHz);

% calculate Sy2:
% (power out of the upper edge)
mX = x;
mY = repmat(top, 1, cols);
mHz = EvalField(vSol, pos, oGd, k0, beta, 'Hz', mX, mY, 1);
mHx = EvalField(vSol, pos, oGd, k0, beta, 'Hx', mX, mY, 1);
mEz = EvalField(vSol, pos, oGd, k0, beta, 'Ez', mX, mY, 1);
mEx = EvalField(vSol, pos, oGd, k0, beta, 'Ex', mX, mY, 1);
mSy2 = mEz.*conj(mHx)-mEx.*conj(mHz);

% calculate Sx1:
% (power into left edge)
mY = y(:);
mX = repmat(left, rows, 1);
mHz = EvalField(vSol, pos, oGd, k0, beta, 'Hz', mX, mY, 1);
mHy = EvalField(vSol, pos, oGd, k0, beta, 'Hy', mX, mY, 1);
mEz = EvalField(vSol, pos, oGd, k0, beta, 'Ez', mX, mY, 1);
mEy = EvalField(vSol, pos, oGd, k0, beta, 'Ey', mX, mY, 1);
mSx1 = mEy.*conj(mHz)-mEz.*conj(mHy);

% calculate Sx2:
% (power out of right edge)
mY = y(:);
mX = repmat(right, rows, 1);
mHz = EvalField(vSol, pos, oGd, k0, beta, 'Hz', mX, mY, 1);
mHy = EvalField(vSol, pos, oGd, k0, beta, 'Hy', mX, mY, 1);
mEz = EvalField(vSol, pos, oGd, k0, beta, 'Ez', mX, mY, 1);
mEy = EvalField(vSol, pos, oGd, k0, beta, 'Ey', mX, mY, 1);
mSx2 = mEy.*conj(mHz)-mEz.*conj(mHy);

dx = abs(x(2)-x(1)); 
dy = abs(y(2)-y(1)); 
Ptrans = -sum(real(mSx1))*dx+sum(real(mSx2))*dx-sum(real(mSy1))*dy+sum(real(mSy2))*dy;
Pz = sum(sum(real(mSz)))*dx*dy;
Neff_i = -Ptrans/(2*Pz)/k0; %should be negative