function SurCharge(Re1, Re2, Type)
%   SurCharge shows the surface charge distribution on the metal walls
%   of the rectangular waveguide.

%   Xiaowei He & Markus Zahn, 2004
%   All Rights Reserved   
helpinfo=strvcat(['The Demo is not ready'... 
      .'],[' '],...
      [' '],[' @2004 Xiaowei He & Markus Zahn. All Rights Reserved.']);
  
if nargin==0,
   diagHndl=helpdlg(helpinfo,'Information');
   
elseif nargin==1,
   diagHndl=helpdlg(helpinfo,'Information');
   
elseif nargin==2,
   diagHndl=helpdlg(helpinfo,'Information');
   
elseif nargin==3,
   % Define constants
   step=61;
   % Geometry and Mode numbers
   a=1;
   b=1;
   c=1;
   m=Re1;
   n=Re2;
   p=1;
   M = moviein(10);
   figmovie=figure( ...
      'Name','Surface Charge Distribution on Rectangular Metal Waveguide', ...
      'NumberTitle','off', ...
      'Visible','on', ...
      'BackingStore','off');
   set(figmovie,'Position',[40,93,955,625]);
   axis([0 a 0 b 0 c]);
   hold on;
   box on;
   % Set the view 
   campos([2 -1 -1]);
   camtarget([1 0 0]);
   camup([0 -1 0]);
   colorbar;
   % Set the illustration position
   xslice = [0,1]; 
   yslice = [0,1]; 
   zslice = 1;
   xinter=a/(2*n);
   yinter=b/(2*m);
   % Wave numbers
   kx=m*pi/a;
   ky=n*pi/b;
   kz=p*pi/c;
   x1=linspace(0,a,step);
   y1=linspace(0,b,step);
   z=linspace(0,c,step);
   [X1,Y1,Z]=meshgrid(x1,y1,z);
   kk=0;
      
   for tt=0:30:360;
       t = 2*pi*tt/360;
       kk=kk+1;
       % TM Mode
       if Type == 2
           U1=-kx*(cos(kx*X1)).*(sin(ky*Y1)).*(cos(kz*Z-t)); % Ex
           V1=-ky*(sin(kx*X1)).*(cos(ky*Y1)).*(cos(kz*Z-t)); % Ey
           W1=(sin(kx*X1)).*(sin(ky*Y1)).*(cos(kz*Z-t)); % Ez
           U2=ky*(sin(kx*X1)).*(cos(ky*Y1)).*(cos(kz*Z-t)); % Hx
           V2=-kx*(cos(kx*X1)).*(sin(ky*Y1)).*(cos(kz*Z-t)); % Hy
           W2=0*V2; % Hz
       % TE Mode
       else
           U1=ky*(cos(kx*X1)).*(sin(ky*Y1)).*(cos(kz*Z-t));
           V1=-kx*(sin(kx*X1)).*(cos(ky*Y1)).*(cos(kz*Z-t));
           W1=0*V1;
           U2=kx*(sin(kx*X1)).*(cos(ky*Y1)).*(cos(kz*Z-t));
           V2=ky*(cos(kx*X1)).*(sin(ky*Y1)).*(cos(kz*Z-t));
           W2=(cos(kx*X1)).*(cos(ky*Y1)).*(cos(kz*Z-t));;
       end
       Etotal = U1 + V1 + W1; 
       hsurfaces = slice(X1,Y1,Z,Etotal,xslice,yslice,zslice);
	   set(hsurfaces,'FaceColor','interp','EdgeColor','none');
%        hcurrentsx0 = streamslice(X1,Y1,Z,U2*0,-W2,V2,[0],[],[],0.5); % At x=0 plane normal is -x direction, Ix=0, Iy=-Hz, Iz=Hy;
%        hcurrentsx1 = streamslice(X1,Y1,Z,U2*0,W2,-V2,[1],[],[],0.5); % At x=1 plane normal is x direction, Ix=0, Iy=Hz, Iz=-Hy;
%        hcurrentsy0 = streamslice(X1,Y1,Z,W2,V2*0,-U2,[],[0],[],0.5); % At y=0 plane normal is -y direction, Ix=Hz, Iy=0, Iz=-Hx;
%        hcurrentsy1 = streamslice(X1,Y1,Z,-W2,V2*0,U2,[],[1],[],0.5); % At y=1 plane normal is y direction, Ix=-Hz, Iy=0 Iz=Hx;
% %        set(hcurrents,'FaceColor','interp','EdgeColor','none');
       
          
       % Define and make nice notation
       if kk==1
           set(gca,'XTick', [0 1], 'YTick', [0 1], 'ZTick',[0 1])
           daspect([2,2,1]);  
           axis tight
           text(0.19, 0.00, 'x/a', 'Units', 'normalized');
           text(-0.09, 0.34, 'y/b', 'Units', 'normalized');
           text(0.18, 0.86, 'z/c', 'Units', 'normalized');
           Tmode = strvcat(['TE'], ['TM']);
           text(1.16, 1.04,'Relative Charge Density','Units','normalized');
           colorbar;
       end  
       if Type == 1
           title(['NonDimensional Surface Charge Distribution \sigma_fc/H_0 of ', Tmode(Type,:), num2str(Re1), num2str(Re2), ' Mode, \omegat = ', num2str(tt),' degree.']);
       else
           title(['NonDimensional Surface Charge Distribution \sigma_f/\epsilonE_0 of ', Tmode(Type,:), num2str(Re1), num2str(Re2), ' Mode, \omegat = ', num2str(tt),' degree.']);
       end
       
       M(kk) = getframe(gcf);
       
	end
    
end

