% 6.014 Parallel plate waveguide simulations
% J. Pacheco (4/6/01)

% Send all bug reports to pacheco@mit.edu

% clear workpace
more off; clear all; close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
minFigWidth        = 620; % pixels. For full display we need a big
minFigHeight       = 520; % figure window.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% different movie types available.
TIMEVAR=1;   %time variation
MODEVAR=2;   %integer mode variation
ANGVAR=3;    %bouncing angle variations
LAMVAR=4;    %lambda variation

% Physical constants
speedlight= 2.9979247917e8;
mu = 4*pi*1e-7; 
epsilon = 1/mu/speedlight/speedlight;
eta = sqrt(mu/epsilon);

disp(' ');
disp('1. Store frames as a matlab movie');
disp('2. Print frames as a series of jpegs');
disp('3. Don''t store frames at all');
storetype=input('Enter storage type: ');

if (storetype == 2)
  basename=input('Enter base name of jpeg images: ','s');
end

% set the movie type
disp('OPTIONS for MOVIE');
disp('1. Vary time');
disp('2. Vary mode (always satisfies guidance condition)');
disp('3. Vary bouncing waves angle (does not always satisfy G.C.)');
disp('4. Vary lambda');
movietype = input('Enter movie type that you wish to view: ');

disp('');
d = input('Enter the distance between the plates (m): ');
len = input('Enter length of plates to show (along z) (m): ');

if (movietype ~= LAMVAR) 
  lam = input('Enter the wavelength (m): ');
  lammin = lam;
else
  lammin = input('Enter smallest wavelength (m): ');
  lammax = input('Enter largest wavelength (m): ');
  nn = input('Enter number of lambda samples (m): ');
  nn = max(nn,2);
  dlam = (lammax-lammin)/(nn-1);
  lam = lammin:dlam:lammax;
  varparmlist = lam;
end

% ensure that only positive nonzero lambda are used
lam = lam(find(lam>0));

kk = 2*pi./lam;  % wave number
omega = kk.*speedlight;  % angular frequency

tem = input('TE (enter 1) or TM (enter 2): ');

if (tem == 1) 
  ctitle= 'V/m';
  ctype = 'TE';
  ftype = 'Electric';
else
  ctitle= 'A/m';
  ctype = 'TM';
  ftype = 'Magnetic';
end

if (movietype ~= ANGVAR)
  if (movietype ~= MODEVAR)
    modenum = input('Enter the mode number (m): ');
  else
    modelow = input('Enter the lowest mode to show: ');
    modehigh = input ('Enter the highest mode to show: ');
    modenum = modelow:modehigh;
    varparmlist = modenum;
  end
else
  modenum = 0.0;
end

if (movietype ~= TIMEVAR)
  t = 0;
else
  % calculate fields for one cycle
  nn = input('Enter number of time samples to plot: ');
  nn = max(nn,2);
  dt = 2*pi./omega/(nn-1);
  times = 0:dt:((nn-1)*dt);
  varparmlist = times;
end

if (movietype == ANGVAR) 
  % define all angles to check if BC are satisfied.
  % next two lines ensure exact modes are included.
  nn = input('Enter number of angles to check: ');
  nn = max(2,nn);
  hm = floor(2*d/lam)+1;
  tmodes = asin([0:hm]*pi/d/kk);
  tmodes = sort([0:(pi/2/(nn-1)):(pi/2) tmodes]);

  % remove any duplicates
  dtmodes = diff(tmodes);
  varparmlist = tmodes([find(abs(diff(tmodes))>1e-10) length(tmodes)]);
end

disp('  ');
disp('1. No scaling');
disp('2. Scale (x,z) axes by d');
if (movietype ~= LAMVAR)
  disp('3. Scale (x,z) axes by lambda');
end
scaleopt = input('Enter your scaling option: ');

if (scaleopt==2)
  scale = d;
  scalet = '/d';
elseif ((scaleopt==3)*(movietype ~= LAMVAR))
  scale = lam;
  scalet = '/\lambda';
else
  scale = 1.0;
  scalet = '';
end

disp(' ');
disp(['1. ' ftype ' Field distribution along x at z=0']);
disp(['2. ' ftype ' Field distribution along z']);
disp('3. Poynting power density (z comp)');
disp('4. K-surface plot');
plot4type=input('Enter 4th plot option: ');

% determine # of grid points to plot based on lambda
dx = d/ceil(20*d/lammin);
dz = len/ceil(20*len/lammin);
x = 0:dx:d;
z = 0:dz:len;
xs = x; zs = z;
Nx = length(x);
Nz = length(z);
[x,z] = meshgrid(x,z);
phi = 0:pi/20:(pi/2);

% find maximum value for axis when doing k-surface plot 
kkm = max(2*pi/lammin,max(modenum)*pi/d);

% reset the movie counter
count = 0;

% Make sure the figure window is big enough.  If not, grow
% it towards the bottom left.  
pos = get(gcf,'Position');
if pos(3) < minFigWidth
  pos(1) = pos(1) - (minFigWidth-pos(3)-1);
  if pos(1) < 1
    pos(1) = 1;
  end;
  pos(3) = minFigWidth;
end;
if pos(4) < minFigHeight
  pos(2) = pos(2) - (minFigHeight-pos(4)-1);
  if pos(2) < 1
    pos(2) = 1;
  end;
  pos(4) = minFigHeight;
end;
set(gcf,'Position',pos);
	
% Define the area to be recorded for movie
rect = get(gcf,'Position');
rect(1:2) = [0 0];

% use guidance condition/dispersion relation to find:
kx = modenum*pi/d;
kz = sqrt(kk.*kk - kx.*kx);
kzr = real(kz); kzi = imag(kz);

for varparm = varparmlist,
  if (movietype==TIMEVAR)
    t=varparm;
  elseif (movietype==MODEVAR)
    modenum = varparm;
    kx = modenum*pi/d;
    kz = sqrt(kk*kk - kx*kx);
    kzr = real(kz); kzi = imag(kz);
  elseif (movietype==ANGVAR)
    kx = real(kk*sin(varparm));
    kz = sqrt(kk*kk - kx*kx);
    kzr = real(kz); kzi = imag(kz);
    modenum = kx*d/pi;
  elseif (movietype==LAMVAR)
    lam = varparm;
    kk = 2*pi./lam;  % wave number
    kx = modenum*pi/d;
    kz = sqrt(kk*kk - kx*kx);
    kzr = real(kz); kzi = imag(kz);
    omega = kk.*speedlight;  % angular frequency
    clear x z xs zs;
    % determine # of grid points to plot based on lambda
    dx = d/ceil(20*d/lam);
    dz = len/ceil(20*len/lam);
    x = 0:dx:d;
    z = 0:dz:len;
    xs = x; zs = z;
    Nx = length(x);
    Nz = length(z);
    [x,z] = meshgrid(x,z);
  else
    error('Unknown movietype!!!!');
  end
    
  % choose different offset phases (cos OR sin) so that plot of 
  % field distribution along x at z=0 is maximum
  if (tem==2)
    % TM case (F=Hy, G=Ex)
    F1 = exp(-kzi*z).*cos(kx*x + kzr*z - omega*t);
    F2 = exp(-kzi*z).*cos(-kx*x + kzr*z - omega*t);
    if (plot4type==3)
      Gtx = (kzr/omega/epsilon)*cos(kx*x).*cos(kzr*z - omega*t) ... 
	  +(kzi/omega/epsilon)*cos(kx*x).*sin(omega*t);
    end
  else
    % TE Case (F=Ey, Gt=Hx_tot)
    F1 = exp(-kzi*z).*sin(kx*x + kzr*z - omega*t);
    F2 = -exp(-kzi*z).*sin(-kx*x + kzr*z - omega*t);
    if (plot4type==3)
      Gtx = -(kzr/omega/mu).*sin(kx*x).*cos(kzr*z-omega*t) ...
	  -(kzi/omega/mu).*sin(kx*x).*exp(-kzi*z).*sin(omega*t);
      Gtx = -Gtx; % multiply by (-1) b/c E x H = -z direction
    end
  end    
   
  Ft = F1 + F2;

  start = [d/2 len/2]/scale;
  stop1 = start + [min(kx/kk,1) kzr/kk]/scale;
  stop2 = start + [-min(kx/kk,1) kzr/kk]/scale;
  stop3 = start + [0 1.0]/scale;
  
  subplot(221)
  pcolor(xs/scale,zs/scale,F1); shading interp; colorbar;
  caxis([-2 2]); colorbar; 
  set(get(colorbar,'Title'),'String',ctitle) 
  axis([-0.1*d/scale 1.1*d/scale 0 len/scale])
  axis equal
  axis([-0.1*d/scale 1.1*d/scale 0 len/scale])
  xlabel(['x' scalet]); ylabel(['z' scalet]);
  arrow(start,stop1)
  line([0 0],[0 len/scale],'linewidth',3,'color','k');
  line([d d]/scale,[0 len/scale],'linewidth',3,'color','k');
  title([ftype ' Field 1']);

  subplot(222)
  pcolor(xs/scale,zs/scale,F2); shading interp; colorbar;
  caxis([-2 2]); colorbar; 
  set(get(colorbar,'Title'),'String',ctitle) 
  axis([-0.1*d/scale 1.1*d/scale 0 len/scale])
  axis equal
  axis([-0.1*d/scale 1.1*d/scale 0 len/scale])
  xlabel(['x' scalet]); ylabel(['z' scalet]);
  arrow(start,stop2)
  line([0 0],[0 len/scale],'linewidth',3,'color','k');
  line([d d]/scale,[0 len/scale],'linewidth',3,'color','k');
  title([ftype ' Field 2'])
  
  subplot(223)
  pcolor(xs/scale,zs/scale,Ft); shading interp; 
  caxis([-2 2]); colorbar; 
  set(get(colorbar,'Title'),'String',ctitle) 
  axis([-0.1*d/scale 1.1*d/scale 0 len/scale])
  axis equal
  axis([-0.1*d/scale 1.1*d/scale 0 len/scale])
  xlabel(['x' scalet]); ylabel(['z' scalet]);
  arrow(start,stop3);
  line([0 0],[0 len/scale],'linewidth',3,'color','k');
  line([d d]/scale,[0 len/scale],'linewidth',3,'color','k');
  title(['Total ' ftype ' Field']);
  
  subplot(224)
  if (plot4type==2)
    PP = round(Nx/2);
    plot(zs.'/scale,Ft(:,PP));
    axis([0 len/scale -2.1 2.1]);
    grid on;
    xlabel(['z' scalet]); ylabel(ctitle);
    title([ftype ' Field distribution along z']);
  elseif (plot4type==1)
    plot(xs.'/scale,Ft(1,:));
    axis([0 d/scale -2.1 2.1])
    xlabel(['x' scalet]); ylabel(ctitle);
    title([ftype ' Field distribution along x at z=0']);
  elseif (plot4type==3)
    if (tem==1)
      PP = round(Nx/2/modenum);
    else
      PP = 2;
    end
    if (tem == 1)
      if (movietype == LAMVAR)
	plot(xs./scale,2*eta*(kzr/omega/mu)*(sin(kx*x(1,:))).^2)
      else
	plot(xs./scale,eta*Ft(1,:).*Gtx(1,:));
      end
      ylabel('\eta Watts / m^2');
    else
      if (movietype == LAMVAR)
	plot(xs./scale,2*(1/eta)*(kzr/omega/epsilon)*(cos(kx*x(1,:))).^2)
      else
	plot(xs./scale,(1/eta)*Ft(1,:).*Gtx(1,:));
      end
      ylabel('Watts / \eta m^2');
    end
    xlabel(['x' scalet]); 
    
    axis([0 d/scale -2 2])
    if (movietype ~= LAMVAR) 
      title('Poynting power density, (z comp) ');
    else
      title('Time averaged S_z');
    end
  elseif (plot4type==4)
    plot(kk*cos(phi),kk*sin(phi));
    axis([0 kkm+kkm/10 0 kkm+kkm/10]);axis equal;
    axis([0 kkm+kkm/10 0 kkm+kkm/10]);grid on;
    %axis([0 kk+kk/10 0 kk+kk/10]);axis equal;
    %axis([0 kk+kk/10 0 kk+kk/10]);grid on;
    arrow([0 0],[kx kzr]);
    hold on;plot([kx kx],[kzr 0],'k--');hold off
    title('K-surface');xlabel('k_x');ylabel('k_z');
  end 

  grid on
  if (movietype==TIMEVAR)
    suptitle([ctype ' ' sprintf('%1.0f',modenum) ' Mode  \omega t =' ...
	  sprintf('%5.2f',omega*t/pi) ' \pi']);
  elseif (movietype==LAMVAR)
    suptitle([ctype ' ' sprintf('%1.0f',modenum) ' Mode  \lambda =' ...
	  sprintf('%5.2f',lam/d) ' d']) 
  else 
    theta = atan2(kx,kzr);
    suptitle([ctype ' Mode, m=' sprintf('%6.3f',kx*d/pi) ... 
	  ', Angle = ' sprintf('%4.1f',theta*180/pi) ' degs'])
  end

  % record movie frame
  if (storetype==1) 
    count = count + 1;
    Mv(count) = getframe(gcf,rect);
    %if (count == 1)
    %print('-djpeg75',strcat(basename,'.jpg'));
    %end
    
  elseif (storetype==2)
    if (count > 9)
      sn = num2str(count);
    else
      sn = strcat('0',num2str(count));
    end
     
    print('-djpeg75',strcat(basename,sn,'.jpg'));
    count = count + 1;
  else
    pause(2);
  end
     	  
end

close all;

if (storetype == 1)
  % execute script to replay movie.
  FPS = 1; N = 2; % set frame rate and number of times to replay movie
  replay;
  disp('*****************************')
  disp('Movie stored in variable: Mv');
  disp('To display movie again type: replay')
end

more on