cut_pos.m

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

function posOut = cut_pos(oGd, pos, strClass)
% outPos = CUT_POS(oGd, pos, sClass)
% Cuts all sources and testing points in pos, outside of the minimum
% sector defined by the class in sClass. 
% Argument sClass is a two letter string. The first letter is either
% 'p' - for McIsaac's classes or 'f' for Fini's.
% The second character is a number which corresponds to the numbering
% of mode type in these two schemes. Possible strings are:
% 'p1', 'p2' ... 'p8', and 'f1' ... 'f6'.
vGreen = green_fun_vector(strClass);
ang = 0;
bAbove = 1;
vN = [pos.testing.n];
vNS = [pos.sources.n];
posOut = cut_plane(oGd, pos, ang, bAbove, vGreen);
vNCut = [posOut.testing.n];
vNCutS = [posOut.sources.n];
iCutFirst = find(vNCut < vN & vNCut > 0);
iCutFirstS = find(vNCutS < vNS & vNCutS > 0);
ang = sym_angle(strClass);
bAbove = 0;
vN = [posOut.testing.n];
vNS = [posOut.sources.n];
posOut = cut_plane(oGd, posOut, ang, bAbove, vGreen);
if lower(strClass(1)) ~= 'f'
        return
end
% find pairs of curves that are now one
vNCut = [posOut.testing.n];
vNCutS = [posOut.sources.n];
iCutSecond = find(vNCut < vN & vNCut > 0);
iCutSecondS = find(vNCutS < vNS & vNCutS > 0);
iCutFirst = intersect(iCutFirst, find(vNCut~=0));
iCutFirstS = intersect(iCutFirstS, find(vNCutS~=0));
iCutTwice = intersect(iCutFirst, iCutSecond);
iCutTwiceS = intersect(iCutFirstS, iCutSecondS);
iCutFirst = setdiff(iCutFirst, iCutTwice);
iCutFirstS = setdiff(iCutFirstS, iCutTwiceS);
iCutSecond = setdiff(iCutSecond, iCutTwice);
iCutSecondS = setdiff(iCutSecondS, iCutTwiceS);
iCutAll = [iCutFirst iCutSecond];
iCutAllS = [iCutFirstS iCutSecondS];
if isempty(iCutAll) 
        return
end
% sort sd's and curves in order of increasing distance from center
for iCenter = 1:length(iCutAll)
        curve = oGd.curveArray(iCutAll(iCenter));
        vR(iCenter) = abs(oGd.offsetX - curve.xc + j*(oGd.offsetY - curve.yc));
end
[dummy, ind] = sort(vR);
iCutAll = iCutAll(ind);
clear vR
for iCenter = 1:length(iCutAllS)
        curve = oGd.sdArray(iCutAllS(iCenter)).insideCurve;
        vR(iCenter) = abs(oGd.offsetX - curve.xc + j*(oGd.offsetY - curve.yc));
end
[dummy, ind] = sort(vR);
iCutAllS = iCutAllS(ind);
for ind = 1:2:length(iCutAll)-1
        testing1 = posOut.testing(iCutAll(ind));
        testing2 = posOut.testing(iCutAll(ind+1));
        unified.vX = [testing1.vX testing2.vX];
        unified.vY = [testing1.vY testing2.vY];
        unified.vNx = [testing1.vNx testing2.vNx];
        unified.vNy = [testing1.vNy testing2.vNy];
        unified.n = testing1.n + testing2.n;
        unified.iCurve = iCutAll(ind);
        posOut.testing(iCutAll(ind)) = unified;
        posOut.testing(iCutAll(ind+1)).vX = [];
        posOut.testing(iCutAll(ind+1)).vY = [];
        posOut.testing(iCutAll(ind+1)).vNx = [];
        posOut.testing(iCutAll(ind+1)).vNy = [];
        posOut.testing(iCutAll(ind+1)).n = 0;
        % iCurve is a reference to the testing points
        posOut.testing(iCutAll(ind+1)).iCurve = iCutAll(ind); 
end

for ind = 1:2:length(iCutAllS)-1
        % iSd is a reference to the sources
        vAll = [iCutAllS(ind) iCutAllS(ind+1)]; 
        posOut.sources(iCutAllS(ind)).iSd = vAll;
        posOut.sources(iCutAllS(ind+1)).iSd = vAll;
end