In [1]:
from IPython.display import HTML
from IPython.display import Image
import os
%pylab
%matplotlib inline
%run ../../utils/load_notebook.py
In [2]:
from photometry import *
In [3]:
from instabilities import *
In [4]:
from utils import *
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]:
Ставим классику:
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']
Вычисление $Q$:
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
Области звездообразования:
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()
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']
Трехкомпонентная версия с одинаковой дисперсией равной $c$:
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()
Видно, что отличия существенные - например n2985_model36max.
Двухкомпонентное приближение:
In [24]:
# def romeo_Qinv2(r=None, epicycl=None, sound_vel=None, sigma_R=None, star_density=None,
# HI_density=None, CO_density=None, alpha=None, verbose=False, thin=True, He_corr=False, **kwargs):
# '''Возвращает приближение 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
# else:
# CO_density_ = CO_density
# HI_density_ = HI_density
# gas = CO_density_ + HI_density_
# Q_star = kappa*sigma_R/(pi*G*star_density) # не 3.36
# Q_g = kappa*sound_vel/(pi*G*gas)
# if not thin:
# T_g = 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_g, T_star = 1., 1.
# dispersions = [sigma_R, sound_vel]
# QTs = [Q_star*T_star, Q_g*T_g]
# components = ['star', 'gas']
# 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):
# 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]))
# return W_i(sig_m, dispersions[0])/QTs[0] + W_i(sig_m, dispersions[1])/QTs[1], components[mindex]
In [25]:
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_Qinv2(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=verbose, 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_Qinv2(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=verbose, 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 [26]:
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()
Двухкомпонентная приближает хорошо, ее в данном случае и будем использовать.
In [27]:
def calc_romeo_Q(r_g_dens=None, HI_gas_dens=None, CO_gas_dens=None, epicycl=None, sound_vel=None, sound_vel_CO=6., sound_vel_HI=11., 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, color=None, **kwargs):
totgas = zip(r_g_dens, [He_coeff*(l[0]+l[1]) for l in zip(HI_gas_dens, CO_gas_dens)])
romeo_min = []
for ind, (r, g, co) in enumerate(zip(r_g_dens, HI_gas_dens, CO_gas_dens)):
rom, _ = romeo_Qinv2(r=r, epicycl=epicycl[ind], sound_vel=sound_vel, sound_vel_CO=sound_vel, sound_vel_HI=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, thin=thin, He_corr=True)
romeo_min.append(rom)
romeo_max = []
for ind, (r, g, co) in enumerate(zip(r_g_dens, HI_gas_dens, CO_gas_dens)):
rom, _ = romeo_Qinv2(r=r, epicycl=epicycl[ind], sound_vel=sound_vel, sound_vel_CO=sound_vel, sound_vel_HI=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, thin=thin, He_corr=True)
romeo_max.append(rom)
return romeo_min, romeo_max
for ind, key in enumerate(models.keys()):
print key
model = models[key]
romeo_min, romeo_max = calc_romeo_Q(alpha_max=0.7, alpha_min=0.3, **model)
model['romeo_min'] = romeo_min
model['romeo_max'] = romeo_max
In [29]:
inverse = lambda l: 1./l
fig, (ax, ax2) = plt.subplots(ncols=2, nrows=1, figsize=[20, 8], sharey=True)
ax.plot([0., 7.0], [0., 7.0], '--', alpha=0.9, color='gray', lw=5)
ax2.plot([0., 7.0], [0., 7.0], '--', alpha=0.9, color='gray', lw=5)
WS_ratio = []
RW_ratio = []
for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
for key in models.keys():
if name in key:
print key
model = models[key]
color = model['color']
romeo_min = model['romeo_min']
romeo_max = model['romeo_max']
invQwff_WS_min = model['invQeff_WS_min']
invQwff_WS_max = model['invQeff_WS_max']
invQeff_min = model['invQeff_min']
invQeff_max = model['invQeff_max']
alpha_ = 0.9
ax2.scatter(map(inverse,romeo_min), map(inverse, invQeff_min),
map(lambda l: 4000./l[0], model['total_gas_data']), marker='o', color=color, alpha=alpha_)
ax2.scatter(map(inverse,romeo_max), map(inverse,invQeff_max),
map(lambda l: 4000./l[0], model['total_gas_data']), marker='o', color=color, alpha=alpha_)
ax.scatter(map(inverse,invQwff_WS_min), map(inverse,invQeff_min),
map(lambda l: 4000./l[0], model['total_gas_data']), marker='o', color=color, alpha=alpha_)
ax.scatter(map(inverse,invQwff_WS_max), map(inverse,invQeff_max),
map(lambda l: 4000./l[0], model['total_gas_data']), marker='o', color=color, alpha=alpha_)
RW_ratio.extend([l[0]/l[1] for l in zip(map(inverse,romeo_min), map(inverse, invQeff_min))])
RW_ratio.extend([l[0]/l[1] for l in zip(map(inverse,romeo_max), map(inverse,invQeff_max))])
WS_ratio.extend([l[0]/l[1] for l in zip(map(inverse,invQwff_WS_min), map(inverse,invQeff_min))])
WS_ratio.extend([l[0]/l[1] for l in zip(map(inverse,invQwff_WS_max), map(inverse,invQeff_max))])
if name in ['n338', 'n1167', 'n2985', 'n3898']:
ax.plot([-1, -2], [-1, -2], '-', color=color,
label=r'$\mathrm{NGC\, '+name[1:]+'}\:\: ' + model['band'] + '$', lw=9, alpha=alpha_)
else:
ax2.plot([-1, -2], [-1, -2], '-', color=color,
label=r'$\mathrm{NGC\, '+name[1:]+'}\:\: ' + model['band'] + '$', lw=9, alpha=alpha_)
# ax.legend(loc='lower right')
ax.set_xlim(0, 7.0)
ax.set_ylim(0, 7.0)
ax.set_ylabel(r'$Q_\mathrm{eff}$', fontsize=32)
ax.set_xlabel(r'$Q_\mathrm{WS}$', fontsize=32)
ax.grid()
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(25)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(25)
ax2.legend(loc='upper left', fontsize=19.5)
ax.legend(loc='lower right', fontsize=19.5)
ax2.set_xlim(0, 7.)
ax2.set_ylim(0, 7.)
# ax2.set_ylabel(r'$Q_{eff}^{-1}$', fontsize=30)
ax2.set_xlabel(r'$Q_\mathrm{RW}$', fontsize=32)
ax2.grid()
for tick in ax2.yaxis.get_major_ticks():
tick.label.set_fontsize(25)
for tick in ax2.xaxis.get_major_ticks():
tick.label.set_fontsize(25)
# for _ in np.linspace(1.5, 2., 2):
for _ in [2.]:
# line, = ax2.plot([0., _*2], [0., 2.0], '--', alpha=0.3, color='gray', lw=2.)
# label_line(line, 'y={:2.2f}x'.format(1./_), 0.8*_, 0.8, color='black', size=18)
line2, = ax.plot([0., 7.0], [0., _*7], '--', alpha=0.3, color='gray', lw=2.)
label_line(line2, 'y={:2.1f}x'.format(_), 4./_, 4., color='black', size=25)
# line, = ax2.plot([0., 0.88*2], [0., 2.0], '--', alpha=0.5, color='gray', lw=2.)
# label_line(line, 'y={:2.2f}x'.format(1./0.88), 1.2*0.88, 1.2, 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.05, hspace=0.02)
plt.savefig(paper_imgs_dir+'Qeff_vs_QWS_vs_QRF.eps', format='eps', bbox_inches='tight')
plt.savefig(paper_imgs_dir+'Qeff_vs_QWS_vs_QRF.png', format='png', bbox_inches='tight')
plt.savefig(paper_imgs_dir+'Qeff_vs_QWS_vs_QRF.pdf', format='pdf', dpi=150, bbox_inches='tight')
plt.show()
In [43]:
plt.plot(range(len(WS_ratio)), WS_ratio, 'o')
plt.plot(range(len(RW_ratio)), RW_ratio, 's')
plt.ylim(0.49, 1.1)
plt.axhline(y=1.065, alpha=0.3)
Out[43]:
In [36]:
max(RW_ratio), min(WS_ratio)
Out[36]:
In [38]:
from sklearn.metrics import mean_squared_error
from math import sqrt
rms = sqrt(mean_squared_error([1. for _ in range(len(RW_ratio))], RW_ratio))
print rms
rms = sqrt(mean_squared_error([1. for _ in range(len(WS_ratio))], WS_ratio))
print rms
In [ ]:
Одинаковые маркеры для верхней и нижней оценки (EDIT: теперь выше тоже):
In [29]:
fig, (ax, ax2) = plt.subplots(ncols=2, nrows=1, figsize=[20, 8], sharey=True)
ax.plot([0., 2.0], [0., 2.0], '--', alpha=0.9, color='gray', lw=5)
ax2.plot([0., 2.0], [0., 2.0], '--', alpha=0.9, color='gray', lw=5)
for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
for key in models.keys():
if name in key:
print key
model = models[key]
color = model['color']
romeo_min = model['romeo_min']
romeo_max = model['romeo_max']
invQwff_WS_min = model['invQeff_WS_min']
invQwff_WS_max = model['invQeff_WS_max']
invQeff_min = model['invQeff_min']
invQeff_max = model['invQeff_max']
alpha_ = 0.9
ax2.scatter(romeo_min, invQeff_min, map(lambda l: 4000./l[0], model['total_gas_data']), marker='o', color=color, alpha=alpha_)
ax2.scatter(romeo_max, invQeff_max, map(lambda l: 4000./l[0], model['total_gas_data']), marker='o', color=color, alpha=alpha_)
ax.scatter(invQwff_WS_min, invQeff_min, map(lambda l: 4000./l[0], model['total_gas_data']), marker='o', color=color, alpha=alpha_)
ax.scatter(invQwff_WS_max, invQeff_max, map(lambda l: 4000./l[0], model['total_gas_data']), marker='o', color=color, alpha=alpha_)
ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\mathrm{NGC\, '+name[1:]+'}\:\: ' + model['band'] + '$', lw=9, alpha=alpha_)
ax2.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\mathrm{NGC\, '+name[1:]+'}\:\: ' + model['band'] + '$', lw=9, alpha=alpha_ )
# ax.legend(loc='lower right')
ax.set_xlim(0, 2.0)
ax.set_ylim(0, 1.6)
ax.set_ylabel(r'$Q_{2F}^{-1}$', fontsize=30)
ax.set_xlabel(r'$Q_{WS}^{-1}$', fontsize=30)
ax.grid()
ax2.legend(loc='lower right', fontsize=17)
ax2.set_xlim(0, 2.0)
ax2.set_ylim(0, 1.6)
# ax2.set_ylabel(r'$Q_{eff}^{-1}$', fontsize=30)
ax2.set_xlabel(r'$Q_{RW}^{-1}$', fontsize=30)
ax2.grid()
for _ in np.linspace(1.5, 2., 2):
line, = ax.plot([0., _*2], [0., 2.0], '--', alpha=0.3, color='gray', lw=2.)
label_line(line, 'y={:2.2f}x'.format(1./_), 0.8*_, 0.8, color='black')
line2, = ax2.plot([0., 2.0], [0., _*2], '--', alpha=0.3, color='gray', lw=2.)
label_line(line2, 'y={:2.2f}x'.format(_), 1.2/_, 1.2, color='black')
line, = ax2.plot([0., 0.88*2], [0., 2.0], '--', alpha=0.5, color='gray', lw=2.)
label_line(line, 'y={:2.2f}x'.format(1./0.88), 1.2*0.88, 1.2, 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+'Qeff_vs_QWS_vs_QRF.eps', format='eps', bbox_inches='tight')
# plt.savefig(paper_imgs_dir+'Qeff_vs_QWS_vs_QRF.png', format='png', bbox_inches='tight')
# plt.savefig(paper_imgs_dir+'Qeff_vs_QWS_vs_QRF.pdf', format='pdf', dpi=150, bbox_inches='tight')
plt.show()
Неполные маркеры:
In [30]:
fig, (ax, ax2) = plt.subplots(ncols=2, nrows=1, figsize=[20, 8], sharey=True)
ax.plot([0., 2.0], [0., 2.0], '--', alpha=0.9, color='gray', lw=5)
ax2.plot([0., 2.0], [0., 2.0], '--', alpha=0.9, color='gray', lw=5)
for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
for key in models.keys():
if name in key:
print key
model = models[key]
color = model['color']
romeo_min = model['romeo_min']
romeo_max = model['romeo_max']
invQwff_WS_min = model['invQeff_WS_min']
invQwff_WS_max = model['invQeff_WS_max']
invQeff_min = model['invQeff_min']
invQeff_max = model['invQeff_max']
alpha_ = 0.6
ax2.scatter(romeo_min, invQeff_min, map(lambda l: 4000./l[0], model['total_gas_data']), marker='s', color='w', alpha=alpha_, linewidth=4,
edgecolors=matplotlib.colors.colorConverter.to_rgba(color, alpha=alpha_))
ax2.scatter(romeo_max, invQeff_max, map(lambda l: 4000./l[0], model['total_gas_data']), marker='o', color='w', alpha=alpha_, linewidth=4,
edgecolors=matplotlib.colors.colorConverter.to_rgba(color, alpha=alpha_))
ax.scatter(invQwff_WS_min, invQeff_min, map(lambda l: 4000./l[0], model['total_gas_data']), marker='^', color='w', alpha=alpha_, linewidth=4,
edgecolors=matplotlib.colors.colorConverter.to_rgba(color, alpha=alpha_))
ax.scatter(invQwff_WS_max, invQeff_max, map(lambda l: 4000./l[0], model['total_gas_data']), marker='D', color='w', alpha=alpha_, linewidth=4,
edgecolors=matplotlib.colors.colorConverter.to_rgba(color, alpha=alpha_))
# ax2.loglog(romeo_min, invQeff_min, 's', color=color, alpha=alpha_)
# ax2.loglog(romeo_max, invQeff_max, 'o', color=color, alpha=alpha_)
# ax.loglog(invQwff_WS_min, invQeff_min, '^', color=color, alpha=alpha_)
# ax.loglog(invQwff_WS_max, invQeff_max, 'D', color=color, alpha=alpha_)
ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\mathrm{NGC\, '+name[1:]+'}\:\: ' + model['band'] + '$', lw=9, alpha=alpha_)
ax2.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\mathrm{NGC\, '+name[1:]+'}\:\: ' + model['band'] + '$', lw=9, alpha=alpha_ )
# ax.legend(loc='lower right')
ax.set_xlim(0, 2.0)
ax.set_ylim(0, 1.6)
ax.set_ylabel(r'$Q_{2F}^{-1}$', fontsize=30)
ax.set_xlabel(r'$Q_{WS}^{-1}$', fontsize=30)
ax.grid()
ax2.legend(loc='lower right', fontsize=17)
ax2.set_xlim(0, 2.0)
ax2.set_ylim(0, 1.6)
# ax2.set_ylabel(r'$Q_{eff}^{-1}$', fontsize=30)
ax2.set_xlabel(r'$Q_{RW}^{-1}$', fontsize=30)
ax2.grid()
for _ in np.linspace(1.5, 2., 2):
line, = ax.plot([0., _*2], [0., 2.0], '--', alpha=0.3, color='gray', lw=2.)
label_line(line, 'y={:2.2f}x'.format(1./_), 0.8*_, 0.8, color='black')
line2, = ax2.plot([0., 2.0], [0., _*2], '--', alpha=0.3, color='gray', lw=2.)
label_line(line2, 'y={:2.2f}x'.format(_), 1.2/_, 1.2, color='black')
line, = ax2.plot([0., 0.873*2], [0., 2.0], '--', alpha=0.5, color='gray', lw=2.)
label_line(line, 'y={:2.2f}x'.format(1./0.873), 1.2*0.873, 1.2, 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.show()
Влияние толщины и разных дисперсий - двухкомпонентная с толщиной и трех с дисперсиями, поскольку трех с дисперсиями и толщиной уже представлена в fiducial:
In [31]:
for ind, key in enumerate(models.keys()):
model = models[key]
romeo_min, romeo_max = calc_romeo_Q(alpha_max=0.7, alpha_min=0.3, thin=False, **model)
model['romeo_min_h'] = romeo_min
model['romeo_max_h'] = romeo_max
In [32]:
def calc_romeo_Q_(r_g_dens=None, HI_gas_dens=None, CO_gas_dens=None, epicycl=None, sound_vel=None, sound_vel_CO=6., sound_vel_HI=11., 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, color=None, **kwargs):
totgas = zip(r_g_dens, [He_coeff*(l[0]+l[1]) for l in zip(HI_gas_dens, CO_gas_dens)])
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_CO=sound_vel, sound_vel_HI=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, thin=thin, He_corr=True)
romeo_min.append(rom)
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_CO=sound_vel, sound_vel_HI=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, thin=thin, He_corr=True)
romeo_max.append(rom)
return romeo_min, romeo_max
for ind, key in enumerate(models.keys()):
model = models[key]
model['sound_vel_CO'] = 6.
model['sound_vel_HI'] = 11.
romeo_min, romeo_max = calc_romeo_Q_(alpha_max=0.7, alpha_min=0.3, thin=True, **model)
model['romeo_min_11'] = romeo_min
model['romeo_max_11'] = romeo_max
In [33]:
fig, (ax, ax2) = plt.subplots(ncols=2, nrows=1, figsize=[20, 8], sharey=True)
ax.plot([0., 2.0], [0., 2.0], '--', alpha=0.9, color='gray', lw=5)
ax2.plot([0., 2.0], [0., 2.0], '--', alpha=0.9, color='gray', lw=5)
for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
# for ind, name in enumerate(['n338']):
for key in models.keys():
if name in key:
print key
model = models[key]
color = model['color']
romeo_min_h = model['romeo_min_h']
romeo_max_h = model['romeo_max_h']
romeo_min_11 = model['romeo_min_11']
romeo_max_11 = model['romeo_max_11']
invQeff_min = model['invQeff_min']
invQeff_max = model['invQeff_max']
alpha_ = 0.9
ax.scatter(romeo_min_h, invQeff_min, map(lambda l: 4000./l[0], model['total_gas_data']), marker='s', color=color, alpha=alpha_)
ax.scatter(romeo_max_h, invQeff_max, map(lambda l: 4000./l[0], model['total_gas_data']), marker='o', color=color, alpha=alpha_)
ax2.scatter(romeo_min_11, invQeff_min, map(lambda l: 4000./l[0], model['total_gas_data']), marker='^', color=color, alpha=alpha_)
ax2.scatter(romeo_max_11, invQeff_max, map(lambda l: 4000./l[0], model['total_gas_data']), marker='D', color=color, alpha=alpha_)
ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\mathrm{NGC\, '+name[1:]+'}\:\: ' + model['band'] + '$', lw=9, alpha=alpha_)
ax2.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\mathrm{NGC\, '+name[1:]+'}\:\: ' + model['band'] + '$', lw=9, alpha=alpha_ )
# ax.legend(loc='lower right')
ax.set_xlim(0, 2.0)
ax.set_ylim(0, 1.6)
ax.set_ylabel(r'$Q_{2F}^{-1}$', fontsize=30)
ax.set_xlabel(r'$Q_{RW, thick}^{-1}$', fontsize=30)
ax.grid()
ax2.legend(loc='lower right', fontsize=17)
ax2.set_xlim(0, 2.0)
ax2.set_ylim(0, 1.6)
# ax2.set_ylabel(r'$Q_{eff}^{-1}$', fontsize=30)
ax2.set_xlabel(r'$Q_{RW}^{-1},c_{\rm{HI}}=11\,km\,s^{-1}$', fontsize=30)
ax2.grid()
for _ in np.linspace(1.5, 2., 2):
# line, = ax.plot([0., _*2], [0., 2.0], '--', alpha=0.3, color='gray', lw=2.)
# label_line(line, 'y={:2.2f}x'.format(1./_), 0.8*_, 0.8, color='black')
line2, = ax2.plot([0., 2.0], [0., _*2], '--', alpha=0.3, color='gray', lw=2.)
label_line(line2, 'y={:2.2f}x'.format(_), 1.2/_, 1.2, color='black')
line1, = ax.plot([0., 2.0], [0., _*2], '--', alpha=0.3, color='gray', lw=2.)
label_line(line1, 'y={:2.2f}x'.format(_), 1.2/_, 1.2, color='black')
line, = ax2.plot([0., 0.88*2], [0., 2.0], '--', alpha=0.5, color='gray', lw=2.)
label_line(line, 'y={:2.2f}x'.format(1./0.88), 1.2*0.88, 1.2, color='black')
T_07 = 1.05
line, = ax.plot([0., 1/T_07 * 0.88*2], [0., 2.0], '--', alpha=0.5, color='gray', lw=2.)
label_line(line, 'y={:2.2f}x'.format(T_07/0.88), 1.2*0.88/T_07, 1.2, 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+'Qeff_vs_QRFh_vs_QRF11.eps', format='eps', bbox_inches='tight')
plt.savefig(paper_imgs_dir+'Qeff_vs_QRFh_vs_QRF11.png', format='png', bbox_inches='tight')
plt.savefig(paper_imgs_dir+'Qeff_vs_QRFh_vs_QRF11.pdf', format='pdf', dpi=150, bbox_inches='tight')
plt.show()
Проверим вообще приближенно как у нас распределены параметры по пространству q-s как Romeo любит:
In [36]:
fig = plt.figure(figsize=[9,9])
ax = plt.gca()
for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
for key in models.keys():
if name in key:
# print key
model = models[key]
# plt.plot([model['sound_vel']/l for l in model['sigma_R_min']], [l[0]/l[1] for l in zip(model['Qgs'], model['Qss_min'])], '.', color=model['color'])
plt.plot([model['sound_vel']/l for l in model['sigma_R_max']], [l[0]/l[1] for l in zip(model['Qgs'], model['Qss_max'])], 's', color=model['color'])
# print model.keys()
plt.plot([0.013, 0.21], [10., 0.88], '-')
plt.plot([0.017, 0.21], [0.1, 0.88], '-')
plt.ylim(0.1, 10.)
plt.xlim(0.01, 1.)
plt.xlabel('s')
plt.ylabel('q')
ax.set_xscale("log", nonposx='clip')
ax.set_yscale("log", nonposy='clip');
В качестве реперной модели возьмем RW с коррекцией за толщину, c $\sigma_{CO}=6$ km/s, $\sigma_{HI}=11$ km/s, $\alpha=0.5=(0.7+0.3)/2$ ну и $\sigma_R$ соответствующее принятому $\alpha$.
Наклоны только добавим сперва:
In [37]:
models['n338_modelB']['incl'] = 64.
models['n338_modelR']['incl'] = 64.
models['n1167_modelRsubmax']['incl'] = 38.
models['n1167_modelRmax']['incl'] = 38.
models['n2985_modelKmax']['incl'] = 36.
models['n2985_model36max']['incl'] = 36.
models['n3898_modelRmax']['incl'] = 61.
models['n3898_modelR2dmax']['incl'] = 61.
models['n4258_model36max']['incl'] = 65.
models['n4258_modelImax']['incl'] = 65.
models['n4725_model36max']['incl'] = 50.
models['n4725_modelHmax']['incl'] = 50.
models['n5533_modelRzeroH2']['incl'] = 52.
models['n5533_modelr']['incl'] = 52.
models['n5533_modelRmax']['incl'] = 52.
In [ ]:
In [38]:
def sig_R_min_(incl):
sin_i, cos_i = np.sin(incl*np.pi/180.), np.cos(incl*np.pi/180.)
return 1./sqrt(sin_i**2 + 0.49*cos_i**2)
def sig_R_max_(incl):
sin_i, cos_i = np.sin(incl*np.pi/180.), np.cos(incl*np.pi/180.)
return 1./sqrt(0.5*sin_i**2 + 0.09*cos_i**2)
In [94]:
for i in np.arange(15., 85.):
plt.scatter(i, 1.15*(sig_R_min_(i)+sig_R_max_(i))/2.)
plt.scatter(i, sig_R_max_(i))
plt.scatter(i, sig_R_max_(i)-1.15*(sig_R_min_(i)+sig_R_max_(i))/2., marker='s', color='r')
plt.scatter(i, (1.15*(sig_R_min_(i)+sig_R_max_(i))/2.)/sig_R_max_(i), marker='+', color='g')
for ind, key in enumerate(models.keys()):
plt.scatter(models[key]['incl'], 0., marker='x', color='b')
plt.xlabel(r'$i$')
plt.axvline(x=35.)
plt.axvline(x=65.)
plt.axhline(y=0.95, alpha=0.3, color='m')
plt.ylim(0., 3.)
plt.grid()
In [39]:
factors=[]
for ind, key in enumerate(models.keys()):
model = models[key]
factors.append([sig_R_max_(model['incl']), sig_R_min_(model['incl'])])
plt.plot(range(len(factors)), zip(*factors)[0], 'o-')
plt.plot(range(len(factors)), zip(*factors)[1], 'o-')
plt.plot(range(len(factors)), [(l[0]+l[1])/2. for l in factors], 'o-')
plt.plot(range(len(factors)), [1.15*(l[0]+l[1])/2. for l in factors], 'o-')
Out[39]:
In [97]:
ratios = []
factors_ = []
borders = []
border = 0
for ind2, key in enumerate(models.keys()):
model = models[key]
# model['sound_vel_CO'] = 6.
# model['sound_vel_HI'] = 11.
# print '============='
# print key
# print '============='
romeo_min = []
fiducial = []
mincomp = []
for ind, (r, g, co) in enumerate(zip(model['r_g_dens'], model['HI_gas_dens'], model['CO_gas_dens'])):
# print 'r={:5.3f}'.format(r)
rom, _ = romeo_Qinv(r=r, epicycl=model['epicycl'][ind], sound_vel_CO=6., sound_vel_HI=11.,
sigma_R=(model['sigma_R_max'][ind]+model['sigma_R_min'][ind])/2.,
star_density=model['star_density'][ind], HI_density=He_coeff*g, CO_density=He_coeff*co,
alpha=0.5, thin=False, verbose=False, He_corr=True)
# print 'QRW={:5.3f}'.format(1./rom)
fiducial.append(rom)
mincomp.append(_)
# print 'Qs={:5.3f}//{:5.3f} Qg={:5.3f} Qeff={:5.3f}//{:5.3f}'.format(model['Qss_min'][ind], model['Qss_max'][ind],
# model['Qgs'][ind], 1./model['invQeff_min'][ind],
# 1./model['invQeff_max'][ind])
# print 'QRW < Qg : {}'.format(1./rom < model['Qgs'][ind])
# print 'sR_max={:5.2f} sR_min={:5.2f} sR_RW={:5.2f}, sR_RW*1.15={:5.2f}'.format(model['sigma_R_max'][ind],
# model['sigma_R_min'][ind],
# (model['sigma_R_max'][ind]+model['sigma_R_min'][ind])/2.,
# 1.15*(model['sigma_R_max'][ind]+model['sigma_R_min'][ind])/2.)
# print 'sR_max_factor={:4.2f} sR_min_factor={:4.2f}'.format(sig_R_max_(model['incl']), sig_R_min_(model['incl']))
# if 1./rom < model['Qgs'][ind]:
# ratios.append(rom/model['invQeff_min'][ind])
# factors_.append([sig_R_max_(model['incl']), sig_R_min_(model['incl'])])
ratios.append(rom/model['invQeff_min'][ind])
factors_.append([sig_R_max_(model['incl']), sig_R_min_(model['incl'])])
border+=1
borders.append(border)
model['fiducialQ'] = fiducial
model['mincomp'] = mincomp
In [98]:
fig = plt.figure(figsize=[16, 7])
plt.plot(range(len(ratios)), ratios, 'd', label='Qeff/QRW')
plt.plot(range(len(factors_)), [l[0]/(1.15*(l[0]+l[1])/2.) for l in factors_], 'o-', label='sigRmax/sigR_RW')
plt.plot(range(len(factors_)), [l[0][0]/(1.15*(l[0][0]+l[0][1])/2.)-l[1] for l in zip(factors_, ratios)], 'o', label='diff')
for border in borders:
plt.axvline(x=border, alpha=0.3)
plt.legend(loc='lower right')
plt.xticks(np.array(borders)-2, models.keys(), rotation=45)
plt.axhline(y=1.0, alpha=0.4, color='gray')
# plt.ylim(0, 1.1);
plt.show()
In [101]:
ratios = []
factors_ = []
for ind2, key in enumerate(models.keys()):
# if key == 'n1167_modelRmax':
model = models[key]
plt.plot(model['r_g_dens'], [abs(1./l[0]-1./l[1]) if 1./l[0] < l[2] else 20 for l in zip(model['fiducialQ'],
model['invQeff_min'], model['Qgs'])], '.', color=model['color'])
i = model['incl']
plt.plot(model['r_g_dens'], [(1.-(1.15*(sig_R_min_(i)+sig_R_max_(i))/2.)/sig_R_max_(i))* 1./l[1] if 1./l[0] < l[2] else 20 for l in zip(model['fiducialQ'],
model['invQeff_min'], model['Qgs'])], '+', color=model['color'])
# plt.plot(model['r_g_dens'], [1./l[0]-1./l[1] for l in zip(model['fiducialQ'], model['invQeff_min'], model['Qgs'])], '.', color=model['color'])
plt.ylim(0., 0.6)
plt.xlim(0, 130)
Out[101]:
In [ ]:
In [ ]:
In [75]:
for ind, key in enumerate(models.keys()):
model = models[key]
# model['sound_vel_CO'] = 6.
# model['sound_vel_HI'] = 11.
print '============='
print key
print '============='
romeo_min = []
fiducial = []
mincomp = []
for ind, (r, g, co) in enumerate(zip(model['r_g_dens'], model['HI_gas_dens'], model['CO_gas_dens'])):
print 'r={:5.3f}'.format(r)
rom, _ = romeo_Qinv(r=r, epicycl=model['epicycl'][ind], sound_vel_CO=6., sound_vel_HI=11.,
sigma_R=(model['sigma_R_max'][ind]+model['sigma_R_min'][ind])/2.,
star_density=model['star_density'][ind], HI_density=He_coeff*g, CO_density=He_coeff*co,
alpha=0.5, thin=False, verbose=True, He_corr=True)
print 'QRW={:5.3f}'.format(1./rom)
fiducial.append(rom)
mincomp.append(_)
print 'Qs={:5.3f}//{:5.3f} Qg={:5.3f} Qeff={:5.3f}//{:5.3f}'.format(model['Qss_min'][ind], model['Qss_max'][ind],
model['Qgs'][ind], 1./model['invQeff_min'][ind],
1./model['invQeff_max'][ind])
print 'QRW < Qg : {}'.format(1./rom < model['Qgs'][ind])
print 'sR_max={:5.2f} sR_min={:5.2f} sR_RW={:5.2f}, sR_RW*1.15={:5.2f}'.format(model['sigma_R_max'][ind],
model['sigma_R_min'][ind],
(model['sigma_R_max'][ind]+model['sigma_R_min'][ind])/2.,
1.15*(model['sigma_R_max'][ind]+model['sigma_R_min'][ind])/2.)
print 'sR_max_factor={:4.2f} sR_min_factor={:4.2f}'.format(sig_R_max_(model['incl']), sig_R_min_(model['incl']))
model['fiducialQ'] = fiducial
model['mincomp'] = mincomp
Посмотрим на нее и на грубо уменьшенные по толщине для 338, а также что где доминирует:
In [76]:
model = models['n338_modelR']
plt.plot(model['r_g_dens'], model['invQeff_min'], '-')
plt.plot(model['r_g_dens'], np.array(model['invQeff_min'])/1.15, '--')
plt.plot(model['r_g_dens'], np.array(model['invQeff_min'])/3., '--')
plt.plot(model['r_g_dens'], model['invQeff_max'], '-')
plt.plot(model['r_g_dens'], model['fiducialQ'], '--')
for ind, comp in enumerate(model['mincomp']):
if comp == 'star':
marker = '*'
elif comp == 'HI':
marker = 's'
elif comp == 'H2':
marker = 'd'
plt.scatter(model['r_g_dens'][ind], model['fiducialQ'][ind], 30., marker=marker)
Видно, что почти в три раза уменьшилось - ибо газ делится на два, больше дисперсия у части и еще и толщина.
Проверим также на примере одной точки из работы Romeo что все правильно работает:
In [77]:
# https://arxiv.org/pdf/1602.03049v2.pdf, R=~1 kpc
rom, _ = romeo_Qinv(r=r, epicycl=290., sound_vel_CO=28., sound_vel_HI=11., sigma_R=112.,
star_density=5426., HI_density=0.1, CO_density=730., alpha=0.6, thin=False, verbose=True, He_corr=True)
print(1./rom, _)
Масштабы
In [78]:
models['n338_modelB']['disk_scales'] = 17.7
models['n338_modelR']['disk_scales'] = 18.3
models['n1167_modelRsubmax']['disk_scales'] = -10.
models['n1167_modelRmax']['disk_scales'] = 24.2
models['n2985_modelKmax']['disk_scales'] = 31.1
models['n2985_model36max']['disk_scales'] = [12.8, 48.9]
models['n3898_modelRmax']['disk_scales'] = 36.2
models['n3898_modelR2dmax']['disk_scales'] = [19.11, 59.9]
models['n4258_model36max']['disk_scales'] = 80.7
models['n4258_modelImax']['disk_scales'] = 74.2
models['n4725_model36max']['disk_scales'] = 73.2
models['n4725_modelHmax']['disk_scales'] = 50.28
models['n5533_modelRzeroH2']['disk_scales'] = -10.
models['n5533_modelr']['disk_scales'] = 28.0
models['n5533_modelRmax']['disk_scales'] = 34.4
In [79]:
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, band=None,
star_density_min=None, data_lim=None, color=None, alpha=0.4, disk_scales=[], label=None, name=None, invQeff_min=None,
invQeff_max = None, Qgs = None, hatch=None,**kwargs):
rr = zip(*total_gas_data)[0]
invQg = [1/_ for _ in Qgs]
ax.fill_between(rr, invQeff_min, invQeff_max, color=color, alpha=alpha, label=r'$' + band + '$', hatch=hatch, edgecolor='grey')
ax.plot(rr, invQeff_min, 'd-', color=color, alpha=0.6, lw=2)
ax.plot(rr, invQeff_max, 'd-', color=color, alpha=0.6, lw=2)
ax.plot(rr, invQg, 'o-', color='b', lw=2)
ax.set_ylim(0., 1.3)
ax.set_xlim(0., data_lim+50.)
ax.axvline(x=data_lim, ls=':', color='black', alpha=0.5, lw=4)
if type(disk_scales) == list:
for h in disk_scales:
ax.fill_between([h-0.3, h+0.3], [0., 0.], [0.1, 0.1], color=color, alpha=0.6, hatch=hatch, edgecolor='grey')
else:
ax.fill_between([disk_scales-0.3, disk_scales+0.3], [0.,0.], [0.1, 0.1], color=color, alpha=0.6, hatch=hatch, edgecolor='grey')
if name in ['NGC 5533', 'NGC 4725']:
ax.legend(fontsize=30, loc='lower right')
elif name == 'NGC 5533':
ax.legend(fontsize=30, ncol=2)
else:
ax.legend(fontsize=30)
ax.xaxis.grid(True)
fi, axes = plt.subplots(ncols=1, nrows=7, figsize=[20., 30.], sharex=True)
for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
ax = axes[ind]
count = 0
hatch = None
for key in models.keys():
if name in key:
model = models[key]
if count > 0:
hatch = [None, '.', 'o'][count]
print 'plot {}'.format(key)
plot_2f_vs_1f_(ax, hatch=hatch, **model)
SF[name](ax)
count+=1
ax.text(5, 1.1, r'$\rm{NGC\:'+name[1:]+'}}$', fontsize=40)
if name == 'n5533':
ax.legend(fontsize=30, ncol=2, loc='lower right')
if ind == 6:
ax.set_xlabel(r'$R, arcsec$', fontsize=30)
ax.set_ylabel(r'$Q_{eff}^{-1}$', fontsize=25)
for q_ in [1., 1.5, 2., 3.]:
ax.axhline(y=1./q_, lw=10, alpha=0.15, color='grey')
# ax.set_yticks(['', '0.2', '0.4', '0.6', '0.8', '1.0', ''])
plt.setp(ax.get_yticklabels()[0], visible=False)
plt.setp(ax.get_yticklabels()[-1], visible=False)
plt.tight_layout(pad=0.1, w_pad=0.1, h_pad=0.0)
# fig.subplots_adjust(wspace=0.01, hspace=0.02)
plt.savefig(paper_imgs_dir+'allQeff.eps', format='eps', bbox_inches='tight')
plt.savefig(paper_imgs_dir+'allQeff.png', format='png', bbox_inches='tight')
plt.savefig(paper_imgs_dir+'allQeff.pdf', format='pdf', dpi=150, bbox_inches='tight')
plt.show();