In [1]:
%run ../../../utils/load_notebook.py
In [2]:
return_back_path = '../../notebooks/2f/test_short/'
In [3]:
import datetime
start = datetime.datetime.now()
start
Out[3]:
In [4]:
%%time
from n338 import *
In [5]:
name
Out[5]:
In [6]:
n338dict = dict(globals(), **locals())
In [7]:
n338dict['name']
Out[7]:
In [8]:
os.chdir(return_back_path)
In [9]:
%%time
from n1167 import *
In [10]:
name
Out[10]:
In [11]:
n1167dict = dict(globals(), **locals())
In [12]:
n1167dict['name']
Out[12]:
In [13]:
del n1167dict['n338dict']
In [14]:
os.chdir(return_back_path)
In [15]:
%%time
from n2985 import *
In [16]:
n2985dict = dict(globals(), **locals())
In [17]:
n2985dict['name']
Out[17]:
In [18]:
del n2985dict['n338dict']
del n2985dict['n1167dict']
In [19]:
os.chdir(return_back_path)
In [20]:
%%time
from n3898 import *
In [21]:
n3898dict = dict(globals(), **locals())
In [22]:
n3898dict['name']
Out[22]:
In [23]:
del n3898dict['n338dict']
del n3898dict['n1167dict']
del n3898dict['n2985dict']
In [24]:
os.chdir(return_back_path)
In [25]:
%%time
from n4258 import *
In [26]:
n4258dict = dict(globals(), **locals())
In [27]:
n4258dict['name']
Out[27]:
In [28]:
del n4258dict['n338dict']
del n4258dict['n1167dict']
del n4258dict['n2985dict']
del n4258dict['n3898dict']
In [29]:
os.chdir(return_back_path)
In [30]:
%%time
from n4725 import *
In [31]:
n4725dict = dict(globals(), **locals())
In [32]:
n4725dict['name']
Out[32]:
In [33]:
del n4725dict['n338dict']
del n4725dict['n1167dict']
del n4725dict['n2985dict']
del n4725dict['n3898dict']
del n4725dict['n4258dict']
In [34]:
os.chdir(return_back_path)
In [35]:
%%time
from n5533 import *
In [36]:
n5533dict = dict(globals(), **locals())
In [37]:
n5533dict['name']
Out[37]:
In [38]:
del n5533dict['n338dict']
del n5533dict['n1167dict']
del n5533dict['n2985dict']
del n5533dict['n3898dict']
del n5533dict['n4258dict']
del n5533dict['n4725dict']
In [39]:
end = datetime.datetime.now()
end
Out[39]:
In [40]:
finish = end-start
In [41]:
finish.total_seconds()/3600.
Out[41]:
In [ ]:
In [42]:
# %%time
# for ind, name1 in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
# for name2 in ['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533'][ind+1:]:
# dict1 = locals()[name1+'dict']
# dict2 = locals()[name2+'dict']
# for key in dict1.keys():
# if key in dict2.keys():
# try:
# if type(dict1[key]) == type(dict2[key]) == np.ndarray:
# if (dict1[key] == dict2[key]).all():
# del dict2[key]
# elif dict1[key] == dict2[key]:
# del dict2[key]
# except Exception:
# print key, name1, name2
In [ ]:
In [43]:
n338dict['star_approx'] == n2985dict['star_approx']
Out[43]:
In [ ]:
In [ ]:
In [ ]:
In [44]:
tex_imgs_dir = 'C:\\Users\\root\\Dropbox\\RotationCurves\\PhD\\paper2\\imgs'
In [45]:
os.chdir(tex_imgs_dir)
fig, axs = plt.subplots(nrows=7, ncols=3, sharex=False, sharey=False, figsize=[16,30])
names = [r'$\rm{NGC\, 338}$', r'$\rm{NGC\, 1167}$', r'$\rm{NGC\, 2985}$', r'$\rm{NGC\, 3898}$', r'$\rm{NGC\, 4258}$', r'$\rm{NGC\, 4725}$', r'$\rm{NGC\, 5533}$',]
ax1 = axs[0,0]
ax2 = axs[1,0]
ax3 = axs[2,0]
ax4 = axs[3,0]
ax5 = axs[4,0]
ax6 = axs[5,0]
ax7 = axs[6,0]
axes = [ax1, ax2, ax3, ax4, ax5, ax6, ax7]
# ===================================================
# Дисперии
# ===================================================
for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
ax = axes[ind]
try:
ax.errorbar(map(abs, locals()[name+'dict']['r_sig_ma']), locals()[name+'dict']['sig_ma'], yerr=locals()[name+'dict']['e_sig_ma'],
fmt='.', marker='.', mew=0, color='red', label='$\sigma_{los}^{maj}$')
ax.plot(points, map(locals()[name+'dict']['sig_R_maj_min'], locals()[name+'dict']['points']))
ax.plot(points, map(locals()[name+'dict']['sig_R_maj_max'], locals()[name+'dict']['points']))
# ax.plot(points, map(locals()[name+'dict']['sig_R_maj_maxmaxtrue'], locals()[name+'dict']['points']))
except Exception:
print 'WARNING S:{}'.format(name)
ax.xaxis.set_major_locator(MultipleLocator(10))
ax.xaxis.set_minor_locator(MultipleLocator(2))
# ax.plot([reb, reb],[0, (50 - 30*(ind/2))], color='black', lw=2)
# ax.axvline(x=reb, ls='--', color='black')
ax.text(0.8, 0.9, names[ind], fontsize=20, ha='center', va='center', transform=ax.transAxes)
ax.yaxis.set_label_coords(-0.055, 0.5)
ax.set_ylabel(r'$\sigma_{\rm{los}},\, \rm{km/s}$', fontsize=20)
ax.set_xlabel(r'$R,\, \rm{arcsec}$', fontsize=20)
ax.yaxis.set_major_locator(MultipleLocator(100))
# ===================================================
# Кривая вращения
# ===================================================
ax1 = axs[0,1]
ax2 = axs[1,1]
ax3 = axs[2,1]
ax4 = axs[3,1]
ax5 = axs[4,1]
ax6 = axs[5,1]
ax7 = axs[6,1]
axes = [ax1, ax2, ax3, ax4, ax5, ax6, ax7]
for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
ax = axes[ind]
try:
ax.plot(points, map(locals()[name+'dict']['spl_gas'], locals()[name+'dict']['points']), '--')
ax.plot(locals()[name+'dict']['r_g_b'], locals()[name+'dict']['vel_g_b'], '.')
except Exception:
print 'WARNING V:{}'.format(name)
try:
ax.plot(locals()[name+'dict']['r_wsrt'], locals()[name+'dict']['vel_wsrt'], '.')
except Exception:
print 'WARNING V:{}'.format(name)
try:
ax.plot(locals()[name+'dict']['r_noord'], locals()[name+'dict']['vel_noord'], '.')
ax.plot(locals()[name+'dict']['r_ma_n'], locals()[name+'dict']['vel_ma_n'], 's')
except Exception:
print 'WARNING V:{}'.format(name)
ax.text(0.8, 0.2, names[ind], fontsize=20, ha='center', va='center', transform=ax.transAxes)
ax.yaxis.set_label_coords(-0.055, 0.5)
ax.set_ylabel(r'$V,\, \rm{km/s}$', fontsize=20)
ax.set_xlabel(r'$R,\, \rm{arcsec}$', fontsize=20)
ax.set_ylim(0, 400)
ax.xaxis.set_major_locator(MultipleLocator(10))
# ax.xaxis.set_minor_locator(MultipleLocator(2))
ax.yaxis.set_major_locator(MultipleLocator(100))
# ax2.set_ylim(0, 400)
# ===================================================
# Поверхностная плотность
# ===================================================
ax1 = axs[0,2]
ax2 = axs[1,2]
ax3 = axs[2,2]
ax4 = axs[3,2]
ax5 = axs[4,2]
ax6 = axs[5,2]
ax7 = axs[6,2]
axes = [ax1, ax2, ax3, ax4, ax5, ax6, ax7]
# axes[3].plot(r_g_dens, gas_dens, 'd-')
# axes[3].plot(r_g_dens, [He_coeff*(y_interp_(l[0], h_disc_R) + l[1]) for l in zip(r_g_dens, gas_dens)], '*-')
# axes[3].plot(r_g_dens, [y_interp_(l, h_disc_R) for l in r_g_dens], '--', label='H2 (R-photom)')
# axes[3].set_title('Gas')
# axes[3].grid()
# axes[3].set_xlim(0, 200)
# axes[3].legend()
# axes[3].plot(r_g_dens, gas_dens, 'd-')
# axes[3].plot(r_g_dens, [He_coeff*(y_interp_(l[0], h_disc_I) + l[1]) for l in zip(r_g_dens, gas_dens)], '*-')
# axes[3].plot(r_g_dens, [y_interp_(l, h_disc_I) for l in r_g_dens], '--', label='H2 (I-photom)')
# axes[3].set_title('Gas')
# axes[3].grid()
# axes[3].set_xlim(0, 200)
# axes[3].legend()
# axes[3].plot(r_HI_dens, HI_dens, '--', label='HI')
# axes[3].plot(zip(*total_gas_data)[0], zip(*total_gas_data)[1], '*-')
# axes[3].plot(r_mol_dens, mol_dens, '--', label='mol')
# axes[3].set_title('Gas')
# axes[3].grid()
# axes[3].set_xlim(0, 200)
# axes[3].legend()
for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
ax = axes[ind]
try:
ax.plot(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'], 'o-')
except Exception:
print 'WARNING G:{}'.format(name)
try:
ax.plot(locals()[name+'dict']['r_HI_dens'], locals()[name+'dict']['HI_dens'], 'o-')
except Exception:
print 'WARNING G:{}'.format(name)
try:
ax.plot(zip(*locals()[name+'dict']['total_gas_data'])[0], zip(*locals()[name+'dict']['total_gas_data'])[1], 'o-')
except Exception:
print 'WARNING G:{}'.format(name)
try:
ax.plot(zip(*locals()[name+'dict']['total_gas_data_'])[0], zip(*locals()[name+'dict']['total_gas_data_'])[1], 'o-')
except Exception:
print 'WARNING G:{}'.format(name)
ax.xaxis.set_major_locator(MultipleLocator(100))
ax.xaxis.set_minor_locator(MultipleLocator(50))
# ax.plot([reb, reb],[0, (50 - 30*(ind/2))], color='black', lw=2)
# ax.axvline(x=reb, ls='--', color='black')
ax.text(0.8, 0.9, names[ind], fontsize=20, ha='center', va='center', transform=ax.transAxes)
ax.yaxis.set_label_coords(-0.055, 0.5)
ax.set_ylabel(r'$\Sigma_{\rm{los}},\, \rm{M_{sun}/{pc}^2}$', fontsize=20)
ax.set_xlabel(r'$R,\, \rm{arcsec}$', fontsize=20)
ax.yaxis.set_major_locator(MultipleLocator(10))
ax.yaxis.set_minor_locator(MultipleLocator(1))
if name == 'n4258':
ax.set_ylim(0, 30)
# ax1.set_xlim(0, 45)
# ax1.set_ylim(0)
# ax2.set_ylim(0, 240)
# ax2.set_xlim(0, 45)
# # ax2.set_ylabel(r'$\sigma_{\rm{los}},\, \rm{km/s}$', fontsize=15)
# ax3.set_ylim(0, 225)
# ax3.set_xlim(0, 60)
# # ax3.axes.yaxis.set_ticklabels([])
# # ax4.set_ylim(0, 150)
# # ax4.set_xlim(0, 72)
# # ax4.set_ylabel(r'$\sigma_{\rm{los}},\, \rm{km/s}$', fontsize=10)
# # ax4.set_xlabel(r'$R,\, \rm{arcsec}$', fontsize=15)
# ax4.set_ylim(0, 110)
# # ax5.set_xlim(0, 50)
# # ax5.axes.yaxis.set_ticklabels([])
# # ax5.set_ylabel(r'$\sigma_{\rm{gas}},\, \rm{km/s}$', fontsize=20)
# # ax5.set_xlabel(r'$R,\, \rm{arcsec}$', fontsize=20)
fig.subplots_adjust(wspace=0.15, hspace=0.25)
# plt.savefig('observ_data.eps', format='eps')
# plt.savefig('observ_data.png', format='png')
# plt.savefig('observ_data.pdf', format='pdf', dpi=150)
plt.show()
In [ ]:
In [46]:
def plot_kennicutt(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, sfrange=None):
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))
rr = zip(*total_gas_data)[0]
# ax.fill_between(rr, invQeff_min, invQeff_max, color=color, alpha=alpha, label=label)
ax.plot(np.array(rr)/sfrange, invQeff_min, '.-', color=color, alpha=0.8)
ax.plot(np.array(rr)/sfrange, invQeff_max, '--', color=color, alpha=0.8)
# 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()
In [47]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}
fig = plt.figure(figsize=[16, 6])
for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
ax = plt.gca()
color = cm.rainbow(np.linspace(0, 1, 7))[ind]
phot = kenn_photom[name]
mud = locals()[name+'dict'][phot[0]]
h = locals()[name+'dict'][phot[1]]
# print name
# print mud
# print h
if name in ['n4258', 'n4725']:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
elif name == 'n2985':
total_gas_data=locals()[name+'dict']['total_gas_data_']
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in
zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
if total_gas_data[0][0] < 0.1:
total_gas_data = total_gas_data[1:]
# print total_gas_data
plot_kennicutt(ax=ax, total_gas_data=total_gas_data,
epicycl=locals()[name+'dict']['epicyclicFreq_real'],
gas_approx=locals()[name+'dict']['spl_gas'],
sound_vel=locals()[name+'dict']['sound_vel'],
scale=locals()[name+'dict']['scale'],
sigma_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_min=locals()[name+'dict']['sig_R_maj_min'],
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
data_lim=locals()[name+'dict']['data_lim'],
color=color, alpha=0.2,
label=name+' '+phot[3],
sfrange=phot[4])
plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
plt.axvline(x=1, alpha=0.5)
plot_Q_levels(ax, [1.5, 2., 3.])
plt.legend()
plt.xlim(0, 3)
plt.ylim(0, 1.3)
plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
plt.xlabel(r'$R/R_{SF}$', fontsize=20)
plt.show()
Среднее значение
In [54]:
def plot_kennicutt_half(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, sfrange=None):
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))
rr = zip(*total_gas_data)[0]
ax.plot(np.array(rr)/sfrange, (np.array(invQeff_min)+np.array(invQeff_max))/2., '.-', color=color, alpha=0.8)
In [55]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}
fig = plt.figure(figsize=[16, 6])
for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
ax = plt.gca()
color = cm.rainbow(np.linspace(0, 1, 7))[ind]
phot = kenn_photom[name]
mud = locals()[name+'dict'][phot[0]]
h = locals()[name+'dict'][phot[1]]
# print name
# print mud
# print h
if name in ['n4258', 'n4725']:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
elif name == 'n2985':
total_gas_data=locals()[name+'dict']['total_gas_data_']
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in
zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
if total_gas_data[0][0] < 0.1:
total_gas_data = total_gas_data[1:]
# print total_gas_data
plot_kennicutt_half(ax=ax, total_gas_data=total_gas_data,
epicycl=locals()[name+'dict']['epicyclicFreq_real'],
gas_approx=locals()[name+'dict']['spl_gas'],
sound_vel=locals()[name+'dict']['sound_vel'],
scale=locals()[name+'dict']['scale'],
sigma_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_min=locals()[name+'dict']['sig_R_maj_min'],
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
data_lim=locals()[name+'dict']['data_lim'],
color=color, alpha=0.2,
label=name+' '+phot[3],
sfrange=phot[4])
plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
plt.axvline(x=1, alpha=0.5)
plot_Q_levels(ax, [1.5, 2., 3.])
plt.legend()
plt.xlim(0, 3)
plt.ylim(0, 1.3)
plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
plt.xlabel(r'$R/R_{SF}$', fontsize=20)
plt.show()
Только для галактик с молек. газом:
In [58]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}
fig = plt.figure(figsize=[16, 6])
for ind, name in enumerate(['n2985', 'n4258', 'n4725']):
ax = plt.gca()
color = cm.rainbow(np.linspace(0, 1, 7))[ind]
phot = kenn_photom[name]
mud = locals()[name+'dict'][phot[0]]
h = locals()[name+'dict'][phot[1]]
# print name
# print mud
# print h
if name in ['n4258', 'n4725']:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
elif name == 'n2985':
total_gas_data=locals()[name+'dict']['total_gas_data_']
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in
zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
if total_gas_data[0][0] < 0.1:
total_gas_data = total_gas_data[1:]
# print total_gas_data
plot_kennicutt_half(ax=ax, total_gas_data=total_gas_data,
epicycl=locals()[name+'dict']['epicyclicFreq_real'],
gas_approx=locals()[name+'dict']['spl_gas'],
sound_vel=locals()[name+'dict']['sound_vel'],
scale=locals()[name+'dict']['scale'],
sigma_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_min=locals()[name+'dict']['sig_R_maj_min'],
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
data_lim=locals()[name+'dict']['data_lim'],
color=color, alpha=0.2,
label=name+' '+phot[3],
sfrange=phot[4])
plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
plt.axvline(x=1, alpha=0.5)
plot_Q_levels(ax, [1.5, 2., 3.])
plt.legend()
plt.xlim(0, 3)
plt.ylim(0, 1.3)
plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
plt.xlabel(r'$R/R_{SF}$', fontsize=20)
plt.show()
ВАУ.
Если изменить масштаб для 4725 (там непонятно), то тоже хорошо, но хуже.
In [59]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 100.],
'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}
fig = plt.figure(figsize=[16, 6])
for ind, name in enumerate(['n2985', 'n4258', 'n4725']):
ax = plt.gca()
color = cm.rainbow(np.linspace(0, 1, 7))[ind]
phot = kenn_photom[name]
mud = locals()[name+'dict'][phot[0]]
h = locals()[name+'dict'][phot[1]]
# print name
# print mud
# print h
if name in ['n4258', 'n4725']:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
elif name == 'n2985':
total_gas_data=locals()[name+'dict']['total_gas_data_']
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in
zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
if total_gas_data[0][0] < 0.1:
total_gas_data = total_gas_data[1:]
# print total_gas_data
plot_kennicutt_half(ax=ax, total_gas_data=total_gas_data,
epicycl=locals()[name+'dict']['epicyclicFreq_real'],
gas_approx=locals()[name+'dict']['spl_gas'],
sound_vel=locals()[name+'dict']['sound_vel'],
scale=locals()[name+'dict']['scale'],
sigma_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_min=locals()[name+'dict']['sig_R_maj_min'],
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
data_lim=locals()[name+'dict']['data_lim'],
color=color, alpha=0.2,
label=name+' '+phot[3],
sfrange=phot[4])
plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
plt.axvline(x=1, alpha=0.5)
plot_Q_levels(ax, [1.5, 2., 3.])
plt.legend()
plt.xlim(0, 3)
plt.ylim(0, 1.3)
plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
plt.xlabel(r'$R/R_{SF}$', fontsize=20)
plt.show()
Разброс для молек.
In [60]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}
fig = plt.figure(figsize=[16, 6])
for ind, name in enumerate(['n2985', 'n4258', 'n4725']):
ax = plt.gca()
color = cm.rainbow(np.linspace(0, 1, 7))[ind]
phot = kenn_photom[name]
mud = locals()[name+'dict'][phot[0]]
h = locals()[name+'dict'][phot[1]]
# print name
# print mud
# print h
if name in ['n4258', 'n4725']:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
elif name == 'n2985':
total_gas_data=locals()[name+'dict']['total_gas_data_']
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in
zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
if total_gas_data[0][0] < 0.1:
total_gas_data = total_gas_data[1:]
# print total_gas_data
plot_kennicutt(ax=ax, total_gas_data=total_gas_data,
epicycl=locals()[name+'dict']['epicyclicFreq_real'],
gas_approx=locals()[name+'dict']['spl_gas'],
sound_vel=locals()[name+'dict']['sound_vel'],
scale=locals()[name+'dict']['scale'],
sigma_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_min=locals()[name+'dict']['sig_R_maj_min'],
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
data_lim=locals()[name+'dict']['data_lim'],
color=color, alpha=0.2,
label=name+' '+phot[3],
sfrange=phot[4])
plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
plt.axvline(x=1, alpha=0.5)
plot_Q_levels(ax, [1.5, 2., 3.])
plt.legend()
plt.xlim(0, 3)
plt.ylim(0, 1.3)
plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
plt.xlabel(r'$R/R_{SF}$', fontsize=20)
plt.show()
Одножидкостный
In [64]:
def plot_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, sfrange=None):
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))
rr = zip(*total_gas_data)[0]
ax.plot(np.array(rr)/sfrange, invQg, '--', color=color, alpha=0.8)
In [65]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}
fig = plt.figure(figsize=[16, 6])
for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
ax = plt.gca()
color = cm.rainbow(np.linspace(0, 1, 7))[ind]
phot = kenn_photom[name]
mud = locals()[name+'dict'][phot[0]]
h = locals()[name+'dict'][phot[1]]
# print name
# print mud
# print h
if name in ['n4258', 'n4725']:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
elif name == 'n2985':
total_gas_data=locals()[name+'dict']['total_gas_data_']
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in
zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
if total_gas_data[0][0] < 0.1:
total_gas_data = total_gas_data[1:]
# print total_gas_data
plot_1f(ax=ax, total_gas_data=total_gas_data,
epicycl=locals()[name+'dict']['epicyclicFreq_real'],
gas_approx=locals()[name+'dict']['spl_gas'],
sound_vel=locals()[name+'dict']['sound_vel'],
scale=locals()[name+'dict']['scale'],
sigma_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_min=locals()[name+'dict']['sig_R_maj_min'],
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
data_lim=locals()[name+'dict']['data_lim'],
color=color, alpha=0.2,
label=name+' '+phot[3],
sfrange=phot[4])
plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
plt.axvline(x=1, alpha=0.5)
plot_Q_levels(ax, [1.5, 2., 3.])
plt.legend()
plt.xlim(0, 3)
plt.ylim(0, 1.3)
plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
plt.xlabel(r'$R/R_{SF}$', fontsize=20)
plt.show()
Только с молек.:
In [66]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}
fig = plt.figure(figsize=[16, 6])
for ind, name in enumerate(['n2985', 'n4258', 'n4725']):
ax = plt.gca()
color = cm.rainbow(np.linspace(0, 1, 7))[ind]
phot = kenn_photom[name]
mud = locals()[name+'dict'][phot[0]]
h = locals()[name+'dict'][phot[1]]
# print name
# print mud
# print h
if name in ['n4258', 'n4725']:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
elif name == 'n2985':
total_gas_data=locals()[name+'dict']['total_gas_data_']
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in
zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
if total_gas_data[0][0] < 0.1:
total_gas_data = total_gas_data[1:]
# print total_gas_data
plot_1f(ax=ax, total_gas_data=total_gas_data,
epicycl=locals()[name+'dict']['epicyclicFreq_real'],
gas_approx=locals()[name+'dict']['spl_gas'],
sound_vel=locals()[name+'dict']['sound_vel'],
scale=locals()[name+'dict']['scale'],
sigma_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_min=locals()[name+'dict']['sig_R_maj_min'],
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
data_lim=locals()[name+'dict']['data_lim'],
color=color, alpha=0.2,
label=name+' '+phot[3],
sfrange=phot[4])
plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
plt.axvline(x=1, alpha=0.5)
plot_Q_levels(ax, [1.5, 2., 3.])
plt.legend()
plt.xlim(0, 3)
plt.ylim(0, 1.3)
plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
plt.xlabel(r'$R/R_{SF}$', fontsize=20)
plt.show()
Vs 2F:
In [67]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}
fig = plt.figure(figsize=[16, 6])
for ind, name in enumerate(['n2985', 'n4258', 'n4725']):
ax = plt.gca()
color = cm.rainbow(np.linspace(0, 1, 7))[ind]
phot = kenn_photom[name]
mud = locals()[name+'dict'][phot[0]]
h = locals()[name+'dict'][phot[1]]
# print name
# print mud
# print h
if name in ['n4258', 'n4725']:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
elif name == 'n2985':
total_gas_data=locals()[name+'dict']['total_gas_data_']
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in
zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
if total_gas_data[0][0] < 0.1:
total_gas_data = total_gas_data[1:]
# print total_gas_data
plot_1f(ax=ax, total_gas_data=total_gas_data,
epicycl=locals()[name+'dict']['epicyclicFreq_real'],
gas_approx=locals()[name+'dict']['spl_gas'],
sound_vel=locals()[name+'dict']['sound_vel'],
scale=locals()[name+'dict']['scale'],
sigma_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_min=locals()[name+'dict']['sig_R_maj_min'],
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
data_lim=locals()[name+'dict']['data_lim'],
color=color, alpha=0.2,
label=name+' '+phot[3],
sfrange=phot[4])
plot_kennicutt_half(ax=ax, total_gas_data=total_gas_data,
epicycl=locals()[name+'dict']['epicyclicFreq_real'],
gas_approx=locals()[name+'dict']['spl_gas'],
sound_vel=locals()[name+'dict']['sound_vel'],
scale=locals()[name+'dict']['scale'],
sigma_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_min=locals()[name+'dict']['sig_R_maj_min'],
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
data_lim=locals()[name+'dict']['data_lim'],
color=color, alpha=0.2,
label=name+' '+phot[3],
sfrange=phot[4])
plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
plt.axvline(x=1, alpha=0.5)
plot_Q_levels(ax, [1.5, 2., 3.])
plt.legend()
plt.xlim(0, 3)
plt.ylim(0, 1.3)
plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
plt.xlabel(r'$R/R_{SF}$', fontsize=20)
plt.show()
1F vs 2F (среднее):
In [68]:
def plot_1f_vs_2f(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, sfrange=None):
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))
rr = zip(*total_gas_data)[0]
ax.plot(invQg, (np.array(invQeff_min)+np.array(invQeff_max))/2., '.', color=color, alpha=0.8)
In [81]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}
fig = plt.figure(figsize=[8, 8])
for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
ax = plt.gca()
color = cm.rainbow(np.linspace(0, 1, 7))[ind]
phot = kenn_photom[name]
mud = locals()[name+'dict'][phot[0]]
h = locals()[name+'dict'][phot[1]]
# print name
# print mud
# print h
if name in ['n4258', 'n4725']:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
elif name == 'n2985':
total_gas_data=locals()[name+'dict']['total_gas_data_']
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in
zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
if total_gas_data[0][0] < 0.1:
total_gas_data = total_gas_data[1:]
# print total_gas_data
plot_1f_vs_2f(ax=ax, total_gas_data=total_gas_data,
epicycl=locals()[name+'dict']['epicyclicFreq_real'],
gas_approx=locals()[name+'dict']['spl_gas'],
sound_vel=locals()[name+'dict']['sound_vel'],
scale=locals()[name+'dict']['scale'],
sigma_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_min=locals()[name+'dict']['sig_R_maj_min'],
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
data_lim=locals()[name+'dict']['data_lim'],
color=color, alpha=0.2,
label=name+' '+phot[3],
sfrange=phot[4])
plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
plt.plot([0., 1.], [0., 1.], '--', alpha=0.1)
plt.legend(loc='lower right')
plt.xlim(0, 1.)
plt.ylim(0, 1.)
plt.ylabel(r'$Q_{eff}^{-1}$', fontsize=20)
plt.xlabel(r'$Q_g^{-1}$', fontsize=20)
plt.show()
Только молек
In [82]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}
fig = plt.figure(figsize=[8, 8])
for ind, name in enumerate(['n2985', 'n4258', 'n4725']):
ax = plt.gca()
color = cm.rainbow(np.linspace(0, 1, 7))[ind]
phot = kenn_photom[name]
mud = locals()[name+'dict'][phot[0]]
h = locals()[name+'dict'][phot[1]]
# print name
# print mud
# print h
if name in ['n4258', 'n4725']:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
elif name == 'n2985':
total_gas_data=locals()[name+'dict']['total_gas_data_']
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in
zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
if total_gas_data[0][0] < 0.1:
total_gas_data = total_gas_data[1:]
# print total_gas_data
plot_1f_vs_2f(ax=ax, total_gas_data=total_gas_data,
epicycl=locals()[name+'dict']['epicyclicFreq_real'],
gas_approx=locals()[name+'dict']['spl_gas'],
sound_vel=locals()[name+'dict']['sound_vel'],
scale=locals()[name+'dict']['scale'],
sigma_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_min=locals()[name+'dict']['sig_R_maj_min'],
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
data_lim=locals()[name+'dict']['data_lim'],
color=color, alpha=0.2,
label=name+' '+phot[3],
sfrange=phot[4])
plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
plt.plot([0., 1.], [0., 1.], '--', alpha=0.1)
plt.legend(loc='lower right')
plt.xlim(0, 1.)
plt.ylim(0, 1.)
plt.ylabel(r'$Q_{eff}^{-1}$', fontsize=20)
plt.xlabel(r'$Q_g^{-1}$', fontsize=20)
plt.show()
Верхнее и нижнее значения:
In [83]:
def plot_1f_vs_2f_full(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, sfrange=None):
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))
rr = zip(*total_gas_data)[0]
ax.plot(invQg, invQeff_min, 'o', color=color, alpha=0.8)
ax.plot(invQg, invQeff_max, '*', color=color, alpha=0.8)
In [84]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}
fig = plt.figure(figsize=[8, 8])
for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
ax = plt.gca()
color = cm.rainbow(np.linspace(0, 1, 7))[ind]
phot = kenn_photom[name]
mud = locals()[name+'dict'][phot[0]]
h = locals()[name+'dict'][phot[1]]
# print name
# print mud
# print h
if name in ['n4258', 'n4725']:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
elif name == 'n2985':
total_gas_data=locals()[name+'dict']['total_gas_data_']
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in
zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
if total_gas_data[0][0] < 0.1:
total_gas_data = total_gas_data[1:]
# print total_gas_data
plot_1f_vs_2f_full(ax=ax, total_gas_data=total_gas_data,
epicycl=locals()[name+'dict']['epicyclicFreq_real'],
gas_approx=locals()[name+'dict']['spl_gas'],
sound_vel=locals()[name+'dict']['sound_vel'],
scale=locals()[name+'dict']['scale'],
sigma_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_min=locals()[name+'dict']['sig_R_maj_min'],
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
data_lim=locals()[name+'dict']['data_lim'],
color=color, alpha=0.2,
label=name+' '+phot[3],
sfrange=phot[4])
plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
plt.plot([0., 1.], [0., 1.], '--', alpha=0.1)
plt.legend(loc='lower right')
plt.xlim(0, 1.)
plt.ylim(0, 1.)
plt.ylabel(r'$Q_{eff}^{-1}$', fontsize=20)
plt.xlabel(r'$Q_g^{-1}$', fontsize=20)
plt.show()
Только молек
In [85]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}
fig = plt.figure(figsize=[8, 8])
for ind, name in enumerate(['n2985', 'n4258', 'n4725']):
ax = plt.gca()
color = cm.rainbow(np.linspace(0, 1, 7))[ind]
phot = kenn_photom[name]
mud = locals()[name+'dict'][phot[0]]
h = locals()[name+'dict'][phot[1]]
# print name
# print mud
# print h
if name in ['n4258', 'n4725']:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
elif name == 'n2985':
total_gas_data=locals()[name+'dict']['total_gas_data_']
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in
zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
if total_gas_data[0][0] < 0.1:
total_gas_data = total_gas_data[1:]
# print total_gas_data
plot_1f_vs_2f_full(ax=ax, total_gas_data=total_gas_data,
epicycl=locals()[name+'dict']['epicyclicFreq_real'],
gas_approx=locals()[name+'dict']['spl_gas'],
sound_vel=locals()[name+'dict']['sound_vel'],
scale=locals()[name+'dict']['scale'],
sigma_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_min=locals()[name+'dict']['sig_R_maj_min'],
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
data_lim=locals()[name+'dict']['data_lim'],
color=color, alpha=0.2,
label=name+' '+phot[3],
sfrange=phot[4])
plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
plt.plot([0., 1.], [0., 1.], '--', alpha=0.1)
plt.legend(loc='lower right')
plt.xlim(0, 1.)
plt.ylim(0, 1.)
plt.ylabel(r'$Q_{eff}^{-1}$', fontsize=20)
plt.xlabel(r'$Q_g^{-1}$', fontsize=20)
plt.show()
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [48]:
print locals()[name+'dict']['total_gas_data_']
print locals()[name+'dict']['epicyclicFreq_real']
print locals()[name+'dict']['spl_gas'](15.)
print locals()[name+'dict']['sound_vel']
print locals()[name+'dict']['scale']
print locals()[name+'dict']['sig_R_maj_max'](15.)
print locals()[name+'dict']['sig_R_maj_min'](15.)
print surf_density(mu=mu_disc(15., mu0=locals()[name+'dict']['mudK'], h=locals()[name+'dict']['hK']), M_to_L=1.42, band='K')
print surf_density(mu=mu_disc(15., mu0=locals()[name+'dict']['mudK'], h=locals()[name+'dict']['hK']), M_to_L=1.42, band='K')
print locals()[name+'dict']['data_lim']
print locals()[name+'dict']['disk_scales']
In [ ]:
In [ ]:
In [49]:
# 338
plot_2f_vs_1f(ax=axes[4], total_gas_data=zip(r_g_dens, [He_coeff*(y_interp_(l[0], h_disc_I) + l[1]) for l in zip(r_g_dens, gas_dens)])[1:],
epicycl=epicyclicFreq_real_,
gas_approx=spl_gas,
sound_vel=sound_vel,
scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: surf_density(mu_disc(l, mu0=mu0d_Ic, h=h_disc_I), 6.23, 'I'),
star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_Ic, h=h_disc_I), 6.23, 'I'),
data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='I maxdisc')
plot_2f_vs_1f(ax=axes[4], total_gas_data=zip(r_g_dens, [He_coeff*(y_interp_(l[0], h_disc_B) + l[1]) for l in zip(r_g_dens, gas_dens)])[1:],
epicycl=epicyclicFreq_real_,
gas_approx=spl_gas,
sound_vel=sound_vel,
scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: surf_density(mu_disc(l, mu0=mu0d_Bc, h=h_disc_B), 6.99, 'B'),
star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_Bc, h=h_disc_B), 6.99, 'B'),
data_lim=data_lim, color='y', alpha=0.2, disk_scales=disk_scales, label='B submaxdisc')
# 1167
plot_2f_vs_1f(ax=axes[4], total_gas_data=zip(r_g_dens, [He_coeff*(y_interp_(l[0], h_disc_R) + l[1]) for l in zip(r_g_dens, gas_dens)])[1:],
epicycl=epicyclicFreq_real,
gas_approx=spl_gas,
sound_vel=sound_vel,
scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), 5.53, 'R'),
star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), 5.53, 'R'),
data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='R maxdisc')
# 2985
plot_2f_vs_1f(ax=ax, total_gas_data=total_gas_data_,
epicycl=epicyclicFreq_real,
gas_approx=spl_gas,
sound_vel=sound_vel,
scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.42, band='K'),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.42, band='K'),
data_lim=data_lim, color='g', alpha=0.2, disk_scales=disk_scales2, label='K Heidt maxdisc')
plot_2f_vs_1f(ax=ax, total_gas_data=total_gas_data_,
epicycl=epicyclicFreq_real,
gas_approx=spl_gas,
sound_vel=sound_vel,
scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu0d_s4g, h=h_disc_s4g), 1.35),
lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu0d_s4g_2, h=h_disc_s4g_2), 1.35)))(l),
star_density_min=lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu0d_s4g, h=h_disc_s4g), 1.35),
lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu0d_s4g_2, h=h_disc_s4g_2), 1.35)))(l),
data_lim=data_lim, color='y', alpha=0.2, disk_scales=disk_scales2, label='S4G 2d maxdisc')
#3898
plot_2f_vs_1f(ax=axes[4], total_gas_data=zip(r_g_dens, [He_coeff*(y_interp_(l[0], h_disc_R) + l[1]) for l in zip(r_g_dens, gas_dens)])[1:7],
epicycl=epicyclicFreq_real,
gas_approx=spl_gas, sound_vel=sound_vel, scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: surf_density(mu_disc(l, mu0=mu0d_R, h=h_disc_R), 7.80, 'R'),
star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_R, h=h_disc_R), 7.80, 'R'),
data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='R Noord maxdisc')
plot_2f_vs_1f(ax=axes[4], total_gas_data=zip(r_g_dens, [He_coeff*(y_interp_(l[0], h_approx) + l[1]) for l in zip(r_g_dens, gas_dens)])[1:7],
epicycl=epicyclicFreq_real,
gas_approx=spl_gas, sound_vel=sound_vel, scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: tot_dens((lambda l1: surf_density(mu=mu_disc(l1, mu0=inner_approx[0], h=h_approx), M_to_L=2.93, band='R'),
lambda l2: surf_density(mu=mu_disc(l2, mu0=mu0_out, h=h_out), M_to_L=2.93, band='R')))(l),
star_density_min=lambda l: tot_dens((lambda l1: surf_density(mu=mu_disc(l1, mu0=inner_approx[0], h=h_approx), M_to_L=2.93, band='R'),
lambda l2: surf_density(mu=mu_disc(l2, mu0=mu0_out, h=h_out), M_to_L=2.93, band='R')))(l),
data_lim=data_lim, color='m', alpha=0.2, disk_scales=disk_scales, label='R Gutierrez 2d maxdisc')
#4258
plot_2f_vs_1f(ax=axes[4], total_gas_data=zip(r_g_dens, map(lambda l: l, gas_dens))[:15], epicycl=epicyclicFreq_real,
gas_approx=spl_gas, sound_vel=sound_vel, scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_I, h=h_disc_I), M_to_L=0.82, band='I'),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_I, h=h_disc_I), M_to_L=0.82, band='I'),
data_lim=data_lim, color='g', alpha=0.2, disk_scales=disk_scales, label='I Yoshino maxdisc')
plot_2f_vs_1f(ax=axes[4], total_gas_data=zip(r_g_dens, map(lambda l: l, gas_dens))[:15], epicycl=epicyclicFreq_real,
gas_approx=spl_gas, sound_vel=sound_vel, scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: s4g_surf_density(mu_disc(l, mu0=mu0d_36, h=h_disc_36), M_to_L=0.38),
star_density_min=lambda l: s4g_surf_density(mu_disc(l, mu0=mu0d_36, h=h_disc_36), M_to_L=0.38),
data_lim=data_lim, color='m', alpha=0.2, disk_scales=disk_scales, label='SPITZER 3.6 maxdisc')
# 4725
plot_2f_vs_1f(ax=axes[4], total_gas_data=zip(r_g_dens, map(lambda l: l, gas_dens))[:15], epicycl=epicyclicFreq_real,
gas_approx=spl_gas, sound_vel=sound_vel, scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_H, h=h_disc_H), M_to_L=0.68, band='H'),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_H, h=h_disc_H), M_to_L=0.68, band='H'),
data_lim=data_lim, color='y', alpha=0.3, disk_scales=disk_scales, label='H Heidt maxdisc')
plot_2f_vs_1f(ax=axes[4], total_gas_data=zip(r_g_dens, map(lambda l: l, gas_dens))[:15], epicycl=epicyclicFreq_real,
gas_approx=spl_gas, sound_vel=sound_vel, scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: s4g_surf_density(mu_disc(l, mu0=mu0d_s4g, h=h_disc_s4g), 1.61),
star_density_min=lambda l: s4g_surf_density(mu_disc(l, mu0=mu0d_s4g, h=h_disc_s4g), 1.61),
data_lim=data_lim, color='g', alpha=0.2, disk_scales=disk_scales, label='S4G maxdisc')
# SF 55/100/175
# 5533
#change this
total_gas_data_ = zip(r_g_dens, [He_coeff*(y_interp_(l[0], hi_r) + l[1]) for l in zip(r_g_dens, gas_dens)])[1:10]
total_gas_data_2 = zip(r_g_dens, [He_coeff*(y_interp_2(l[0], hi_r) + l[1]) for l in zip(r_g_dens, gas_dens)])[1:10]
plot_2f_vs_1f(ax=axes[4], total_gas_data=total_gas_data_,
epicycl=epicyclicFreq_real,
gas_approx=spl_gas,
sound_vel=sound_vel,
scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=MU0_r, h=hi_r), M_to_L=4.41, band='r'),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=MU0_r, h=hi_r), M_to_L=4.41, band='r'),
data_lim=data_lim, color='g', alpha=0.2, disk_scales=disk_scales, label='r(g-i) maxdisc')
total_gas_data_ = zip(r_g_dens, [He_coeff*(y_interp_(l[0], h_disc_R) + l[1]) for l in zip(r_g_dens, gas_dens)])[1:10]
total_gas_data_2 = zip(r_g_dens, [He_coeff*(y_interp_2(l[0], h_disc_R) + l[1]) for l in zip(r_g_dens, gas_dens)])[1:10]
plot_2f_vs_1f(ax=axes[4], total_gas_data=total_gas_data_2,
epicycl=epicyclicFreq_real,
gas_approx=spl_gas,
sound_vel=sound_vel,
scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L=7.15, band='R'),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L=7.15, band='R'),
data_lim=data_lim, color='r', alpha=0.2, disk_scales=disk_scales, label='R Noord maxdisc 2')
In [ ]:
In [16]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}
fig = plt.figure(figsize=[16, 6])
for ind, name in enumerate(['n1167', 'n3898']):
ax = plt.gca()
color = cm.rainbow(np.linspace(0, 1, 7))[ind]
phot = kenn_photom[name]
mud = locals()[name+'dict'][phot[0]]
h = locals()[name+'dict'][phot[1]]
# print name
# print mud
# print h
if name in ['n4258', 'n4725']:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
elif name == 'n2985':
total_gas_data=locals()[name+'dict']['total_gas_data_']
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in
zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
if total_gas_data[0][0] < 0.1:
total_gas_data = total_gas_data[1:]
# print total_gas_data
plot_2f_vs_1f(ax=ax, total_gas_data=total_gas_data,
epicycl=locals()[name+'dict']['epicyclicFreq_real'],
gas_approx=locals()[name+'dict']['spl_gas'],
sound_vel=locals()[name+'dict']['sound_vel'],
scale=locals()[name+'dict']['scale'],
sigma_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_min=locals()[name+'dict']['sig_R_maj_min'],
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
data_lim=locals()[name+'dict']['data_lim'],
color=color, alpha=0.2,
label=name+' '+phot[3])
plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
# plt.axvline(x=1, alpha=0.5)
plot_Q_levels(ax, [1.5, 2., 3.])
plt.legend()
plt.xlim(0, 130.)
plt.ylim(0, 1.3)
plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
plt.xlabel(r'$R$', fontsize=20)
plt.show()
Wang & Silk 2011 approximation:
In [45]:
def plot_WS_vs_2f(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, sfrange=None):
rr = zip(*total_gas_data)[0]
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))
invQwff_WS_min = [1./Qg(epicycl=epicycl(gas_approx, r, scale), sound_vel=sound_vel, gas_density=gd) + 1./Qs(
epicycl=epicycl(gas_approx, r, scale), sigma=sigma_max(r), star_density=star_density_min(r)) for r, gd in total_gas_data]
ax.plot(invQwff_WS_min, invQeff_min, 'd', color=color, alpha=0.8)
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))
invQwff_WS_max = [1./Qg(epicycl=epicycl(gas_approx, r, scale), sound_vel=sound_vel, gas_density=gd) + 1./Qs(
epicycl=epicycl(gas_approx, r, scale), sigma=sigma_min(r), star_density=star_density_max(r)) for r, gd in total_gas_data]
ax.plot(invQwff_WS_max, invQeff_max, 'o', color=color, alpha=0.8)
In [46]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}
fig = plt.figure(figsize=[8, 8])
for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
ax = plt.gca()
color = cm.rainbow(np.linspace(0, 1, 7))[ind]
phot = kenn_photom[name]
mud = locals()[name+'dict'][phot[0]]
h = locals()[name+'dict'][phot[1]]
# print name
# print mud
# print h
if name in ['n4258', 'n4725']:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
elif name == 'n2985':
total_gas_data=locals()[name+'dict']['total_gas_data_']
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in
zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
if total_gas_data[0][0] < 0.1:
total_gas_data = total_gas_data[1:]
# print total_gas_data
plot_WS_vs_2f(ax=ax, total_gas_data=total_gas_data,
epicycl=locals()[name+'dict']['epicyclicFreq_real'],
gas_approx=locals()[name+'dict']['spl_gas'],
sound_vel=locals()[name+'dict']['sound_vel'],
scale=locals()[name+'dict']['scale'],
sigma_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_min=locals()[name+'dict']['sig_R_maj_min'],
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
data_lim=locals()[name+'dict']['data_lim'],
color=color, alpha=0.2,
label=name+' '+phot[3],
sfrange=phot[4])
plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
plt.plot([0., 1.], [0., 1.], '--', alpha=0.1)
plt.legend(loc='lower right')
plt.xlim(0, 1.)
plt.ylim(0, 1.)
plt.ylabel(r'$Q_{eff}^{-1}$', fontsize=20)
plt.xlabel(r'$Q_{WS}^{-1}$', fontsize=20)
plt.grid()
plt.show()
In [42]:
def plot_WS_vs_2f_dist(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, sfrange=None):
rr = zip(*total_gas_data)[0]
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))
invQwff_WS_min = [1./Qg(epicycl=epicycl(gas_approx, r, scale), sound_vel=sound_vel, gas_density=gd) + 1./Qs(
epicycl=epicycl(gas_approx, r, scale), sigma=sigma_max(r), star_density=star_density_min(r)) for r, gd in total_gas_data]
ax.plot(rr, [l[0]/l[1] for l in zip(invQwff_WS_min, invQeff_min)], 'd', color=color, alpha=0.8)
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))
invQwff_WS_max = [1./Qg(epicycl=epicycl(gas_approx, r, scale), sound_vel=sound_vel, gas_density=gd) + 1./Qs(
epicycl=epicycl(gas_approx, r, scale), sigma=sigma_min(r), star_density=star_density_max(r)) for r, gd in total_gas_data]
ax.plot(rr, [l[0]/l[1] for l in zip(invQwff_WS_max, invQeff_max)], 'o', color=color, alpha=0.8)
In [43]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}
fig = plt.figure(figsize=[16, 6])
for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
ax = plt.gca()
color = cm.rainbow(np.linspace(0, 1, 7))[ind]
phot = kenn_photom[name]
mud = locals()[name+'dict'][phot[0]]
h = locals()[name+'dict'][phot[1]]
# print name
# print mud
# print h
if name in ['n4258', 'n4725']:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
elif name == 'n2985':
total_gas_data=locals()[name+'dict']['total_gas_data_']
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in
zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
if total_gas_data[0][0] < 0.1:
total_gas_data = total_gas_data[1:]
# print total_gas_data
plot_WS_vs_2f_dist(ax=ax, total_gas_data=total_gas_data,
epicycl=locals()[name+'dict']['epicyclicFreq_real'],
gas_approx=locals()[name+'dict']['spl_gas'],
sound_vel=locals()[name+'dict']['sound_vel'],
scale=locals()[name+'dict']['scale'],
sigma_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_min=locals()[name+'dict']['sig_R_maj_min'],
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
data_lim=locals()[name+'dict']['data_lim'],
color=color, alpha=0.2,
label=name+' '+phot[3],
sfrange=phot[4])
plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
# plt.plot([0., 1.], [0., 1.], '--', alpha=0.1)
plt.legend(loc='upper right')
plt.xlim(0)
plt.ylim(1.0)
plt.ylabel(r'$Q_{WS}^{-1}/Q_{eff}^{-1}$', fontsize=20)
plt.xlabel(r'$R$', fontsize=20)
plt.grid()
plt.show()
In [44]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}
fi, axes = plt.subplots(ncols=1, nrows=7, figsize=[20., 26.], sharex=True)nbv
for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
# for ind, name in enumerate(['n338']):
ax = axes[ind]
color = cm.rainbow(np.linspace(0, 1, 7))[ind]
phot = kenn_photom[name]
mud = locals()[name+'dict'][phot[0]]
h = locals()[name+'dict'][phot[1]]
# print name
# print mud
# print h
if name in ['n4258', 'n4725']:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
elif name == 'n2985':
total_gas_data=locals()[name+'dict']['total_gas_data_']
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in
zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
if total_gas_data[0][0] < 0.1:
total_gas_data = total_gas_data[1:]
# print total_gas_data
plot_2f_vs_1f(ax=ax, total_gas_data=total_gas_data,
epicycl=locals()[name+'dict']['epicyclicFreq_real'],
gas_approx=locals()[name+'dict']['spl_gas'],
sound_vel=locals()[name+'dict']['sound_vel'],
scale=locals()[name+'dict']['scale'],
sigma_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_min=locals()[name+'dict']['sig_R_maj_min'],
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
data_lim=locals()[name+'dict']['data_lim'],
color=color, alpha=0.2,
label='')
ax.set_xlim(0, 130)
ax.set_ylim(0, 1.2)
ax.xaxis.grid(True)
if ind == 6:
ax.set_xlabel(r'$R, arcsec$', fontsize=20)
ax.set_ylabel(r'$Q_{eff}^{-1}$', fontsize=20)
for q_ in [1., 1.5, 2., 3.]:
ax.axhline(y=1./q_, lw=10, alpha=0.15, color='grey')
ax.text(95, 0.85, 'NGC '+name[1:]+' $'+phot[3]+'$', fontsize=30)
# 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.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
# plt.axvline(x=1, alpha=0.5)
# plot_Q_levels(ax, [1.5, 2., 3.])
# plt.legend()
# plt.xlim(0, 130.)
# plt.ylim(0, 1.3)
# plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
# plt.xlabel(r'$R$', fontsize=20)
plt.subplots_adjust(wspace=0, hspace=0)
plt.show()
In [51]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}
fi, axes = plt.subplots(ncols=1, nrows=7, figsize=[20., 26.], sharex=True)
for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
# for ind, name in enumerate(['n338']):
ax = axes[ind]
color = cm.rainbow(np.linspace(0, 1, 7))[ind]
phot = kenn_photom[name]
mud = locals()[name+'dict'][phot[0]]
h = locals()[name+'dict'][phot[1]]
# print name
# print mud
# print h
if name in ['n4258', 'n4725']:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
elif name == 'n2985':
total_gas_data=locals()[name+'dict']['total_gas_data_']
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in
zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
if total_gas_data[0][0] < 0.1:
total_gas_data = total_gas_data[1:]
# print total_gas_data
plot_2f_vs_1f(ax=ax, total_gas_data=total_gas_data,
epicycl=locals()[name+'dict']['epicyclicFreq_real'],
gas_approx=locals()[name+'dict']['spl_gas'],
sound_vel=locals()[name+'dict']['sound_vel'],
scale=locals()[name+'dict']['scale'],
sigma_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_min=locals()[name+'dict']['sig_R_maj_min'],
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]),
data_lim=locals()[name+'dict']['data_lim'],
color=color, alpha=0.2,
label='')
ax.set_xlim(0., 130)
ax.set_ylim(0.01, 1.5)
ax.xaxis.grid(True)
if ind == 6:
ax.set_xlabel(r'$R, arcsec$', fontsize=20)
ax.set_ylabel(r'$Q_{eff}^{-1}$', fontsize=20)
for q_ in [1., 1.5, 2., 3.]:
ax.axhline(y=1./q_, lw=10, alpha=0.15, color='grey')
ax.text(95, 0.85, 'NGC '+name[1:]+' $'+phot[3]+'$', fontsize=30)
# 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)
ax.set_yscale('log')
# plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
# plt.axvline(x=1, alpha=0.5)
# plot_Q_levels(ax, [1.5, 2., 3.])
# plt.legend()
# plt.xlim(0, 130.)
# plt.ylim(0, 1.3)
# plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
# plt.xlabel(r'$R$', fontsize=20)
plt.subplots_adjust(wspace=0, hspace=0)
plt.show()
In [42]:
from random import shuffle
random.seed(42)
kenn_photom={'n338_I' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60., 15.],
'n338_B' : ['mu0d_Bc', 'h_disc_B', 6.99, 'B', 60., 15.],
'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40., 6.7],
'n2985_K' : ['mudK', 'hK', 1.42, 'K', 70., 14.],
'n2985_S4G' : [('mu0d_s4g', 'mu0d_s4g_2'), ('h_disc_s4g', 'h_disc_s4g_2'), 1.35, 'S4G', 70., 6.3],
'n3898_R' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80., 8.8],
'n3898_R2' : [('mu_inner_approx', 'mu_out'), ('h_approx', 'h_out'), 2.93, 'R', 80., 8.8],
'n4258_I' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150., 15.],
'n4258_S4G' : ['mu0d_36', 'h_disc_36', 0.38, 'S4G', 150., 15.],
'n4725_H' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55., 10.14],
'n4725_S4G' : ['mu0d_s4g', 'h_disc_s4g', 1.61, 'S4G', 55., 10.14],
'n5533_r' : ['MU0_r', 'hi_r', 4.41, 'r', 100., 8.9],
'n5533_R' : ['mu0d_Rc', 'h_disc_R', 7.15, 'R', 100., 9.9]}
color_shuffle = range(13)
shuffle(color_shuffle)
In [43]:
%%time
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, '')
plot_Q_levels(ax, [1., 1.5, 2., 3.])
ax.legend(fontsize=20)
fi, axes = plt.subplots(ncols=1, nrows=7, figsize=[20., 26.], sharex=True)
color_counter = 0
for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
# for ind, name in enumerate(['n3898']):
ax = axes[ind]
for key in kenn_photom.keys():
if name in key:
print key
phot = kenn_photom[key]
color = cm.rainbow(np.linspace(0, 1, 13))[color_shuffle[color_counter]]
try:
h = locals()[name+'dict'][phot[1][0]]
except KeyError:
h = locals()[name+'dict'][phot[1]]
# print name
# print mud
# print h
if name in ['n4258', 'n4725']:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
elif name == 'n2985':
total_gas_data=locals()[name+'dict']['total_gas_data_']
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in
zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
while total_gas_data[0][0] < phot[5]:
total_gas_data = total_gas_data[1:]
# print total_gas_data
ML = phot[2]
band = phot[3]
if type(phot[0]) == tuple and phot[3] == 'S4G':
h1 = globals()[name+'dict'][phot[1][0]]
h2 = globals()[name+'dict'][phot[1][1]]
mu01 = globals()[name+'dict'][phot[0][0]]
mu02 = globals()[name+'dict'][phot[0][1]]
star_density = lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu01, h=h1), ML),
lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu02, h=h2), ML)))(l)
elif type(phot[0]) == tuple:
h1 = globals()[name+'dict'][phot[1][0]]
h2 = globals()[name+'dict'][phot[1][1]]
mu01 = globals()[name+'dict'][phot[0][0]]
mu02 = globals()[name+'dict'][phot[0][1]]
star_density = lambda l: tot_dens((lambda l1: surf_density(mu=mu_disc(l1, mu0=mu01, h=h1), M_to_L=ML, band=band),
lambda l2: surf_density(mu=mu_disc(l2, mu0=mu02, h=h2), M_to_L=ML, band=band)))(l)
elif phot[3] == 'S4G':
star_density = lambda l: s4g_surf_density(mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), ML)
else:
star_density = lambda l: surf_density(mu=mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), M_to_L=ML, band=band)
if name == 'n338':
epicycl = locals()[name+'dict']['epicyclicFreq_real_']
else:
epicycl = locals()[name+'dict']['epicyclicFreq_real']
plot_2f_vs_1f(ax=ax, total_gas_data=total_gas_data,
epicycl=epicycl,
gas_approx=locals()[name+'dict']['spl_gas'],
sound_vel=locals()[name+'dict']['sound_vel'],
scale=locals()[name+'dict']['scale'],
sigma_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_min=locals()[name+'dict']['sig_R_maj_min'],
star_density_max=star_density,
star_density_min=star_density,
data_lim=locals()[name+'dict']['data_lim'],
color=color, alpha=0.3,
label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$',
disk_scales=[(h, band)])
ax.set_xlim(0, 130)
if name == 'n338':
ax.set_ylim(0, 1.4)
else:
ax.set_ylim(0, 1.2)
ax.xaxis.grid(True)
locals()[name+'dict']['plot_SF'](ax)
if ind == 6:
ax.set_xlabel(r'$R, arcsec$', fontsize=20)
ax.set_ylabel(r'$Q_{eff}^{-1}$', fontsize=20)
for q_ in [1., 1.5, 2., 3.]:
ax.axhline(y=1./q_, lw=10, alpha=0.15, color='grey')
# ax.text(95, 0.85, 'NGC '+name[1:]+' $'+phot[3]+'$', fontsize=30)
# 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)
color_counter+=1
# plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
# plt.axvline(x=1, alpha=0.5)
# plot_Q_levels(ax, [1.5, 2., 3.])
# plt.legend()
# plt.xlim(0, 130.)
# plt.ylim(0, 1.3)
# plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
# plt.xlabel(r'$R$', fontsize=20)
plt.subplots_adjust(wspace=0, hspace=0)
plt.show()
In [ ]:
In [49]:
%%time
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_mi, 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_ma, 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, invQs_mi, '--', color=color, alpha=0.75)
ax.plot(rr, invQs_ma, '--', color=color, alpha=0.75)
ax.plot(rr, invQg, '--', color='b')
fi, axes = plt.subplots(ncols=1, nrows=13, figsize=[20., 50.], sharex=True)
for ind, key in enumerate(kenn_photom.keys()):
ax = axes[ind]
name = key.split('_')[0]
print key, name
phot = kenn_photom[key]
color = cm.rainbow(np.linspace(0, 1, 13))[color_shuffle[ind]]
try:
h = locals()[name+'dict'][phot[1][0]]
except KeyError:
h = locals()[name+'dict'][phot[1]]
# print name
# print mud
# print h
if name in ['n4258', 'n4725']:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
elif name == 'n2985':
total_gas_data=locals()[name+'dict']['total_gas_data_']
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in
zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
while total_gas_data[0][0] < phot[5]:
total_gas_data = total_gas_data[1:]
# print total_gas_data
ML = phot[2]
band = phot[3]
if type(phot[0]) == tuple and phot[3] == 'S4G':
h1 = globals()[name+'dict'][phot[1][0]]
h2 = globals()[name+'dict'][phot[1][1]]
mu01 = globals()[name+'dict'][phot[0][0]]
mu02 = globals()[name+'dict'][phot[0][1]]
star_density = lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu01, h=h1), ML),
lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu02, h=h2), ML)))(l)
elif type(phot[0]) == tuple:
h1 = globals()[name+'dict'][phot[1][0]]
h2 = globals()[name+'dict'][phot[1][1]]
mu01 = globals()[name+'dict'][phot[0][0]]
mu02 = globals()[name+'dict'][phot[0][1]]
star_density = lambda l: tot_dens((lambda l1: surf_density(mu=mu_disc(l1, mu0=mu01, h=h1), M_to_L=ML, band=band),
lambda l2: surf_density(mu=mu_disc(l2, mu0=mu02, h=h2), M_to_L=ML, band=band)))(l)
elif phot[3] == 'S4G':
star_density = lambda l: s4g_surf_density(mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), ML)
else:
star_density = lambda l: surf_density(mu=mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), M_to_L=ML, band=band)
if name == 'n338':
epicycl = locals()[name+'dict']['epicyclicFreq_real_']
else:
epicycl = locals()[name+'dict']['epicyclicFreq_real']
plot_param_depend(ax=ax,
epicycl=epicycl,
gas_approx=locals()[name+'dict']['spl_gas'],
sound_vel=list(np.linspace(4., 20., 50)),
N=50,
scale=locals()[name+'dict']['scale'],
sigma_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_min=locals()[name+'dict']['sig_R_maj_min'],
star_density_max=star_density,
star_density_min=star_density,
data_lim=locals()[name+'dict']['data_lim'],
color=color, alpha=0.3,
label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$',
total_gas_data=total_gas_data)
plot_2f_vs_1f_(ax=ax, total_gas_data=total_gas_data,
epicycl=epicycl,
gas_approx=locals()[name+'dict']['spl_gas'],
sound_vel=locals()[name+'dict']['sound_vel'],
scale=locals()[name+'dict']['scale'],
sigma_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_min=locals()[name+'dict']['sig_R_maj_min'],
star_density_max=star_density,
star_density_min=star_density,
data_lim=locals()[name+'dict']['data_lim'],
color=color, alpha=0.3,
label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$',
disk_scales=[(h, band)])
ax.set_xlim(0, 130)
if name == 'n338' or name == 'n5533':
ax.set_ylim(0, 1.4)
else:
ax.set_ylim(0, 1.25)
ax.xaxis.grid(True)
locals()[name+'dict']['plot_SF'](ax)
if ind == 6:
ax.set_xlabel(r'$R, arcsec$', fontsize=20)
ax.set_ylabel(r'$Q_{eff}^{-1}$', fontsize=20)
for q_ in [1., 1.5, 2., 3.]:
ax.axhline(y=1./q_, lw=10, alpha=0.15, color='grey')
ax.text(95, 0.95, 'NGC '+name[1:]+' $'+phot[3]+'$', fontsize=30)
# 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.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
# plt.axvline(x=1, alpha=0.5)
# plot_Q_levels(ax, [1.5, 2., 3.])
# plt.legend()
# plt.xlim(0, 130.)
# plt.ylim(0, 1.3)
# plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
# plt.xlabel(r'$R$', fontsize=20)
plt.subplots_adjust(wspace=0, hspace=0)
plt.show()
In [ ]:
In [43]:
%%time
import scipy.interpolate
def plot_WS_vs_2f(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, sfrange=None):
rr = zip(*total_gas_data)[0]
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))
invQwff_WS_min = [1./Qg(epicycl=epicycl(gas_approx, r, scale), sound_vel=sound_vel, gas_density=gd) + 1./Qs(
epicycl=epicycl(gas_approx, r, scale), sigma=sigma_max(r), star_density=star_density_min(r)) for r, gd in total_gas_data]
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))
invQwff_WS_max = [1./Qg(epicycl=epicycl(gas_approx, r, scale), sound_vel=sound_vel, gas_density=gd) + 1./Qs(
epicycl=epicycl(gas_approx, r, scale), sigma=sigma_min(r), star_density=star_density_max(r)) for r, gd in total_gas_data]
ax.plot(invQwff_WS_min, invQeff_min, '^', color=color, alpha=0.4, ms=8, lw=1)
ax.plot(invQwff_WS_max, invQeff_max, 'D', color=color, alpha=0.4, ms=8, lw=1)
return invQwff_WS_min, invQwff_WS_max
def plot_RF13_vs_2F(ax=None, 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, color=None):
'''Плотности газа передаются не скорр. за гелий.'''
if thin:
romeo_Q = romeo_Qinv_thin
else:
romeo_Q = romeo_Qinv
totgas = zip(r_g_dens, [He_coeff*(l[0]+l[1]) for l in zip(HI_gas_dens, CO_gas_dens)])
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)
romeo_min.append(rom)
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)
romeo_max.append(rom)
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(romeo_min, invQeff_min, 's', color=color, alpha=0.4, ms=8)
ax.plot(romeo_max, invQeff_max, 'o', color=color, alpha=0.4, ms=8)
return invQeff_min, invQeff_max, romeo_min, romeo_max
markers = {'.': 'point', ',': 'pixel', 'o': 'circle', 'v': 'triangle_down', '^': 'triangle_up', '<': 'triangle_left', '>': 'triangle_right', '1': 'tri_down', '2':
'tri_up', '3': 'tri_left', '4': 'tri_right', '8': 'octagon', 's': 'square', 'p': 'pentagon', '*': 'star', 'h': 'hexagon1', 'H': 'hexagon2', '+': 'plus',
'x': 'x', 'D': 'diamond', 'd': 'thin_diamond', '|': 'vline', '_': 'hline', 'P': 'plus_filled', 'X': 'x_filled', 0: 'tickleft', 1: 'tickright', 2: 'tickup',
3: 'tickdown', 4: 'caretleft', 5: 'caretright', 6: 'caretup', 7: 'caretdown', 8: 'caretleftbase', 9: 'caretrightbase', 10: 'caretupbase', 11: 'caretdownbase'}
kenn_photom={
'n338_I' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60., 15.],
'n338_B' : ['mu0d_Bc', 'h_disc_B', 6.99, 'B', 60., 15.],
'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40., 6.7],
'n2985_K' : ['mudK', 'hK', 1.42, 'K', 70., 14.],
'n2985_S4G' : [('mu0d_s4g', 'mu0d_s4g_2'), ('h_disc_s4g', 'h_disc_s4g_2'), 1.35, 'S4G', 70., 6.3],
'n3898_R' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80., 8.8],
'n3898_R2' : [('mu_inner_approx', 'mu_out'), ('h_approx', 'h_out'), 2.93, 'R', 80., 8.8],
'n4258_I' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150., 15.],
'n4258_S4G' : ['mu0d_36', 'h_disc_36', 0.38, 'S4G', 150., 15.],
'n4725_H' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55., 10.14],
'n4725_S4G' : ['mu0d_s4g', 'h_disc_s4g', 1.61, 'S4G', 55., 10.14],
'n5533_r' : ['MU0_r', 'hi_r', 4.41, 'r', 100., 8.9],
'n5533_R' : ['mu0d_Rc', 'h_disc_R', 7.15, 'R', 100., 9.9]
}
fig, (ax, ax2) = plt.subplots(ncols=2, nrows=1, figsize=[16, 8])
color_counter = 0
color_shuffle = range(13)
shuffle(color_shuffle)
saved_data_ = []
for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
# for ind, name in enumerate(['n4258']):
for key in kenn_photom.keys():
if name in key:
print key
phot = kenn_photom[key]
color = cm.rainbow(np.linspace(0, 1, 13))[color_shuffle[color_counter]]
try:
h = locals()[name+'dict'][phot[1][0]]
except KeyError:
h = locals()[name+'dict'][phot[1]]
# print name
# print mud
# print h
if name in ['n4258', 'n4725']:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
elif name == 'n2985':
total_gas_data=locals()[name+'dict']['total_gas_data_']
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in
zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
while total_gas_data[0][0] < phot[5]:
total_gas_data = total_gas_data[1:]
while total_gas_data[0][-1] > 130.:
total_gas_data = total_gas_data[:-1]
# print total_gas_data
ML = phot[2]
band = phot[3]
if type(phot[0]) == tuple and phot[3] == 'S4G':
h1 = globals()[name+'dict'][phot[1][0]]
h2 = globals()[name+'dict'][phot[1][1]]
mu01 = globals()[name+'dict'][phot[0][0]]
mu02 = globals()[name+'dict'][phot[0][1]]
star_density = lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu01, h=h1), ML),
lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu02, h=h2), ML)))(l)
elif type(phot[0]) == tuple:
h1 = globals()[name+'dict'][phot[1][0]]
h2 = globals()[name+'dict'][phot[1][1]]
mu01 = globals()[name+'dict'][phot[0][0]]
mu02 = globals()[name+'dict'][phot[0][1]]
star_density = lambda l: tot_dens((lambda l1: surf_density(mu=mu_disc(l1, mu0=mu01, h=h1), M_to_L=ML, band=band),
lambda l2: surf_density(mu=mu_disc(l2, mu0=mu02, h=h2), M_to_L=ML, band=band)))(l)
elif phot[3] == 'S4G':
star_density = lambda l: s4g_surf_density(mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), ML)
else:
star_density = lambda l: surf_density(mu=mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), M_to_L=ML, band=band)
if name == 'n338':
epicycl = locals()[name+'dict']['epicyclicFreq_real_']
else:
epicycl = locals()[name+'dict']['epicyclicFreq_real']
invQwff_WS_min, invQwff_WS_max = plot_WS_vs_2f(ax=ax, total_gas_data=total_gas_data,
epicycl=epicycl,
gas_approx=locals()[name+'dict']['spl_gas'],
sound_vel=locals()[name+'dict']['sound_vel'],
scale=locals()[name+'dict']['scale'],
sigma_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_min=locals()[name+'dict']['sig_R_maj_min'],
star_density_max=star_density,
star_density_min=star_density,
data_lim=locals()[name+'dict']['data_lim'],
color=color, alpha=0.5,
label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$',
sfrange=phot[4])
ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$')
if name in ['n4258', 'n4725']:
r_mol_dens = locals()[name+'dict']['r_mol_dens']
mol_dens = locals()[name+'dict']['mol_dens']
y_interp = scipy.interpolate.interp1d(list(r_mol_dens), list(mol_dens))
def y_interp_(r):
if r <= min(r_mol_dens):
return y_interp(min(r_mol_dens))
elif r < max(r_mol_dens):
return y_interp(r)
else:
return 0.
# не совсем точно так, но все же
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [l[1]/He_coeff - y_interp_(l[0]) for l in zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])],
[y_interp_(l) for l in locals()[name+'dict']['r_g_dens']])[1:15]
elif name == 'n2985':
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'],
[locals()[name+'dict']['y_interp_'](l) for l in locals()[name+'dict']['r_g_dens']])[1:15]
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'],
[locals()[name+'dict']['y_interp_'](l, h) for l in locals()[name+'dict']['r_g_dens']])[1:15]
while total_gas_data[0][0] < phot[5]:
total_gas_data = total_gas_data[1:]
while total_gas_data[0][-1] > 130.:
total_gas_data = total_gas_data[:-1]
invQeff_min, invQeff_max, romeo_min, romeo_max = plot_RF13_vs_2F(ax=ax, r_g_dens=zip(*total_gas_data)[0],
HI_gas_dens=zip(*total_gas_data)[1],
CO_gas_dens=zip(*total_gas_data)[2],
epicycl=epicycl,
sound_vel=6.,
sigma_R_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_R_min=locals()[name+'dict']['sig_R_maj_min'],
star_density=star_density,
alpha_max=0.7,
alpha_min=0.3,
scale=locals()[name+'dict']['scale'],
gas_approx=locals()[name+'dict']['spl_gas'],
thin=True,
color=color)
saved_data_.append(zip(zip(*total_gas_data)[0], invQeff_min, invQeff_max, romeo_min, romeo_max, invQwff_WS_min, invQwff_WS_max))
ax2.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$')
color_counter+=1
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
ax.legend(loc='lower right')
ax.set_xlim(0, 1.75)
ax.set_ylim(0, 1.05)
ax.set_ylabel(r'$Q_{2F}^{-1}$', fontsize=20)
ax.set_xlabel(r'$Q_{WS,RF}^{-1}$', fontsize=20)
ax.grid()
ax2.legend(loc='lower right')
ax2.set_xlim(0, 1.05)
ax2.set_ylim(0, 1.05)
ax2.set_ylabel(r'$Q_{eff}^{-1}$', fontsize=20)
ax2.set_xlabel(r'$Q_{RF}^{-1}$', fontsize=20)
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)
plt.show()
In [44]:
os.chdir(return_back_path)
In [45]:
np.save('2f_appr_data.npy', saved_data_)
In [46]:
approx_data = np.load('2f_appr_data.npy')
In [47]:
fig = plt.figure(figsize=[12,12])
ax = plt.gca()
color_counter = 0
color_shuffle = range(13)
shuffle(color_shuffle)
for ind2, key in enumerate(['n338_I', 'n338_B','n1167', 'n2985_S4G','n2985_K','n3898_R2','n3898_R','n4258_S4G','n4258_I','n4725_S4G','n4725_H','n5533_R','n5533_r']):
name = key.split('_')[0]
phot = kenn_photom[key]
rr, invQeff_min, invQeff_max, romeo_min, romeo_max, invQwff_WS_min, invQwff_WS_max = zip(*approx_data[ind2])
color = cm.rainbow(np.linspace(0, 1, 13))[color_shuffle[color_counter]]
ax.plot(romeo_min, invQeff_min, 's', color=color, alpha=0.4, ms=8)
ax.plot(romeo_max, invQeff_max, 'o', color=color, alpha=0.4, ms=8)
ax.plot(invQwff_WS_min, invQeff_min, '^', color=color, alpha=0.4, ms=8, lw=1)
ax.plot(invQwff_WS_max, invQeff_max, 'D', color=color, alpha=0.4, ms=8, lw=1)
ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$', lw=5)
color_counter += 1
ax.legend(loc='lower right')
ax.set_xlim(0, 1.75)
ax.set_ylim(0, 1.05)
ax.set_ylabel(r'$Q_{2F}^{-1}$', fontsize=20)
ax.set_xlabel(r'$Q_{WS,RF}^{-1}$', fontsize=20)
ax.grid()
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
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')
line, = ax.plot([0., 0.89*2], [0., 2.0], '--', alpha=0.3, color='gray')
label_line(line, 'y={:2.2f}x'.format(1./0.89), 0.8*0.89, 0.8, color='black')
ax.plot([0., 2.0], [0., 2.0], '--', alpha=0.9, color='gray', lw=3)
plt.show()
In [48]:
fig = plt.figure(figsize=[24,8])
ax = plt.gca()
color_counter = 0
color_shuffle = range(13)
shuffle(color_shuffle)
for ind2, key in enumerate(['n338_I', 'n338_B','n1167', 'n2985_S4G','n2985_K','n3898_R2','n3898_R','n4258_S4G','n4258_I','n4725_S4G','n4725_H','n5533_R','n5533_r']):
name = key.split('_')[0]
phot = kenn_photom[key]
rr, invQeff_min, invQeff_max, romeo_min, romeo_max, invQwff_WS_min, invQwff_WS_max = zip(*approx_data[ind2])
color = cm.rainbow(np.linspace(0, 1, 13))[color_shuffle[color_counter]]
ax.plot(rr, [l[0]/l[1] for l in zip(invQeff_min, romeo_min)], 's-', color=color, alpha=0.4, ms=8)
ax.plot(rr, [l[0]/l[1] for l in zip(invQeff_max, romeo_max)], 'o-', color=color, alpha=0.4, ms=8)
# ax.plot(rr, romeo_min, '^-', color=color, alpha=0.4, ms=8)
# ax.plot(rr, romeo_max, 'D-', color=color, alpha=0.4, ms=8)
# ax.plot(invQwff_WS_min, invQeff_min, '^', color=color, alpha=0.4, ms=8, lw=1)
# ax.plot(invQwff_WS_max, invQeff_max, 'D', color=color, alpha=0.4, ms=8, lw=1)
ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$', lw=5)
color_counter += 1
ax.legend(loc='lower right')
ax.set_xlim(0)
ax.set_ylim(0.9, 1.6)
# ax.set_ylabel(r'$Q_{2F}^{-1}$', fontsize=20)
# ax.set_xlabel(r'$Q_{WS,RF}^{-1}$', fontsize=20)
ax.grid()
plt.show()
Плохо видно - разнесем:
In [49]:
fig = plt.figure(figsize=[24,8])
ax = plt.gca()
color_counter = 0
color_shuffle = range(13)
shuffle(color_shuffle)
for ind2, key in enumerate(['n338_I', 'n338_B','n1167', 'n2985_S4G','n2985_K','n3898_R2','n3898_R','n4258_S4G','n4258_I','n4725_S4G','n4725_H','n5533_R','n5533_r']):
name = key.split('_')[0]
phot = kenn_photom[key]
rr, invQeff_min, invQeff_max, romeo_min, romeo_max, invQwff_WS_min, invQwff_WS_max = zip(*approx_data[ind2])
color = cm.rainbow(np.linspace(0, 1, 13))[color_shuffle[color_counter]]
a_ = filter(lambda l: l[1] > 1.1 and l[1] < 1.2, zip(rr, [l[0]/l[1] for l in zip(invQeff_min, romeo_min)]))
try:
ax.plot(zip(*a_)[0], map(lambda l: l+ind2*0.1, zip(*a_)[1]), 's-', color=color, alpha=0.4, ms=8)
# ax.plot(rr, [ind2 + l[0]/l[1] for l in zip(invQeff_max, romeo_max)], 'o-', color=color, alpha=0.4, ms=8)
# ax.plot(rr, romeo_min, '^-', color=color, alpha=0.4, ms=8)
# ax.plot(rr, romeo_max, 'D-', color=color, alpha=0.4, ms=8)
# ax.plot(invQwff_WS_min, invQeff_min, '^', color=color, alpha=0.4, ms=8, lw=1)
# ax.plot(invQwff_WS_max, invQeff_max, 'D', color=color, alpha=0.4, ms=8, lw=1)
ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$', lw=5)
except Exception:
print name
color_counter += 1
ax.legend(loc='lower right')
ax.set_xlim(0)
ax.set_ylim(1.2, 2.5)
# ax.set_ylabel(r'$Q_{2F}^{-1}$', fontsize=20)
# ax.set_xlabel(r'$Q_{WS,RF}^{-1}$', fontsize=20)
ax.grid()
plt.show()
In [77]:
import scipy.interpolate
from math import pi
def plot_RF13_vs_2F(ax=None, r_g_dens=None, HI_gas_dens=None, CO_gas_dens=None, epicycl=None, sound_vel=6., 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, verbose=True):
'''Плотности газа передаются не скорр. за гелий.'''
if thin:
romeo_Q = romeo_Qinv_thin
else:
romeo_Q = romeo_Qinv
totgas = zip(r_g_dens, [He_coeff*(l[0]+l[1]) for l in zip(HI_gas_dens, CO_gas_dens)])
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, verbose=verbose)
romeo_min.append(rom)
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, verbose=verbose)
romeo_max.append(rom)
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(romeo_min, invQeff_min, 's', color=color, alpha=0.4, ms=8)
# ax.plot(romeo_max, invQeff_max, 'o', color=color, alpha=0.4, ms=8)
return invQeff_min, invQeff_max, romeo_min, romeo_max
kenn_photom={
'n338_I' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60., 15.],
'n338_B' : ['mu0d_Bc', 'h_disc_B', 6.99, 'B', 60., 15.],
'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40., 6.7],
'n2985_K' : ['mudK', 'hK', 1.42, 'K', 70., 14.],
'n2985_S4G' : [('mu0d_s4g', 'mu0d_s4g_2'), ('h_disc_s4g', 'h_disc_s4g_2'), 1.35, 'S4G', 70., 6.3],
'n3898_R' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80., 8.8],
'n3898_R2' : [('mu_inner_approx', 'mu_out'), ('h_approx', 'h_out'), 2.93, 'R', 80., 8.8],
'n4258_I' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150., 15.],
'n4258_S4G' : ['mu0d_36', 'h_disc_36', 0.38, 'S4G', 150., 15.],
'n4725_H' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55., 10.14],
'n4725_S4G' : ['mu0d_s4g', 'h_disc_s4g', 1.61, 'S4G', 55., 10.14],
'n5533_r' : ['MU0_r', 'hi_r', 4.41, 'r', 100., 8.9],
'n5533_R' : ['mu0d_Rc', 'h_disc_R', 7.15, 'R', 100., 9.9]
}
fig = plt.figure(figsize=[16, 8])
color_counter = 0
color_shuffle = range(13)
shuffle(color_shuffle)
ax = plt.gca()
for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
# for ind, name in enumerate(['n1167']):
for key in kenn_photom.keys():
if name in key:
print key
phot = kenn_photom[key]
color = cm.rainbow(np.linspace(0, 1, 13))[color_shuffle[color_counter]]
try:
h = locals()[name+'dict'][phot[1][0]]
except KeyError:
h = locals()[name+'dict'][phot[1]]
# print name
# print mud
# print h
if name in ['n4258', 'n4725']:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
elif name == 'n2985':
total_gas_data=locals()[name+'dict']['total_gas_data_']
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in
zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
while total_gas_data[0][0] < phot[5]:
total_gas_data = total_gas_data[1:]
while total_gas_data[0][-1] > 130.:
total_gas_data = total_gas_data[:-1]
# print total_gas_data
ML = phot[2]
band = phot[3]
if type(phot[0]) == tuple and phot[3] == 'S4G':
h1 = globals()[name+'dict'][phot[1][0]]
h2 = globals()[name+'dict'][phot[1][1]]
mu01 = globals()[name+'dict'][phot[0][0]]
mu02 = globals()[name+'dict'][phot[0][1]]
star_density = lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu01, h=h1), ML),
lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu02, h=h2), ML)))(l)
elif type(phot[0]) == tuple:
h1 = globals()[name+'dict'][phot[1][0]]
h2 = globals()[name+'dict'][phot[1][1]]
mu01 = globals()[name+'dict'][phot[0][0]]
mu02 = globals()[name+'dict'][phot[0][1]]
star_density = lambda l: tot_dens((lambda l1: surf_density(mu=mu_disc(l1, mu0=mu01, h=h1), M_to_L=ML, band=band),
lambda l2: surf_density(mu=mu_disc(l2, mu0=mu02, h=h2), M_to_L=ML, band=band)))(l)
elif phot[3] == 'S4G':
star_density = lambda l: s4g_surf_density(mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), ML)
else:
star_density = lambda l: surf_density(mu=mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), M_to_L=ML, band=band)
if name == 'n338':
epicycl = locals()[name+'dict']['epicyclicFreq_real_']
else:
epicycl = locals()[name+'dict']['epicyclicFreq_real']
# invQwff_WS_min, invQwff_WS_max = plot_WS_vs_2f(ax=ax, total_gas_data=total_gas_data,
# epicycl=epicycl,
# gas_approx=locals()[name+'dict']['spl_gas'],
# sound_vel=locals()[name+'dict']['sound_vel'],
# scale=locals()[name+'dict']['scale'],
# sigma_max=locals()[name+'dict']['sig_R_maj_max'],
# sigma_min=locals()[name+'dict']['sig_R_maj_min'],
# star_density_max=star_density,
# star_density_min=star_density,
# data_lim=locals()[name+'dict']['data_lim'],
# color=color, alpha=0.5,
# label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$',
# sfrange=phot[4])
# ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$')
if name in ['n4258', 'n4725']:
r_mol_dens = locals()[name+'dict']['r_mol_dens']
mol_dens = locals()[name+'dict']['mol_dens']
y_interp = scipy.interpolate.interp1d(list(r_mol_dens), list(mol_dens))
def y_interp_(r):
if r <= min(r_mol_dens):
return y_interp(min(r_mol_dens))
elif r < max(r_mol_dens):
return y_interp(r)
else:
return 0.
# не совсем точно так, но все же
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [l[1]/He_coeff - y_interp_(l[0]) for l in zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])],
[y_interp_(l) for l in locals()[name+'dict']['r_g_dens']])[1:15]
elif name == 'n2985':
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'],
[locals()[name+'dict']['y_interp_'](l) for l in locals()[name+'dict']['r_g_dens']])[1:15]
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'],
[locals()[name+'dict']['y_interp_'](l, h) for l in locals()[name+'dict']['r_g_dens']])[1:15]
while total_gas_data[0][0] < phot[5]:
total_gas_data = total_gas_data[1:]
while total_gas_data[0][-1] > 130.:
total_gas_data = total_gas_data[:-1]
invQeff_min, invQeff_max, romeo_min, romeo_max = plot_RF13_vs_2F(ax=ax, r_g_dens=zip(*total_gas_data)[0],
HI_gas_dens=zip(*total_gas_data)[1],
CO_gas_dens=zip(*total_gas_data)[2],
epicycl=epicycl,
sound_vel=sound_vel,
sigma_R_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_R_min=locals()[name+'dict']['sig_R_maj_min'],
star_density=star_density,
alpha_max=0.7,
alpha_min=0.3,
scale=locals()[name+'dict']['scale'],
gas_approx=locals()[name+'dict']['spl_gas'],
thin=True,
color=color, verbose=False)
romeo_min, components = zip(*romeo_min)
a_ = filter(lambda l: l[1] > 1.1 and l[1] < 1.2, zip(rr, [l[0]/l[1] for l in zip(invQeff_min, romeo_min)], components))
# print a_
try:
for _ in a_:
if _[2] == 'star':
marker = '*'
elif _[2] == 'HI':
marker = 'o'
else:
marker = 's'
# ax.plot([_[0], _[1]], [_[0]+0.01, _[1]+0.01])
ax.scatter(_[0], _[1]+color_counter*0.1, 30, marker=marker, color=color)
# ax.plot(zip(*a_)[0], zip(*a_)[1],'-', lw=10)
# print zip(*a_)[0]
except Exception:
print 'Exception for ' + name
# ax.plot(romeo_min, romeo_min2, '.', ms=9, color=color)
# ax.plot(romeo_max, romeo_max2, '.', ms=9, color=color)
ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$')
color_counter+=1
ax.legend(loc='lower right')
ax.set_xlim(0, 100)
ax.set_ylim(1.0, 2.5)
ax.set_ylabel(r'$Q_{RF, thick}$', fontsize=20)
ax.set_xlabel(r'$Q_{RF}^{-1}$', fontsize=20)
ax.grid()
plt.show()
Влияние разных дисперсий H2 и HI:
In [50]:
import scipy.interpolate
from math import pi
def romeo_Qinv_thin(r=None, epicycl=None, sound_vel_CO=11., sound_vel_HI=6., 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_CO/(pi*G*CO_density)
Q_HI = kappa*sound_vel_HI/(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_HI, sound_vel_CO]
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]
def plot_RF13_vs_2F(ax=None, r_g_dens=None, HI_gas_dens=None, CO_gas_dens=None, epicycl=None, sound_vel_CO=11., sound_vel_HI=6., 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):
'''Плотности газа передаются не скорр. за гелий.'''
if thin:
romeo_Q = romeo_Qinv_thin
else:
romeo_Q = romeo_Qinv
totgas = zip(r_g_dens, [He_coeff*(l[0]+l[1]) for l in zip(HI_gas_dens, CO_gas_dens)])
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_CO=sound_vel_CO, sound_vel_HI=sound_vel_HI, 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)
romeo_min.append(rom)
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_CO=sound_vel_CO, sound_vel_HI=sound_vel_HI, 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)
romeo_max.append(rom)
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(romeo_min, invQeff_min, 's', color=color, alpha=0.4, ms=8)
# ax.plot(romeo_max, invQeff_max, 'o', color=color, alpha=0.4, ms=8)
return invQeff_min, invQeff_max, romeo_min, romeo_max
kenn_photom={
'n338_I' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60., 15.],
'n338_B' : ['mu0d_Bc', 'h_disc_B', 6.99, 'B', 60., 15.],
'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40., 6.7],
'n2985_K' : ['mudK', 'hK', 1.42, 'K', 70., 14.],
'n2985_S4G' : [('mu0d_s4g', 'mu0d_s4g_2'), ('h_disc_s4g', 'h_disc_s4g_2'), 1.35, 'S4G', 70., 6.3],
'n3898_R' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80., 8.8],
'n3898_R2' : [('mu_inner_approx', 'mu_out'), ('h_approx', 'h_out'), 2.93, 'R', 80., 8.8],
'n4258_I' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150., 15.],
'n4258_S4G' : ['mu0d_36', 'h_disc_36', 0.38, 'S4G', 150., 15.],
'n4725_H' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55., 10.14],
'n4725_S4G' : ['mu0d_s4g', 'h_disc_s4g', 1.61, 'S4G', 55., 10.14],
'n5533_r' : ['MU0_r', 'hi_r', 4.41, 'r', 100., 8.9],
'n5533_R' : ['mu0d_Rc', 'h_disc_R', 7.15, 'R', 100., 9.9]
}
fig = plt.figure(figsize=[8, 8])
color_counter = 0
color_shuffle = range(13)
shuffle(color_shuffle)
ax = plt.gca()
for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
for key in kenn_photom.keys():
if name in key:
print key
phot = kenn_photom[key]
color = cm.rainbow(np.linspace(0, 1, 13))[color_shuffle[color_counter]]
try:
h = locals()[name+'dict'][phot[1][0]]
except KeyError:
h = locals()[name+'dict'][phot[1]]
# print name
# print mud
# print h
if name in ['n4258', 'n4725']:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
elif name == 'n2985':
total_gas_data=locals()[name+'dict']['total_gas_data_']
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in
zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
while total_gas_data[0][0] < phot[5]:
total_gas_data = total_gas_data[1:]
while total_gas_data[0][-1] > 130.:
total_gas_data = total_gas_data[:-1]
# print total_gas_data
ML = phot[2]
band = phot[3]
if type(phot[0]) == tuple and phot[3] == 'S4G':
h1 = globals()[name+'dict'][phot[1][0]]
h2 = globals()[name+'dict'][phot[1][1]]
mu01 = globals()[name+'dict'][phot[0][0]]
mu02 = globals()[name+'dict'][phot[0][1]]
star_density = lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu01, h=h1), ML),
lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu02, h=h2), ML)))(l)
elif type(phot[0]) == tuple:
h1 = globals()[name+'dict'][phot[1][0]]
h2 = globals()[name+'dict'][phot[1][1]]
mu01 = globals()[name+'dict'][phot[0][0]]
mu02 = globals()[name+'dict'][phot[0][1]]
star_density = lambda l: tot_dens((lambda l1: surf_density(mu=mu_disc(l1, mu0=mu01, h=h1), M_to_L=ML, band=band),
lambda l2: surf_density(mu=mu_disc(l2, mu0=mu02, h=h2), M_to_L=ML, band=band)))(l)
elif phot[3] == 'S4G':
star_density = lambda l: s4g_surf_density(mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), ML)
else:
star_density = lambda l: surf_density(mu=mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), M_to_L=ML, band=band)
if name == 'n338':
epicycl = locals()[name+'dict']['epicyclicFreq_real_']
else:
epicycl = locals()[name+'dict']['epicyclicFreq_real']
# invQwff_WS_min, invQwff_WS_max = plot_WS_vs_2f(ax=ax, total_gas_data=total_gas_data,
# epicycl=epicycl,
# gas_approx=locals()[name+'dict']['spl_gas'],
# sound_vel=locals()[name+'dict']['sound_vel'],
# scale=locals()[name+'dict']['scale'],
# sigma_max=locals()[name+'dict']['sig_R_maj_max'],
# sigma_min=locals()[name+'dict']['sig_R_maj_min'],
# star_density_max=star_density,
# star_density_min=star_density,
# data_lim=locals()[name+'dict']['data_lim'],
# color=color, alpha=0.5,
# label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$',
# sfrange=phot[4])
# ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$')
if name in ['n4258', 'n4725']:
r_mol_dens = locals()[name+'dict']['r_mol_dens']
mol_dens = locals()[name+'dict']['mol_dens']
y_interp = scipy.interpolate.interp1d(list(r_mol_dens), list(mol_dens))
def y_interp_(r):
if r <= min(r_mol_dens):
return y_interp(min(r_mol_dens))
elif r < max(r_mol_dens):
return y_interp(r)
else:
return 0.
# не совсем точно так, но все же
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [l[1]/He_coeff - y_interp_(l[0]) for l in zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])],
[y_interp_(l) for l in locals()[name+'dict']['r_g_dens']])[1:15]
elif name == 'n2985':
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'],
[locals()[name+'dict']['y_interp_'](l) for l in locals()[name+'dict']['r_g_dens']])[1:15]
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'],
[locals()[name+'dict']['y_interp_'](l, h) for l in locals()[name+'dict']['r_g_dens']])[1:15]
while total_gas_data[0][0] < phot[5]:
total_gas_data = total_gas_data[1:]
while total_gas_data[0][-1] > 130.:
total_gas_data = total_gas_data[:-1]
invQeff_min, invQeff_max, romeo_min, romeo_max = plot_RF13_vs_2F(ax=ax, r_g_dens=zip(*total_gas_data)[0],
HI_gas_dens=zip(*total_gas_data)[1],
CO_gas_dens=zip(*total_gas_data)[2],
epicycl=epicycl,
sound_vel_HI=6.,
sound_vel_CO=6.,
sigma_R_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_R_min=locals()[name+'dict']['sig_R_maj_min'],
star_density=star_density,
alpha_max=0.7,
alpha_min=0.3,
scale=locals()[name+'dict']['scale'],
gas_approx=locals()[name+'dict']['spl_gas'],
thin=True,
color=color)
invQeff_min, invQeff_max, romeo_min2, romeo_max2 = plot_RF13_vs_2F(ax=ax, r_g_dens=zip(*total_gas_data)[0],
HI_gas_dens=zip(*total_gas_data)[1],
CO_gas_dens=zip(*total_gas_data)[2],
epicycl=epicycl,
sound_vel_HI=6.,
sound_vel_CO=11.,
sigma_R_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_R_min=locals()[name+'dict']['sig_R_maj_min'],
star_density=star_density,
alpha_max=0.7,
alpha_min=0.3,
scale=locals()[name+'dict']['scale'],
gas_approx=locals()[name+'dict']['spl_gas'],
thin=True,
color=color)
ax.plot(romeo_min, romeo_min2, 's', ms=9, color=color)
ax.plot(romeo_max, romeo_max2, 'o', ms=9, color=color)
ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$')
color_counter+=1
ax.legend(loc='lower right')
ax.set_xlim(0, 1.05)
ax.set_ylim(0, 1.05)
ax.set_ylabel(r'$Q_{RF, \sigma_{CO}=11}^{-1}$', fontsize=20)
ax.set_xlabel(r'$Q_{RF}^{-1}$', fontsize=20)
ax.grid()
line, = ax.plot([0., 1*2], [0., 2.0], '--', alpha=0.3, color='gray')
label_line(line, 'y=x', 0.8*1, 0.8, color='black')
plt.show()
Толщины:
In [55]:
import scipy.interpolate
from math import pi
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]
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]
def plot_RF13_vs_2F(ax=None, r_g_dens=None, HI_gas_dens=None, CO_gas_dens=None, epicycl=None, sound_vel=6., 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, verbose=True):
'''Плотности газа передаются не скорр. за гелий.'''
if thin:
romeo_Q = romeo_Qinv_thin
else:
romeo_Q = romeo_Qinv
totgas = zip(r_g_dens, [He_coeff*(l[0]+l[1]) for l in zip(HI_gas_dens, CO_gas_dens)])
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, verbose=verbose)
romeo_min.append(rom)
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, verbose=verbose)
romeo_max.append(rom)
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(romeo_min, invQeff_min, 's', color=color, alpha=0.4, ms=8)
# ax.plot(romeo_max, invQeff_max, 'o', color=color, alpha=0.4, ms=8)
return invQeff_min, invQeff_max, romeo_min, romeo_max
kenn_photom={
'n338_I' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60., 15.],
'n338_B' : ['mu0d_Bc', 'h_disc_B', 6.99, 'B', 60., 15.],
'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40., 6.7],
'n2985_K' : ['mudK', 'hK', 1.42, 'K', 70., 14.],
'n2985_S4G' : [('mu0d_s4g', 'mu0d_s4g_2'), ('h_disc_s4g', 'h_disc_s4g_2'), 1.35, 'S4G', 70., 6.3],
'n3898_R' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80., 8.8],
'n3898_R2' : [('mu_inner_approx', 'mu_out'), ('h_approx', 'h_out'), 2.93, 'R', 80., 8.8],
'n4258_I' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150., 15.],
'n4258_S4G' : ['mu0d_36', 'h_disc_36', 0.38, 'S4G', 150., 15.],
'n4725_H' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55., 10.14],
'n4725_S4G' : ['mu0d_s4g', 'h_disc_s4g', 1.61, 'S4G', 55., 10.14],
'n5533_r' : ['MU0_r', 'hi_r', 4.41, 'r', 100., 8.9],
'n5533_R' : ['mu0d_Rc', 'h_disc_R', 7.15, 'R', 100., 9.9]
}
fig = plt.figure(figsize=[8, 8])
color_counter = 0
color_shuffle = range(13)
shuffle(color_shuffle)
ax = plt.gca()
for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
# for ind, name in enumerate(['n338']):
for key in kenn_photom.keys():
if name in key:
print key
phot = kenn_photom[key]
color = cm.rainbow(np.linspace(0, 1, 13))[color_shuffle[color_counter]]
try:
h = locals()[name+'dict'][phot[1][0]]
except KeyError:
h = locals()[name+'dict'][phot[1]]
# print name
# print mud
# print h
if name in ['n4258', 'n4725']:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
elif name == 'n2985':
total_gas_data=locals()[name+'dict']['total_gas_data_']
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in
zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
while total_gas_data[0][0] < phot[5]:
total_gas_data = total_gas_data[1:]
while total_gas_data[0][-1] > 130.:
total_gas_data = total_gas_data[:-1]
# print total_gas_data
ML = phot[2]
band = phot[3]
if type(phot[0]) == tuple and phot[3] == 'S4G':
h1 = globals()[name+'dict'][phot[1][0]]
h2 = globals()[name+'dict'][phot[1][1]]
mu01 = globals()[name+'dict'][phot[0][0]]
mu02 = globals()[name+'dict'][phot[0][1]]
star_density = lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu01, h=h1), ML),
lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu02, h=h2), ML)))(l)
elif type(phot[0]) == tuple:
h1 = globals()[name+'dict'][phot[1][0]]
h2 = globals()[name+'dict'][phot[1][1]]
mu01 = globals()[name+'dict'][phot[0][0]]
mu02 = globals()[name+'dict'][phot[0][1]]
star_density = lambda l: tot_dens((lambda l1: surf_density(mu=mu_disc(l1, mu0=mu01, h=h1), M_to_L=ML, band=band),
lambda l2: surf_density(mu=mu_disc(l2, mu0=mu02, h=h2), M_to_L=ML, band=band)))(l)
elif phot[3] == 'S4G':
star_density = lambda l: s4g_surf_density(mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), ML)
else:
star_density = lambda l: surf_density(mu=mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), M_to_L=ML, band=band)
if name == 'n338':
epicycl = locals()[name+'dict']['epicyclicFreq_real_']
else:
epicycl = locals()[name+'dict']['epicyclicFreq_real']
# invQwff_WS_min, invQwff_WS_max = plot_WS_vs_2f(ax=ax, total_gas_data=total_gas_data,
# epicycl=epicycl,
# gas_approx=locals()[name+'dict']['spl_gas'],
# sound_vel=locals()[name+'dict']['sound_vel'],
# scale=locals()[name+'dict']['scale'],
# sigma_max=locals()[name+'dict']['sig_R_maj_max'],
# sigma_min=locals()[name+'dict']['sig_R_maj_min'],
# star_density_max=star_density,
# star_density_min=star_density,
# data_lim=locals()[name+'dict']['data_lim'],
# color=color, alpha=0.5,
# label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$',
# sfrange=phot[4])
# ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$')
if name in ['n4258', 'n4725']:
r_mol_dens = locals()[name+'dict']['r_mol_dens']
mol_dens = locals()[name+'dict']['mol_dens']
y_interp = scipy.interpolate.interp1d(list(r_mol_dens), list(mol_dens))
def y_interp_(r):
if r <= min(r_mol_dens):
return y_interp(min(r_mol_dens))
elif r < max(r_mol_dens):
return y_interp(r)
else:
return 0.
# не совсем точно так, но все же
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [l[1]/He_coeff - y_interp_(l[0]) for l in zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])],
[y_interp_(l) for l in locals()[name+'dict']['r_g_dens']])[1:15]
elif name == 'n2985':
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'],
[locals()[name+'dict']['y_interp_'](l) for l in locals()[name+'dict']['r_g_dens']])[1:15]
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'],
[locals()[name+'dict']['y_interp_'](l, h) for l in locals()[name+'dict']['r_g_dens']])[1:15]
while total_gas_data[0][0] < phot[5]:
total_gas_data = total_gas_data[1:]
while total_gas_data[0][-1] > 130.:
total_gas_data = total_gas_data[:-1]
invQeff_min, invQeff_max, romeo_min, romeo_max = plot_RF13_vs_2F(ax=ax, r_g_dens=zip(*total_gas_data)[0],
HI_gas_dens=zip(*total_gas_data)[1],
CO_gas_dens=zip(*total_gas_data)[2],
epicycl=epicycl,
sound_vel=sound_vel,
sigma_R_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_R_min=locals()[name+'dict']['sig_R_maj_min'],
star_density=star_density,
alpha_max=0.7,
alpha_min=0.3,
scale=locals()[name+'dict']['scale'],
gas_approx=locals()[name+'dict']['spl_gas'],
thin=True,
color=color, verbose=False)
invQeff_min, invQeff_max, romeo_min2, romeo_max2 = plot_RF13_vs_2F(ax=ax, r_g_dens=zip(*total_gas_data)[0],
HI_gas_dens=zip(*total_gas_data)[1],
CO_gas_dens=zip(*total_gas_data)[2],
epicycl=epicycl,
sound_vel=sound_vel,
sigma_R_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_R_min=locals()[name+'dict']['sig_R_maj_min'],
star_density=star_density,
alpha_max=0.7,
alpha_min=0.3,
scale=locals()[name+'dict']['scale'],
gas_approx=locals()[name+'dict']['spl_gas'],
thin=False,
color=color, verbose=False)
ax.plot(romeo_min, romeo_min2, '.', ms=9, color=color)
ax.plot(romeo_max, romeo_max2, '.', ms=9, color=color)
ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$')
color_counter+=1
ax.legend(loc='lower right')
ax.set_xlim(0, 1.05)
ax.set_ylim(0, 1.05)
ax.set_ylabel(r'$Q_{RF, thick}$', fontsize=20)
ax.set_xlabel(r'$Q_{RF}^{-1}$', fontsize=20)
ax.grid()
line, = ax.plot([0., 1*2], [0., 2.0], '--', alpha=0.3, color='gray')
label_line(line, 'y=x', 0.8*1, 0.8, color='black')
plt.show()
Хорошо видно проблему, которая была с газом:
In [69]:
%%time
import scipy.interpolate
def plot_WS_vs_2f(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, sfrange=None):
rr = zip(*total_gas_data)[0]
print total_gas_data
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))
# print invQeff_min
invQwff_WS_min = [1./Qg(epicycl=epicycl(gas_approx, r, scale), sound_vel=sound_vel, gas_density=gd) + 1./Qs(
epicycl=epicycl(gas_approx, r, scale), sigma=sigma_max(r), star_density=star_density_min(r)) for r, gd in total_gas_data]
ax.plot(rr, invQeff_min, 'o', color=color, alpha=0.8, ms=10)
# ax.plot(rr, invQwff_WS_min, 's', color=color, alpha=0.8, ms=10)
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))
# print invQeff_max
invQwff_WS_max = [1./Qg(epicycl=epicycl(gas_approx, r, scale), sound_vel=sound_vel, gas_density=gd) + 1./Qs(
epicycl=epicycl(gas_approx, r, scale), sigma=sigma_min(r), star_density=star_density_max(r)) for r, gd in total_gas_data]
# ax.plot(invQwff_WS_max, invQeff_max, 'o', color=color, alpha=0.8, ms=10)
def plot_RF13_vs_2F(ax=None, 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, color=None):
'''Плотности газа передаются не скорр. за гелий.'''
if thin:
romeo_Q = romeo_Qinv_thin
else:
romeo_Q = romeo_Qinv
totgas = zip(r_g_dens, [He_coeff*(l[0]+l[1]) for l in zip(HI_gas_dens, CO_gas_dens)])
print totgas
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)
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))
# print invQeff_min
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)
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))
# print invQeff_max
# ax.plot(r_g_dens, romeo_min, 's', color=color, alpha=0.8, ms=10)
ax.plot(r_g_dens, invQeff_min, 'o', color=color, alpha=0.8, ms=10)
ax.plot(r_g_dens, invQg, '-', color='b', alpha=0.8)
# ax.plot(romeo_min, invQeff_min, 's', color=color, alpha=0.8, ms=10)
# ax.plot(romeo_max, invQeff_max, 'o', color=color, alpha=0.8, ms=10)
markers = {'.': 'point', ',': 'pixel', 'o': 'circle', 'v': 'triangle_down', '^': 'triangle_up', '<': 'triangle_left', '>': 'triangle_right', '1': 'tri_down', '2':
'tri_up', '3': 'tri_left', '4': 'tri_right', '8': 'octagon', 's': 'square', 'p': 'pentagon', '*': 'star', 'h': 'hexagon1', 'H': 'hexagon2', '+': 'plus',
'x': 'x', 'D': 'diamond', 'd': 'thin_diamond', '|': 'vline', '_': 'hline', 'P': 'plus_filled', 'X': 'x_filled', 0: 'tickleft', 1: 'tickright', 2: 'tickup',
3: 'tickdown', 4: 'caretleft', 5: 'caretright', 6: 'caretup', 7: 'caretdown', 8: 'caretleftbase', 9: 'caretrightbase', 10: 'caretupbase', 11: 'caretdownbase'}
kenn_photom={'n338_I' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60., 15.],
'n338_B' : ['mu0d_Bc', 'h_disc_B', 6.99, 'B', 60., 15.],
'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40., 6.7],
'n2985_K' : ['mudK', 'hK', 1.42, 'K', 70., 14.],
'n2985_S4G' : [('mu0d_s4g', 'mu0d_s4g_2'), ('h_disc_s4g', 'h_disc_s4g_2'), 1.35, 'S4G', 70., 6.3],
'n3898_R' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80., 8.8],
'n3898_R2' : [('mu_inner_approx', 'mu_out'), ('h_approx', 'h_out'), 2.93, 'R', 80., 8.8],
'n4258_I' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150., 15.],
# 'n4258_S4G' : ['mu0d_36', 'h_disc_36', 0.38, 'S4G', 150., 15.],
'n4725_H' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55., 10.14],
'n4725_S4G' : ['mu0d_s4g', 'h_disc_s4g', 1.61, 'S4G', 55., 10.14],
'n5533_r' : ['MU0_r', 'hi_r', 4.41, 'r', 100., 8.9],
'n5533_R' : ['mu0d_Rc', 'h_disc_R', 7.15, 'R', 100., 9.9]}
fig, (ax, ax2) = plt.subplots(ncols=2, nrows=1, figsize=[16, 8])
color_counter = 0
color_shuffle = range(13)
shuffle(color_shuffle)
# for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
for ind, name in enumerate(['n4258']):
for key in kenn_photom.keys():
if name in key:
print key
phot = kenn_photom[key]
color = cm.rainbow(np.linspace(0, 1, 13))[color_shuffle[color_counter]]
try:
h = locals()[name+'dict'][phot[1][0]]
except KeyError:
h = locals()[name+'dict'][phot[1]]
# print name
# print mud
# print h
if name in ['n4258', 'n4725']:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
elif name == 'n2985':
total_gas_data=locals()[name+'dict']['total_gas_data_']
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in
zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
while total_gas_data[0][0] < phot[5]:
total_gas_data = total_gas_data[1:]
# print total_gas_data
ML = phot[2]
band = phot[3]
if type(phot[0]) == tuple and phot[3] == 'S4G':
h1 = globals()[name+'dict'][phot[1][0]]
h2 = globals()[name+'dict'][phot[1][1]]
mu01 = globals()[name+'dict'][phot[0][0]]
mu02 = globals()[name+'dict'][phot[0][1]]
star_density = lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu01, h=h1), ML),
lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu02, h=h2), ML)))(l)
elif type(phot[0]) == tuple:
h1 = globals()[name+'dict'][phot[1][0]]
h2 = globals()[name+'dict'][phot[1][1]]
mu01 = globals()[name+'dict'][phot[0][0]]
mu02 = globals()[name+'dict'][phot[0][1]]
star_density = lambda l: tot_dens((lambda l1: surf_density(mu=mu_disc(l1, mu0=mu01, h=h1), M_to_L=ML, band=band),
lambda l2: surf_density(mu=mu_disc(l2, mu0=mu02, h=h2), M_to_L=ML, band=band)))(l)
elif phot[3] == 'S4G':
star_density = lambda l: s4g_surf_density(mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), ML)
else:
star_density = lambda l: surf_density(mu=mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), M_to_L=ML, band=band)
if name == 'n338':
epicycl = locals()[name+'dict']['epicyclicFreq_real_']
else:
epicycl = locals()[name+'dict']['epicyclicFreq_real']
plot_WS_vs_2f(ax=ax, total_gas_data=total_gas_data,
epicycl=epicycl,
gas_approx=locals()[name+'dict']['spl_gas'],
sound_vel=locals()[name+'dict']['sound_vel'],
scale=locals()[name+'dict']['scale'],
sigma_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_min=locals()[name+'dict']['sig_R_maj_min'],
star_density_max=star_density,
star_density_min=star_density,
data_lim=locals()[name+'dict']['data_lim'],
color=color, alpha=0.5,
label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$',
sfrange=phot[4])
ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$')
if name in ['n4258', 'n4725']:
r_mol_dens = locals()[name+'dict']['r_mol_dens']
mol_dens = locals()[name+'dict']['mol_dens']
y_interp = scipy.interpolate.interp1d(list(r_mol_dens), list(mol_dens))
def y_interp_(r):
if r <= min(r_mol_dens):
return y_interp(min(r_mol_dens))
elif r < max(r_mol_dens):
return y_interp(r)
else:
return 0.
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['HI_dens'],
[y_interp_(l) for l in locals()[name+'dict']['r_g_dens']])[1:15]
elif name == 'n2985':
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'],
[locals()[name+'dict']['y_interp_'](l) for l in locals()[name+'dict']['r_g_dens']])[1:15]
else:
total_gas_data=zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'],
[locals()[name+'dict']['y_interp_'](l, h) for l in locals()[name+'dict']['r_g_dens']])[1:15]
while total_gas_data[0][0] < phot[5]:
total_gas_data = total_gas_data[1:]
plot_RF13_vs_2F(ax=ax, r_g_dens=zip(*total_gas_data)[0],
HI_gas_dens=zip(*total_gas_data)[1],
CO_gas_dens=zip(*total_gas_data)[2],
epicycl=epicycl,
sound_vel=6.,
sigma_R_max=locals()[name+'dict']['sig_R_maj_max'],
sigma_R_min=locals()[name+'dict']['sig_R_maj_min'],
star_density=star_density,
alpha_max=0.7,
alpha_min=0.3,
scale=locals()[name+'dict']['scale'],
gas_approx=locals()[name+'dict']['spl_gas'],
thin=True,
color=color)
ax2.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$')
color_counter+=1
for _ in np.linspace(1., 2., 5):
ax.plot([0., _], [0., 1.0], '--', alpha=0.3, color='gray')
ax2.plot([0., 1.0], [0., _], '--', alpha=0.3, color='gray')
# ax.legend(loc='upper left')
# ax.set_xlim(0, 2.05)
ax.set_ylim(0, 1.05)
# ax.set_ylabel(r'$Q_{eff}^{-1}$', fontsize=20)
# ax.set_xlabel(r'$Q_{WS}^{-1}$', fontsize=20)
ax.grid()
# ax2.legend(loc='lower right')
# ax2.set_xlim(0, 1.05)
ax2.set_ylim(0, 1.05)
# ax2.set_ylabel(r'$Q_{eff}^{-1}$', fontsize=20)
# ax2.set_xlabel(r'$Q_{RF}^{-1}$', fontsize=20)
ax2.grid()
plt.show()
In [ ]:
In [ ]: