    
    ## H3O - water system:
    ## ------------------------
    
    p_react_mbar = 9;
    T_react_K = 273+20;
    p_H2O_mbar = 100/1013*p_react_mbar;
    
    # plot:
    w = 600;
    h = 500;
    
    kB = 1.38E-23; % J/K
    eV = 1.602E-19; % eV/J
    u = 1.66E-27; % kg/u
    R = 8.314; % J/mol K
    m_air = 29;
    m_water = 18;
    
    % Primary ion distribution - H3O+:
    enthalpies = -[31.6 19.5 17.9 12.7 11.6 10.7]*4.184; % kJ/mol, Lau et al., 1982
    entropies = -[24.3 21.7 28.4 23.4 25 26.1]*4.184; %J/mol K, Lau et al., 1982
    mobilities = [2.76 2.28 2.13 2 1.9 1.8]; % cm2/Vs, Dotan et al.,1976 and approx. values
    masses = [19 37 55 73 91 109 127]; % u
    E_over_N = 0.1:0.1:120; % Td
    E_120Td = 120*p_react_mbar*100/(1.38*T_react_K); % V/cm at 120 Td
    v_ions_120Td = mobilities*E_120Td*1013/p_react_mbar*T_react_K/273.15/100; %m/s
    v_ions = 3*E_over_N*p_react_mbar*100/(1.38*T_react_K)*1013/p_react_mbar*T_react_K/273.15/100; %m/s
    T_eff = zeros(size(E_over_N,2),size(mobilities,2));
    for i=1:size(mobilities,2)
        %     T_eff(:,i) = 2/(3*kB)*(1/2*30*1.66E-27*(v_ions_120Td(1,i)/120*E_over_N').^2+3/2*kB*T_react);
        T_eff(:,i) = 2/(3*kB)* ( 0.5*(masses(i)+m_air)*m_water*u/(masses(i)+m_water)*(v_ions_120Td(1,i)/120*E_over_N').^2+3/2*kB*T_react_K );
    end
    K01 = exp(-enthalpies(1)./(R/1000*T_eff(:,1)')+entropies(1)/R)/1013; %1/mbar
    K12 = exp(-enthalpies(2)./(R/1000*T_eff(:,2)')+entropies(2)/R)/1013;
    K23 = exp(-enthalpies(3)./(R/1000*T_eff(:,3)')+entropies(3)/R)/1013;
    K34 = exp(-enthalpies(4)./(R/1000*T_eff(:,4)')+entropies(4)/R)/1013;
    K45 = exp(-enthalpies(5)./(R/1000*T_eff(:,5)')+entropies(5)/R)/1013;
    K56 = exp(-enthalpies(6)./(R/1000*T_eff(:,6)')+entropies(6)/R)/1013;
    
    ions = ones(size(T_eff,1),7);
    %m19 = T*0+1;
    ions(:,2) = ions(:,1).*K01'*p_H2O_mbar;
    ions(:,3) = ions(:,2).*K12'*p_H2O_mbar;
    ions(:,4) = ions(:,3).*K23'*p_H2O_mbar;
    ions(:,5) = ions(:,4).*K34'*p_H2O_mbar;
    ions(:,6) = ions(:,5).*K45'*p_H2O_mbar;
    ions(:,7) = ions(:,6).*K56'*p_H2O_mbar;

    for i=1:size(T_eff,1)
        ions(i,:)=ions(i,:)/sum(ions(i,:));
    end
%     ions=ions

    
    
    figure('color','white','Position',[50 50 600 400]);
    #title('H_3O^+ cluster distribution');
    plot(E_over_N,ions(:,1),'LineWidth',1);
    hold on
    plot(E_over_N,ions(:,2),'LineWidth',1);
    plot(E_over_N,ions(:,3),'LineWidth',1);
    plot(E_over_N,ions(:,4),'LineWidth',1);
    plot(E_over_N,ions(:,5),'LineWidth',1);
    plot(E_over_N,ions(:,6),'LineWidth',1);
    plot(E_over_N,ions(:,7),'LineWidth',1);
set(gca,"fontsize",18)
    xlabel('E/N [Td]',"fontsize",20);
    ylabel('Relative Intensity',"fontsize",20);
    h = legend({'H_3O^+','H_3O^+H_2O','H_3O^+(H_2O)_2','H_3O^+(H_2O)_3','H_3O^+(H_2O)_4','H_3O^+(H_2O)_5','H_3O^+(H_2O)_6'},...
        'location','eastoutside');
    set (h, "fontsize", 18);

##set (h, "fontsize", 20);
    
    
    
##    
##    % plot versus reaction time:
##    t = 0.1./v_ions*1000; % ms
##    
##    figure('color','white','Position',[50 50 800 400]);
##    title('H_3O^+ cluster distribution');
##    plot(t,ions(:,1),'LineWidth',2);
##    hold on
##    plot(t,ions(:,2),'LineWidth',2);
##    plot(t,ions(:,3),'LineWidth',2);
##    plot(t,ions(:,4),'LineWidth',2);
##    plot(t,ions(:,5),'LineWidth',2);
##    plot(t,ions(:,6),'LineWidth',2);
##    plot(t,ions(:,7),'LineWidth',2);
##    
##%     set(gca,'XScale','log');
##    xlabel('Reaction time [ms]');
##    ylabel('Relative Intensity');
##    xlim([0.1 100]);
##    legend({'H_3O^+','H_3O^+H_2O','H_3O^+(H_2O)_2','H_3O^+(H_2O)_3','H_3O^+(H_2O)_4','H_3O^+(H_2O)_5','H_3O^+(H_2O)_6'},...
##        'location','eastoutside');
##    


    ## Iodide - water system:
    ## ------------------------
    
    p_react_mbar = 9;
    p_H2O_mbar_0 = 0.48*9/4;
    p_H2O_mbar = p_H2O_mbar_0;

    T_react_K = 273+20;
    
    % Primary ion distribution - I- water:
    enthalpies = -[43 41 38.9 38.5 37.7 32.6]; % kJ/mol, NIST webbook
    entropies = -[64 73.6 85.4 98.7 100 102]; %J/mol K, Lau et al., 1982
    mobilities = [2.13 2 1.9 1.8 1.75 1.72]; % cm2/Vs, Dotan et al.,1976 and approx. values
    masses = [0 18 36 54 72 90]+127; % u
    E_over_N = 0.1:0.1:120; % Td
    E_120Td = 120*p_react_mbar*100/(1.38*T_react_K); % V/cm at 120 Td
    v_ions_120Td = mobilities*E_120Td*1013/p_react_mbar*T_react_K/273.15/100; %m/s
    v_ions = 3*E_over_N*p_react_mbar*100/(1.38*T_react_K)*1013/p_react_mbar*T_react_K/273.15/100; %m/s
    T_eff = zeros(size(E_over_N,2),size(mobilities,2));
    for i=1:size(mobilities,2)
        %     T_eff(:,i) = 2/(3*kB)*(1/2*30*1.66E-27*(v_ions_120Td(1,i)/120*E_over_N').^2+3/2*kB*T_react);
        T_eff(:,i) = 2/(3*kB)* ( 0.5*(masses(i)+m_air)*m_water*u/(masses(i)+m_water)*(v_ions_120Td(1,i)/120*E_over_N').^2+3/2*kB*T_react_K );
    end
    K01 = exp(-enthalpies(1)./(R/1000*T_eff(:,1)')+entropies(1)/R)/1013; %1/mbar
    K12 = exp(-enthalpies(2)./(R/1000*T_eff(:,2)')+entropies(2)/R)/1013;
    K23 = exp(-enthalpies(3)./(R/1000*T_eff(:,3)')+entropies(3)/R)/1013;
    K34 = exp(-enthalpies(4)./(R/1000*T_eff(:,4)')+entropies(4)/R)/1013;
    K45 = exp(-enthalpies(5)./(R/1000*T_eff(:,5)')+entropies(5)/R)/1013;
    K56 = exp(-enthalpies(6)./(R/1000*T_eff(:,6)')+entropies(6)/R)/1013;
    
    ions = ones(size(T_eff,1),7);
    %m19 = T*0+1;
    ions(:,2) = ions(:,1).*K01'*p_H2O_mbar;
    ions(:,3) = ions(:,2).*K12'*p_H2O_mbar;
    ions(:,4) = ions(:,3).*K23'*p_H2O_mbar;
    ions(:,5) = ions(:,4).*K34'*p_H2O_mbar;
    ions(:,6) = ions(:,5).*K45'*p_H2O_mbar;
    ions(:,7) = ions(:,6).*K56'*p_H2O_mbar;

    for i=1:size(T_eff,1)
        ions(i,:)=ions(i,:)/sum(ions(i,:));
    end

    # --- add 100% RH at 25C:
    p_react_mbar = 9;
    p_H2O_mbar = p_H2O_mbar_0+31.7*p_react_mbar/1013/5;
    T_react_K = 273+20;
    
    % Primary ion distribution - I- water:
    E_120Td = 120*p_react_mbar*100/(1.38*T_react_K); % V/cm at 120 Td
    v_ions_120Td = mobilities*E_120Td*1013/p_react_mbar*T_react_K/273.15/100; %m/s
    v_ions = 3*E_over_N*p_react_mbar*100/(1.38*T_react_K)*1013/p_react_mbar*T_react_K/273.15/100; %m/s
    T_eff = zeros(size(E_over_N,2),size(mobilities,2));
    for i=1:size(mobilities,2)
        %     T_eff(:,i) = 2/(3*kB)*(1/2*30*1.66E-27*(v_ions_120Td(1,i)/120*E_over_N').^2+3/2*kB*T_react);
        T_eff(:,i) = 2/(3*kB)* ( 0.5*(masses(i)+m_air)*m_water*u/(masses(i)+m_water)*(v_ions_120Td(1,i)/120*E_over_N').^2+3/2*kB*T_react_K );
    end
    K01 = exp(-enthalpies(1)./(R/1000*T_eff(:,1)')+entropies(1)/R)/1013; %1/mbar
    K12 = exp(-enthalpies(2)./(R/1000*T_eff(:,2)')+entropies(2)/R)/1013;
    K23 = exp(-enthalpies(3)./(R/1000*T_eff(:,3)')+entropies(3)/R)/1013;
    K34 = exp(-enthalpies(4)./(R/1000*T_eff(:,4)')+entropies(4)/R)/1013;
    K45 = exp(-enthalpies(5)./(R/1000*T_eff(:,5)')+entropies(5)/R)/1013;
    K56 = exp(-enthalpies(6)./(R/1000*T_eff(:,6)')+entropies(6)/R)/1013;
    
    ions2 = ones(size(T_eff,1),7);
    %m19 = T*0+1;
    ions2(:,2) = ions2(:,1).*K01'*p_H2O_mbar;
    ions2(:,3) = ions2(:,2).*K12'*p_H2O_mbar;
    ions2(:,4) = ions2(:,3).*K23'*p_H2O_mbar;
    ions2(:,5) = ions2(:,4).*K34'*p_H2O_mbar;
    ions2(:,6) = ions2(:,5).*K45'*p_H2O_mbar;
    ions2(:,7) = ions2(:,6).*K56'*p_H2O_mbar;

    for i=1:size(T_eff,1)
        ions2(i,:)=ions2(i,:)/sum(ions2(i,:));
    end

    # --- 100% RH at 25C:
    p_react_mbar = 9;
    p_H2O_mbar = 31.7*p_react_mbar/1013/5;
    T_react_K = 273+20;
    
    % Primary ion distribution - I- water:
    E_120Td = 120*p_react_mbar*100/(1.38*T_react_K); % V/cm at 120 Td
    v_ions_120Td = mobilities*E_120Td*1013/p_react_mbar*T_react_K/273.15/100; %m/s
    v_ions = 3*E_over_N*p_react_mbar*100/(1.38*T_react_K)*1013/p_react_mbar*T_react_K/273.15/100; %m/s
    T_eff = zeros(size(E_over_N,2),size(mobilities,2));
    for i=1:size(mobilities,2)
        %     T_eff(:,i) = 2/(3*kB)*(1/2*30*1.66E-27*(v_ions_120Td(1,i)/120*E_over_N').^2+3/2*kB*T_react);
        T_eff(:,i) = 2/(3*kB)* ( 0.5*(masses(i)+m_air)*m_water*u/(masses(i)+m_water)*(v_ions_120Td(1,i)/120*E_over_N').^2+3/2*kB*T_react_K );
    end
    K01 = exp(-enthalpies(1)./(R/1000*T_eff(:,1)')+entropies(1)/R)/1013; %1/mbar
    K12 = exp(-enthalpies(2)./(R/1000*T_eff(:,2)')+entropies(2)/R)/1013;
    K23 = exp(-enthalpies(3)./(R/1000*T_eff(:,3)')+entropies(3)/R)/1013;
    K34 = exp(-enthalpies(4)./(R/1000*T_eff(:,4)')+entropies(4)/R)/1013;
    K45 = exp(-enthalpies(5)./(R/1000*T_eff(:,5)')+entropies(5)/R)/1013;
    K56 = exp(-enthalpies(6)./(R/1000*T_eff(:,6)')+entropies(6)/R)/1013;
    
    ions0 = ones(size(T_eff,1),7);
    %m19 = T*0+1;
    ions0(:,2) = ions0(:,1).*K01'*p_H2O_mbar;
    ions0(:,3) = ions0(:,2).*K12'*p_H2O_mbar;
    ions0(:,4) = ions0(:,3).*K23'*p_H2O_mbar;
    ions0(:,5) = ions0(:,4).*K34'*p_H2O_mbar;
    ions0(:,6) = ions0(:,5).*K45'*p_H2O_mbar;
    ions0(:,7) = ions0(:,6).*K56'*p_H2O_mbar;

    for i=1:size(T_eff,1)
        ions0(i,:)=ions0(i,:)/sum(ions0(i,:));
    end


    # --- 10% RH at 25C:
    p_react_mbar = 9;
    p_H2O_mbar = 31.7/10*p_react_mbar/1013/5;
    T_react_K = 273+20;
    
    % Primary ion distribution - I- water:
    E_120Td = 120*p_react_mbar*100/(1.38*T_react_K); % V/cm at 120 Td
    v_ions_120Td = mobilities*E_120Td*1013/p_react_mbar*T_react_K/273.15/100; %m/s
    v_ions = 3*E_over_N*p_react_mbar*100/(1.38*T_react_K)*1013/p_react_mbar*T_react_K/273.15/100; %m/s
    T_eff = zeros(size(E_over_N,2),size(mobilities,2));
    for i=1:size(mobilities,2)
        %     T_eff(:,i) = 2/(3*kB)*(1/2*30*1.66E-27*(v_ions_120Td(1,i)/120*E_over_N').^2+3/2*kB*T_react);
        T_eff(:,i) = 2/(3*kB)* ( 0.5*(masses(i)+m_air)*m_water*u/(masses(i)+m_water)*(v_ions_120Td(1,i)/120*E_over_N').^2+3/2*kB*T_react_K );
    end
    K01 = exp(-enthalpies(1)./(R/1000*T_eff(:,1)')+entropies(1)/R)/1013; %1/mbar
    K12 = exp(-enthalpies(2)./(R/1000*T_eff(:,2)')+entropies(2)/R)/1013;
    K23 = exp(-enthalpies(3)./(R/1000*T_eff(:,3)')+entropies(3)/R)/1013;
    K34 = exp(-enthalpies(4)./(R/1000*T_eff(:,4)')+entropies(4)/R)/1013;
    K45 = exp(-enthalpies(5)./(R/1000*T_eff(:,5)')+entropies(5)/R)/1013;
    K56 = exp(-enthalpies(6)./(R/1000*T_eff(:,6)')+entropies(6)/R)/1013;
    
    ions3 = ones(size(T_eff,1),7);
    %m19 = T*0+1;
    ions3(:,2) = ions3(:,1).*K01'*p_H2O_mbar;
    ions3(:,3) = ions3(:,2).*K12'*p_H2O_mbar;
    ions3(:,4) = ions3(:,3).*K23'*p_H2O_mbar;
    ions3(:,5) = ions3(:,4).*K34'*p_H2O_mbar;
    ions3(:,6) = ions3(:,5).*K45'*p_H2O_mbar;
    ions3(:,7) = ions3(:,6).*K56'*p_H2O_mbar;

    for i=1:size(T_eff,1)
        ions3(i,:)=ions3(i,:)/sum(ions3(i,:));
    end

    
    
    

    
    
    figure('color','white','Position',[50 50 600 400]);
    #title('H_3O^+ cluster distribution');
    plot(E_over_N,ions(:,1),'LineWidth',1);
    hold on
    plot(E_over_N,ions(:,2),'LineWidth',1);
    plot(E_over_N,ions(:,3),'LineWidth',1);
    plot(E_over_N,ions(:,4),'LineWidth',1);
    #plot(E_over_N,ions(:,5),'LineWidth',1);
    #plot(E_over_N,ions(:,6),'LineWidth',1);
    #plot(E_over_N,ions(:,7),'LineWidth',1);
    set(gca,"fontsize",18)
    xlabel('E/N [Td]',"fontsize",20);
    ylabel('Relative Intensity',"fontsize",20);
    h = legend({'I^-','I^-H_2O','I^-(H_2O)_2','I^-(H_2O)_3','I^-(H_2O)_4','I^-(H_2O)_5','I^-(H_2O)_6'},...
        'location','eastoutside');
    set (h, "fontsize", 18);

    ## import excel-data:
    EN_exp = [2 4 6 8 10 12 14 16 18 20 22];
    m127_exp = [2.23E-01 2.66E-01 3.28E-01 4.21E-01 5.55E-01 7.37E-01 ...
    8.92E-01 9.62E-01 9.79E-01 9.83E-01 9.86E-01];
    m145_exp = [6.60E-01 6.41E-01 6.03E-01 5.32E-01 4.17E-01 2.50E-01 ...
    1.04E-01 3.73E-02 2.01E-02 1.62E-02 1.34E-02];
    m163_exp = [1.16E-01 9.30E-02 6.89E-02 4.74E-02 2.81E-02 1.31E-02 ...
    3.80E-03 1.18E-03 7.91E-04 7.57E-04 7.38E-04];
    #scatter(EN_exp,m127_exp);


    
    blue = [32 51 158]/255;
    green = [32 158 51]/255;
    red = [158, 32, 32]/255;
    lw=4;
    
    # -------------------------------- RH, RH+delta, water added  ----------------------------

    figure('color','white','Position',[50 50 600 400]);
    
    # ions: 0.48 mbar H2O
    plot(E_over_N,ions(:,1),'LineWidth',lw,'color',green);
    hold on
    plot(E_over_N,ions(:,2),'LineWidth',lw,'color',blue);
    plot(E_over_N,ions(:,3),'LineWidth',lw,'color',red);

    # ions2: 0.48+delta mbar H2O
    plot(E_over_N,ions2(:,1),'LineWidth',lw,'linestyle','--','color',green);
    plot(E_over_N,ions2(:,2),'LineWidth',lw,'linestyle','--','color',blue);
    plot(E_over_N,ions2(:,3),'LineWidth',lw,'linestyle','--','color',red);

    
    
    set(gca,"fontsize",24)
    xlabel('E/N [Td]',"fontsize",24);
    ylabel('Relative Intensity',"fontsize",24);
    #h = legend({'I^-','I^-H_2O','I^-(H_2O)_2','I^-(H_2O)_3','I^-(H_2O)_4','I^-(H_2O)_5','I^-(H_2O)_6'},...
     #   'location','eastoutside');
    #set (h, "fontsize", 18);
        set(gca, "linewidth", 2, "fontsize", 24)

    xlim([0,60]);

    #scatter(EN_exp*2.5,m127_exp,'marker','d','markeredgecolor',green,'markerfacecolor','w');

    
    
    # -------------------------------- RH, RH+delta, no water added  ----------------------------

    figure('color','white','Position',[50 50 600 400]);
    
    # ions0: without water addition, 10% RH
    plot(E_over_N,ions3(:,1),'LineWidth',lw,'linestyle','-','color',green);
    hold on
    plot(E_over_N,ions3(:,2),'LineWidth',lw,'linestyle','-','color',blue);
    plot(E_over_N,ions3(:,3),'LineWidth',lw,'linestyle','-','color',red);

    # ions2: 0.48+delta mbar H2O
    plot(E_over_N,ions0(:,1),'LineWidth',lw,'linestyle','--','color',green);
    plot(E_over_N,ions0(:,2),'LineWidth',lw,'linestyle','--','color',blue);
    plot(E_over_N,ions0(:,3),'LineWidth',lw,'linestyle','--','color',red);

    
    
    set(gca,"fontsize",24)
    xlabel('E/N [Td]',"fontsize",24);
    ylabel('Relative Intensity',"fontsize",24);
    #h = legend({'I^-','I^-H_2O','I^-(H_2O)_2','I^-','I^-H_2O','I^-(H_2O)_2'},...
    #    'location','east');
    #set (h, "fontsize", 24);
    set(gca, "linewidth", 2, "fontsize", 24)

    xlim([0,60]);
    
    
    # reaction time:
    EN_react = [0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 ...
    1 2 3 4 5 6 7 8 9 10 15 20 25 30 35 40 45 50 55 60];
    t_react = [50.6593 40.5412 33.7920 28.9692 25.3511 22.5365 20.2844 18.4415 16.9055 15.6058 14.4916 8.4552 5.9688 4.6125 3.7584 3.1712 2.7427 2.4162 2.1592 1.9516 1.3180 0.9950 0.7991 0.6677 0.5734 0.5024 0.4471 0.4027 0.3664 0.3361];

    figure('color','white','Position',[50 50 600 400]);
    
    # ions0: without water addition, 10% RH
    semilogy(EN_react,t_react,'LineWidth',lw,'linestyle','-','color','k');
    set(gca,"fontsize",24)
    xlabel('E/N [Td]',"fontsize",24);
    ylabel('Reaction time [ms]',"fontsize",24);
    set(gca, "linewidth", 2, "fontsize", 24)
    set (gca, "ygrid", "on",'gridcolor','k');
    ax = gca;
ax.GridColor = 'k'
    
    
    
    ## Iodide - Acetonitrile:
    ## ------------------------
    p_react_mbar = 5;
    T_react_K = 273+20;
    p_H2O_mbar = 100/1013*p_react_mbar;
    
    kB = 1.38E-23; % J/K
    eV = 1.602E-19; % eV/J
    u = 1.66E-27; % kg/u
    R = 8.314; % J/mol K
    m_air = 29;
    m_water = 41;
    
    
    % Primary ion distribution - I- acetonitrile:
    enthalpies = -[49.8 43.5 38.5 32.6 29.7 10]; % kJ/mol, NIST webbook
    entropies = -[57.7 77 83.3 80.8 79 77]; %J/mol K, Lau et al., 1982
    mobilities = [2.13 2 1.9 1.8 1.75 1.72]; % cm2/Vs, Dotan et al.,1976 and approx. values
    masses = [0 41 82 123 164 205]+127; % u
    E_over_N = 0.1:0.1:120; % Td
    E_120Td = 120*p_react_mbar*100/(1.38*T_react_K); % V/cm at 120 Td
    v_ions_120Td = mobilities*E_120Td*1013/p_react_mbar*T_react_K/273.15/100; %m/s
    v_ions = 3*E_over_N*p_react_mbar*100/(1.38*T_react_K)*1013/p_react_mbar*T_react_K/273.15/100; %m/s
    T_eff = zeros(size(E_over_N,2),size(mobilities,2));
    for i=1:size(mobilities,2)
        %     T_eff(:,i) = 2/(3*kB)*(1/2*30*1.66E-27*(v_ions_120Td(1,i)/120*E_over_N').^2+3/2*kB*T_react);
        T_eff(:,i) = 2/(3*kB)* ( 0.5*(masses(i)+m_air)*m_water*u/(masses(i)+m_water)*(v_ions_120Td(1,i)/120*E_over_N').^2+3/2*kB*T_react_K );
    end
    K01 = exp(-enthalpies(1)./(R/1000*T_eff(:,1)')+entropies(1)/R)/1013; %1/mbar
    K12 = exp(-enthalpies(2)./(R/1000*T_eff(:,2)')+entropies(2)/R)/1013;
    K23 = exp(-enthalpies(3)./(R/1000*T_eff(:,3)')+entropies(3)/R)/1013;
    K34 = exp(-enthalpies(4)./(R/1000*T_eff(:,4)')+entropies(4)/R)/1013;
    K45 = exp(-enthalpies(5)./(R/1000*T_eff(:,5)')+entropies(5)/R)/1013;
    K56 = exp(-enthalpies(6)./(R/1000*T_eff(:,6)')+entropies(6)/R)/1013;
    
    ions = ones(size(T_eff,1),7);
    %m19 = T*0+1;
    ions(:,2) = ions(:,1).*K01'*p_H2O_mbar;
    ions(:,3) = ions(:,2).*K12'*p_H2O_mbar;
    ions(:,4) = ions(:,3).*K23'*p_H2O_mbar;
    ions(:,5) = ions(:,4).*K34'*p_H2O_mbar;
    ions(:,6) = ions(:,5).*K45'*p_H2O_mbar;
    ions(:,7) = ions(:,6).*K56'*p_H2O_mbar;

    for i=1:size(T_eff,1)
        ions(i,:)=ions(i,:)/sum(ions(i,:));
    end
%     ions=ions

    
    
    figure('color','white','Position',[50 50 600 400]);
    #title('H_3O^+ cluster distribution');
    plot(E_over_N,ions(:,1),'LineWidth',1);
    hold on
    plot(E_over_N,ions(:,2),'LineWidth',1);
    plot(E_over_N,ions(:,3),'LineWidth',1);
    plot(E_over_N,ions(:,4),'LineWidth',1);
    #plot(E_over_N,ions(:,5),'LineWidth',1);
    #plot(E_over_N,ions(:,6),'LineWidth',1);
    #plot(E_over_N,ions(:,7),'LineWidth',1);
    set(gca,"fontsize",18)
    xlabel('E/N [Td]',"fontsize",20);
    ylabel('Relative Intensity',"fontsize",20);
    h = legend({'I^-','I^-C_2H_3N','I^-(C_2H_3N)_2','I^-(C_2H_3N)_3','I^-(C_2H_3N)_4','I^-(C_2H_3N)_5','I^-(C_2H_3N)_6'},...
        'location','eastoutside');
    set (h, "fontsize", 18);


    ## Iodide - SO2:
    ## ------------------------
    p_react_mbar = 5;
    T_react_K = 273;
    p_H2O_mbar = 100/1013*p_react_mbar;
    
    kB = 1.38E-23; % J/K
    eV = 1.602E-19; % eV/J
    u = 1.66E-27; % kg/u
    R = 8.314; % J/mol K
    m_air = 29;
    m_water = 64;
    
    
    % Primary ion distribution - I- acetonitrile:
    enthalpies = -[59.8 42.3 38.5 36 34 32]; % kJ/mol, NIST webbook
    entropies = -[73.6 90.4 103 103 103 103]; %J/mol K, Lau et al., 1982
    mobilities = [2.13 2 1.9 1.8 1.75 1.72]; % cm2/Vs, Dotan et al.,1976 and approx. values
    masses = [0 64 128 192 256 320]+127; % u
    E_over_N = 0.1:0.1:120; % Td
    E_120Td = 120*p_react_mbar*100/(1.38*T_react_K); % V/cm at 120 Td
    v_ions_120Td = mobilities*E_120Td*1013/p_react_mbar*T_react_K/273.15/100; %m/s
    v_ions = 3*E_over_N*p_react_mbar*100/(1.38*T_react_K)*1013/p_react_mbar*T_react_K/273.15/100; %m/s
    T_eff = zeros(size(E_over_N,2),size(mobilities,2));
    for i=1:size(mobilities,2)
        %     T_eff(:,i) = 2/(3*kB)*(1/2*30*1.66E-27*(v_ions_120Td(1,i)/120*E_over_N').^2+3/2*kB*T_react);
        T_eff(:,i) = 2/(3*kB)* ( 0.5*(masses(i)+m_air)*m_water*u/(masses(i)+m_water)*(v_ions_120Td(1,i)/120*E_over_N').^2+3/2*kB*T_react_K );
    end
    K01 = exp(-enthalpies(1)./(R/1000*T_eff(:,1)')+entropies(1)/R)/1013; %1/mbar
    K12 = exp(-enthalpies(2)./(R/1000*T_eff(:,2)')+entropies(2)/R)/1013;
    K23 = exp(-enthalpies(3)./(R/1000*T_eff(:,3)')+entropies(3)/R)/1013;
    K34 = exp(-enthalpies(4)./(R/1000*T_eff(:,4)')+entropies(4)/R)/1013;
    K45 = exp(-enthalpies(5)./(R/1000*T_eff(:,5)')+entropies(5)/R)/1013;
    K56 = exp(-enthalpies(6)./(R/1000*T_eff(:,6)')+entropies(6)/R)/1013;
    
    ions = ones(size(T_eff,1),7);
    %m19 = T*0+1;
    ions(:,2) = ions(:,1).*K01'*p_H2O_mbar;
    ions(:,3) = ions(:,2).*K12'*p_H2O_mbar;
    ions(:,4) = ions(:,3).*K23'*p_H2O_mbar;
    ions(:,5) = ions(:,4).*K34'*p_H2O_mbar;
    ions(:,6) = ions(:,5).*K45'*p_H2O_mbar;
    ions(:,7) = ions(:,6).*K56'*p_H2O_mbar;

    for i=1:size(T_eff,1)
        ions(i,:)=ions(i,:)/sum(ions(i,:));
    end
%     ions=ions

    
    
    figure('color','white','Position',[50 50 600 400]);
    #title('H_3O^+ cluster distribution');
    plot(E_over_N,ions(:,1),'LineWidth',1);
    hold on
    plot(E_over_N,ions(:,2),'LineWidth',1);
    plot(E_over_N,ions(:,3),'LineWidth',1);
    plot(E_over_N,ions(:,4),'LineWidth',1);
    #plot(E_over_N,ions(:,5),'LineWidth',1);
    #plot(E_over_N,ions(:,6),'LineWidth',1);
    #plot(E_over_N,ions(:,7),'LineWidth',1);
    set(gca,"fontsize",18)
    xlabel('E/N [Td]',"fontsize",20);
    ylabel('Relative Intensity',"fontsize",20);
    h = legend({'I^-','I^-SO_2','I^-(SO_2)_2','I^-(SO_2)_3','I^-(SO_2)_4','I^-(SO_2)_5','I^-(SO_2)_6'},...
        'location','eastoutside');
    set (h, "fontsize", 18);

    
    
    
    
