%% Kronig-Penney Modell

% Vers. 1.1; 08.2015 
% Copyright W.Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB 
% 2016
%% Visualisierung der transzendenten Gleichung
qa = linspace(-6*pi,6*pi,500);
vorf = 5;%;7.5;%10;%10;
y = cos(qa) + vorf*sin(qa)./qa;  % y=cos(k*a)
plot(qa,y),shg
hold on
plot([qa(1),qa(end)],[1,1],[qa(1),qa(end)],[-1,-1])
ynull = zeros(size(y));
ynull(abs(y)>=1)=NaN;
plot(qa,ynull)
hold off
%%
qaErlaubt = qa(abs(y)<1);
figure,
plot(qaErlaubt.^2),shg
%%    Bestimmen der Schnittpunkt abs(y) = 1
yg1 = @(x) cos(x)+vorf*sin(x)./x-1;
yg2 = @(x) cos(x)+vorf*sin(x)./x+1;
%yn1 = [-6*pi,-4*pi,-2*pi,+2*pi,+4*pi,+6*pi];
%yn2 = [-5*pi,-3*pi,-pi,pi,3*pi,5*pi];
yn1 = [-6*pi,-4*pi,-2*pi];%,+2*pi,+4*pi,+6*pi];
yn2 = [-5*pi,-3*pi,-pi];%,pi,3*pi,5*pi];
%yg1(yn2)
%yg2(yn1)
for n=1:length(yn1)
    ym1(n) = fzero(yg1,yn2(n));
    ym2(n) = fzero(yg2,yn1(n));
end
%%    Plotten von (q*a)^2 versus k*a
figure, hold on, shg
ycontrol = [];
for n=1:length(yn1)
    win1 = linspace(yn1(n),ym2(n));
    win2 = linspace(ym1(n),yn2(n));
    E1 = win1.^2;
    y1 = cos(win1) + vorf*sin(win1)./win1;
    ycontrol=[ycontrol;y1];
    ka1 = real(acos(y1));
    plot(ka1,E1,-ka1,E1,[-ka1(end),ka1(end)],[E1(end),E1(end)]),shg

    E2 = win2.^2;
    y2 = cos(win2) + vorf*sin(win2)./win2;
    ycontrol=[ycontrol;y2];
    ka2 = real(acos(y2));
    plot(ka2,E2,-ka2,E2,[-ka2(end),ka2(end)],[E2(end),E2(end)]),shg
end
hold off