TODO: сделать так, чтобы можно было импортировать



In [1]:

    
%run ../../utils/load_notebook.py
from instabilities import *
import numpy as np









    



importing Jupyter notebook from instabilities.ipynb
Using matplotlib backend: Qt4Agg
Populating the interactive namespace from numpy and matplotlib

Коэффициент для учета вклада гелия в массу газа (см. Notes):



In [1]:

    
He_coeff = 1.34



In [ ]:



In [2]:

    
def flat_end(argument):
    '''декоратор для того, чтобы продолжать функцию на уровне последнего значения'''
    def real_decorator(function):
        def wrapper(*args, **kwargs):
            if args[0] < argument:
                return function(*args, **kwargs)
            else:
                return function(argument, *args[1:], **kwargs)
        return wrapper
    return real_decorator

Для большой оси: $\sigma^2_{maj} = \sigma^2_{\varphi}\sin^2 i + \sigma^2_{z}\cos^2 i$, следовательно примерные ограничения $$\sigma_{maj} < \frac{\sigma_{maj}}{\sqrt{\sin^2 i + 0.49\cos^2 i}}< \sigma_R = \frac{\sigma_{maj}}{\sqrt{f\sin^2 i + \alpha^2\cos^2 i}} ~< \frac{\sigma_{maj}}{\sqrt{0.5\sin^2 i + 0.09\cos^2 i}} < \frac{\sqrt{2}\sigma_{maj}}{\sin i} (или \frac{\sigma_{maj}}{\sqrt{f}\sin i}),$$ или можно более точную оценку дать, если построить $f$ (сейчас $0.5 < f < 1$).

Для малой оси: $\sigma^2_{min} = \sigma^2_{R}\sin^2 i + \sigma^2_{z}\cos^2 i$ и ограничения $$\sigma_{min} < \frac{\sigma_{min}}{\sqrt{\sin^2 i + 0.49\cos^2 i}} < \sigma_R = \frac{\sigma_{min}}{\sqrt{\sin^2 i + \alpha^2\cos^2 i}} ~< \frac{\sigma_{min}}{\sqrt{\sin^2 i + 0.09\cos^2 i}} < \frac{\sigma_{min}}{\sin i}$$

Соответственно имеем 5 оценок из maj и 4 оценки из min.



In [3]:

    
# sig_maj_lim=None
# spl_maj=None

# @flat_end(sig_maj_lim)
# def sig_R_maj_minmin(r, spl_maj=spl_maj):
#     return spl_maj(r).item()

# @flat_end(sig_maj_lim)
# def sig_R_maj_min(r, spl_maj=spl_maj):
#     return spl_maj(r).item()/sqrt(sin_i**2 + 0.49*cos_i**2)
    
# @flat_end(sig_maj_lim)
# def sig_R_maj_max(r, spl_maj=spl_maj):
#     return spl_maj(r).item()/sqrt(0.5*sin_i**2 + 0.09*cos_i**2)

# @flat_end(sig_maj_lim)
# def sig_R_maj_maxmax(r, spl_maj=spl_maj):
#     return spl_maj(r)*sqrt(2)/sin_i
    
# @flat_end(sig_maj_lim)
# def sig_R_maj_maxmaxtrue(r, spl_maj=spl_maj):
#     return spl_maj(r)/sin_i/sqrt(sigPhi_to_sigR_real(r))



In [4]:

    
# sig_min_lim=None
# spl_min=None

# @flat_end(sig_min_lim)
# def sig_R_minor_minmin(r, spl_min=spl_min):
#     return spl_min(r).item()

# @flat_end(sig_min_lim)
# def sig_R_minor_min(r, spl_min=spl_min):
#     return spl_min(r).item()/sqrt(sin_i**2 + 0.49*cos_i**2)
    
# @flat_end(sig_min_lim)
# def sig_R_minor_max(r, spl_min=spl_min):
#     return spl_min(r).item()/sqrt(sin_i**2 + 0.09*cos_i**2)

# @flat_end(sig_min_lim)
# def sig_R_minor_maxmax(r, spl_min=spl_min):
#     return spl_min(r)/sin_i



In [6]:

    
# TODO: move to proper place

def plot_data_lim(ax, data_lim):
    '''Вертикальная линия, обозначающая конец данных'''
    ax.axvline(x=data_lim, ls='-.', color='black', alpha=0.5)
    
def plot_disc_scale(scale, ax, text=None):
    '''Обозначает масштаб диска'''
    ax.plot([scale, scale], [0., 0.05], '-', lw=6., color='black')
    if text:
        ax.annotate(text, xy=(scale, 0.025), xytext=(scale, 0.065), textcoords='data', arrowprops=dict(arrowstyle="->"))
        
def plot_Q_levels(ax, Qs, style='--', color='grey', alpha=0.4):
    '''Функция, чтобы рисовать горизонтальные линии различных уровней $Q^{-1}$:'''
    for Q in Qs:
        ax.axhline(y=1./Q, ls=style, color=color, alpha=alpha)
        
def plot_2f_vs_1f(ax=None, total_gas_data=None, epicycl=None, gas_approx=None, sound_vel=None, scale=None, sigma_max=None, sigma_min=None, star_density_max=None, 
                  star_density_min=None, data_lim=None, color=None, alpha=0.3, disk_scales=[], label=None):
    '''Картинка сравнения 2F и 1F критерия для разных фотометрий и величин sig_R, 
    куда подается весь газ, результат НЕ исправляется за осесимметричные возмущения.'''

    invQg, invQs, invQeff_min = zip(*get_invQeff_from_data(gas_data=total_gas_data, 
                                    epicycl=epicycl, 
                                    gas_approx=gas_approx,
                                    sound_vel=sound_vel, 
                                    scale=scale,
                                    sigma=sigma_max,
                                    star_density=star_density_min))

    invQg, invQs, invQeff_max = zip(*get_invQeff_from_data(gas_data=total_gas_data, 
                                    epicycl=epicycl, 
                                    gas_approx=gas_approx,
                                    sound_vel=sound_vel, 
                                    scale=scale,
                                    sigma=sigma_min,
                                    star_density=star_density_max))

#     invQg = map(lambda l: l*1.6, invQg)
#     invQeff_min = map(lambda l: l*1.6, invQeff_min)
#     invQeff_max = map(lambda l: l*1.6, invQeff_max)
    
    rr = zip(*total_gas_data)[0]
    
    ax.fill_between(rr, invQeff_min, invQeff_max, color=color, alpha=alpha, label=label)
    ax.plot(rr, invQeff_min, 'd-', color=color, alpha=0.6)
    ax.plot(rr, invQeff_max, 'd-', color=color, alpha=0.6)
    ax.plot(rr, invQg, 'v-', color='b')

    ax.set_ylim(0., 1.5)
    ax.set_xlim(0., data_lim+50.)
#     plot_SF(ax)
    plot_data_lim(ax, data_lim)
    for h, annot in disk_scales:
        plot_disc_scale(h, ax, annot)
    plot_Q_levels(ax, [1., 1.5, 2., 3.])
    ax.legend()

Для случая бесконечного тонкого диска: $$\kappa=\frac{3}{R}\frac{d\Phi}{dR}+\frac{d^2\Phi}{dR^2}$$ где $\Phi$ - гравпотенциал, однако его знать не надо, т.к. есть проще формула: $$\kappa=\sqrt{2}\frac{\vartheta_c}{R}\sqrt{1+\frac{R}{\vartheta_c}\frac{d\vartheta_c}{dR}}$$



In [1]:

    
def epicyclicFreq_real(poly_gas, R, resolution):
    '''Честное вычисление эпициклической частоты на расстоянии R для сплайна или полинома'''
    try:
        return sqrt(2.0) * poly_gas(R) * sqrt(1 + R * poly_gas.deriv()(R) / poly_gas(R)) / (R * resolution )
    except:
        return sqrt(2.0) * poly_gas(R) * sqrt(1 + R * poly_gas.derivative()(R) / poly_gas(R)) / (R * resolution )



In [ ]:

Два других механизма из http://iopscience.iop.org/article/10.1088/0004-6256/148/4/69/pdf:

Schaye (2004), 'cold gas phase': $$\Sigma_g > 6.1 f_g^{0.3} Z^{-0.3} I^{0.23}$$ или при constant metallicity of 0.1 $Z_{sun}$ and interstellar flux of ionizing photons 10^6 cm−2 s−1: $$\Sigma_g > 6.1 \frac{\Sigma_g}{\Sigma_g + \Sigma_s}$$



In [59]:

    
def Sigma_crit_S04(gas_dens, r_gas, star_surf_dens):
    return 6.1 * gas_dens / (gas_dens + star_surf_dens(r_gas))

Hunter et al (1998), 'competition with shear' according to Leroy: $$\Sigma_A = \alpha_A\frac{\sigma_g A}{\pi G}$$



In [2]:

    
def oort_a(r, gas_vel):
    try:
        return 0.5 * (gas_vel(r)/r - gas_vel.deriv()(r))
    except:
        return 0.5 * (gas_vel(r)/r - gas_vel.derivative()(r))



In [3]:

    
def Sigma_crit_A(r, gas_vel, alpha, sound_vel):
    G = 4.32
    return alpha * (sound_vel*oort_a(r, gas_vel)) / (np.pi*G)



In [ ]:

Кривая вращения тонкого диска: $$\frac{v^2}{r} = 2\pi G \Sigma_0 \frac{r}{2h} \left[I_0(\frac{r}{2h})K_0(\frac{r}{2h}) - I_1(\frac{r}{2h})K_1(\frac{r}{2h})\right]$$



In [4]:

    
from scipy.special import i0, i1, k0, k1

def disc_vel(r, Sigma0, h, scale, Sigma0_2=None, h_2=None):
    G = 4.3
    bessels = i0(0.5*r/h)*k0(0.5*r/h) - i1(0.5*r/h)*k1(0.5*r/h)
    if h_2 is None:
        return np.sqrt(2*np.pi*G*Sigma0*r*scale * 0.5*r/h * bessels)
    else: #двухдисковая модель
        bessels2 = i0(0.5*r/h_2)*k0(0.5*r/h_2) - i1(0.5*r/h_2)*k1(0.5*r/h_2)
        return np.sqrt(2*np.pi*G*Sigma0*r*scale * 0.5*r/h * bessels + 2*np.pi*G*Sigma0_2*r*scale * 0.5*r/h_2 * bessels2)



In [ ]:

Функция для печати статистики по фотометриям, куда добавлена также информация о полной массе диска $M_d = 2\pi h^2 \Sigma(0)$ (только нужно учесть, что там пк в arcsec надо перевести):



In [3]:

    
from tabulate import tabulate
import pandas as pd

def show_all_photometry_table(all_photometry, scale):
    '''scale in kpc/arcsec'''
    
    copy = [list(l) for l in all_photometry]
    
    #все это дальше из-за того, что бывают двухдисковые модели и их нужно по другому использовать
    for entry in copy:
        if type(entry[5]) == tuple:
            entry[5] = (round(entry[5][0], 2), round(entry[5][1], 2))
        else:
            entry[5] = round(entry[5], 2)
            
    for entry in copy:
        if type(entry[4]) == tuple:
            entry[4] = (round(entry[4][0], 2), round(entry[4][1], 2))
        else:
            entry[4] = round(entry[4], 2)
            
    for entry in copy:
        if type(entry[5]) == tuple:
            entry.append(2*math.pi*entry[5][0]**2 * entry[-1][0](0) * (scale * 1000.)**2 + 
                         2*math.pi*entry[5][1]**2 * entry[-1][1](0) * (scale * 1000.)**2)
        else:
            entry.append(2*math.pi*entry[5]**2 * entry[-1](0) * (scale * 1000.)**2)
            
    for entry in copy:
        if type(entry[5]) == tuple:
            entry.append(entry[7][0](0) + entry[7][1](0))
        else:
            entry.append(entry[7](0))
    
    df = pd.DataFrame(data=copy, columns=['Name', 'r_eff', 'mu_eff', 'n', 'mu0_d', 'h_disc', 'M/L', 'surf', 'M_d/M_sun', 'Sigma_0'])
    df['M/L'] = df['M/L'].apply(lambda l: '%2.2f'%l)

#     df['Sigma_0'] = df['surf'].map(lambda l:l(0))
    df['Sigma_0'] = df['Sigma_0'].apply(lambda l: '%2.0f' % l)
       
#     df['M_d/M_sun'] = 2*math.pi*df['h_disc']**2 * df['surf'].map(lambda l:l(0)) * (scale * 1000.)**2
    df['M_d/M_sun'] = df['M_d/M_sun'].apply(lambda l: '%.2E.' % l)
    
    df.drop('surf', axis=1, inplace=True)
    print tabulate(df, headers='keys', tablefmt='psql', floatfmt=".2f")



In [ ]:

Префикс для плохих фотометрий, чтобы их не брать в конце.



In [1]:

    
BAD_MODEL_PREFIX = 'b:'



In [ ]:

Функция, которая возвращает суммарный профиль плотности для двухдисковой модели, если это необходимо:



In [1]:

    
def tot_dens(dens):
    if type(dens) == tuple:
        star_density = lambda l: dens[0](l) + dens[1](l)
    else:
        star_density = lambda l: dens(l)
    return star_density



In [ ]:

Цикл по стилям линий для того, что бы они отличались там, где их много на одной картинке и цвета сливаются.



In [1]:

    
from itertools import cycle
lines = ["-","--","-.",":"]
linecycler = cycle(lines)



In [ ]:

Раскрашиваем задний фон в "зебру" с некоей периодичностью.



In [1]:

    
def foreground_zebra(ax, step, alpha):
    for i in range(int(ax.get_xlim()[1])+1):
        if i%2 == 0:
            ax.axvspan(i*step, (i+1)*step, color='grey', alpha=alpha)



In [ ]:

Сравним с оценкой Romeo & Falstad (2013) https://ui.adsabs.harvard.edu/#abs/2013MNRAS.433.1389R/abstract:

$$Q_N^{-1} = \sum_{i=1}^{N}\frac{W_i}{Q_iT_i}$$

где

$$Q_i = \frac{\kappa\sigma_{R,i}}{\pi G\Sigma_i}$$$$T_i= \begin{cases} 1 + 0.6(\frac{\sigma_z}{\sigma_R})^2_i, & \mbox{if } 0.0 \le \frac{\sigma_z}{\sigma_R} \le 0.5, \\ 0.8 + 0.7(\frac{\sigma_z}{\sigma_R})_i, & \mbox{if } 0.5 \le \frac{\sigma_z}{\sigma_R} \le 1.0 \end{cases}$$$$W_i = \frac{2\sigma_{R,m}\sigma_{R,i}}{\sigma_{R,m}^2 + \sigma_{R,i}^2},$$$$m:\ index\ of\ min(T_iQ_i)$$

В самой развитой модели 3 компонента - HI, CO, stars. В их же модели полагается $(\sigma_z/\sigma_R)_{CO} = (\sigma_z/\sigma_R)_{HI} = 1$, т.е. $T_{CO} = T_{HI} = 1.5$. Скорость звука я буду полагать в обеих средах равной 11 км/c. Звезды у меня в двух крайних случаях бывают с $\alpha$ равным 0.3 и 0.7, соответственно $T_{0.3} = 1.05;\ T_{0.7}=1.29$.



In [3]:

    
from math import pi

def romeo_Qinv(r=None, epicycl=None, sound_vel=11., sigma_R=None, star_density=None, 
               HI_density=None, CO_density=None, alpha=None, scale=None, gas_approx=None, verbose=False, show=False):
    G = 4.32
    kappa = epicycl(gas_approx, r, scale)
    Q_star = kappa*sigma_R(r)/(pi*G*star_density(r))
    Q_CO = kappa*sound_vel/(pi*G*CO_density)
    Q_HI = kappa*sound_vel/(pi*G*HI_density)
    T_CO, T_HI = 1.5, 1.5
    if alpha > 0 and alpha <= 0.5:
        T_star = 1. + 0.6*alpha**2
    else:
        T_star = 0.8 + 0.7*alpha
        
    #     TODO: оставить только show или verbose
    if show:
        print 'r={:7.3f} Qg={:7.3f} Qs={:7.3f} Qg^-1={:7.3f} Qs^-1={:7.3f}'.format(r, Q_HI, Q_star, 1./Q_HI, 1./Q_star)
    
    dispersions = [sigma_R(r), sound_vel, sound_vel]
    QTs = [Q_star*T_star, Q_HI*T_HI, Q_CO*T_CO]
    components = ['star', 'HI', 'H2']
    index = QTs.index(min(QTs))
    
    if verbose:
        print 'QTs: {}'.format(QTs)
        print 'min index: {}'.format(index)
        print 'min component: {}'.format(components[index])
    
    sig_m = dispersions[index]
    def W_i(sig_m, sig_i):
        return 2*sig_m*sig_i/(sig_m**2 + sig_i**2)
    return W_i(sig_m, dispersions[0])/QTs[0] + W_i(sig_m, dispersions[1])/QTs[1] + W_i(sig_m, dispersions[2])/QTs[2], components[index]

Для тонкого диска:



In [1]:

    
def romeo_Qinv_thin(r=None, epicycl=None, sound_vel=11., sigma_R=None, star_density=None, 
               HI_density=None, CO_density=None, alpha=None, scale=None, gas_approx=None, verbose=False, show=False):
    G = 4.32
    kappa = epicycl(gas_approx, r, scale)
    Q_star = kappa*sigma_R(r)/(pi*G*star_density(r))
    Q_CO = kappa*sound_vel/(pi*G*CO_density)
    Q_HI = kappa*sound_vel/(pi*G*HI_density)
    
#     TODO: оставить только show или verbose
    if show:
        print 'r={:7.3f} Qg={:7.3f} Qs={:7.3f} Qg^-1={:7.3f} Qs^-1={:7.3f}'.format(r, Q_HI, Q_star, 1./Q_HI, 1./Q_star)
        
    dispersions = [sigma_R(r), sound_vel, sound_vel]
    QTs = [Q_star, Q_HI, Q_CO]
    components = ['star', 'HI', 'H2']
    index = QTs.index(min(QTs))
    
    if verbose:
        print 'QTs: {}'.format(QTs)
        print 'min index: {}'.format(index)
        print 'min component: {}'.format(components[index])
    
    sig_m = dispersions[index]
    def W_i(sig_m, sig_i):
        return 2*sig_m*sig_i/(sig_m**2 + sig_i**2)
    return W_i(sig_m, dispersions[0])/QTs[0] + W_i(sig_m, dispersions[1])/QTs[1] + W_i(sig_m, dispersions[2])/QTs[2], components[index]

Функция-сравнение с наблюдениями:



In [ ]:

    
def plot_RF13_vs_2F(r_g_dens=None, HI_gas_dens=None, CO_gas_dens=None, epicycl=None, sound_vel=None, sigma_R_max=None, sigma_R_min=None,  
           star_density=None, alpha_max=None, alpha_min=None, scale=None, gas_approx=None, thin=True, show=False):
    '''Плотности газа передаются не скорр. за гелий.'''
    
    if thin:
        romeo_Q = romeo_Qinv_thin
    else:
        romeo_Q = romeo_Qinv
        
    fig = plt.figure(figsize=[20, 5])
    ax = plt.subplot(131)
    
    totgas = zip(r_g_dens, [He_coeff*(l[0]+l[1]) for l in zip(HI_gas_dens, CO_gas_dens)])[1:]
    
    if show:
        print 'sig_R_max case:'
    romeo_min = []
    for r, g, co in zip(r_g_dens, HI_gas_dens, CO_gas_dens):
        rom, _ = romeo_Q(r=r, epicycl=epicycl, sound_vel=sound_vel, sigma_R=sigma_R_max, 
               star_density=star_density, HI_density=He_coeff*g, CO_density=He_coeff*co, 
                         alpha=alpha_min, scale=scale, gas_approx=gas_approx, show=show)
        romeo_min.append(rom)
        if _ == 'star':
            color = 'g' 
        elif _ == 'HI':
            color = 'b'
        else:
            color = 'm'
        ax.scatter(r, rom, 10, marker='o', color=color)


    invQg, invQs, invQeff_min = zip(*get_invQeff_from_data(gas_data=totgas, 
                                    epicycl=epicycl, 
                                    gas_approx=gas_approx,
                                    sound_vel=sound_vel, 
                                    scale=scale,
                                    sigma=sigma_R_max,
                                    star_density=star_density))

    if show:
        print 'sig_R_min case:'
    
    romeo_max = []
    for r, g, co in zip(r_g_dens, HI_gas_dens, CO_gas_dens):
        rom, _ = romeo_Q(r=r, epicycl=epicycl, sound_vel=sound_vel, sigma_R=sigma_R_min, 
               star_density=star_density, HI_density=He_coeff*g, CO_density=He_coeff*co, 
                         alpha=alpha_max, scale=scale, gas_approx=gas_approx, show=show)
        romeo_max.append(rom)
        if _ == 'star':
            color = 'g' 
        elif _ == 'HI':
            color = 'b'
        else:
            color = 'm'
        ax.scatter(r, rom, 10, marker = 's', color=color)


    invQg, invQs, invQeff_max = zip(*get_invQeff_from_data(gas_data=totgas, 
                                    epicycl=epicycl, 
                                    gas_approx=gas_approx,
                                    sound_vel=sound_vel, 
                                    scale=scale,
                                    sigma=sigma_R_min,
                                    star_density=star_density))

    
    ax.plot(r_g_dens[1:], invQeff_min, '-', alpha=0.5, color='r')
    ax.plot(r_g_dens[1:], invQeff_max, '-', alpha=0.5, color='r')
    
    plot_Q_levels(ax, [1., 1.5, 2., 3.])
    ax.set_xlim(0)
    ax.set_ylim(0)
    ax.legend([matplotlib.lines.Line2D([0], [0], linestyle='none', mfc='g', mec='none', marker='o'), 
                      matplotlib.lines.Line2D([0], [0], linestyle='none', mfc='b', mec='none', marker='o'), 
                      matplotlib.lines.Line2D([0], [0], linestyle='none', mfc='m', mec='none', marker='o')], 
                     ['star', 'HI', 'H2'], numpoints=1, markerscale=1, loc='upper right') #add custom legend
    ax.set_title('RF13: major component')
    
    ax = plt.subplot(132)
    ax.plot(romeo_min[1:], invQeff_min, 'o')
    ax.plot(romeo_max[1:], invQeff_max, 'o', color='m', alpha=0.5)
    ax.set_xlabel('Romeo')
    ax.set_ylabel('2F')
    ax.set_xlim(0., 1.)
    ax.set_ylim(0., 1.)
    ax.plot(ax.get_xlim(), ax.get_ylim(), '--')
    
    ax = plt.subplot(133)
    ax.plot(r_g_dens[1:], [l[1]/l[0] for l in zip(romeo_min[1:], invQeff_min)], 'o-')
    ax.plot(r_g_dens[1:], [l[1]/l[0] for l in zip(romeo_max[1:], invQeff_max)], 'o-', color='m', alpha=0.5)
    ax.set_xlabel('R')
    ax.set_ylabel('[2F]/[Romeo]');



In [ ]:

Функция для исправления центральной поверхностной яркости за наклон (приведение диска к виду "плашмя"). Взята из http://www.astronet.ru/db/msg/1166765/node20.html, (61).



In [2]:

    
def mu_face_on(mu0d, cos_i):
    return mu0d + 2.5*np.log10(1./cos_i)



In [ ]:

Функция для анализа влияния параметров. Берет стандартные параметры, делает из них списки и прогоняет для всех в списке, после чего измеряет среднее и std. Можно варьировать несколько параметров одновременно.



In [1]:

    
def plot_param_depend(ax=None, N=None, data_lim=None, color=None, alpha=0.3, disk_scales=[], label=None, **kwargs):

    params = kwargs.copy()
       
    for p in params.keys():
        if p == 'total_gas_data':
            depth = lambda L: isinstance(L, list) and max(map(depth, L))+1 #depth of nested lists
            if depth(params[p]) == 1:
                params[p] = [params[p]]*N
        elif type(params[p]) is not list:
            params[p] = [params[p]]*N
                
    result = []
    for i in range(N):

        invQg, invQs, invQeff_min = zip(*get_invQeff_from_data(gas_data=params['total_gas_data'][i], 
                                        epicycl=params['epicycl'][i], 
                                        gas_approx=params['gas_approx'][i],
                                        sound_vel=params['sound_vel'][i], 
                                        scale=params['scale'][i],
                                        sigma=params['sigma_max'][i],
                                        star_density=params['star_density_min'][i]))

        invQg, invQs, invQeff_max = zip(*get_invQeff_from_data(gas_data=params['total_gas_data'][i], 
                                        epicycl=params['epicycl'][i], 
                                        gas_approx=params['gas_approx'][i],
                                        sound_vel=params['sound_vel'][i], 
                                        scale=params['scale'][i],
                                        sigma=params['sigma_min'][i],
                                        star_density=params['star_density_max'][i]))


        result.append((invQeff_min, invQeff_max))
    
    rr = zip(*params['total_gas_data'][0])[0]
    qmins = []
    qmaxs = []
    for ind, rrr in enumerate(rr):
        qmin = [result[l][0][ind] for l in range(len(result))]
        qmax = [result[l][1][ind] for l in range(len(result))]
        qmins.append((np.mean(qmin), np.std(qmin)))
        qmaxs.append((np.mean(qmax), np.std(qmax)))
    ax.errorbar(rr, zip(*qmins)[0], fmt='o-', yerr=zip(*qmins)[1], elinewidth=6, alpha=0.3);
    ax.errorbar(rr, zip(*qmaxs)[0], fmt='o-', yerr=zip(*qmaxs)[1])
    ax.axhline(y=1., ls='-', color='grey')
    
    ax.set_ylim(0.)
    ax.set_xlim(0.)
    plot_data_lim(ax, data_lim)
    plot_Q_levels(ax, [1., 1.5, 2., 3.]);



In [ ]:



In [ ]:



In [ ]: