%%%% Ionisationsgrad einer SternatmosphŠre aus purem Wasserstoff
%%% Uebung 1, Uebungsblatt 5
%%%% Saha-Gleichung, Boltzmann-Gleichung
%%%% matlab code, L. Baudis, 21.11.06

% Definition einiger Konstanten


k = 1.38e-16; % Boltzmann Konstante,  in erg K^-1
h=6.63e-27; % Planck Konstante, in erg s
Ei=13.6*1.6e-12; % H-Ionisationsenergie, in erg
me=9.11e-28; % Masse des Elektrons, in g
pi=3.14

% WŠhle den Vektor fŸr den Temperaturbereich

T=[5000:100:20000]; % Temperaturbereich, in K

Pe = 200; % Elektronendruck in dyne/cm^2 = 10^-5 N/cm^2

% Berechne die Partitionsfunktionen

ZI = 2*(1)^2;
ZII = 1; % Proton (H-Ion) hat keine Entartung

% Berechne das VerhŠltnis nII/nI

R_II_I =(2*k.*T*ZII)/(Pe*ZI).*((2*pi*me*k.*T/h^2).^(3/2)).*exp(-Ei./(k*T));

% Berechne das VerhŠltnis nII/ntot = (nII/nI)/(1+nII/nI)

R_II_tot =R_II_I./(1+R_II_I);

% plotte R_II_tot als Funktion der Temperatur

figure(1); clf;
plot(T,R_II_tot, 'b-');
grid on;
xlabel('T [K]'); ylabel('n_{II}/n_{tot}');
axis([5000 20000 0 1.05]);
set(1,'PaperPosition',[0.25 1.5 8 6]); 
print -dpng saha_nI_ntot;

%%%%%%%%%%%%%%%%%%%%%

% Berechne das VerhŠltnis n2/n1 aus der Boltzmann Gleichung

n=1; g1=2*(n)^2;
n=2; g2=2*(n)^2;
DE = 10.2*1.6e-12; % Energiedifferenz in erg

R21 = (g2/g1)*exp(-DE./(k*T));

% plotte R21 als Funktion der Temperatur

figure(2); clf;
plot(T,R21, 'r-');
grid on;
xlabel('T [K]'); ylabel('n_{2}/n_{1}');
%axis([5000 20000 0 1.05]);
set(2,'PaperPosition',[0.25 1.5 8 6]); 
print -dpng boltz_n2_n1;

%%%

% Berechne das VerhŠltnis nI/ntot = 1/(1+nII/nI) durch Kombination von Saha- und Boltzmann-Gleichung

R_I_tot =1./(1+R_II_I);

% Berechne das VerhŠltnis n2/(n1+n2)

R22 = R21./(1+R21);

% Berechne das VerhŠltnis n2/(ntot)

R2tot = R22.*R_I_tot;

% plotte R2tot (n2/ntot) als Funktion der Temperatur

figure(3); clf;
plot(T,R2tot, 'r-');
grid on;
xlabel('T [K]'); ylabel('n_{2}/n_{tot}');
%axis([5000 20000 0 1.05]);
set(3,'PaperPosition',[0.25 1.5 8 6]); 
print -dpng saha_boltz_n2_ntot;

% plotte alles zusammen

figure(4); clf;

subplot(3,1,1);
plot(T,R_II_tot, 'b-');
grid on;
xlabel('T [K]'); ylabel('n_{II}/n_{tot}');
axis([5000 20000 0 1.05]);

subplot(3,1,2);
plot(T,R21, 'r-');
grid on;
xlabel('T [K]'); ylabel('n_{2}/n_{1}');
%axis([5000 20000 0 1.05]);

subplot(3,1,3);
plot(T,R2tot, 'r-');
grid on;
xlabel('T [K]'); ylabel('n_{2}/n_{tot}');
%axis([5000 20000 0 1.05]);
set(4,'PaperPosition',[0.25 1.5 6 10]); 
print -dpng saha_boltz_all;