function Xmatrix = ToCoup_simple(junction,varargin)
% Calculation of the coupling matrix of two curved waveguides
% Kuldeep Amarnath
% August, 2002
%
% This program calculates the coupling matrix T using ref [1], 
% as defined below, between two curved waveguides, represented by junction #, 
% at a specific wavelength lambda:
% [a2 b2]' = [t k; -k* t*] [a1 b1]'
% where a and b represent the waveguides and 1 and 2 represent the
% i/p and o/p ports respectively. The following functions accompany
% this m-file:
% a) shape_fun_*(x,varargin) : returns y and local coupling coeffs. c1 & c2, 
%    where y=f(x), x,y being the parametric coords of the two waveguides.
%    * is the junction identifier
%
% [1] Matsuhara, M., Watanabe, A., J. Opt. Soc. Amer., V65, n2, p 163, (1975)


if (nargin ~= 1) & (nargin ~= 2);
    error('Usage: tocoup_simple(junction,<{lambda,R,L,D,H,n1_eff,n2_eff}>)');
end


% get shape function for requested junction
shape_fun = strcat('shape_fun_',int2str(junction));

% call shape_func to initialize
wg_dat = feval(shape_fun,'init',varargin);
x0 = wg_dat(3);
xN = wg_dat(4);
beta_1 = wg_dat(1); % propgn. vector of wg 1 at lambda
beta_2 = wg_dat(2);

% get start and end points
wg_dat = feval(shape_fun,x0,varargin);
y0 = wg_dat(1);
wg_dat = feval(shape_fun,xN,varargin);
yN = wg_dat(1);


% mesh it up!
N = 1000;
Xpos = linspace(x0,xN,N+1);

% start her rollin !
for k=1:N+1
    % get waveguide data
    wg_dat = feval(shape_fun,Xpos(k),varargin);
    Ypos(k) = wg_dat(1);
    c1(k) = wg_dat(2);
    c2(k) = wg_dat(3);
end

% calculate T-matrix components 
T = [1 0; 0 1];
for k=2:N+1
    dx = Xpos(k) - Xpos(k-1);
    dy = Ypos(k) - Ypos(k-1);
    dydx = dy/dx;
    dc1dx = (c1(k)-c1(k-1))/dx;
    dc2dy = (c2(k)-c2(k-1))/dy;
    c1_avg = (c1(k)+c1(k-1))/2;
    c2_avg = (c2(k)+c2(k-1))/2;
    
    B1 = ((beta_2*dydx-beta_1) + j*dc1dx/c1_avg)/2;
    B2 = ((beta_1/dydx-beta_2) + j*dc2dy/c2_avg)/2;
    C1 = sqrt(c1_avg*c2_avg*dydx);
    C2 = sqrt(c1_avg*c2_avg/dydx);
    
    dT(1,1) = (cos(sqrt(B1*B1+C1*C1)*dx) + ...
               j*B1/sqrt(B1*B1+C1*C1)*sin(sqrt(B1*B1+C1*C1)*dx)) * exp(-j*B1*dx);
    dT(2,2) = (cos(sqrt(B2*B2+C2*C2)*dy) + ...
               j*B2/sqrt(B2*B2+C2*C2)*sin(sqrt(B2*B2+C2*C2)*dy)) * exp(-j*B2*dy);
    dT(1,2) = j*c1(k-1) * exp(-j*(beta_2*(Ypos(k-1)-y0)-beta_1*(Xpos(k-1)-x0))) * ...
              1/sqrt(B1*B1+C1*C1) * sin(sqrt(B1*B1+C1*C1)*dx) * exp(-j*B1*dx);
    dT(2,1) = j*c2(k-1) * exp(-j*(beta_1*(Xpos(k-1)-x0)-beta_2*(Ypos(k-1)-y0))) * ...
              1/sqrt(B2*B2+C2*C2) * sin(sqrt(B2*B2+C2*C2)*dy) * exp(-j*B2*dy);

    T = dT * T;
        
end

% now multiply the phase change factor in 
% (not needed for lumped element representation):
% 
% T(1,1) = T(1,1) * exp(-j*beta_1*(xN-x0));
% T(1,2) = T(1,2) * exp(-j*beta_1*(xN-x0));
% T(2,1) = T(2,1) * exp(-j*beta_2*(yN-y0));
% T(2,2) = T(2,2) * exp(-j*beta_2*(yN-y0));

% remember: x-coord is for bus waveguide and y-coord for the ring.
% so interchange along diagonals for conventional structure with ring above bus:
%     
%     0
%     x
%    ---
% for a lumped element coupler:
Xmatrix = [abs(T(2,2)),-j*abs(T(2,1));-j*abs(T(1,2)),abs(T(1,1))];