
	

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%					                    %
	%    FILE NAME: phase_double_barrier	%
	%		 		                       	%
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




% Program to plot wavefunctions and scattering phases
% for a double-barrier.
% Written by Einat Frishman, Sep. 00.
%
%  The program solves the eigenstates of the following potential
%
%
%                           v(2)         v(4)
%                         _______        _______
%                        |       |      |       |    
%            v(1)        |       | v(3) |       |  v(5)
%    ____________________|       |______|       |_________________
%
%    -xf               -x2     -x1      x1      x2                xf
%
%           

%  The main engin of this program is two_wells_boundary.m
%               see comments there.
%



clear all;

xf=6;
x_mid=1.75;
x2=x_mid+0.25;
x1=x_mid-0.25;

k_res=pi/(2*x_mid);
e_res=k_res^2;

%v=[0 8 0 8 0];   % potential
v=[0 12 0 12 0];
n_res=sqrt(8)/k_res;

x=[-xf -x2 -x1 x1 x2 xf];
% x(1) is initial x, x(end) is final x.
% x(2) to x(end-1) are boundaries.

% \psi = a*e(ikx+if) + b*e(-ikx-if)
%      = a*e(if)*e(ikx) + b*e(-if)*e(-ikx)

for ifig=[1 3];
bot=1.02;



figure(ifig);
f=figure(ifig);


clf;hold on;
for iframe=1:5;
switch ifig;
% for a potential barrier of height 12
case 1, a(3)=1; b(3)=1;k0=(0.97+(iframe-3)^3*0.03)*k_res; % symmetric wfn.
case 3, a(3)=-i;b(3)=i;k0=(1.93+(iframe-3)^3*0.05)*k_res; % antisym. wfn.
end;
  e=(k0^2)*ones(1,5);
  k=sqrt(e-v);
%%%  bot=bot-0.18;
  if (iframe ~=3);bot=bot-0.18;else bot=bot-0.22;end;
  if (iframe ~=3);height=0.15;else height=0.19;end;
%%%  pos=[0.3 bot 0.35 0.15];
  pos=[0.3 bot 0.35 height];
  subplot('position',pos);hold on;
% --- calculate coefficients for areas 1,2 (left)
  for j=2:-1:1;
    [a(j) b(j)]=two_wells_boundary(a(j+1),b(j+1),k(j+1),k(j),x(j+1));
  end;
% --- calculate coefficients for areas 4,5 (right)
  for j=4:5;
    [a(j) b(j)]=two_wells_boundary(a(j-1),b(j-1),k(j-1),k(j),x(j));
  end;

% --- norm of asymptotic region when norm of inner region is 1 :
  sin_phi=(a(5)+b(5))/2; cos_phi=i*(a(5)-b(5))/2;
  norm=sqrt(sin_phi^2+cos_phi^2);
  for j=1:5;
    xv=linspace(x(j),x(j+1),200);

    func=( a(j)*exp(i*k(j)*xv)+b(j)*exp(-i*k(j)*xv) ) / norm;
    if (ifig==1); func=func*1.0;end; % sharp resonance
    if (ifig==3); func=func*1.5;end; % shallower resonance

    plot(xv,func);plot([xv(1) xv(end)],[v(j) v(j)]);
  end;
  for j=1:4;
    plot([x(j+1) x(j+1)],[v(j),v(j+1)]);
  end;

  plot([x(1) x(end)],[e(1) e(end)],':');
  text(xf*1.3,0,['k=',num2str(k(3)/k_res),'* \pi/2c']);
  set(gca,'TickDir','out');
  set(gca,'XLim',[x(1) x(end)]);
  if (ifig==1);
  if (iframe ~=3);set(gca,'YLim',[-3 17]);
  else set(gca,'YLim',19/15*[-3 17]); end;end; % center frame is higher.
  if (ifig==3);
  if (iframe ~=3);set(gca,'YLim',[-13 13]);
  else set(gca,'YLim',19/15*[-13 13]); end;end; % center frame is higher.
  set(gca,'XTick',[(x(2)+x(3))/2 (x(4)+x(5))/2]);
  if (iframe ~=5);set(gca,'XTickLabel',[]);
  else set(gca,'XTickLabel',['-c';' c']);end;
   set(gca,'YTickLabel',[]);
end; % --- loop over iframe.

end; % --- loop over ifig.





% ============================ phase as a function of energy:

f3=figure;
clf;
ncalc=200;
nres_max=7.5; % if phas1 is without ramp of kv*x2
nres_max=5.0; % if phas1 is with    ramp of kv*x2

kv=linspace(0.1*k_res,nres_max*k_res,ncalc); % --- momentum for area 3.
for isym=1:2;
if (isym==1);a(3)=1;b(3)=1;else a(3)=-i;b(3)=i;end;
for icalc=1:ncalc;
  e=(kv(icalc)^2)*ones(1,5);
  k=sqrt(e-v);
% --- calculate coefficients for areas 4,5 (right)
  for j=4:5;
    [a(j) b(j)]=two_wells_boundary(a(j-1),b(j-1),k(j-1),k(j),x(j));
  end;

% --- \psi =     a   e(ikx)    +     b   e(-ikx)
% ---      = (a+b)/2 cos kx    + i(a-b)/2 sin kx

  sin_phi=(a(5)+b(5))/2; cos_phi=i*(a(5)-b(5))/2;
  norm=sqrt(sin_phi^2+cos_phi^2);
  sin_phi=sin_phi/norm; cos_phi=cos_phi/norm;
  if (isym==1);phas1(icalc)=phase(cos_phi+i*sin_phi)-pi/2;
  else phas1(icalc)=phase(cos_phi+i*sin_phi);end;
  norm1(icalc)=norm; % --- ratio of inner to asymptotic regions norms.
end;

  subplot('position',[0.2 0.63 0.5 0.32]);hold on;
% plot absolute phase at x=x2: therefore add kv*x2 to phas1.
% without kv*x2 , phas1 is the difference between phases with
% and without the barrier.
  if (isym==1);plot(kv/k_res,phas1+kv*x2);
  else plot(kv/k_res,phas1+kv*x2,'--');end;

% plot the background phase of k*x2.

  if (isym==2);plot(kv/k_res,kv*x2,':');end;

% set(gca,'XLim',[0 nres_max]);set(gca,'YLim',[-3 0]);
  set(gca,'XLim',[0 nres_max]);set(gca,'YLim',[-2 8]);
 % set(gca,'FontSize',18);
  set(gca,'TickDir','out');

% if phas1 is without ramp of kv*x2

% set(gca,'YTick',[-pi/2 0]);
% set(gca,'YTickLabel',['-p/2';'  0 '],'FontName','Symbol');

% if phas1 is with    ramp of kv*x2

set(gca,'YTick',[0 pi 2*pi]);
set(gca,'YTickLabel',[' 0 ';' p ';' 2p'],'FontName','Symbol');

  xlabel('k/k_{res}','FontName','Times');
  ylabel('\chi (rad)','FontName','Times');

% --------------------------------------------
%  RATIO OF INNER TO ASYMPTOTIC REGIONS NORMS
% --------------------------------------------

  subplot('position',[0.2 0.15 0.5 0.32]);hold on;
  if (isym==1);plot(kv/k_res,1./norm1);
  else plot(kv/k_res,1./norm1,'--');end;
  set(gca,'XLim',[0 nres_max]);set(gca,'YLim',[0 7]);
  set(gca,'TickDir','out');
  xlabel('k/k_{res}');
  ylabel('N_{inner}/N_{asymptotic}');

end; % --- loop over isym



f4=figure;

clf;
ncalc=2000;
nres_max=7.5; % if phas1 is without ramp of kv*x2
nres_max=5.0; % if phas1 is with    ramp of kv*x2
nres_max=nres_max+5; 

kv=linspace(0.1*k_res,nres_max*k_res,ncalc); % --- momentum for area 3.
a(5)=1;b(5)=0;
for icalc=1:ncalc;
  e=(kv(icalc)^2)*ones(1,5);
  k=sqrt(e-v);
  for j=4:-1:1;
     [a(j) b(j)]=two_wells_boundary(a(j+1),b(j+1),k(j+1),k(j),x(j+1));
  end;

  R(icalc)=b(1)/a(1);%reflection coefficient
  T(icalc)=a(5)/a(1);%transmission coefficient
end;



   subplot('position',[0.2 0.63 0.5 0.32]);hold on;
   axis([0 10 0 2]);
   plot(kv/k_res,abs(T).^2);
   plot(kv/k_res,abs(R).^2,'r');
   
   set(gca,'XLim',[0 nres_max]);set(gca,'YLim',[0 2]);
   set(gca,'TickDir','out');


  set(gca,'YTick',[0 1 2]);
  xlabel('k/k_{res}','FontSize',20,'FontName','Times');
  ylabel('|T(k)|^2,|R(k)|^2','FontSize',15,'FontName','Times');




