In [1]:

    

from IPython.display import HTML
from IPython.display import Image
import os

%pylab
%matplotlib inline
%run ../../utils/load_notebook.py


    

    




Using matplotlib backend: Qt5Agg
Populating the interactive namespace from numpy and matplotlib

In [2]:

    

from photometry import *


    

    




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

In [3]:

    

from instabilities import *


    

    




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

In [4]:

    

from utils import *


    

    




importing Jupyter notebook from utils.ipynb

In [5]:

    

def label_line(line, label, x, y, color='0.5', size=12):
    """Add a label to a line, at the proper angle.

    Arguments
    ---------
    line : matplotlib.lines.Line2D object,
    label : str
    x : float
        x-position to place center of text (in data coordinated
    y : float
        y-position to place center of text (in data coordinates)
    color : str
    size : float
    """
    xdata, ydata = line.get_data()
    x1 = xdata[0]
    x2 = xdata[-1]
    y1 = ydata[0]
    y2 = ydata[-1]

    ax = line.get_axes()
    text = ax.annotate(label, xy=(x, y), xytext=(-10, 0),
                       textcoords='offset points',
                       size=size, color=color,
                       horizontalalignment='left',
                       verticalalignment='bottom')

    sp1 = ax.transData.transform_point((x1, y1))
    sp2 = ax.transData.transform_point((x2, y2))

    rise = (sp2[1] - sp1[1])
    run = (sp2[0] - sp1[0])

    slope_degrees = np.degrees(np.arctan2(rise, run))
    text.set_rotation(slope_degrees)
    return text


    

In [6]:

    

dictionary = np.load('test_short//models//n1167_modelRmax.npy').tolist()


    

In [7]:

    

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, **kwargs):
    '''Картинка сравнения 2F и 1F критерия для разных фотометрий и величин sig_R, 
    куда подается весь газ, результат НЕ исправляется за осесимметричные возмущения.'''
    
    Qgs = []
    Qss = []
    invQeff_min = []
    for ind, (r, gd) in enumerate(total_gas_data):
        Qgs.append(Qg(epicycl=epicycl[ind], sound_vel=sound_vel, gas_density=gd))
        Qss.append(Qs(epicycl=epicycl[ind], sigma=sigma_max[ind], star_density=star_density_min[ind]))
        qeff = findInvKinemQeffBrentq(Qss[-1], Qgs[-1], sound_vel/sigma_max[ind], np.arange(0.01, 60000., 1.))
        invQeff_min.append(qeff[1])


    Qgs = []
    Qss = []
    invQeff_max = []
    for ind, (r, gd) in enumerate(total_gas_data):
        Qgs.append(Qg(epicycl=epicycl[ind], sound_vel=sound_vel, gas_density=gd))
        Qss.append(Qs(epicycl=epicycl[ind], sigma=sigma_min[ind], star_density=star_density_max[ind]))
        qeff = findInvKinemQeffBrentq(Qss[-1], Qgs[-1], sound_vel/sigma_min[ind], np.arange(0.01, 60000., 1.))
        invQeff_max.append(qeff[1])
           
    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, [1./_ for _ in Qgs], 'v-', color='b')

    ax.set_ylim(0., 1.5)
    ax.set_xlim(0., data_lim+50.)
    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()


    

In [8]:

    

plot_2f_vs_1f_(ax=plt.gca(), **dictionary);


    

In [9]:

    

dictionary = np.load('test_short//models//n1167_modelRsubmax.npy').tolist()


    

In [10]:

    

plot_2f_vs_1f_(ax=plt.gca(), **dictionary);


    

In [11]:

    

models = {}
for model in os.listdir('test_short//models/'):
    if not model[0].startswith('_'):
        models[model[:-4]] = np.load('test_short//models//'+model).tolist()


    

In [12]:

    

models.keys()


    

    
Out[12]:





['n5533_modelRzeroH2',
 'n4725_model36max',
 'n5533_modelr',
 'n1167_modelRsubmax',
 'n338_modelB',
 'n3898_modelRmax',
 'n4725_modelHmax',
 'n3898_modelR2dmax',
 'n2985_modelKmax',
 'n1167_modelRmax',
 'n5533_modelRmax',
 'n2985_model36max',
 'n338_modelR',
 'n4258_model36max',
 'n4258_modelImax']

In [13]:

    

import matplotlib as mpl
mpl.style.use('classic')


    

In [14]:

    

for ind, model in enumerate(models.keys()):
    models[model]['color'] = cm.rainbow(np.linspace(0, 1, 15))[ind]
    
# swap for better look
models['n4258_model36max']['color'], models['n3898_modelRmax']['color'] = models['n3898_modelRmax']['color'], models['n4258_model36max']['color']


    

In [15]:

    

def calc_Qs(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, sfrange=None, **kwargs):
    
    rr = zip(*total_gas_data)[0]

    Qgs = []
    Qss = []
    invQeff_min = []
    for ind, (r, gd) in enumerate(total_gas_data):
#         print 'r={}'.format(r)
        if type(sound_vel) == list:
            sound_vel_ = sound_vel[ind]
        else:
            sound_vel_ = sound_vel
        Qgs.append(Qg(epicycl=epicycl[ind], sound_vel=sound_vel_, gas_density=gd))
        Qss.append(Qs(epicycl=epicycl[ind], sigma=sigma_max[ind], star_density=star_density_min[ind]))
#         print 'gas_density={}'.format(gd)
#         print 'epicycl={}'.format(epicycl[ind])
#         print 'sigma={}'.format(sigma_max[ind])
#         print 'star_density={}'.format(star_density_min[ind])
#         print 'Qgs={}'.format(Qgs[-1])
#         print 'Qss={}'.format(Qss[-1])
        qeff = findInvKinemQeffBrentq(Qss[-1], Qgs[-1], sound_vel_/sigma_max[ind], np.arange(0.01, 60000., 1.))
        invQeff_min.append(qeff[1])
#         print 'invQeff_min={}'.format(invQeff_min[-1])
        
    invQwff_WS_min = [1./Qg_ + 1./Qs_ for Qg_, Qs_ in zip(Qgs, Qss)]

    Qgs = []
    Qss_max = []
    invQeff_max = []
    for ind, (r, gd) in enumerate(total_gas_data):
        if type(sound_vel) == list:
            sound_vel_ = sound_vel[ind]
        else:
            sound_vel_ = sound_vel
        Qgs.append(Qg(epicycl=epicycl[ind], sound_vel=sound_vel_, gas_density=gd))
        Qss_max.append(Qs(epicycl=epicycl[ind], sigma=sigma_min[ind], star_density=star_density_max[ind]))
#         print 'sigma={}'.format(sigma_min[ind])
#         print 'star_density={}'.format(star_density_max[ind])
#         print 'Qss={}'.format(Qss_max[-1])
        qeff = findInvKinemQeffBrentq(Qss_max[-1], Qgs[-1], sound_vel_/sigma_min[ind], np.arange(0.01, 60000., 1.))
        invQeff_max.append(qeff[1])
#         print 'invQeff_max={}'.format(invQeff_max[-1])
        
    invQwff_WS_max = [1./Qg_ + 1./Qs_ for Qg_, Qs_ in zip(Qgs, Qss_max)]
    
    return Qgs, Qss, Qss_max, invQeff_min, invQeff_max, invQwff_WS_min, invQwff_WS_max


for ind, key in enumerate(models.keys()):
    print key
    model = models[key]
    Qgs, Qss, Qss_max, invQeff_min, invQeff_max, invQeff_WS_min, invQeff_WS_max = calc_Qs(**model)
    model['Qgs'] = Qgs
    model['Qss_min'] = Qss
    model['Qss_max'] = Qss_max
    model['invQeff_min'] = invQeff_min
    model['invQeff_max'] = invQeff_max
    model['invQeff_WS_min'] = invQeff_WS_min
    model['invQeff_WS_max'] = invQeff_WS_max


    

    




n5533_modelRzeroH2
n4725_model36max
n5533_modelr
n1167_modelRsubmax
n338_modelB
n3898_modelRmax
n4725_modelHmax
n3898_modelR2dmax
n2985_modelKmax
n1167_modelRmax
n5533_modelRmax
n2985_model36max
n338_modelR
n4258_model36max
n4258_modelImax

In [16]:

    

def plot_SF_338(ax, y=0):
    ax.plot([5., 60.], [y, y], '-', lw=7., color='red')
    
def plot_SF_1167(ax, y=0):
    ax.plot([16., 40.], [y, y], '-', lw=7., color='b') #GALEX
#     ax.plot([55., 80.], [y, y], '-', lw=7., color='b') #GALEX
    ax.plot([49., 80.], [y, y], '-', lw=7., color='b') #GALEX
#     ax.plot([22., 45.], [y, y], '-', lw=7., color='red') #Halpha
    ax.plot([22., 49.], [y, y], '-', lw=7., color='red') #Halpha
    
def plot_SF_2985(ax, y=0):
    ax.plot([10., 70.], [y, y], '-', lw=7., color='b')
    ax.plot([10., 7.2/(0.102*22.4/21.1)], [y, y], '-', lw=7., color='r') #TODO: исправить менее грубо
    
def plot_SF_3898(ax, y=0):
    ax.plot([70., 80.], [y, y], '-', lw=7., color='red')
    ax.plot([190., 210.], [y, y], '-', lw=7., color='red')
    ax.plot([0., 70.], [y, y], '-', lw=7., color='blue')
    ax.plot([40., 50.], [y, y], '-', lw=7., color='blue') #GALEX
    ax.plot([60., 75.], [y, y], '-', lw=7., color='blue') #GALEX
    
def plot_SF_4258(ax, y=0):
    ax.plot([0., 86./2/(175./600.)], [y, y], '-', lw=7., color='b')
    ax.plot([0., 150*58./180.], [y, y], '-', lw=7., color='r')
    ax.plot([250./2/(175./600.), 260./2/(175./600.)], [y, y], '-', lw=7., color='b') #внешняя спираль
    ax.plot([220./2/(175./600.), 310./2/(175./600.)], [y, y], '-', lw=7., color='b') #внешняя спираль, на глаз
    
def plot_SF_4725(ax, y=0):
    ax.plot([0., 300./238 * 161./4.], [y, y], '-', lw=7., color='red') #очень примерно, потому что в SDSS не видно, но оно там есть
    ax.plot([300./238 * 161./4., 300./238 *268./2.], [y, y], '--', lw=6., color='red', alpha=0.5) #в спитцере видно слабое
    ax.plot([300./238 * 161./2., 300./238 *268./2.], [y, y], '-', lw=7., color='red')
    ax.plot([290./238 *220., 300./238 *240], [y, y], '-', lw=7., color='b') #внешняя спираль, ширина условна
    
def plot_SF_5533(ax, y=0):
    ax.plot([10., 95.], [y, y], '-', lw=7., color='red')
    ax.plot([58.*(60/149.), 78.*(60/149.)], [y, y], '-', lw=7., color='b') #спирали оценил из картинки SDSS с масштабом
    ax.plot([126.*(60/149.), 155.*(60/149.)], [y, y], '-', lw=7., color='b')
    ax.plot([204.*(60/149.), 212.*(60/149.)], [y, y], '-', lw=7., color='b')
    ax.plot([235.*(60/149.), 250.*(60/149.)], [y, y], '-', lw=7., color='b')

SF = {'n338' : plot_SF_338, 'n1167' : plot_SF_1167, 'n2985' : plot_SF_2985, 'n3898' : plot_SF_3898, 'n4258': plot_SF_4258, 'n4725' : plot_SF_4725, 'n5533' : plot_SF_5533}


    

In [17]:

    

for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
    for key in models.keys():
        if name in key:
            models[key]['name'] = 'NGC '+name[1:]
            models[key]['band'] = models[key]['label'].split(' ')[0]
            if name == 'n338':
                if key == 'n338_modelR':
                    models[key]['band'] = 'R^{sub}'
                else:
                    models[key]['band'] = 'B^{sub}'
            if name == 'n1167':
                if key == 'n1167_modelRsubmax':
                    models[key]['band'] = 'R^{sub}'
                else:
                    models[key]['band'] = 'R^{max}'
            if key == 'n5533_modelr':
                models[key]['band'] = 'r'
            if key == 'n3898_modelRmax':
                models[key]['band'] = 'R^{N}'
            if key == 'n3898_modelR2dmax':
                models[key]['band'] = 'R^{G}'
            if models[key]['band'] in ['S4G', 'SPITZER']:
                models[key]['band'] = '3.6\mu m'
            if key == 'n5533_modelRzeroH2':
                models[key]['band'] = 'R\,(H_2=0)'


    

In [18]:

    

fig, (ax, ax2) = plt.subplots(ncols=2, nrows=1, figsize=[20, 10], sharey=True)

for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
    for key in models.keys():
        if name in key:
            color = models[key]['color']
            model = models[key]
            
            alpha_ = 0.5
            
            ax.scatter(map(lambda l: 1/l, model['Qgs']), map(lambda l: 1/l, model['Qss_min']), map(lambda l: 2000./l[0], model['total_gas_data']), marker='o', color=color, alpha=alpha_)
            ax.scatter(map(lambda l: 1/l, model['Qgs']), map(lambda l: 1/l, model['Qss_max']), map(lambda l: 2000./l[0], model['total_gas_data']), marker='s', color=color, alpha=alpha_)
            ax2.scatter(model['invQeff_min'], map(lambda l: 1/l, model['Qss_min']), map(lambda l: l[0]/2., model['total_gas_data']), marker='o', color=color, alpha=alpha_)
            ax2.scatter(model['invQeff_max'], map(lambda l: 1/l, model['Qss_max']), map(lambda l: l[0]/2., model['total_gas_data']), marker='s', color=color, alpha=alpha_)
            
#             model['Qgs'] = Qgs
#             model['Qss_min'] = Qss
#             model['Qss_max'] = Qss_max
#             model['invQeff_min'] = invQeff_min
#             model['invQeff_max'] = invQeff_max
#             model['invQwff_WS_min'] = invQwff_WS_min
#             model['invQwff_WS_max'] = invQwff_WS_max
            
            ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: ' + model['band'] + '$', lw=6, alpha=alpha_)
            ax2.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: ' + model['band'] + '$', lw=6, alpha=alpha_)
            
# ax.legend(loc='lower right')
ax.set_xlim(0, 1.6)
ax.set_ylim(0, 1.05)
ax.set_ylabel(r'$Q_{s}^{-1}$', fontsize=30)
ax.set_xlabel(r'$Q_{g}^{-1}$', fontsize=30)
ax.grid()

ax2.legend(loc='upper left', fontsize=16)
ax2.set_xlim(0, 1.6)
# ax2.set_ylim(0, 1.05)
ax2.set_xlabel(r'$Q_{eff}^{-1}$', fontsize=30)
# ax2.set_ylabel(r'$Q_{s}^{-1}$', fontsize=30)
ax2.grid()

# for _ in np.linspace(1., 2., 5):
       
#     line, = ax.plot([0., _*2], [0., 2.0], '--', alpha=0.3, color='gray')
#     label_line(line, 'y={:2.2f}x'.format(1./_), 0.8*_, 0.8, color='black')
    
#     line2, = ax.plot([0., 1.0], [0., _], '--', alpha=0.3, color='gray')
#     label_line(line2, 'y={:2.2f}x'.format(_), 0.8/_, 0.8, color='black')
# #     print np.arctan(1./_)*(180./np.pi)
# #     ax.text(0.8*_, 0.8, 'y={}x'.format(_), rotation=np.arctan(1./_)*(180./np.pi))

ax.plot([0., 2.0], [0., 2.0], '--', alpha=0.9, color='gray', lw=3)
ax2.plot([0., 2.0], [0., 2.0], '--', alpha=0.9, color='gray', lw=3)

line, = ax2.plot([0., 2/0.89], [0., 2.0], '--', alpha=0.5, color='gray', lw=2.)
label_line(line, 'y={:2.2f}x'.format(0.935), 0.9/0.935, 0.9, color='black')

plt.setp(ax2.get_xticklabels()[0], visible=False)

plt.tight_layout(pad=0.1, w_pad=0.1, h_pad=0.1)
fig.subplots_adjust(wspace=0.01, hspace=0.02)
plt.savefig(paper_imgs_dir+'Qs_vs_Qg_vs_Qeff.eps', format='eps', bbox_inches='tight')
plt.savefig(paper_imgs_dir+'Qs_vs_Qg_vs_Qeff.png', format='png', bbox_inches='tight')
plt.savefig(paper_imgs_dir+'Qs_vs_Qg_vs_Qeff.pdf', format='pdf', dpi=150, bbox_inches='tight')

plt.show()


    

    




C:\Anaconda\lib\site-packages\matplotlib\artist.py:233: MatplotlibDeprecationWarning: get_axes has been deprecated in mpl 1.5, please use the
axes property.  A removal date has not been set.
  stacklevel=1)






    

In [19]:

    

fig, (ax, ax2) = plt.subplots(ncols=2, nrows=1, figsize=[20, 10], sharey=True)

for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
    for key in models.keys():
        if name in key:
            color = models[key]['color']
            model = models[key]
            
            alpha_ = 0.5
            
            ax.scatter(map(lambda l: 1/l, model['Qss_min']), map(lambda l: 1/l, model['Qgs']), map(lambda l: 2000./l[0], model['total_gas_data']), marker='o', color=color, alpha=alpha_)
            ax.scatter(map(lambda l: 1/l, model['Qss_max']), map(lambda l: 1/l, model['Qgs']), map(lambda l: 2000./l[0], model['total_gas_data']), marker='s', color=color, alpha=alpha_)
            ax2.scatter(model['invQeff_min'], map(lambda l: 1/l, model['Qgs']), map(lambda l: l[0]/2., model['total_gas_data']), marker='o', color=color, alpha=alpha_)
            ax2.scatter(model['invQeff_max'], map(lambda l: 1/l, model['Qgs']), map(lambda l: l[0]/2., model['total_gas_data']), marker='s', color=color, alpha=alpha_)
            
#             model['Qgs'] = Qgs
#             model['Qss_min'] = Qss
#             model['Qss_max'] = Qss_max
#             model['invQeff_min'] = invQeff_min
#             model['invQeff_max'] = invQeff_max
#             model['invQwff_WS_min'] = invQwff_WS_min
#             model['invQwff_WS_max'] = invQwff_WS_max
            
            ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: ' + model['band'] + '$', lw=6, alpha=alpha_)
            ax2.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: ' + model['band'] + '$', lw=6, alpha=alpha_)
            
# ax.legend(loc='lower right')
ax.set_xlim(0, 1.05)
ax.set_ylim(0, 1.6)
ax.set_ylabel(r'$Q_{g}^{-1}$', fontsize=30)
ax.set_xlabel(r'$Q_{s}^{-1}$', fontsize=30)
ax.grid()

ax2.legend(loc='upper left', fontsize=16)
ax2.set_xlim(0, 1.6)
# ax2.set_ylim(0, 1.05)
ax2.set_xlabel(r'$Q_{eff}^{-1}$', fontsize=30)
# ax2.set_ylabel(r'$Q_{s}^{-1}$', fontsize=30)
ax2.grid()

# for _ in np.linspace(1., 2., 5):
       
#     line, = ax.plot([0., _*2], [0., 2.0], '--', alpha=0.3, color='gray')
#     label_line(line, 'y={:2.2f}x'.format(1./_), 0.8*_, 0.8, color='black')
    
#     line2, = ax.plot([0., 1.0], [0., _], '--', alpha=0.3, color='gray')
#     label_line(line2, 'y={:2.2f}x'.format(_), 0.8/_, 0.8, color='black')
# #     print np.arctan(1./_)*(180./np.pi)
# #     ax.text(0.8*_, 0.8, 'y={}x'.format(_), rotation=np.arctan(1./_)*(180./np.pi))

ax.plot([0., 2.0], [0., 2.0], '--', alpha=0.9, color='gray', lw=3)
ax2.plot([0., 2.0], [0., 2.0], '--', alpha=0.9, color='gray', lw=3)

# line, = ax2.plot([0., 2/0.89], [0., 2.0], '--', alpha=0.5, color='gray', lw=2.)
# label_line(line, 'y={:2.2f}x'.format(0.935), 0.9/0.935, 0.9, color='black')

plt.setp(ax2.get_xticklabels()[0], visible=False)

plt.tight_layout(pad=0.1, w_pad=0.1, h_pad=0.1)
fig.subplots_adjust(wspace=0.01, hspace=0.02)
# plt.savefig(paper_imgs_dir+'Qs_vs_Qg_vs_Qeff.eps', format='eps', bbox_inches='tight')
# plt.savefig(paper_imgs_dir+'Qs_vs_Qg_vs_Qeff.png', format='png', bbox_inches='tight')
# plt.savefig(paper_imgs_dir+'Qs_vs_Qg_vs_Qeff.pdf', format='pdf', dpi=150, bbox_inches='tight')

plt.show()


    

In [20]:

    

for ind, key in enumerate(models.keys()):
    print key
    model = models[key]
    model['r_g_dens'] = zip(*model['total_gas_data'])[0]
    model['HI_gas_dens'] = [l[0][1]/He_coeff - l[1] for l in zip(model['total_gas_data'], model['CO'])]
    model['CO_gas_dens'] = model['CO']
    model['star_density'] = model['star_density_min']
    model['sigma_R_max'] = model['sigma_max']
    model['sigma_R_min'] = model['sigma_min']
    model['sound_vel_CO'] = model['sound_vel']
    model['sound_vel_HI'] = model['sound_vel']


    

    




n5533_modelRzeroH2
n4725_model36max
n5533_modelr
n1167_modelRsubmax
n338_modelB
n3898_modelRmax
n4725_modelHmax
n3898_modelR2dmax
n2985_modelKmax
n1167_modelRmax
n5533_modelRmax
n2985_model36max
n338_modelR
n4258_model36max
n4258_modelImax

In [21]:

    

# from math import pi

# def romeo_Qinv(r=None, epicycl=None, sound_vel=None, sound_vel_CO=6., sound_vel_HI=11., sigma_R=None, star_density=None, 
#                HI_density=None, CO_density=None, alpha=None, verbose=False, thin=True, He_corr=False):
#     '''Возвращает приближение Q из Romeo&Falstad(2013), оно же трехкомпонентное Romeo&Wiegert(2011) и наиболее неустойчивую компоненту.'''
#     G = 4.32
#     kappa = epicycl
    
#     if not He_corr:
#         CO_density = He_coeff*CO_density
#         HI_density = He_coeff*HI_density
        
#     if sound_vel is not None:
#         sound_vel_CO=sound_vel
#         sound_vel_HI=sound_vel
        
#     Q_star = kappa*sigma_R/(pi*G*star_density) # не 3.36
#     Q_CO = kappa*sound_vel_CO/(pi*G*CO_density)
#     Q_HI = kappa*sound_vel_HI/(pi*G*HI_density)
#     if not thin:
#         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
#     else:
#         T_CO, T_HI, T_star = 1., 1., 1.
#     dispersions = [sigma_R, sound_vel_HI, sound_vel_CO]
# #     print dispersions
#     QTs = [Q_star*T_star, Q_HI*T_HI, Q_CO*T_CO]
#     components = ['star', 'HI', 'H2']
#     mindex = QTs.index(min(QTs))
    
#     if verbose:
#         print 'QTs: {}'.format(QTs)
#         print 'min index: {}'.format(mindex)
#         print 'min component: {}'.format(components[mindex])
    
#     sig_m = dispersions[mindex]
#     def W_i(sig_m, sig_i):
# #         print sig_m, sig_i
#         return 2*sig_m*sig_i/(sig_m**2 + sig_i**2)
    
#     if verbose:
#         print 'Ws/TQs={:5.3f} WHI/TQHI={:5.3f} WCO/TQCO={:5.3f}'.format(W_i(sig_m, dispersions[0])/QTs[0], W_i(sig_m, dispersions[1])/QTs[1], W_i(sig_m, dispersions[2])/QTs[2])
#         print 'Ws={:5.3f} WHI={:5.3f} WCO={:5.3f}'.format(W_i(sig_m, dispersions[0]), W_i(sig_m, dispersions[1]), W_i(sig_m, dispersions[2]))
#         print 'THI = {:5.3f} QHI={:5.3f}'.format(T_HI, Q_HI)
        
#     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[mindex]


    

In [22]:

    

def plot_RF13_vs_2F_(r_g_dens=None, HI_gas_dens=None, CO_gas_dens=None, epicycl=None, sound_vel_CO=6., sound_vel_HI=11., sound_vel=None, sigma_R_max=None, sigma_R_min=None,  
           star_density=None, alpha_max=None, alpha_min=None, thin=True, verbose=False, scale=None, gas_approx=None, invQeff_min=None, invQeff_max=None, **kwargs):
    '''Плотности газа передаются не скорр. за гелий.'''
           
    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 sound_vel is not None:
        sound_vel_CO=sound_vel
        sound_vel_HI=sound_vel
    
    if verbose:
        print 'sig_R_max case:'
    romeo_min = []
    for ind, (r, g, co) in enumerate(zip(r_g_dens, HI_gas_dens, CO_gas_dens)):
        rom, _ = romeo_Qinv(r=r, epicycl=epicycl[ind], sound_vel=sound_vel, sigma_R=sigma_R_max[ind], 
               star_density=star_density[ind], HI_density=He_coeff*g, CO_density=He_coeff*co, 
                         alpha=alpha_min, verbose=False, He_corr=True, thin=thin)
        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[ind], 
#                                     gas_approx=gas_approx,
#                                     sound_vel=sound_vel, 
#                                     scale=scale,
#                                     sigma=sigma_R_max,
#                                     star_density=star_density[ind]))

    if verbose:
        print 'sig_R_min case:'
    
    romeo_max = []
    for ind, (r, g, co) in enumerate(zip(r_g_dens, HI_gas_dens, CO_gas_dens)):
        rom, _ = romeo_Qinv(r=r, epicycl=epicycl[ind], sound_vel=sound_vel, sigma_R=sigma_R_min[ind], 
               star_density=star_density[ind], HI_density=He_coeff*g, CO_density=He_coeff*co, 
                         alpha=alpha_max, verbose=False, He_corr=True, thin=thin)
        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, invQeff_min, '-', alpha=0.5, color='r')
    ax.plot(r_g_dens, 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.set_xlim(0., 130.)
    
    ax = plt.subplot(132)
    ax.plot(romeo_min, invQeff_min, 'o')
    ax.plot(romeo_max, 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, [l[1]/l[0] for l in zip(romeo_min, invQeff_min)], 'o-')
    ax.plot(r_g_dens, [l[1]/l[0] for l in zip(romeo_max, invQeff_max)], 'o-', color='m', alpha=0.5)
    ax.set_xlabel('R')
    ax.set_ylabel('[2F]/[Romeo]')
    ax.set_xlim(0., 130.)
    ax.set_ylim(1., 1.5);


    

In [23]:

    

for ind, key in enumerate(models.keys()):
    print key
    model = models[key]   
    plot_RF13_vs_2F_(r_g_dens=model['r_g_dens'], HI_gas_dens=model['HI_gas_dens'], CO_gas_dens=model['CO_gas_dens'], 
                epicycl=model['epicycl'], sound_vel=model['sound_vel'], sigma_R_max=model['sigma_R_max'], sigma_R_min=model['sigma_R_min'],  
                star_density=model['star_density'], 
                alpha_max=0.7, alpha_min=0.3, scale=model['scale'], gas_approx=model['gas_approx'], thin=True, invQeff_min=model['invQeff_min'], invQeff_max=model['invQeff_max'])
#     plot_RF13_vs_2F(alpha_max=0.7, alpha_min=0.3, thin=True, **model)
    plt.show()


    

    




n5533_modelRzeroH2






    













    




n4725_model36max






    













    




n5533_modelr