function [x,y,xc,yc,nx,ny,eps] = coupler(n,h,rh,d,rw,side,dx,dy);

% This function creates an index mesh for the finite-difference
% mode solver.  The function will accommodate a generalized 
% N-layer rib waveguide structure.  (Note: channel waveguides can
% also be treated by selecting the parameters appropriately.) 
% 
% USAGE:
% 
% [x,y,xc,yc,nx,ny,eps] = waveguidemesh(n,h,rh,d,rw,side,dx,dy)
%
% INPUT
%
% n - indices of refraction for layers in waveguide
% h - height of each layer in waveguide
% rh - height of waveguide feature
% d - center-to-center separation between waveguides
% rw - full-width of waveguide
% side - excess space to the right of waveguide
% dx - horizontal grid spacing
% dy - vertical grid spacing
% 
% OUTPUT
% 
% x,y - vectors specifying mesh coordinates
% xc,yc - vectors specifying grid-center coordinates
% nx,ny - size of index mesh
% eps - index mesh (n^2)
%
% AUTHOR:  Thomas E. Murphy (tem@umd.edu)

ih = round(h/dy);
irh = round (rh/dy);
irw = round (rw/dx);
iside = round (side/dx);
isideleft = round((d-rw)/(2*dx));
iside = round (side/dx);
nlayers = length(h);

nx = isideleft+irw+iside+1;
ny = sum(ih)+1;

xc = (1:(nx-1))'*dx - dx/2;
yc = (1:(ny-1))*dy - dy/2;
x = (0:(nx-1))'*dx;
y = (0:(ny-1))*dy;

eps = zeros(nx-1,ny-1);

iy = 1;

for jj = 1:nlayers,
  for i = 1:ih(jj),
	eps(:,iy) = n(jj)^2*ones(nx-1,1);
	iy = iy+1;
  end
end

iy = sum(ih)-ih(nlayers);
for i = 1:irh,
   eps(1:isideleft,iy) = n(nlayers)^2*ones(isideleft,1);
   eps(isideleft+irw+1:isideleft+irw+iside,iy) = n(nlayers)^2*ones(iside,1);
   iy = iy-1;
end

nx = length(xc);
ny = length(yc);