fullvector_all_fields.m

强大的计算电磁场本征函数与本征模的程序

% This example demonstrates how to use the auxilliary routine
% (postprocess.m) to calculate the four remaining field
% components (Hz, Ex, Ey, Ez) once the transverse magnetic
% field components (Hx, Hy) are known.  The waveguide
% considered here is a uniaxial channel waveguide with the
% c-axis oriented at pi/4 relative to the x-y axes.

n1 = 1.55;
n2x = 2.156;
n2y = 2.232;
n2z = 2.232;
theta = pi/4;

e2xx = n2x^2*cos(theta)^2 + n2y^2*sin(theta)^2;
e2yy = n2y^2*cos(theta)^2 + n2x^2*sin(theta)^2;
e2xy = cos(theta)*sin(theta)*(n2x^2-n2y^2);
e2yx = e2xy;

Rx = 0.30;
Ry = 0.20;
side = 0.20;

dx = 0.0025;        % grid size (x)
dy = dx;            % grid size (y)

lambda = 1.00;      % wavelength
nmodes = 1;         % number of modes to compute

fprintf (1,'generating index mesh...\n');

[x,y,xc,yc,nx,ny,epsxx] = ...
    waveguidemeshfull([n1,sqrt(e2xx),n1],[side,2*Ry,side],2*Ry,Rx, ...
                  side,dx,dy); 
[x,y,xc,yc,nx,ny,epsxy] = ...
    waveguidemeshfull([0,sqrt(e2xy),0],[side,2*Ry,side],2*Ry,Rx, ...
                  side,dx,dy); 
[x,y,xc,yc,nx,ny,epsyx] = ...
    waveguidemeshfull([0,sqrt(e2yx),0],[side,2*Ry,side],2*Ry,Rx, ...
                  side,dx,dy); 
[x,y,xc,yc,nx,ny,epsyy] = ...
    waveguidemeshfull([n1,sqrt(e2yy),n1],[side,2*Ry,side],2*Ry,Rx, ...
                  side,dx,dy); 
[x,y,xc,yc,nx,ny,epszz] = ...
    waveguidemeshfull([n1,n2z,n1],[side,2*Ry,side],2*Ry,Rx, ...
                  side,dx,dy); 

% Now we stretch out the mesh at the boundaries:
[x,y,xc,yc,dx,dy] = stretchmesh(x,y,[80,80,80,80],[4,4,4,4]);

fprintf (1,'solving for eigenmodes...'); t = cputime;
[Hx,Hy,neff] = wgmodes (lambda, n2y, nmodes, dx, dy, epsxx, epsxy, epsyx, epsyy, epszz, '0000');
fprintf (1,'done (cputime = %7.3f)\n', cputime-t);

fprintf(1,'neff = %7.5f\n',neff);

fprintf (1,'post-processing...'); t = cputime;
[Hz,Ex,Ey,Ez] = postprocess (lambda, neff, Hx, Hy, dx, dy, epsxx, epsxy, epsyx, epsyy, epszz, '0000');
fprintf (1,'done (cputime = %7.3f)\n', cputime-t);

figure(1);

ii = 1;

subplot(231);
contourmode(x,y,Hx,1,3,60,sprintf('Hx (mode %d)',ii));
line([-Rx,-Rx,+Rx,+Rx,-Rx], ...
     [side,side+2*Ry,side+2*Ry,side,side]);

subplot(232);
contourmode(x,y,Hy,1,3,60,sprintf('Hy (mode %d)',ii));
line([-Rx,-Rx,+Rx,+Rx,-Rx], ...
     [side,side+2*Ry,side+2*Ry,side,side]);

subplot(233);
contourmode(x,y,Hz,1,3,60,sprintf('Hz (mode %d)',ii));
line([-Rx,-Rx,+Rx,+Rx,-Rx], ...
     [side,side+2*Ry,side+2*Ry,side,side]);

subplot(234);
contourmode(xc,yc,Ex,1,3,60,sprintf('Ex (mode %d)',ii));
line([-Rx,-Rx,+Rx,+Rx,-Rx], ...
     [side,side+2*Ry,side+2*Ry,side,side]);

subplot(235);
contourmode(xc,yc,Ey,1,3,60,sprintf('Ey (mode %d)',ii));
line([-Rx,-Rx,+Rx,+Rx,-Rx], ...
     [side,side+2*Ry,side+2*Ry,side,side]);

subplot(236);
contourmode(xc,yc,Ez,1,3,60,sprintf('Ez (mode %d)',ii));
line([-Rx,-Rx,+Rx,+Rx,-Rx], ...
     [side,side+2*Ry,side+2*Ry,side,side]);

figure(2);

ii = 1;
colormap(jet(256));
hn = abs(interp2(y,x,Hy,side+Ry,0));
en = abs(interp2(yc,xc,Ex,side+Ry,0));

subplot(231);
imagemode(x,y,Hx/hn,sprintf('Hx (mode %d)',ii));
line([-Rx,-Rx,+Rx,+Rx,-Rx], ...
     [side,side+2*Ry,side+2*Ry,side,side]);

subplot(232);
imagemode(x,y,Hy/hn,sprintf('Hy (mode %d)',ii));
line([-Rx,-Rx,+Rx,+Rx,-Rx], ...
     [side,side+2*Ry,side+2*Ry,side,side]);

subplot(233);
imagemode(x,y,Hz/hn,sprintf('Hz (mode %d)',ii));
line([-Rx,-Rx,+Rx,+Rx,-Rx], ...
     [side,side+2*Ry,side+2*Ry,side,side]);

subplot(234);
imagemode(xc,yc,Ex/en,sprintf('Ex (mode %d)',ii));
line([-Rx,-Rx,+Rx,+Rx,-Rx], ...
     [side,side+2*Ry,side+2*Ry,side,side]);

subplot(235);
imagemode(xc,yc,Ey/en,sprintf('Ey (mode %d)',ii));
line([-Rx,-Rx,+Rx,+Rx,-Rx], ...
     [side,side+2*Ry,side+2*Ry,side,side]);

subplot(236);
imagemode(xc,yc,Ez/en,sprintf('Ez (mode %d)',ii));
line([-Rx,-Rx,+Rx,+Rx,-Rx], ...
     [side,side+2*Ry,side+2*Ry,side,side]);