
	

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%					                    %
	%      FILE NAME: phase_wall_well	    %
	%		 			                    %
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% Program to plot wavefunctions and scattering phases
% for a wall and a barrier.
% Written by Einat Frishman, Sep. 00.
% Modified on mar. 00.
%
%  The program solves the eigenstates of the following potential
%
%
%                           v(2)        |
%                         _______       |
%                        |       |      |
%            v(1)        |       | v(3) |
%    ____________________|       |______|
%
%    0                   x1     x2     xf
%
%           
%  The main engine of this program is two_well_boudary see comments there.
%
%

clear all;


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

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

v=[0 12 0];   % potential
%v=[0 v2 0];   % potential
x=[0 x1 x2 xf];

bot=1.02;


f=figure(1);

clf;hold on;
for iframe=1:5;
% for a potential barrier of height 12
  a(1)=-i;b(1)=i;k0=(0.965+(iframe-3)^3*0.025)*k_res; % antisym. wfn.
  e=(k0^2)*ones(1,3);
  k=sqrt(e-v);
  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 height];
  subplot('position',pos);hold on;
% --- calculate coefficients for areas 2,3 (right)
  for j=2:3;
    [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(3)+b(3))/2; cos_phi=i*(a(3)-b(3))/2;
  norm=sqrt(sin_phi^2+cos_phi^2);
  for j=1:3;
    xv=linspace(x(j),x(j+1),200);
    func=( a(j)*exp(i*k(j)*xv)+b(j)*exp(-i*k(j)*xv) ) / norm;
    func=func*1.5; % shallower resonance
    plot(xv,func);plot([xv(1) xv(end)],[v(j) v(j)]);
  end;
  for j=1:2;
    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/c']);
  
  %set(gca,'FontSize',20);
  set(gca,'TickDir','out');
  set(gca,'XLim',[x(1) x(end)]);
  if (iframe ~=3);set(gca,'YLim',[-13 13]);
  else set(gca,'YLim',19/15*[-13 13]); end; % center frame is higher.
  set(gca,'XTick',[0 (x(2)+x(3))/2]);
  if (iframe ~=5);set(gca,'XTickLabel',[]);
  else set(gca,'XTickLabel',['0';'c']);end;
   set(gca,'YTickLabel',[]);
   
 
end; % --- loop over iframe.


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


figure;
subplot(1,2,2);%--------------------------------------------------
clf;     
ncalc=200;
nres_max=4.0; % if phas1 is with    ramp of kv*x2

kv=linspace(0.1*k_res,nres_max*k_res,ncalc); % --- momentum for area 3.
a(1)=-i;b(1)=i;
for icalc=1:ncalc;
  e=(kv(icalc)^2)*ones(1,3);
  k=sqrt(e-v);
% --- calculate coefficients for areas 2,3 (right)
  for j=2:3;
    [a(j) b(j)]=two_wells_boundary(a(j-1),b(j-1),k(j-1),k(j),x(j));
  end;
  sin_phi=(a(3)+b(3))/2; cos_phi=i*(a(3)-b(3))/2;
  norm=sqrt(sin_phi^2+cos_phi^2);
  sin_phi=sin_phi/norm; cos_phi=cos_phi/norm;
  phas1(icalc)=phase(cos_phi+i*sin_phi);
  norm1(icalc)=norm; % --- ratio of inner to asymptotic regions norms.
end;
  subplot('position',[0.2 0.63 0.5 0.32]);hold on;
  plot(kv/k_res,phas1+kv*x2,'-');
  plot(kv/k_res,kv*x2,':');
  set(gca,'XLim',[0 nres_max]);set(gca,'YLim',[-2 14]);
  set(gca,'XTick',[0:nres_max]);
 % set(gca,'FontSize',18);
  set(gca,'TickDir','out');
  set(gca,'YTick',[0 pi 2*pi 3*pi]);
  set(gca,'YTickLabel',[' 0 ';' p ';' 2p';' 3p'],'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;
  plot(kv/k_res,1./norm1,'-');
  set(gca,'XLim',[0 nres_max]);set(gca,'YLim',[0 7]);
  set(gca,'XTick',[0:nres_max]);
  %set(gca,'FontSize',18);
  set(gca,'TickDir','out');
  xlabel('k/k_{res}');
  ylabel('N_{inner}/N_{asymptotic}');
  set(get(gca,'YLabel'),'VerticalAlignment','bot');
	
  return;
  

