dispersion.m

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

clear 
clf
smt_init
load -mat step_index.gd
positions = distribute(oGd);

vLambda = 1.5:0.05:2.5;
steps = length(vLambda);

mNeff = []; % matrix for effective indices
mErr = []; % matrix for continuity condition error
tol = 1e-6;
sm = 2.5; % smoothness criterion
nStart = 7;
bPlot = 1;

min_beta = 1.5;
max_beta = 1.75;
tic
for ind = 1:steps 
    lambda = vLambda(ind);
    k0 = 2*pi/lambda;
    figure(1)
    [vNeff, vErr] = find_modes(oGd, positions, k0, min_beta, ...
        max_beta, tol, sm, bPlot, nStart);
    if vNeff == -1
        break;
    end
    figure(2)
    vNeff = vNeff(:);
    vErr = vErr(:);
    if isempty(mNeff)
        mNeff = [mNeff vNeff];
        mErr = [mErr vErr];
    else
        PrevBeta = mNeff(:,end);
        % deal with the possibility of modes going in or out of the 
        % range of Neff and lambda, by filling their places in the 
        % mNeff with NaN
        PrevNaN = ones(length(vNeff)-length(PrevBeta), size(mNeff,2))*NaN;
        CurrNaN = ones(length(PrevBeta)-length(vNeff), 1)*NaN;
        if ~isempty(PrevNaN)
            mNeff = [mNeff; PrevNaN];
        end
        if ~isempty(CurrNaN)
            vNeff = [vNeff; CurrNaN];
        end
        mNeff = [mNeff vNeff];
    end        
    fraction_done = ind/steps;
    fprintf(1, 'percent done: %d\n', round(fraction_done*100))
    t = toc;
    fprintf(1,'minutes to go: %6.2f\n', (1-fraction_done)/fraction_done*t/60)
    plot(2*pi./vLambda(1:ind), mNeff', 'b.')
    xlabel('k_0')
    ylabel('n_{eff}')
end