%% Aufgabe 1
clearvars
A = [1  2
     3 -4];
B = [ 5
     -5];
A1 = [B A(:,2)];
A2 = [A(:,1) B];
%Teil a
x0 = det(A1)/det(A);
y0 = det(A2)/det(A);
%Teil b
x = linspace(-5,5,101); % oder x=-5:0.1:5   
%zum x-Wert dazugehörigen y-Wert berechnen
%erste Möglichkeit
y1 = (5-1*x)/2;  
y2 = (-5-3*x)/(-4);
plot(x,y1,x,y2,x0,y0,'o'); %graphisch x- und y-Werte dastellen
%zweite Möglichkeit ist etwas kompakter, aber auch komplizierter
y = (B-A(:,1).*[x;x])./A(:,2);
plot(x,y,x0,y0,'o');
%% Aufgabe 2
clearvars
%Teil a
n = 3;%Anzahl Winkel des Polygons
a = 2*pi/n; %Winkel alpha
u = [1
     0]; %Einheitsvektor
R = [cos(a) -sin(a)
     sin(a)  cos(a)]; %Rotationsmatrix um den Winkel alpha
M = [u R*u R*R*u R*R*R*u]; %[P0 P1 P2 P3] %M wird erzeugt durch wiederholte Rotation des Einheitsvektors um alpha
plot(M(1,:),M(2,:),'-o'); %Plot des Polygons
axis equal
%Teil b
n = 5; %n wurde hier erneut definiert!
a = 2*pi/n; %a wurde hier erneut definiert!
R = [cos(a) -sin(a)
     sin(a)  cos(a)]; %R wurde hier erneut definiert!
%erste Möglichkeit
for i = 1:(n+1) %n+1 -> für n Linien braucht man (n+1) Punkte. In diesem Fall ist M(:,1) = M(:,n+1)
    M(:,i) = R^(i-1)*u; %nach jedem Durchlauf i, wird eine weitere Spalte an M hinzugefügt
end
%zweite Möglichkeit
M = u; 
for i = 2:(n+1)%erste Spalte in M ist u, fülle ab Spalte 2 auf
    M(:,i) = R*M(:,i-1);%nehme die vorherrige Spalte in M und rotiere sie mit R um alpha
end
%Teil c
a = 5*pi/180; %a wurde hier erneut definiert!
R = [cos(a) -sin(a)
     sin(a)  cos(a)]; %R wurde hier erneut definiert!
N = R*M;
%Teil d
plot(M(1,:),M(2,:),'-o',N(1,:),N(2,:),'-o');
axis equal
legend('Ursprüngliches Polygon','Gedrehtes Polygon');
%% Aufgabe 3
clearvars
G0 = 214.5; %1/cm
xe = 0.0028;
r0 = 2.666; %Å
h = 6.626070040E-34; %J·s
hq = h/(2*pi); %J·s
c = 299792458; %m/s
kB = 1.38064852E-23; %J/K
u = 1.660539040E-27; %Kg
%Teil a
r = linspace(2.0,4.0,101)*1e-10; %m 
x = r-(r0*1e-10);
w0 = 2*pi*c*(G0*100); %rad/s
mu = (127*127)/(127+127)*u; %Kg
k = w0^2*mu; %Kg/s²=N/m
Vh = (1/2)*k*(x.^2); %J 
De = (G0*100)/(4*xe); %1/m
a = sqrt(k/(2*h*c*De)); %1/m
Va = h*c*De*(1-exp(-a*x)).^2; %J, komponentenweise quadrieren 
plot(x*1e10,Vh/(h*c)/100,x*1e10,Va/(h*c)/100);
xlabel('Kernabstandsänderungen (Å)');
ylabel('Wellenzahl (cm^{-1})');
legend('Harmonisches Potential','Morse-Potential');
%Teil b
T = linspace(1,400,400); %K
vmax = floor((h*c*De)/(hq*(w0/2))-1/2); %für vmax siehe Atkins-Friedman, S. 355, Gl. (10.50) Die Funktion 'floor' rundet ab. 
v = 0:vmax;
Ea = (v+1/2)*hq*w0-(v+1/2).^2*hq*w0*xe; %J, Energie anharmonischer Oszillator
for i = 1:length(T) %für alle Temperaturen
    Za(i) = sum(exp(-(Ea-Ea(1))/(kB*T(i)))); %Zustandssumme anharmonischer Oszillator
end
%erste Möglichkeit
Eh = (v+1/2)*hq*w0; %J, Energie harmonischer Oszillator
for i = 1:length(T)
    Zh(i) = sum(exp(-(Eh-Eh(1))/(kB*T(i)))); %Zustandssumme harmonischer Oszillator
end
plot(T,Zh,T,Za)
xlabel('Temperatur (K)');
ylabel('Zustandssumme');
legend('Harmonisches Potential','Morse-Potential');
%zweite Möglichkeit
Zh = 1./(1-exp(-(hq*w0)./(kB*T))); %siehe PC-II...
plot(T,Zh,T,Za)
xlabel('Temperatur (K)');
ylabel('Zustandssumme');
legend('Harmonisches Potential','Morse-Potential');
%% Aufgabe 4
clearvars
va = [0 1 2 3 4]; %Schwingungsenergieniveau Anfang
vz = va+1; %Schwingungsenergieniveau Ende
DG = [213.3 212.1 210.9 209.7 208.5]; %1/cm
Ir = [100 35.7 12.8 4.6 1.7]/100; 
h = 6.626070040E-34; %J·s
c = 299792458; %m/s
kB = 1.38064852E-23; %J/K
%Teil a
f = polyfit(2*vz,DG,1); %E(v+1)-E(v) = (h/(2*pi))*w0 - 2*(v+1)*(h/(2*pi))*w0*xe
G0 = f(2); %y = m·x+q: m = f(1), q = f(2); %vgl. mit Aufgabe 3...
xe = -f(1)/G0; %vgl. mit Aufgabe 3...
%erste Möglichkeit
xf = linspace(2,10,101);
yf = f(1)*xf+f(2);
plot(2*vz,DG,'o',xf,yf);
xlabel('2·(v+1)');
ylabel('Wellenzahl (cm^{-1})')
legend('Messwerte','Lineare Regression');
%zweite Möglichkeit
xf = linspace(2,10,101);
yf = polyval(f,xf);
plot(2*vz,DG,'o',xf,yf);
xlabel('2·(v+1)');
ylabel('Wellenzahl (cm^{-1})')
legend('Messwerte','Lineare Regression');
%Teil b
G = 0;
for i = 1:length(va)-1
    G(i+1) = G(i)+DG(i); %0-->0, 0-->1, 0-->2, 0-->3, 0-->4
end
f = polyfit(G,log(Ir),1); %I(n+1<--n)/I(1<--0)=exp(-(En-E0)/(kBT))
                          %ln[I(n+1<--n)/I(1<--0)]= -(En-E0)/(kBT)
xf = linspace(0,1000,101); %1/cm
yf = f(1)*xf+f(2);
plot(G,log(Ir),'o',xf,yf)
xlabel('G(n)-G(0) (cm^{-1})');
ylabel('ln[I(n+1<--n)/I(1<--0)]');
legend('Messwerte','Lineare Regression');
kBT = (-1/f(1)); %1/cm
T = (kBT*100)*h*c/kB; %K
fprintf('T = %.1f K = %.1f °C\n',T,T-273.15);