a=(1.8+6/18)*10^-10;
b=0.8*10^-10;
a0=sqrt((2.*9.1.*10^-31).*(9.*1.6.*10^-19)/(1.055*10^-34)^2);
beta=5;
u0=9.*1.6e-19;
hb = 1.055E-34;


E_f=(2.75+(beta/36)).*1.6e-19.*9;
x = 0:0.0001:4;

g=((1-2*x)./(2*sqrt(x.*(1-x))));

m=sin(a0.*a.*sqrt(x));

n=sinh(a0.*b.*sqrt(1-x));

p=cos(a0*a*sqrt(x));

q=cosh(a0.*b.*sqrt(1-x));

y=(g.*m.*n)+(p.*q);

plot(x,y);

rectangle('Position',[0.2148 -1 0.3019 2], 'FaceColor', [1 0 0 0.3], 'EdgeColor','r');
rectangle('Position',[0.9545 -1 1.1015 2], 'FaceColor', [1 0 1 0.3], 'EdgeColor','m');
rectangle('Position',[2.412 -1 1.588 2], 'FaceColor', [0 1 0 0.3], 'EdgeColor','g');

annotation('textarrow',[0.199853587115666 0.199121522693997],...
    [0.597333333333333 0.464],'String',{'Band 1'});

annotation('textarrow',[0.425329428989751 0.424597364568082],...
    [0.605333333333333 0.471999999999999],'String',{'Band 2'});

annotation('textarrow',[0.759882869692533 0.759150805270864],...
    [0.602666666666666 0.469333333333333],'String',{'Band 3'});

%% Code useful for Task 2
%Values:
formatSpec = 'afirst=%1.11f, bfirst=%1.11f\n';
fprintf(formatSpec,a,b) %Print the



M = 100; %number of iterations
pi=3.141592654;
kstart = -pi./(a+b); %Starting point
kstop = pi./(a+b);   %Stop point

%E=1;
%x= (E/U0);
delta_k = (kstop-kstart)/M; %Discretization step
%%%%E1 = 0.1012

for n=0:1:M           %begin cycle

    k = kstart + n.*delta_k;
    k1(n+1) = k;
    v = cos(k.*a+k.*b);
    v1(n+1) = v;
    funp = @(x) ((1-(2*x))/(2*sqrt(x.*(1-x))) .* (sin(alpha0.*a.*sqrt(x))) .* sinh(alpha0 .* b .* sqrt(1-x)) + (cos(alpha0.*a.*sqrt(x))).*cosh(alpha0.*b.*sqrt(1-x)))-(cos(k*(a+b)));
    x1 = 0.1;
    Solx1 = fzero(funp,x1);     %Solve a generic equation Sin(2x) = Cos(k) numerically
    Sol1(n+1) = Solx1;
    x2 = 2.0;
    Solx2 = fzero(funp,x2);
    Sol2(n+1) = Solx2;
    x3 = 4.0;
    Solx3 = fzero(funp,x3);
    Sol3(n+1) = Solx3;

end
figure
plot(k1,Sol1,'x')                %Plot the graph using the vectors of numerical values as above
axis([kstart kstop 0 6])
hold on;
plot(k1,Sol2,'x')                %Plot the graph using the vectors of numerical values as above
hold on;
plot(k1,Sol3,'x')                %Plot the graph using the vectors of numerical values as above
hold on

xlabel('wavenumber')
ylabel('Energy eV')

plot( [kstart, kstop] , [E_f/u0 , E_f/u0] );

%Task 3

%klst;

E0 = 2.412;
Eneg = 2.416;
Epos = 2.416;

dEdk = (Eneg - 2*E0 + Epos)/((delta_k)^2);

Mc = (hb^(2))/dEdk/u0

E01 = 2.059;
Eneg1 = 2.052;
Epos1 = 2.052;

dEdk1 = (Eneg1 - 2*E01 + Epos1)/((delta_k)^2);

Mv = (hb^(2))/dEdk1/u0