function [Conc,t_,Molec] = SimAtmos(dt0,nt)
% SimAtmos simmuliert einfache atmosphärische Reaktionen, die durch die
% Ozonphotolyse initiiert werden und durch Wasser ihren Fortgang im OH finden.
%
% Input:
%   nt    - maximale Zeit (Ende der Simulation); Default: 1 Sekunde
%   dt0   - kleinster initialer Zeitschritt; Default: 10^-7 Sekunden (100ns).
%
% Output:
%   Conc  - Konzentration der Moleküle (1.Index: Konz. z.Zt. t; 2.Index: Molekül) 
%   t_    - Zeit / Sekunden
%   Molec - Spezies (Name des Moleküls)

% Die folgenden Reaktionen sollen behandelt werden:
%   O2 + hv  -> O + O       k0 
%   O3 + hv  -> O2+O1D      k1a = k1*f        f = 0.9     O1D= O(1D) 
%            -> O2+O        k1b = k1*f        f = 0.1     O = O(3P)
%   O + O3   -> 2 O2        k2  = 1.5*10-14
%   O+O2+M   -> O3+M        k3  = 5.4*10-34   für N2 und O2
%   O1D+H2O  -> 2 OH        k4a = 2.3*10-10
%            -> O + H2O     k4b = 4*10-12
%   O1D+O2/N2-> O + O2/N2   k5  = 5*10-11     ähnlich groß für N2 &O2
%   OH + OH  -> H2O + O     k6  = 2*10-12
%   OH+OH+M  -> H2O2+M      k7  = 2.5*10-31
%   OH + O3  -> HO2 + O2    k8  = 5*10-14
%   OH + HO2 -> H2O+ O2     k9  = 5*10-11
%   HO2+ HO2 -> H2O2+O2     k10 = 3*10-12
%   HO2 + O3 -> OH+ 2O2     k11 = 2*10-15
%
% Die Photolyse des Ozons wird über den aktinischen Fluß als unimolekularer
% Zerfall gerechnet.
%

% Zunächst nur die Reaktionen mit Ox:
% A) Blitzlichtphotolyse, d.h. zum Zeitpunkt t=0 sind neben Ozon auch die
% Photolyseprodukte mit einer bestimmten Anfangskonzentration vorhanden.
% Anfangskonzentrationen:

O3   = 1.79E13;                             % Ozonanfangskonzentration.
O2   = 0.21*2.4E19;                         % 21% O2-Anfangskonz. 
O    = 0;                                   % O -Anfangsk. nach Photolyse>0; s.u.
O1D  = 0;
H2O  = 4.8E17;                              % Mit Wasser kommt H ins Spiel.
OH   = 0;
HO2  = 0;
H2O2 = 0;
M    = 2.4E19;                              % Gesamtteilchenzahl.
t    = 0;                                   % Beginn zum Zeitpunkt t=0
% Falls Photolyse instantan ist (Blitzlicht), dann setzt man k0 und k1= 0 und 
% setzt eine Anfangskonzentration für O fest. 
% Vorraussetzung: Photolysedauer klein gegen schnellste Reaktionszeit
k0 = 1E-8;                                  % "stetige" Photolyserate, sonst s.o.
                                            %   O2 + hv  -> O + O
k1 = 5E-6;
k1a= 0.9*k1;                                %   O3 + hv  -> O2+ O1D
k1b= 0.1*k1;                                %            -> O2+ O
k2 = 1.5E-14;                               %   O + O3   -> 2 O2
k3 = 5.4E-34;                               %   O+O2+M   -> O3+ M    für N2 & O2 
k4a = 2.3E-10;                              %   O1D+H2O  -> 2 OH        
k4b = 4E-12;                                %            -> O + H2O     
k5  = 5E-11;                                %   O1D+O2/N2-> O + O2/N2    O2/N2=M, da  ähnlich groß für N2&O2
k6  = 2E-12;                                %   OH + OH  -> H2O + O     
k7  = 2.5E-31;                              %   OH+OH+M  -> H2O2+ M      
k8  = 5E-14;                                %   OH + O3  -> HO2 + O2    
k9  = 5E-11;                                %   OH + HO2 -> H2O + O2     
k10 = 3E-12;                                %   HO2+ HO2 -> H2O2+ O2     
k11 = 2E-15;                                %   HO2 + O3 -> OH + 2O2     
if nargin<2, nt = 1; end                    % maximale Zeit.
if nargin<1, dt0 = 1E-10; end               % kleinster Zeitschritt.
nmax = 100000000;                           % maximale Zahl der Schritte.
s = 10000;
ns = nmax/s;
t_(ns)  = 0;                                % Preallocate sizes of arrays. 
O3_(ns) = 0;
O2_(ns) = 0;
O_(ns)  = 0;
O1D_(ns)= 0;
OH_(ns) = 0;
HO2_(ns)= 0;
H2O2_(ns)= 0;

tic                                         % Beginn der Zeitmessung f. Rechendauer.
i = 1;
while i<nmax
  dO3 = -k1*O3 - k2*O*O3 + k3*O*O2*M;                                      % nur "O"s.
  dO3 = dO3 - k8*OH*O3 - k11*HO2*O3;                                          % mit "H"s.
  dO2 = -k0*O2 + k1*O3 + 2*k2*O*O3 - k3*O*O2*M;                         % nur "O"s.
  dO2 = dO2 + k8*OH*O3 + k9*OH*HO2 + k10*HO2*HO2 + 2*k11*HO2*O3;  % mit "H"s.
  dO  = +2*k0*O2 + k1b*O3 - k2*O*O3 - k3*O*O2*M;                        % nur "O"s.
  dO  = dO + k4b*O1D*H2O + k5*M*O1D + k6*OH*OH;                            % mit "H"s.
  dO1D= +k1a*O3 - k5*M*O1D;
  dO1D= dO1D - (k4a+k4b)*H2O*O1D;
  dOH = 2*k4a*O1D*H2O - ((k6+k7*M)*2*OH + k8*O3 + k9*HO2)*OH + k11*HO2*O3; 
  dHO2= +k8*O3*OH - (k9*OH+2*k10*HO2+k11*O3)*HO2;
  dH2O2=+k7*OH*OH*M + k10*HO2*HO2;
  dt = 0.01*min(abs([O3/dO3,O2/dO2,O/dO,O1D/dO1D,OH/dOH,HO2/dHO2,H2O2/dH2O2]))+dt0;
  O3  = O3  + dO3*dt;                       % Neue Ozo-Konzentration nach t+dt
  O2  = O2  + dO2*dt;
  O   = O   + dO*dt;
  O1D = O1D + dO1D*dt;
  OH  = OH  + dOH*dt;
  HO2 = HO2 + dHO2*dt;
  H2O2= H2O2+ dH2O2*dt;
  t   = t   + dt;
  i=i+1;
  if ~mod(i,s)                              % Speichere jeden s. Datenpunkt
    k = i/s;
    t_(k)  = t;
    O3_(k) = O3;
    O2_(k) = O2;
    O_(k)  = O;
    O1D_(k)= O1D;
    OH_(k) = OH;
    HO2_(k)= HO2;
    H2O2_(k) = H2O2;
  end
  if t>=nt, break, end
end
toc                                         % Ausgabe der Rechendauer.

i,t
figure(1); plot(t_,O3_);   title('Ozon');   % Zeige Ergebnisse.
figure(2); plot(t_,O2_);   title('O2');
figure(3); plot(t_,O_);    title('O');
figure(4); plot(t_,O1D_);  title('O1D');
figure(5); plot(t_,OH_);   title('OH');
figure(6); plot(t_,HO2_);  title('HO2');
figure(7); plot(t_t,H2O2_);title('H2O2');

Conc(:,1) = O3_' ; Molec(1) = {'Ozon'};     % Präpariere Ausgabe der Ergebnisse.
Conc(:,2) = O2_' ; Molec(2) = {'O2'};
Conc(:,3) = O_'  ; Molec(3) = {'O'};
Conc(:,4) = O1D_'; Molec(4) = {'O1D'};
Conc(:,5) = OH_' ; Molec(5) = {'OH'};
Conc(:,6) = HO2_'; Molec(6) = {'HO2'};
