find_modes.m

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

function [vBeta, vErr_at_peaks] = find_modes(oGd, pos, k0, min_beta, max_beta, tol, sm, bPlot, nStart)
% [Neff, error] = FIND_MODES(oGd, positions, k0, minNeff, maxNeff, tol, smooth, bShow, nStart)
% Finds modes for the geometry in oGd.
% Sources and testing points are distributed according to positions.
% Free-space wave number: k0.
% Range of effective index: (minNeff, maxNeff)
% Tolerance for effective index: tol
% Show search progress: bShow.
% Number of evenly-spaced samples: nStart.
global FontSize
%#function smt_solve
bSkip = 0;
bCancel = 0;
if bPlot
    if bPlot == 2
        hAx = subplot(20,2,3:2:30);
    else
        hAx = subplot(20,2,3:30);
    end
    cla(hAx);
    xlabel('n_{eff}');
    ylabel('\Delta E^{-1}');
    set(gcf, 'DoubleBuffer','on')
end
hWaitbar = waitbar(0,'searching for modes...','CreateCancelBtn', @cancel_callback);
hButton = uicontrol(hWaitbar, 'style','pushbutton');
set(hButton,'string', 'Jump to Fine Search...')
set(hButton,'position', [130 10 140 23])
set(hButton, 'callback', @skip_callback)
max_iter  = 1000;
max_diff_tol = 10^(-sm);
%vBeta = min_beta:peak_width:max_beta;
vBeta = linspace(min_beta, max_beta, nStart);
vBeta_found = [];
vBeta_full = [];
vErr_full = [];
counter = 1;
num_of_min = 0;
while counter < max_iter
    vErr = [];
    for ind = 1:length(vBeta)
        err = smt_solve(vBeta(ind), oGd, pos, k0);
        [vBeta_full, s_ind] = sort([vBeta_full vBeta(ind)]);
        vErr_full = [vErr_full err];
        vErr_full = vErr_full(s_ind);
        if bPlot
            hold(hAx, 'off')
            plot(hAx, vBeta_full, -2*db(vErr_full),'b', vBeta_full, -2*db(vErr_full),'r.')
            axis(hAx, 'tight')
            fixfonts(hAx, FontSize)
            xlabel(hAx, 'n_{eff}')
            ylabel(hAx, '1/\DeltaE [dB]')
            drawnow
            pause(0.1)
        end
        if bSkip
            break
        end
        if bCancel
            delete(hWaitbar);
            vBeta = -1;
            vErr_at_peaks = [];
            return
        end
    end
    if bSkip
        break
    end
    [dummy, mj] = local_maxima(1./vErr_full);
    prev_num_of_min = num_of_min;
    num_of_min = length(mj);
    vBeta = [];
    counter = counter+1;
    xx = vBeta_full(1:end-1) + diff(vBeta_full)/2;
    yy_lin = interp1(vBeta_full, -db(vErr_full), xx, 'linear');
    yy_sp = interp1(vBeta_full, -db(vErr_full), xx, 'spline');
    [max_diff, new_beta_ind] = max(abs(yy_lin-yy_sp));
    max_diff = max_diff/rms(yy_sp);
    if num_of_min == 1
        sTitle = [num2str(num_of_min) ' mode found...'];
    else
        sTitle = [num2str(num_of_min) ' modes found...'];
    end
    waitbar(max_diff_tol/max_diff, hWaitbar, sTitle)
    if max_diff < max_diff_tol
        break
    end
    vBeta = xx(new_beta_ind);
    counter = counter + 1;
end
pause(0.1)
if counter == max_iter
    disp('maximum number of iterations reached')
end
[dummy, mj] = local_maxima(1./vErr_full);
if isempty(mj)
    vBeta = [];
    vErr_at_peaks = [];
    delete(hWaitbar)
    return
end

delete(hButton)
for ind = 1:length(mj)
    if bCancel
        delete(hWaitbar);
        vBeta = -1;
        vErr_at_peaks = [];
        return
    end
    prev_neff = 0;
    sTitle = ['Searching in a small neighborhood around mode No. ' num2str(ind) '...'];
    waitbar((ind-1)/length(mj), hWaitbar, sTitle)
    min_local = vBeta_full(mj(ind)-1);
    max_local = vBeta_full(mj(ind)+1);
    [vBeta_final(ind), vErr_final(ind)] = fminbnd('smt_solve', min_local, max_local,...
        optimset('TolX', tol/2, 'OutputFcn', @fine_search_cancel), oGd, pos, k0);
    %[vBeta_final(ind), vErr_final(ind)] = fminbnd('smt_solve', min_local, max_local,...
    %    optimset('TolX', tol/2), oGd, pos, k0);
    waitbar((ind-1)/length(mj), hWaitbar, sTitle)
end
delete(hWaitbar);
vBeta = vBeta_final;
vErr_at_peaks = vErr_final;

    function skip_callback(varargin)
        bSkip = 1;
    end

    function cancel_callback(varargin)
        bCancel = 1;
    end

    function stop = fine_search_cancel(beta, val, state, varargin)
        if ~isempty(beta)
            curr_neff = beta/k0;
            if abs(curr_neff - prev_neff) < tol
                this_mode_prog = 1;
            else
                this_mode_prog = tol/abs(curr_neff-prev_neff);
            end
            waitbar((ind-1)/length(mj)+min(1, this_mode_prog)/length(mj), hWaitbar, sTitle)
            prev_neff = curr_neff;
        end
        drawnow
        stop = bCancel;
    end
end