% Written by Einat Frishman, Sep. 00.
function [a2,b2] = two_wells_boundary(a1,b1,k1,k2,x);
% This function matches the coefficients of the incoming and outgoing 
% plain waves at a boundary of a potential step.
% the input are the coefficients of the k1 and k2 for both sides of the
% step and a1 a2 are the coefficients of the exp(+-ik1x). the output are the 
% coefficients on the other side of the step of exp(+-ik2x). This is done by 
% matching the wave function and its derivative at the boundary.
%
%
%                           ___________________________
%                          |
%  a1 exp(ik1x)            |       a2 exp(ik2x)  
%         + a2 exp(-ik1x)  |              + a2 exp(-ik2x)   
%    ______________________|
%
%
% The method utiliz
%  matching the wave packet:
% a_1 exp(i k_1 x) + b_1 exp(-i k_1 x) =
%                      a_2 exp(i k_2 x) + b_2 exp(-i k_2 x)
%  matching the derivative:
% k_1 a_1 exp(i k_1 x) - k_1 b_1 exp(-i k_1 x) =
%                      k_2 a_2 exp(i k_2 x) - k_2 b_2 exp(-i k_2 x)
% This can be written in matrix form as:
%   m1 c1= m2 c2        where  c1=[a1,b1] and c2=[a2, b2] 
%    or 
%   c2=inv(m2)m1 c1
% 
%
%

c1    = [a1;b1];
m1    = [exp(i*k1*x) exp(-i*k1*x) ; k1*exp(i*k1*x) -k1*exp(-i*k1*x)];
m2inv = [k2*exp(-i*k2*x) exp(-i*k2*x) ; k2*exp(i*k2*x) -exp(i*k2*x)]*0.5/k2;
c2=m2inv*m1*c1;
a2=c2(1);b2=c2(2);

