clear all
close all


%%
S=0.2; %distance between wells
% other values hardcoded into the function f2

%%
%root finding and ploting parameters.
Emax=10;%Maximum energy When E>V0 we do not have a bound state
npoints=5000;%number of points=npoints+1 in plot of f2
dE=Emax/npoints;%energy division in plot of f2
Epoints=[0:dE:Emax];%x variable for plot
f2plot=0*Epoints;%yvarable for plot of f2
f1plot=0*Epoints;%In case you want to plot f1 as well
nregions=1000; % of regions which the x axis is divided into for root search
%Must have only one root of f2 in each region.
% need nregions  large enough so that the region containg the
%lowest root is isolated from 0 
regionends=[0:Emax/nregions:Emax];% regions for root search

%%
% finding roots of single well
Esinglew=rootsoff(0.6,regionends); % single well width 0.6 nm
Esingle2w=rootsoff(2*0.6,regionends);% single well width 1.2 nm





%% ploting f1(E) f2(E) to get an idea of what is going on

for j=1:length(f2plot)%generating f1,f2 for ploting
    f2plot(j)=real(f2(dE*j,S));
    f1plot(j)=real(f1(dE*j,0.6));
end

% % plotting single well problem
%  figure 
%  plot(Epoints,f1plot)
%  hold on
%  plot([0 Emax],[0 0],'k')
%  plot(Esinglew,0*Esinglew,'ok','LineWidth',2,'MarkerFaceColor','r')
% plot(Epoints,f2plot,'r:','LineWidth',2)
% axis([0,10,-50,200])

%% plotting double well function
figure
plot(Epoints,f2plot,'LineWidth',2)
hold on
plot([0 Emax],[0 0],'k')%zero line
%set(gca,'XTick',regionends)% x ticks mark out the regions
%plot(X,0.1*dP(X),'r')
xlabel('$E$ (eV)','interpreter','latex','FontSize', 18)
ylabel('$f_1(E)$ nm$^{-1/2}$','interpreter','latex', 'FontSize', 18)
axis([0,10,-50,1e3])
Elevels=rootsoff2(S,regionends); %find roots at S
plot(Elevels,0*Elevels,'ok','LineWidth',2,'MarkerFaceColor','r') %put roots on plot
print( '-depsc2', 'Part2f2.eps' )% save plot
Elevels%print to screen


Region12=[0.2:5e-4:0.9];% from inspection of above roots 1&2 are always in here

% %%
 %how do the roots depend on S?
 
%ElevelsL=rootsoff2(0.0,regionends);% test using these to decide on number
%of regions
%ElevelsR=rootsoff2(0.5,regionends);%
 ns=10;% numbe of different S values
 slist=[0.00:0.45/ns:0.5];% list of S values
 E1=0*slist;% blank lists to load E1 into
 E2=0*slist;% blank lists to load E2 into
for kk=1:length(E1)% Generating E1&E2 
    
    ElofS=rootsoff2(slist(kk),Region12);
    E1(kk)=ElofS(1);
    E2(kk)=ElofS(2);
end
 

%% Ploting E1 vs S
 figure
 plot(slist,E1,'o:','LineWidth',2)
 hold on
 plot([0,0.5],Esinglew(1)*[1,1],'r--')%annotating graph with single well energies 
  plot([0,0.5],Esingle2w(1)*[1,1],'r-.')
 xlabel('$S\,$ (nm)','interpreter','latex','FontSize', 18)
 ylabel('$E_1\,$ (eV)','interpreter','latex', 'FontSize', 18)
 print( '-depsc2', 'Part2E1.eps' )% save plot
 
 %% E2 vs S
  figure
  plot(slist,E2,'o:','LineWidth',2)
  hold on
   plot([0,0.5],Esinglew(1)*[1,1],'r--')%annotating graph with single well energies 
  plot([0,0.5],Esingle2w(2)*[1,1],'r-.')
   xlabel('$S\,$ (nm)','interpreter','latex','FontSize', 18)
 ylabel('$E_2\,$ (eV)','interpreter','latex', 'FontSize', 18)
  print( '-depsc2', 'Part2E2.eps' )% save plot
  
  
  figure %ploting the energy difference (not included in solutions pdf)
 semilogy(slist,E2-E1,'o:','LineWidth',2)
  hold on
   plot([0,0.5],(Esingle2w(2)-Esingle2w(1))*[1,1],'r--')
   xlabel('$S\,$ (nm)','interpreter','latex','FontSize', 18)
 ylabel('$\Delta E_{12}\,$ (eV)','interpreter','latex', 'FontSize', 18)
  print( '-depsc2', 'Part2DE2.eps' )% save plot

 
function roots=rootsoff2(s,regionends)% returns roots of f2 
%using hybrid bisection/secant solver when
%looks for roots between regionends list, Em is the maximum energy
%s is the well seperation
nr=length(regionends)-1;%number of regions
roots=nan*[1:1:nr];%Blank list to load roots into 
 %(will use isnan to remove unfilled elements)

maxiter=100;%max itterations
toleranceb=1e-5;%tolerence of bisection
tolerances=1e-9;%tolerence of secant method
 for i=1:nr % loop over regions
     
         
       
    if f2(regionends(i),s)*f2(regionends(i+1),s)<0% if f2 has a difrent sign 
%         %on each end of a region then that region contains a root.
         
      
        %%bisection method
        accuracy=100*toleranceb;%set acuracy bigger than tollerence
        xL=regionends(i);
        xR=regionends(i+1);
        f2L=f2(xL,s);% store f2 values so we can save time later by not recalculating them
        iter=1;
        while accuracy> toleranceb&iter<maxiter% bisection steps
            iter=iter+1;
            xmiddle=(xL+xR)/2;
            f2m=f2(xmiddle,s);
            if f2L*f2m<0
                xR=xmiddle;
            else
            xL=xmiddle;
            f2L=f2m;
            end
            
            accuracy=abs(f2m);
        end
        
        %% secant method
        f2R=f2(xR,s);
        if abs(f2R)<abs(f2L)%make x the end with smallest |f2(x)| x is x_i
            % and xo the end with |f2(xo)| biggest (xo will be x_{i-1})
            x=xR;
            f2x=f2R;% trying not to wast our calulations of f2
            xo=xL;
            f2xo=f2L;
            
        else
            x=xL;
            f2x=f2L;
            xo=xR;
            f2xo=f2R;
        end
        
        while accuracy> toleranceb&iter<maxiter% secant method iteration
            iter=iter+1;
            xn=x-(x-xo)/(f2x-f2xo);
            f2xn=f2(xn,s);
            if abs(f2xn)<abs(f2x)
                xo=x;
                f2xo=f2x;
                x=xn;
                f2x=f2xn;
            else
                xo=xn;
                f2xo=f2xn;
                
            end
            accuracy=f2x;
        
        end
        
        
        roots(i)=x;
         
         
        
     end

          
        
     end
      



roots=roots(isnan(roots)==0);% removes NaNs from the vector roots
% and print roots to command window.

end

function roots=rootsoff(s,regionends)% returns roots of f1 
%using hybrid bisection/secant solver when (could use function handles here
%but Im to lazy)
%looks for roots between regionends list, Em is the maximum energy
%s is the well width
nr=length(regionends)-1;%number of regions
roots=nan*[1:1:nr];%Blank list to load roots into 
 %(will use isnan to remove unfilled elements)

maxiter=100;%max itterations
toleranceb=1e-5;%tolerence of bisection
tolerances=1e-9;%tolerence of secant method
 for i=1:nr % loop over regions
     
         
       
    if f1(regionends(i),s)*f1(regionends(i+1),s)<0% if f1 has a difrent sign 
%         %on each end of a region then that region contains a root.
         
      
        %%bisection method
        accuracy=100*toleranceb;%set acuracy bigger than tollerence
        xL=regionends(i);
        xR=regionends(i+1);
        f1L=f1(xL,s);% store f1 values so we can save time later by not recalculating them
        iter=1;
        while accuracy> toleranceb&&iter<maxiter% bisection steps
            iter=iter+1;
            xmiddle=(xL+xR)/2;
            f1m=f1(xmiddle,s);
            if f1L*f1m<0
                xR=xmiddle;
            else
            xL=xmiddle;
            f1L=f1m;
            end
            
            accuracy=abs(f1m);
        end
        
        %% secant method
        f1R=f1(xR,s);
        if abs(f1R)<abs(f1L)%make x the end with smallest |f1(x)| x is x_i
            % and xo the end with |f1(xo)| biggest (xo will be x_{i-1})
            x=xR;
            f1x=f1R;% trying not to wast our calulations of f1
            xo=xL;
            f1xo=f1L;
            
        else
            x=xL;
            f1x=f1L;
            xo=xR;
            f1xo=f1R;
        end
        
        while accuracy> toleranceb&iter<maxiter% secant method iteration
            iter=iter+1;
            xn=x-(x-xo)/(f1x-f1xo);
            f1xn=f1(xn,s);
            if abs(f1xn)<abs(f1x)
                xo=x;
                f1xo=f1x;
                x=xn;
                f1x=f1xn;
            else
                xo=xn;
                f1xo=f1xn;
                
            end
            accuracy=f1x;
        
        end
        
        
        roots(i)=x;
         
         
        
     end

          
        
     end
      



roots=roots(isnan(roots)==0);% removes NaNs from the vector roots
% and print roots to command window.

end



function value=f2(X,s)%[ X= energy (eV) s= width (nm)] produces f2(E,S,)

Vo=10; %potential height in the forbidden region (in eV)
W=0.6; %width of well (in nm)
C=5.123167;% \sqrt{2m.(1 eV)/\hbar^2} (in nm^{-1})
alpha=sqrt(X)*C;
beta=sqrt(Vo-X)*C;
A=[cos(alpha*W),sin(alpha*W)/alpha;-alpha*sin(alpha*W),cos(alpha*W)];
F=[cosh(beta*s),sinh(beta*s)/beta;beta*sinh(beta*s),cosh(beta*s)];
value=[beta,1]*A*F*A*[1;beta];

end

function value=f1(X,W)% X=energy (eV), w= width (nm)


Vo=10; %potential height in the forbidden region (in eV)

C=5.123167;% \sqrt{2m.(1 eV)/\hbar^2} (in nm^{-1})
alpha=sqrt(X)*C;
beta=sqrt(Vo-X)*C;
A=[cos(alpha*W),sin(alpha*W)/alpha;-alpha*sin(alpha*W),cos(alpha*W)];

value=[beta,1]*A*[1;beta];

end

