shape.m

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

% equidistant points
function [vSx, vSy, varargout] = shape(x0, y0, r, n, vein)
N = 2500;
cs = vein.cs;
s = linspace(cs.breaks(1), cs.breaks(end), N);
mX = ppval(cs, s);
vSx = mX(1,:) + x0;
vSy = mX(2,:) + y0;
ds = (s(2)-s(1));
dcs = fnder(cs, 1);
mX = ppval(dcs, s);
dxds = mX(1,:);
dyds = mX(2,:);
vLen = cumsum(sqrt(dyds.^2 + dxds.^2))*ds;
L = vLen(end);
dL = L/n;
vLenTestPoints = linspace(0, L-dL, n);
[mLen, mLenT] = ndgrid(vLen, vLenTestPoints);
mD = mLen-mLenT;
[y, ind] = min(abs(mD));

warning('off', 'MATLAB:divideByZero')
m = dyds./dxds;
warning('on', 'MATLAB:divideByZero')
d2cs = fnder(dcs, 1);
mX = ppval(d2cs, s);
dx2ds2 = mX(1,:);
dy2ds2 = mX(2,:);
r_curv  = (dxds.^2+dyds.^2).^(1.5)./abs(dxds.*dy2ds2-dyds.*dx2ds2);

sgn = sign(vSx-x0);
vNx = -m./sqrt(1+m.^2);
vNy = 1./sqrt(1+m.^2);
[dummy, nan_ind] = find(isnan(vNx)+isnan(vNy));
vNx(nan_ind) = 1;
vNy(nan_ind) = 0;

vSx = vSx(ind);
vSy = vSy(ind);
vNx = vNx(ind);
vNy = vNy(ind);
%r_curv = r_curv(ind);
min_r_curv = min(r_curv);

tips_x = vSx+min_r_curv*vNx/3;
tips_y = vSy+min_r_curv*vNy/3;
inv_tips_x = vSx-min_r_curv*vNx/3;
inv_tips_y = vSy-min_r_curv*vNy/3;
tips_r = sqrt((tips_x-x0).^2+(tips_y-y0).^2);
inv_tips_r = sqrt((inv_tips_x-x0).^2+(inv_tips_y-y0).^2);
vSign = (tips_r>inv_tips_r)*2-1;
vNx = vSign.*vNx;
vNy = vSign.*vNy;
if nargout > 2 % calc normals
    varargout(1) = {vNx};
    varargout(2) = {vNy};
end

if r ~= 1
    d = min_r_curv*(1-r);
    if (r < 0) | (r > 2)
        cs = cs_dist(vein, -d);
    else
        xi = vSx-vNx.*d;
        yi = vSy-vNy.*d;
        xy = [[xi xi(1)]; [yi yi(1)]];
        cs = cscvn(xy);
    end
    N = 2500;
    s = linspace(cs.breaks(1), cs.breaks(end), N);
    f = fnval(cs, s);
    xs = f(1,:);
    ys = f(2,:);
    ds = s(2)-s(1);
    vLen = cumsum(sqrt((per_diff(xs)).^2 + (per_diff(ys)).^2));
    L = vLen(end);
    dL = L/n;
    vLenTestPoints = linspace(0, L-dL, n);
    [mLen, mLenT] = ndgrid(vLen, vLenTestPoints);
    mD = mLen-mLenT;
    [y, ind] = min(abs(mD));
    vSx = xs(ind);
    vSy = ys(ind);
end



%uncomment to check
% x = vSx;
% y = vSy;
% r_curv = 1;
% plot(x,y,'.','markersize',1)
% nx = vNx;
% ny = vNy;
% hold on
% plot(x,y,'r.')
% for k = 1:n
% line([x(k),x(k)+nx(k)*r_curv],[y(k),y(k)+ny(k)*r_curv])
% end
% axis equal
% hold off