%% Anharmonischer Oszillator - komplexe Koordinatenrotation
% Vers. 1.2 03/2016
%
% Simulation physikalischer Systeme - Computational Physics mit MATLAB 
% 2016

%%
clear
ndim = 100;   % Anzahl der Basiselemente 
wieviel=3;    % Anzahl der zu berechnenden Zustaende

kl = [0.05,0.08,0.12,0.15,0.17,0.185, 0.2,0.225,0.25,0.275,0.3,0.325,0.35,0.375,0.4,0.425,0.45,0.475,0.5,0.55,0.6,0.65,0.7,0.8,0.9,1];
%    kl: kappa-Werte

% Die Ergebnisse liegen in erg1, erg2, erg3
% Der Aufabu der Arrays ist:
% ergx = [kappa,real(E),abweichung,imag(E),abweichung,winkel min, max]
erg3 = [0,2.5,0,0,0,0,0];
erg2 = [0,1.5,0,0,0,0,0];
erg1 = [0,0.5,0,0,0,0,0];

% Zusaetzlich finden sich die Ergebnisse in einer Tabelle: tergx

for kappa = kl

Nsqrt=zeros(ndim);  % Allokation
Emat = Nsqrt;
x3 = Nsqrt;

qz = [0:ndim-1].';  % Eigenloesungen
N = diag(qz);
Nsqrt = sqrt(N);

%E0 = N + 1/2*eye(ndim);

ad = [Nsqrt(:,2:end),zeros(ndim,1)];  % Erzeugungsoperator
a = ad';

x = sqrt(1/2)*(ad + a);               % Ortsoperator
p = 1i*sqrt(1/2)*(ad - a);            % Impulsoperator
xq0 = x^2;
pq0 = p^2;
xc0 = x*xq0;

nlauf=0;
Enlauf = [];

wl = linspace(1,30,75);              % Komplexe Koordinatenrotation

for winkel=wl
    nlauf=nlauf+1;
    wr = winkel*pi/180;
    xq = exp(2*1i*wr)*xq0;
    pq = exp(-2*1i*wr)*pq0;
    xc = exp(3*1i*wr)*xc0;
    Emat = 1/2*(xq+pq)+1/2*kappa*xc;
    E=eigs(Emat,wieviel,'sm');      % Eigenwertberechnung
    Enlauf = [Enlauf,E];
end
%
E3 = Enlauf(3,:);                   % Grundzustand
auswahl3 = imag(E3)<=0;             % H"aufungspunkt berechnen
[ra3,rb3]=min(abs(diff(real(E3(auswahl3)))));
index = 1:length(auswahl3);
index = index(auswahl3);
%index(rb3)
[ia3,ib3]=min(abs(diff(imag(E3(auswahl3)))));
%index(ib3)
vonwo =min([index(rb3),index(ib3)])-2;
biswo =max([index(rb3),index(ib3)])+2;
if vonwo<=0
    vonwo=1;
end
if biswo >max(index)
    biswo=max([index(rb2),index(ib2)]);
end
abr3=max(real(E3(vonwo:biswo)))-min(real(E3(vonwo:biswo)));
abi3=max(imag(E3(vonwo:biswo)))-min(imag(E3(vonwo:biswo)));
gamma3=mean(imag(E3(vonwo:biswo)));
eig3=mean(real(E3(vonwo:biswo)));
zuerg1 = [kappa,eig3,abr3,gamma3,abi3,wl(vonwo),wl(biswo)];
erg1= [erg1;zuerg1];
%
E2 = Enlauf(2,:);                  % 1. angeregter Zustand
auswahl2 = imag(E2)<=0;
[ra2,rb2]=min(abs(diff(real(E2(auswahl2)))));
index = 1:length(auswahl2);
index = index(auswahl2);
%index(rb2)
[ia2,ib2]=min(abs(diff(imag(E2(auswahl2)))));
%index(ib2)
vonwo =min([index(rb2),index(ib2)])-2;
biswo =max([index(rb2),index(ib2)])+2;
if vonwo<=0
    vonwo=1;
end
if biswo >max(index)
    biswo=max([index(rb2),index(ib2)]);
end
%E2(vonwo:biswo)
abr2=max(real(E2(vonwo:biswo)))-min(real(E2(vonwo:biswo)));
abi2=max(imag(E2(vonwo:biswo)))-min(imag(E2(vonwo:biswo)));
gamma2=mean(imag(E2(vonwo:biswo)));
eig2=mean(real(E2(vonwo:biswo)));
zuerg2 = [kappa,eig2,abr2,gamma2,abi2,wl(vonwo),wl(biswo)];
erg2= [erg2;zuerg2];
%
E1 = Enlauf(1,:);             % 2. angeregter Zustand
auswahl1 = imag(E1)<=0;
%plot(wl(auswahl1),real(E1(auswahl1)),'*',wl(auswahl1),imag(E1(auswahl1)),'p'),shg
[ra1,rb1]=min(abs(diff(real(E1(auswahl1)))));
index = 1:length(auswahl1);
index = index(auswahl1);
%index(rb1)
[ia1,ib1]=min(abs(diff(imag(E1(auswahl1)))));
%index(ib1)
vonwo =min([index(rb1),index(ib1)])-2;
biswo =max([index(rb1),index(ib1)])+2;
if vonwo<=0
    vonwo=1;
end
if biswo >max(index)
    biswo=max([index(rb1),index(ib1)]);
end
%E1(vonwo:biswo)
abr1=max(real(E1(vonwo:biswo)))-min(real(E1(vonwo:biswo)));
abi1=max(imag(E1(vonwo:biswo)))-min(imag(E1(vonwo:biswo)));
gamma1=mean(imag(E1(vonwo:biswo)));
eig1=mean(real(E1(vonwo:biswo)));
zuerg3 = [kappa,eig1,abr1,gamma1,abi1,wl(vonwo),wl(biswo)];
erg3= [erg3;zuerg3];
%


% x = linspace(-4,4);
% potential =[1/2*kappa,1/2,0,0];
% plot(x,polyval(potential,x))
% potabl = polyder(potential);
% polymax = roots(potabl);
% polymax = polymax(polymax<0)

end
%% Tabellen erstellen
Variablennamen = {'kappa','E_Realteil','AbweichungReal','E_Imaginaerteil','AbweichungImag','Winkel_min','und_max'};
terg1 = array2table(erg1,'VariableNames',Variablennamen);
terg2 = array2table(erg2,'VariableNames',Variablennamen);
terg3 = array2table(erg3,'VariableNames',Variablennamen);