Галактика из диплома.
In [37]:
from IPython.display import HTML
from IPython.display import Image
import os
%pylab
%matplotlib inline
%run ../../utils/load_notebook.py
In [38]:
from photometry import *
In [39]:
from instabilities import *
In [40]:
from utils import *
In [41]:
name = 'N3898'
gtype = 'SA(s)ab'
incl = 61. #Mean value from amny sources
distance = 18.9 #Noordermeer
scale = 0.092 #kpc/arcsec according to Noordermeer and the same in ApJ 142 145(31pp) 2011
data_path = '../../data/ngc3898'
sin_i, cos_i = np.sin(incl*np.pi/180.), np.cos(incl*np.pi/180.)
In [42]:
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')
In [43]:
os.chdir(data_path)
# Данные из NED
HTML('<iframe src=http://ned.ipac.caltech.edu/cgi-bin/objsearch?objname=ngc+3898&extend=no&hconst=\
73&omegam=0.27&omegav=0.73&corr_z=1&out_csys=Equatorial&out_equinox=J2000.0&obj_sort=RA+or+Longitude&of=pre_text&zv_breaker=\
30000.0&list_limit=5&img_stamp=YES width=1000 height=350></iframe>')
Out[43]:
In [44]:
# Данные из HYPERLEDA
HTML('<iframe src=http://leda.univ-lyon1.fr/ledacat.cgi?o=ngc3898 width=1000 height=350></iframe>')
Out[44]:
In [45]:
#SDSS
Image('ngc3898_SDSS.jpeg', width=300)
Out[45]:
Картинка из выборки http://cosmo.nyu.edu/hogg/rc3/ и там есть с маштабом изображение:
In [46]:
Image('ngc3898_SDSS_labeled.jpeg', width=500)
Out[46]:
In [47]:
#JHK
Image('ngc3898_JHK.jpg', width=300)
Out[47]:
Есть пара картинок от HST, например вот, но не очень хороших.
Noordermeer, thesis, p. 113
UGC 6787 (NGC 3898) suffers from the same problem as UGC 6787: the luminous bulge distorts the isophotes of the disk out to large radii and the standard method to derive the inclination (section 3.5.1) cannot be used. Only after the model image of the bulge was subtracted could the disk ellipticity be determined. The final disk photometric profile and the intrinsic axis ratio of the bulge were derived using an inclination of $61^\circ$ . Note that this is still somewhat lower than the kinematic inclination derived from the HI velocity field (chapter 4)
Noordermeer, thesis, pp. 188-189
The rotation curve of UGC 6787 (NGC 3898) is well resolved and shows some characteristic ‘wiggles’ with an amplitude of 30 – 50 km/s. The kinematics in the central parts are only barely resolved in the optical spectrum, and due to the high inclination angle and resulting line-of-sight integration effects, the central line-profiles are strongly broadened. After the initial rise of the rotation curve to the peak velocity of 270 km/s, the rotation velocities drop to approximately 220 km/s at a radius of 3000 (∼ 2.75 kpc), after which they gradually rise again to 250 km/s at R ≈ 10000 (∼9 kpc). The rotation velocities then drop again to 220 km/s, after which they rise again to reach a more or less flat plateau at 250 km/s. Although there are clear indications that the gas disk of UGC 6787 is warped, the locations of the ‘wiggles’ in the rotation curve do not coincide with the radii where the position angle and the inclination change and the variations in the rotation velocity seem to be real. This is further confirmed by the fact that the variations occur symmetrically at all position angles over the velocity field. The discrepancy between the kinematic inclination angle derived here and the optical inclination from chapter 3 can be explained by the dominance of the bulge in the optical image. As was noted in chapter 3, the optical image of this galaxy is dominated by the spheroidal bulge out to large radii, which makes it impossible to obtain a reliable estimate for the inclination. The kinematical inclination derived here is free of such effects and thus reflects the true orientation of this galaxy more accurately.
In [48]:
Image('noord_pics/HI_data.png')
Out[48]:
In [49]:
Image('noord_pics/photometry.png')
Out[49]:
In [50]:
Image('noord_pics/rot_vel_all.png')
Out[50]:
In [51]:
#точки, снятые с предыдущей картинки с помощью WebPlotDigitizer
r_ma, vel_ma = zip(*np.loadtxt("v_gas_from_webplotdigitizer.dat", float, delimiter=','))
plt.plot(r_ma, vel_ma, '--')
plt.xlim(0, 100)
plt.ylim(0, 500);
In [52]:
Image('noord_pics/rot_vel_HI.png')
Out[52]:
In [53]:
Image('noord_pics/rot_vel_center.png', width=300)
Out[53]:
In [54]:
Image('noord_pics/rot_vel_large.png', width=300)
Out[54]:
Снимем отсюда точки - уже по наблюдательным данным (отдельно в центре и отдельно на периферии):
In [55]:
r_ma2, vel_ma2 = zip(*np.loadtxt("v_gas_from_webplotdigitizer2.dat", float, delimiter=','))
plt.plot(r_ma2, vel_ma2, '--')
plt.plot(map(lambda l: l/scale, r_ma), vel_ma, '--')
plt.xlim(0, 400)
plt.ylim(0, 320);
Видно, что снял я достаточно точно (с двух картинок данные совпали).
In [56]:
Image('noord_pics/isothermal_halo.png')
Out[56]:
In [57]:
Image('noord_pics/NFW.png')
Out[57]:
Ценное - одножидкостный критерий на сравнение:
In [58]:
Image('noord_pics/instab.png')
Out[58]:
In [59]:
r_q_n, q_n = zip(*np.loadtxt("noord_pics/Q1F.dat", float, delimiter=','))
fig = plt.figure(figsize=[8,2])
plt.plot(r_q_n, q_n, 'o-')
plt.xlim(0, 30)
plt.ylim(0, 1.2);
Наконец последнее - все табличные данные из диссертации:
In [60]:
Image('noord_pics/tables.png')
Out[60]:
TODO: подчистить данные из Noordermeera, убрать ненужное
TODO: добавить pignatelli
In [61]:
# Данные по звездной кинематике Noordermeer вдоль большей полуоси (не исправленные за наклон?)
r_ma3, vel_ma3, e_vel_ma3, sig_ma3, e_sig_ma3 = zip(*np.loadtxt("v_stars_ma.dat", float))
# Данные по звездной кинематике Heraudeau+1999 вдоль большой полуоси (не исправленные за наклон?)
r_ma4, vel_ma4, e_vel_ma4, sig_ma4, e_sig_ma4 = zip(*np.loadtxt("v_stars_HSM.dat", float))
fig = plt.figure(figsize=[10,8])
plt.errorbar(r_ma3, vel_ma3, e_vel_ma3, fmt='-.', marker='.', mew=0, label="Noord stars maj")
plt.errorbar(r_ma4, vel_ma4, e_vel_ma4, fmt='-.', marker='.', mew=0, label="Heraudeau stars maj")
plt.legend()
plt.show()
TODO: проверить, как хорошо были сняты звездные данные (обе кривые)
Определенно надо было минусы проставить в Heraudeau. Перегнем и сравним со снятыми данными:
In [62]:
r_ma3_, vel_ma3_, e_vel_ma3_ = zip(*sorted(zip(np.abs(r_ma3), np.abs(vel_ma3), e_vel_ma3)))
r_ma4_, vel_ma4_, e_vel_ma4_ = zip(*sorted(zip(np.abs(r_ma4), np.abs(vel_ma4), e_vel_ma4)))
In [63]:
fig = plt.figure(figsize=[10,8])
plt.errorbar(r_ma3_, vel_ma3_, e_vel_ma3_, fmt='-.', marker='.', mew=0, label="Noord stars maj")
plt.errorbar(r_ma4_, vel_ma4_, e_vel_ma4_, fmt='-.', marker='.', mew=0, label="Heraudeau stars maj")
plt.plot(r_ma2, vel_ma2, '--')
plt.plot(map(lambda l: l/scale, r_ma), vel_ma, '--')
plt.xlim(0, 400)
plt.ylim(0, 320)
plt.legend()
plt.xlabel('R, arcsec')
plt.show()
Видно, что в статье 2008 года кривая протянута приемлимо далеко, однако газ имеет сильно большие значения, чем звездная кривая вращения.
In [64]:
incl_corr = lambda l: l/sin_i
r_ma_b, vel_ma_b, e_vel_b = r_ma3_, map(incl_corr, vel_ma3_), map(incl_corr, e_vel_ma3_)
In [65]:
fig = plt.figure(figsize=[10,8])
plt.errorbar(r_ma_b, vel_ma_b, yerr=e_vel_b, fmt='.', marker='.', mew=0, color='blue', label = 'Noord star maj')
test_points = np.linspace(0.0, max(r_ma_b), 100)
poly_star = poly1d(polyfit(r_ma_b, vel_ma_b, deg=3))
plt.plot(test_points, poly_star(test_points), '-', label='poly deg=3')
def w(arr):
return map(lambda l: 1/(1. + l**2), arr)
import scipy.interpolate as inter
spl = inter.UnivariateSpline(r_ma_b[::2], vel_ma_b[::2], k=3, s=500.)
plt.plot(test_points, spl(test_points), '-', label='spl s=500')
# spl1 = inter.UnivariateSpline(r_ma_b[::2], vel_ma_b[::2], k=3, s=10000., w=w(e_vel_b))
# plt.plot(test_points, spl1(test_points), '-', label='spl s=500 w^2')
plt.legend(loc='lower right')
plt.ylim(0, 250)
plt.show()
C весами плохо получается, берем обычный достаточно гладкий сплайн.
In [66]:
star_approx = spl
Дисперсии из Pignatelli http://mnras.oxfordjournals.org/content/323/1/188.full.pdf:
In [67]:
Image('pignatelli_2001_kinem.png', width=400)
Out[67]:
In [68]:
r_p, v_p, dv_p, sig_p, dsig_p, _, _, _, _ = zip(*np.loadtxt("pignatelli_2001_kinem.dat", float))
fig = plt.figure(figsize=[8, 3])
plt.errorbar(r_p, sig_p, yerr=dsig_p, fmt='.', marker='.', mew=0, color='blue', label=r'$\sigma_{los}^{maj} Pignatelli$')
plt.xlabel('$R$')
plt.ylabel('$\sigma$')
plt.xlim(-40., 40.)
plt.ylim(0, 320)
plt.legend();
In [69]:
# Исправляем значения вдоль малой оси на синус угла:
def correct_min(R):
return R / cos_i
r_sig_ma, sig_ma, e_sig_ma = zip(*np.loadtxt("s_stars_maN.dat", float))
r_sig_mi, sig_mi, e_sig_mi = zip(*np.loadtxt("s_stars_miN.dat", float))
fig = plt.figure(figsize=[8, 8])
plt.errorbar(r_sig_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='.', mew=0, color='blue', label=r'$\sigma_{los}^{maj} Noord$')
plt.errorbar(map(abs, r_ma4), sig_ma4, yerr=e_sig_ma4, fmt='.', marker='.', mew=0, color='red', label=r'$maj Heraudeau $')
plt.errorbar(map(abs, r_p), sig_p, yerr=dsig_p, fmt='.', marker='.', mew=0, color='m', label=r'$\sigma_{los}^{maj} Pignatelli$')
r_sig_mi = map(correct_min, r_sig_mi)
plt.errorbar(r_sig_mi, sig_mi, yerr=e_sig_mi, fmt='.', marker='.', mew=0, color='black', label='$\sigma_{los}^{min} Noord$')
plt.xlabel('$R$')
plt.ylabel('$\sigma$')
plt.legend();
TODO: с малой осью явно что-то не то
In [70]:
spl_min = inter.UnivariateSpline(r_sig_mi, sig_mi, k=3, s=10000.)
sig_min_lim = max(r_sig_mi)
spl_maj = inter.UnivariateSpline(r_sig_ma, sig_ma, k=3, s=100.)
sig_maj_lim = max(r_sig_ma)
points = np.linspace(0.1, max(r_sig_ma)+15., 100)
In [71]:
@flat_end(sig_maj_lim)
def sig_R_maj_minmin(r, spl_maj=spl_maj):
return spl_maj(r).item()
@flat_end(sig_maj_lim)
def sig_R_maj_min(r, spl_maj=spl_maj):
return spl_maj(r).item()/sqrt(sin_i**2 + 0.49*cos_i**2)
@flat_end(sig_maj_lim)
def sig_R_maj_max(r, spl_maj=spl_maj):
return spl_maj(r).item()/sqrt(0.5*sin_i**2 + 0.09*cos_i**2)
@flat_end(sig_maj_lim)
def sig_R_maj_maxmax(r, spl_maj=spl_maj):
return spl_maj(r)*sqrt(2)/sin_i
@flat_end(sig_maj_lim)
def sig_R_maj_maxmaxtrue(r, spl_maj=spl_maj):
return spl_maj(r)/sin_i/sqrt(sigPhi_to_sigR_real(r))
@flat_end(sig_min_lim)
def sig_R_minor_minmin(r, spl_min=spl_min):
return spl_min(r).item()
@flat_end(sig_min_lim)
def sig_R_minor_min(r, spl_min=spl_min):
return spl_min(r).item()/sqrt(sin_i**2 + 0.49*cos_i**2)
@flat_end(sig_min_lim)
def sig_R_minor_max(r, spl_min=spl_min):
return spl_min(r).item()/sqrt(sin_i**2 + 0.09*cos_i**2)
@flat_end(sig_min_lim)
def sig_R_minor_maxmax(r, spl_min=spl_min):
return spl_min(r)/sin_i
Для малой оси:
In [72]:
plt.errorbar(r_sig_mi, sig_mi, yerr=e_sig_mi, fmt='.', marker='.', mew=0, color='red', label='$\sigma_{los}^{min}$')
plt.plot(points, map(sig_R_minor_minmin, points), label = 'minmin')
plt.plot(points, map(sig_R_minor_min, points), label = 'min')
plt.plot(points, map(sig_R_minor_max, points), label = 'max')
plt.plot(points, map(sig_R_minor_maxmax, points), label = 'maxmax')
plt.legend()
plt.ylim(0,400)
plt.xlim(0,100);
Используем соотношение $\sigma_{\varphi}^{2}/\sigma_{R}^{2}$, которое описывается уравнением ${\displaystyle \sigma_{\varphi}^{2}/\sigma_{R}^{2}=0.5\left(1+\frac{R}{\bar{v}_{\varphi}}\frac{d\bar{v}_{\varphi}}{dR}\right)}$ (Binney & Tremaine, 1987)
In [73]:
def sigPhi_to_sigR_real(R):
return 0.5 * (1 + R*star_approx.derivative()(R) / star_approx(R))
plt.plot(test_points, [sigPhi_to_sigR_real(R) for R in test_points], 'd-', color='green')
plt.axhline(y=0.5)
plt.axhline(y=0.0)
plt.xlabel('$R$')
plt.ylabel(r"$\sigma_{\varphi}^2/\sigma_{R}^2$")
plt.ylim(0);
По большой:
In [74]:
plt.errorbar(r_sig_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='.', mew=0, color='red', label='$\sigma_{los}^{maj}$')
plt.plot(points, map(sig_R_maj_minmin, points), label = 'minmin')
plt.plot(points, map(sig_R_maj_min, points), label = 'min')
plt.plot(points, map(sig_R_maj_max, points), label = 'max')
plt.plot(points, map(sig_R_maj_maxmax, points), label = 'maxmax')
plt.plot(points, map(sig_R_maj_maxmaxtrue, points), label = 'maxmaxtrue')
plt.legend()
plt.ylim(0,400)
plt.xlim(0,100);
Сравним major vs minor оценки:
In [75]:
fig = plt.figure(figsize=[10, 7])
plt.plot(points, map(sig_R_maj_minmin, points), label = 'maj minmin')
plt.plot(points, map(sig_R_maj_min, points), label = 'maj min')
plt.plot(points, map(sig_R_maj_max, points), label = 'maj max')
plt.plot(points, map(sig_R_maj_maxmax, points), label = 'maj maxmax')
plt.plot(points, map(sig_R_maj_maxmaxtrue, points), label = 'maj maxmaxtrue')
plt.plot(points, map(sig_R_minor_minmin, points), '--', label = 'minor minmin')
plt.plot(points, map(sig_R_minor_min, points), '--', label = 'minor min')
plt.plot(points, map(sig_R_minor_max, points), '--', label = 'minor max')
plt.plot(points, map(sig_R_minor_maxmax, points), '--', label = 'minor maxmax')
plt.legend()
plt.ylim(50,300)
plt.xlim(0,100);
"Настоящая" оценка сверху на этот раз имеет какую-то особенность.
Газовая кривая, нужна для эпициклического приближения:
In [76]:
fig = plt.figure(figsize=[10,8])
# Данные по кинематике WSRT
r_wsrt, vel_wsrt, e_vel_wsrt = zip(*np.loadtxt("v_gas_WSRT.dat", float))
# Данные по кинематике газа Epinat+2008 в Halpha
r_ha, _,_,_, vel_ha, e_vel_ha, _,_ = zip(*np.loadtxt("v_gasHa.dat", str))
r_ha, vel_ha, e_vel_ha = np.array(r_ha, dtype='float'), np.array(vel_ha, dtype='float'), np.array(e_vel_ha, dtype='float')
plt.plot(r_ma2, vel_ma2, '--', label='Noord thesis #1')
plt.plot(map(lambda l: l/scale, r_ma), vel_ma, '--', label='Noord thesis #2')
plt.errorbar(r_wsrt, vel_wsrt, yerr=e_vel_wsrt, fmt='.', marker='.', mew=0, label = 'WSRT')
plt.errorbar(r_ha, vel_ha, yerr=e_vel_ha, fmt='.', marker='d', mew=0, label = 'Halpha')
plt.ylim(0, 300)
plt.legend(loc='lower right');
TODO: проверить, не надо ли исправлять за угол
Согласие достаточно хорошее, приблизим:
In [77]:
fig = plt.figure(figsize=[10,6])
gas_approx = poly1d(polyfit(r_wsrt, vel_wsrt, deg=15))
test_points = np.linspace(min(r_wsrt), max(r_wsrt), 100)
plt.plot(test_points, gas_approx(test_points), '--', label='poly approx')
spl_gas = inter.UnivariateSpline(r_wsrt, vel_wsrt, k=3, s=2400.)
plt.plot(test_points, spl_gas(test_points), '-', label='spline')
plt.errorbar(r_wsrt, vel_wsrt, yerr=e_vel_wsrt, fmt='.', marker='.', mew=0, label = 'WSRT')
plt.ylim(0, 300)
plt.legend(loc='lower right');
Для случая бесконечного тонкого диска: $$\kappa=\frac{3}{R}\frac{d\Phi}{dR}+\frac{d^2\Phi}{dR^2}$$ где $\Phi$ - гравпотенциал, однако его знать не надо, т.к. есть проще формула: $$\kappa=\sqrt{2}\frac{\vartheta_c}{R}\sqrt{1+\frac{R}{\vartheta_c}\frac{d\vartheta_c}{dR}}$$
In [78]:
fig = plt.figure(figsize=[12, 8])
plt.plot(test_points, [epicyclicFreq_real(gas_approx, x, scale) for x in test_points], '-', label='poly')
plt.plot(test_points, [epicyclicFreq_real(spl_gas, x, scale) for x in test_points], '-', label='spline')
plt.xlabel('$R, arcsec$')
plt.ylabel('$\kappa,\, km/s/kpc$', fontsize=15)
plt.ylim(0, 1000)
plt.legend()
plt.show()
Плотность HI:
In [79]:
r_g_dens, gas_dens = zip(*np.loadtxt("gas_density.dat", float))
plt.plot(r_g_dens, gas_dens, 's-');
У Ноордермеера масса атомарного газа $3.96\times10^9$, и у нас похоже
In [80]:
import scipy.integrate as integrate
import scipy.interpolate
tmp_ = scipy.interpolate.interp1d(r_g_dens, gas_dens)
result = integrate.quad(lambda l: 2*np.pi*l*tmp_(l), r_g_dens[0], r_g_dens[-1])
print (scale * 1000.)**2 * result[0]/1e9
В этой давней работе http://adsabs.harvard.edu/abs/1987ApJ...322...64S есть указание что плотность постоянна, равна $CO$ 1.74 на протяжении примерно 30 секунд от центра.
В этой работе http://adsabs.harvard.edu/abs/2014A%26A...564A..65B (из трех частей) вроде как есть в последней части точка, в первой работе написано что она измерялась, в таблице интенсивности для нее не даны и надо вычислять по формуле $I(CO) = 5\sigma(W_{HI}\delta V_{CO})^{1/2}K km\, s−1$
Насколько я нашел все параметры, интенсивность $I(CO)$ получается
In [81]:
5*2.00*sqrt(463*15.7)/1000. #почему делим на 1000?
Out[81]:
А центральные плотности
In [82]:
5*2.00*sqrt(463*15.7)/1000. * 3.2
Out[82]:
По порядкам адекватно.
Однако все не так простто - опять же врядли это центральная точка. Массы в работе оценены как $logM_{HI} = 9.49$, $logM_{H_2} = 8.28$ (там две массы для разных фаторов $X_{CO}$, тут первая для постоянного фактора, вторая - 8.08)
In [160]:
np.power(10., 8.28)/1e9, np.power(10., 8.08)/1e9
Out[160]:
В работе $X_{CO}=2.3$ для этой массы и расстояние 16.73. Исправим:
In [84]:
(distance/16.73)**2 *(X_CO/2.3) * np.power(10., 8.28)/1e9
Out[84]:
Вообще-то по порядкам схоже и с Ноордермеером согласуется. Соответственно у нас были бы центральные плотности:
In [85]:
0.20*1e9/(2*math.pi*36.2**2 * (scale * 1000.)**2), 0.20*1e9/(2*math.pi*20.2**2 * (scale * 1000.)**2)
Out[85]:
Разницы прям сильно большой не должно быть.
TODO: добавить
In [86]:
all_photometry = []
Фотометрий четыре - одна в J, две в R и одна в B.
In [87]:
Image('diplom_cite_p38.png')
Out[87]:
In [89]:
# from wand.image import Image as WImage
# img = WImage(filename='ngc3898.pdf', resolution=200)
# img[:, 150:1800]
In [90]:
# img = WImage(filename='ngc3898.pdf[1]', resolution=200)
# img[:, 150:1300]
На этой странице ошибка - фотометрия не в J, она в R!
In [91]:
Image('noord_pics/photom_table.png')
Out[91]:
In [92]:
Image('photom_table2.png')
Out[92]:
In [93]:
Image('noord_pics/n3898_photom.png')
Out[93]:
In [94]:
r_phot, mu_phot = zip(*np.loadtxt("noord_pics/n3898_noord_photoRB.dat", float, delimiter=','))
plt.plot(r_phot, mu_phot, 's')
plt.xlim(-10, 250)
plt.ylim(30, 15);
Посмотрим, как ведет себя ноордермееровская фотометрия в $R$:
In [95]:
# for R band
r_eff_R = 8.8
mu_eff_R = 18.37
n_R = 2.3
mu0d_R = 20.49
h_disc_R = 36.2
In [96]:
p_ = np.arange(0.1, 220., 0.1)
bulge = [mu_bulge(l, mu_eff=mu_eff_R, r_eff=r_eff_R, n=n_R) for l in p_]
disc = [mu_disc(l, mu0=mu0d_R, h=h_disc_R) for l in p_]
total_profile = total_mu_profile(bulge, disc)
fig = plt.figure(figsize=[8, 5])
plt.plot(r_phot, mu_phot, 's')
plt.plot(p_, bulge, '-', label='bulge', color='red')
plt.plot(p_, disc, '-', label='disk', color='green')
plt.plot(p_, total_profile, '-', label='sum', color='m')
plt.xlim(-10, 175)
plt.ylim(30, 15)
plt.legend();
В $B$:
In [97]:
# for B band
r_eff_B = 8.8
mu_eff_B = 19.80
n_B = 2.3
mu0d_B = 22.0
h_disc_B = 42.9
In [98]:
bulge_B = [mu_bulge(l, mu_eff=mu_eff_B, r_eff=r_eff_B, n=n_B) for l in p_]
disc_B = [mu_disc(l, mu0=mu0d_B, h=h_disc_B) for l in p_]
total_profile_B = total_mu_profile(bulge_B, disc_B)
fig = plt.figure(figsize=[8, 5])
plt.plot(r_phot, mu_phot, 's')
plt.plot(p_, bulge_B, '-', label='bulge', color='red')
plt.plot(p_, disc_B, '-', label='disk', color='green')
plt.plot(p_, total_profile_B, '-', label='sum', color='m')
plt.xlim(-10, 175)
plt.ylim(30, 15)
plt.legend();
In [99]:
# TODO: для калибровки нужны светимости диска видимо?
M_R = -20.74
M_B = -19.55 # а у Гутиереза -20.5, что достаточно критично
M_to_L_R = bell_mass_to_light(M_B - M_R, 'R', 'B-R')
M_to_L_R
Out[99]:
In [100]:
# TODO: почему балдж не учитывается?
surf_R = [surf_density(mu=m, M_to_L=M_to_L_R, band='R') for m in disc]
In [101]:
M_to_L_B = bell_mass_to_light(M_B - M_R, 'B', 'B-R')
surf_B = [surf_density(mu=m, M_to_L=M_to_L_B, band='B') for m in disc_B]
plt.plot(p_, surf_R, '-', label='Noord R [M/L={:2.2f}]'.format(M_to_L_R))
plt.plot(p_, surf_B, '-', label='Noord B [M/L={:2.2f}]'.format(M_to_L_B))
plt.legend()
Out[101]:
In [102]:
all_photometry.append(('Noorder R', r_eff_R, mu_eff_R, n_R, mu0d_R, h_disc_R, M_to_L_R,
lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_R, h=h_disc_R), M_to_L=M_to_L_R, band='R')))
all_photometry.append(('Noorder B', r_eff_B, mu_eff_B, n_B, mu0d_B, h_disc_B, M_to_L_B,
lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_B, h=h_disc_B), M_to_L=M_to_L_B, band='B')))
$J$ из Mendez-Abreu:
In [103]:
# for J band
r_eff_J = 11.9
mu_eff_J = 18.13
n_J = 3.75
mu0d_J = 19.07
h_disc_J = 29.2
In [104]:
disc_J = [mu_disc(l, mu0=mu0d_J, h=h_disc_J) for l in p_]
M_to_L_J = bell_mass_to_light(M_B - M_R, 'J', 'B-R')
surf_J = [surf_density(mu=m, M_to_L=M_to_L_J, band='J') for m in disc_J]
plt.plot(p_, surf_R, '-', label='Noord R [M/L={:2.2f}]'.format(M_to_L_R))
plt.plot(p_, surf_B, '-', label='Noord B [M/L={:2.2f}]'.format(M_to_L_B))
plt.plot(p_, surf_J, '-', label='Mendez-Abreu J [M/L={:2.2f}]'.format(M_to_L_J))
plt.legend()
Out[104]:
Тоже весьма похоже.
In [105]:
all_photometry.append(('Mendez-Abreu J', r_eff_J, mu_eff_J, n_J, mu0d_J, h_disc_J, M_to_L_J,
lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_J, h=h_disc_J), M_to_L=M_to_L_J, band='J')))
Последнее - Gutierrez в R:
In [106]:
Image('gutierrez_cite_p24.png')
Out[106]:
Тут два диска, видимо надо брать оба:
In [107]:
# данные из базы в Vizier
_, _, r_phot2, mu_phot2 = zip(*np.loadtxt("gutierrez_photomR.dat", str, delimiter='\t'))
r_phot2, mu_phot2 = map(float, r_phot2), map(float, mu_phot2)
plt.plot(r_phot2, mu_phot2, 's')
plt.ylim(28, 15);
In [108]:
h_inner = 30.
h_out = 59.9
mu0_inner = 19.54
mu0_out = 21.53
In [109]:
Image('gutierrez_photomR.png' , width=500)
Out[109]:
In [110]:
fig = plt.figure(figsize=[5, 5])
plt.plot(p_, [mu_disc(l, mu0=mu0_inner, h=h_inner) for l in p_], '--', label='in', color='#007700')
plt.plot(p_, [mu_disc(l, mu0=mu0_out, h=h_out) for l in p_], '--', label='out', color='#000077')
plt.plot(p_, total_mu_profile([mu_disc(l, mu0=mu0_inner, h=h_inner) for l in p_],
[mu_disc(l, mu0=mu0_out, h=h_out) for l in p_]), '-', label='sum', color='#770000')
plt.plot(r_phot2, mu_phot2, '.', ms=2)
plt.xlim(0., 400)
plt.ylim(28, 15)
plt.legend()
plt.axhline(y=23.3, ls='--', alpha=0.2)
plt.axvline(x = 111., linestyle='-.', alpha=0.2);
In [111]:
disc_1 = [mu_disc(l, mu0=mu0_inner, h=h_inner) for l in p_]
disc_2 = [mu_disc(l, mu0=mu0_out, h=h_out) for l in p_]
total_profile_2 = total_mu_profile(disc_1, disc_2)
fig = plt.figure(figsize=[8, 5])
plt.plot(r_phot2, mu_phot2, '.')
plt.plot(p_, disc_1, '-', label='disk inner', color='green')
plt.plot(p_, disc_2, '-', label='disk out', color='green')
plt.plot(p_, total_profile_2, '-', label='sum', color='m')
plt.ylim(28, 15)
plt.legend();
Сумма выглядит достаточно странно, ибо она больше очевидно чем сам профиль (и в самой работе такого много на самом деле).
In [112]:
surf_R2 = [surf_density(mu=m, M_to_L=M_to_L_R, band='R') for m in total_profile_2]
plt.plot(p_, surf_R2, '-', label='R2 [M/L={:2.2f}]'.format(M_to_L_R))
surf_R2_ = [surf_density(mu=mu_disc(l, mu0=mu0_inner, h=h_inner), M_to_L=M_to_L_R, band='R') for l in p_]
plt.plot(p_, surf_R2_, '-', label='inner only'.format(M_to_L_R)) #только внутренний диск
gutierrez_surf = lambda l: surf_density(total_mu_profile([mu_disc(l, mu0_inner, h_inner)], [mu_disc(l, mu0_out, h_out)]), M_to_L_R, 'R')
plt.legend();
Ожидаемо в принципе, слишком массивно получилось.
Добавлять не будем, явно переоцененная модель.
In [113]:
# all_photometry.append(('Gutierrez R (in)', None, None, None, mu0_inner, h_inner, M_to_L_R,
# lambda l: surf_density(mu=mu_disc(l, mu0=mu0_inner, h=h_inner), M_to_L=M_to_L_R, band='R')))
# all_photometry.append(('Gutierrez R (out)', None, None, None, mu0_out, h_out, M_to_L_R,
# lambda l: surf_density(mu=mu_disc(l, mu0=mu0_out, h=h_out), M_to_L=M_to_L_R, band='R')))
Попробую решить проблему с выступающим профилем:
In [114]:
residual_profile = [-2.5*np.log10(np.power(10, -np.array(l[1])/2.5) - np.power(10, -np.array(mu_disc(l[0], mu0=mu0_out, h=h_out))/2.5)) for l in zip(r_phot2, mu_phot2)]
fig = plt.figure(figsize=[8, 5])
plt.plot(r_phot2, residual_profile, '.')
plt.plot(r_phot2[:161], residual_profile[:161], '.')
plt.ylim(32, 15)
plt.legend();
In [115]:
fig = plt.figure(figsize=[8, 5])
plt.plot(r_phot2, residual_profile, '.')
plt.plot(r_phot2[100:161], residual_profile[100:161], '.')
inner_approx = poly1d(polyfit(r_phot2[100:161], residual_profile[100:161], deg=1))
plt.plot(p_, inner_approx(p_), '--')
plt.ylim(32, 15)
plt.legend();
In [116]:
inner_approx
Out[116]:
Т.к. для диска mu = mu0 + 1.0857*r/h, то
In [117]:
h_approx = 1.0857/inner_approx[1]
h_approx
Out[117]:
In [118]:
mu_inner_approx = inner_approx[0]
Т.е. у меня получилась короче шкала и массивнее внутренний диск.
Проверяем:
In [119]:
fig = plt.figure(figsize=[8, 5])
plt.plot(p_, [mu_disc(l, mu0=inner_approx[0], h=h_approx) for l in p_], '--', label='in', color='#007700')
plt.plot(p_, [mu_disc(l, mu0=mu0_out, h=h_out) for l in p_], '--', label='out', color='#000077')
plt.plot(p_, total_mu_profile([mu_disc(l, mu0=inner_approx[0], h=h_approx) for l in p_],
[mu_disc(l, mu0=mu0_out, h=h_out) for l in p_]), '-', label='sum', color='#770000')
plt.plot(r_phot2, mu_phot2, '.', ms=2)
plt.xlim(0., 400)
plt.ylim(28, 15)
plt.legend();
In [120]:
disc_1_ = [mu_disc(l, mu0=inner_approx[0], h=h_approx) for l in p_]
disc_2 = [mu_disc(l, mu0=mu0_out, h=h_out) for l in p_]
total_profile_2_ = total_mu_profile(disc_1_, disc_2)
surf_R2_ = [surf_density(mu=m, M_to_L=M_to_L_R, band='R') for m in total_profile_2_]
plt.plot(p_, surf_R2_, '-', label='R2 [M/L={:2.2f}]'.format(M_to_L_R))
gutierrez_surf_ = lambda l: surf_density(total_mu_profile([mu_disc(l, inner_approx[0], h_approx)], [mu_disc(l, mu0_out, h_out)]), M_to_L_R, 'R')
plt.legend();
Массивнее стало за счет того, что было $19.5^m$, а стало $19^m$, т.е. примерно в полтора раза увеличилось.
In [121]:
all_photometry.append(('Gutierrez R (new)', None, None, None, (inner_approx[0],mu0_out), (h_approx,h_out), M_to_L_R,
(lambda l: surf_density(mu=mu_disc(l, mu0=inner_approx[0], h=h_approx), M_to_L=M_to_L_R, band='R'),
lambda l: surf_density(mu=mu_disc(l, mu0=mu0_out, h=h_out), M_to_L=M_to_L_R, band='R'))))
Я нашел еще две работы в $V$ - Driel 1994 и Pignatelli 2013, можно использовать другой цвет в калибровке. В первой нет нормальной декомпозиции, возбмем более новую (цвет $B-V$ совпадает примерно 0.85 в первой и 0.90 во второй). Еще ценно, что там есть ссылки на кучу другой фотометрии:
TODO: почему тут написано Pignatelli 2013, хотя картинки из работы 2001 года?
In [122]:
Image('pignatelli_cite_p3.png', width=500)
Out[122]:
In [123]:
# Pignatelli 2013
# у них там две модели, непонятно, какую лучше брать
mu_eff_V = 20.9
r_eff_V = 25.6
n_V = 4 #они использовали Вокулера
mu0d_V = 21.6
h_disc_V = 45.9
In [124]:
Image('pignatelli_photomV.png')
Out[124]:
In [125]:
# данные из базы в Vizier
r_photV, mu_photV = zip(*np.loadtxt("pignatelli_photomV.dat", float, delimiter=','))
plt.plot(r_photV, mu_photV, 's')
plt.ylim(26, 16.5);
In [126]:
bulge_V = [mu_bulge(l, mu_eff=mu_eff_V, r_eff=r_eff_V, n=n_V) for l in p_]
disc_V = [mu_disc(l, mu0=mu0d_V, h=h_disc_V) for l in p_]
total_profile_V = total_mu_profile(bulge_V, disc_V)
fig = plt.figure(figsize=[8, 5])
plt.plot(r_photV, mu_photV, 'o')
plt.plot(p_, bulge_V, '-', label='bulge', color='red')
plt.plot(p_, disc_V, '-', label='disk', color='green')
plt.plot(p_, total_profile_V, '-', label='sum', color='m')
plt.xlim(-10, 175)
plt.ylim(30, 15)
plt.legend();
Вторая модель (2D):
TODO: какую модель брать?
In [127]:
# Pignatelli 2013
# у них там две модели, непонятно, какую лучше брать
mu_eff_V = 20.6
r_eff_V = 18.9
n_V = 4 #они использовали Вокулера
mu0d_V = 20.4
h_disc_V = 29.0
In [128]:
bulge_V = [mu_bulge(l, mu_eff=mu_eff_V, r_eff=r_eff_V, n=n_V) for l in p_]
disc_V = [mu_disc(l, mu0=mu0d_V, h=h_disc_V) for l in p_]
total_profile_V = total_mu_profile(bulge_V, disc_V)
fig = plt.figure(figsize=[8, 5])
plt.plot(r_photV, mu_photV, 'o')
plt.plot(p_, bulge_V, '-', label='bulge', color='red')
plt.plot(p_, disc_V, '-', label='disk', color='green')
plt.plot(p_, total_profile_V, '-', label='sum', color='m')
plt.xlim(-10, 175)
plt.ylim(30, 15)
plt.legend();
У них также есть $M/L$, для диска это 4.2
In [129]:
M_to_L_V = bell_mass_to_light(0.90, 'V', 'B-V')
print M_to_L_V
Достаточно близко.
In [130]:
surf_V = [surf_density(mu=m, M_to_L=M_to_L_V, band='V') for m in disc_V]
plt.plot(p_, surf_R, '-', label='Noord R [M/L={:2.2f}]'.format(M_to_L_R))
plt.plot(p_, surf_B, '-', label='Noord B [M/L={:2.2f}]'.format(M_to_L_B))
plt.plot(p_, surf_J, '-', label='Mendez-Abreu J [M/L={:2.2f}]'.format(M_to_L_J))
plt.plot(p_, surf_R2, '-', label='Gutierrez R [M/L={:2.2f}]'.format(M_to_L_R))
plt.plot(p_, surf_R2_, '-', label='Gutierrez R (new) [M/L={:2.2f}]'.format(M_to_L_R))
plt.plot(p_, surf_V, '-', label='Pignatelli V [M/L={:2.2f}]'.format(M_to_L_V))
plt.legend()
Out[130]:
Хм, занимательно - похоже имеем две конкурирующие модели фотометрии. Причем максимальный диск по ноордермееру даст что-то среднее.
In [131]:
all_photometry.append(('Pignatelli V', r_eff_V, mu_eff_V, n_V, mu0d_V, h_disc_V, M_to_L_V,
lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_V, h=h_disc_V), M_to_L=M_to_L_V, band='V')))
Для $V$ можно вычислить калибровки McGaugh:
In [132]:
mcgaugh_mass_to_light(0.90, 'V')
Out[132]:
Не сильно отличается от используемой, среднее почти такое же.
In [133]:
show_all_photometry_table(all_photometry, scale)
Можно провести тест-сравнение с кривой вращения тонкого диска при заданной фотометрии, если она слишком массивная - то не брать ее (это ограничение сверху).
In [134]:
fig = plt.figure(figsize=[10,6])
plt.plot(test_points, spl_gas(test_points), '-', label='spline')
plt.errorbar(r_wsrt, vel_wsrt, yerr=e_vel_wsrt, fmt='.', marker='.', mew=0, label = 'WSRT')
for ind, photom in enumerate(all_photometry):
if type(photom[5]) == tuple:
plt.plot(test_points, map(lambda l: disc_vel(l, photom[7][0](0), photom[5][0], scale,
Sigma0_2=photom[7][1](0), h_2=photom[5][1]), test_points), next(linecycler), label=photom[0])
else:
plt.plot(test_points, map(lambda l: disc_vel(l, photom[7](0), photom[5], scale), test_points), next(linecycler), label=photom[0])
plt.ylim(0, 300)
plt.xlim(0, 300)
plt.legend(loc='best');
Тут Гутиеррез имеет массивнее хвост, чем Пигнателли, да и в целом стоит его проверить. По максимумам все похожи.
TODO: понять и добавить
TODO: $H_{\alpha}$
Изображение в $H_{\alpha}$ из Hameed 2005 http://iopscience.iop.org/article/10.1086/430211/pdf: (масштабная полоска - 1 кпк)
In [135]:
Image('hameed2005_halpha.png')
Out[135]:
Есть каринка в $\rm{H}_{\alpha}$ из Pignatelli http://mnras.oxfordjournals.org/content/323/1/188.full.pdf:
In [136]:
Image('pignatelli_halpha.png', width=700)
Out[136]:
Кривая вращения и карта из GHASP 2008 http://mnras.oxfordjournals.org/content/388/2/500.full.pdf:
In [137]:
Image('ghasp_2008_halpha_vel.png')
Out[137]:
In [138]:
Image('ghasp_2008_map.png')
Out[138]:
Данные GALEX с расстояниями до образований: http://galex.stsci.edu/GR6/?page=explore&photo=true&objid=6374399074898027266
In [139]:
def plot_SF(ax):
ax.plot([70., 80.], [0., 0.], '-', lw=7., color='red')
ax.plot([190., 210.], [0., 0.], '-', lw=7., color='red')
ax.plot([0., 70.], [0., 0.], '-', lw=7., color='blue')
ax.plot([40., 50.], [0., 0.], '-', lw=7., color='blue') #GALEX
ax.plot([60., 75.], [0., 0.], '-', lw=7., color='blue') #GALEX
plot_SF(plt.gca())
plt.xlim(0, 350)
plt.ylim(0, 200)
Out[139]:
In [140]:
Image('diplom_results.png')
Out[140]:
Устойчиво, когда > 1: $$Q_g = \frac{\Sigma_g^{cr}}{\Sigma_g}=\frac{\kappa c_g}{\pi G \Sigma_g}$$ $$Q_s = \frac{\Sigma_s^{cr}}{\Sigma_s}=\frac{\sigma_R}{\sigma_R^{min}}=\frac{\kappa \sigma_R}{3.36 G \Sigma_s}$$
In [141]:
sound_vel = 6. #скорость звука в газе, км/с
data_lim = max(r_sig_ma)+10. #где заканчиваются данные
In [142]:
fig = plt.figure(figsize=[12, 6])
gd_data = zip(r_g_dens, gas_dens)
#кроме первой точки, т.к. там деление на ноль
plt.plot(r_g_dens, [1./Qg(epicycl=l[0], sound_vel=sound_vel, gas_density=l[1]) for l in
zip([epicyclicFreq_real(gas_approx, x, scale) for x in r_g_dens], gas_dens)], 's-', label='$Q_{gas}$')
plt.plot(r_g_dens, [1./Qg(epicycl=l[0], sound_vel=sound_vel, gas_density=l[1]*He_coeff) for l in
zip([epicyclicFreq_real(gas_approx, x, scale) for x in r_g_dens], gas_dens)], 's-', label='$Q_{gas*He_coeff}$')
plt.plot(r_g_dens, [1./Qg(epicycl=l[0], sound_vel=4., gas_density=l[1]*He_coeff) for l in
zip([epicyclicFreq_real(gas_approx, x, scale) for x in r_g_dens], gas_dens)], 's-', label='$Q_{gas*He_coeff}$ & c=4')
plt.plot(map(lambda l: l/scale, r_q_n), q_n, 'o', label='noordermeer 1F')
plt.plot(r_g_dens, [1./Qs(epicycl=l[0], sigma=l[1], star_density=l[2]) for l in
zip([epicyclicFreq_real(gas_approx, x, scale) for x in r_g_dens],
map(sig_R_maj_max, r_g_dens),
[surf_density(mu_disc(l_, mu0=mu0d_R, h=h_disc_R), M_to_L_R, 'R') for l_ in r_g_dens])], 's-', label='$Q_{star}^{max}$')
plt.plot(r_g_dens, [1./Qs(epicycl=l[0], sigma=l[1], star_density=l[2]) for l in
zip([epicyclicFreq_real(gas_approx, x, scale) for x in r_g_dens],
map(sig_R_maj_min, r_g_dens),
[surf_density(mu_disc(l_, mu0=mu0d_R, h=h_disc_R), M_to_L_R, 'R') for l_ in r_g_dens])], 's-', label='$Q_{star}^{max}$')
plt.axhline(y=1, ls='--')
plt.legend()
plt.ylabel('$Q^{-1}$', fontsize=15);
Одножидкостный уже не похож на диплом, но зато похоже на Ноордермеера (хоть и не совпадает с точностью).
НЕ ИСПРАВЛЕНО ЗА 1.6! Что верно и не надо.
Кинетическое приближение: $$\frac{1}{Q_{\mathrm{eff}}}=\frac{2}{Q_{\mathrm{s}}}\frac{1}{\bar{k}}\left[1-e^{-\bar{k}^{2}}I_{0}(\bar{k}^{2})\right]+\frac{2}{Q_{\mathrm{g}}}s\frac{\bar{k}}{1+\bar{k}^{2}s^{2}}>1\,$$
Гидродинамическое приближение: $$\frac{2\,\pi\, G\, k\,\Sigma_{\mathrm{s}}}{\kappa+k^{2}\sigma_{\mathrm{s}}}+\frac{2\,\pi\, G\, k\,\Sigma_{\mathrm{g}}}{\kappa+k^{2}c_{\mathrm{g}}}>1$$ или $$\frac{1}{Q_{\mathrm{eff}}}=\frac{2}{Q_{\mathrm{s}}}\frac{\bar{k}}{1+\bar{k}^{2}}+\frac{2}{Q_{\mathrm{g}}}s\frac{\bar{k}}{1+\bar{k}^{2}s^{2}}>1$$ для безразмерного волнового числа ${\displaystyle \bar{k}\equiv\frac{k\,\sigma_{\mathrm{s}}}{\kappa}},\, s=c/\sigma$
In [143]:
total_gas_data = zip(r_g_dens, map(lambda l: He_coeff*l, gas_dens))[1:7]
disk_scales = []
for l in all_photometry:
try:
disk_scales.append((l[5][0], l[0].split(' ')[1])) #внутренний диск только
except TypeError:
disk_scales.append((l[5], l[0].split(' ')[1]))
fig = plt.figure(figsize=[10, 6])
ax = plt.gca()
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_disc(l, mu0=mu0d_V, h=h_disc_V), M_to_L_V, 'V'),
star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_V, h=h_disc_V), M_to_L_V, 'V'),
data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='V maj max/min')
plt.ylim(0., 2.5)
plt.axhline(y=1., ls='-', color='grey')
plot_SF(ax)
plt.grid();
In [144]:
total_gas_data = zip(r_g_dens, map(lambda l: He_coeff*l, gas_dens))[1:7]
fig = plt.figure(figsize=[10, 6])
ax = plt.gca()
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=gutierrez_surf,
star_density_min=gutierrez_surf,
data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='Gutierrez R maj max/min')
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=gutierrez_surf_,
star_density_min=gutierrez_surf_,
data_lim=data_lim, color='y', alpha=0.3, disk_scales=disk_scales, label='Gutierrez R (new) maj max/min')
plt.ylim(0., 2.5)
plt.axhline(y=1., ls='-', color='grey')
plot_SF(ax)
plt.grid();
In [145]:
from matplotlib.animation import FuncAnimation
fig = plt.gcf()
plt.figure(figsize=(10,6))
ax = plt.gca()
def animate(i):
ax.cla()
if all_photometry[i][-1] is not None:
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=tot_dens(all_photometry[i][-1]),
star_density_min=tot_dens(all_photometry[i][-1]),
data_lim=data_lim, color=np.random.rand(3), alpha=0.3, disk_scales=[(all_photometry[i][5], '')])
ax.axhline(y=1., ls='-', color='grey')
ax.set_title(all_photometry[i][0])
ax.set_ylim(0., 2.5)
return ax
anim = FuncAnimation(plt.gcf(), animate, frames=len(all_photometry), interval=1000)
In [146]:
# anim.save('..\\..\pics\\'+name+'.gif', writer='imagemagick', fps=1)
In [147]:
# from IPython.display import HTML
# HTML(anim.to_html5_video())
Существует ограничение на максимальный диск в ~0.85 (изотермическое гало) и на субмаксимальный в 0.55-0.6 (NFW гало). Попробуем дотянуть фотметрию до максимальных дисков и посмотрим, как изменятся M/L (скорость зависит как корень из M/L):
In [148]:
fig = plt.figure(figsize=[12,6])
plt.plot(test_points, spl_gas(test_points), '-', label='spline')
plt.plot(test_points, 0.85*spl_gas(test_points), '--', label='max disc')
# plt.plot(test_points, 0.6*spl_gas(test_points), '--', label='submax disc')
plt.errorbar(r_wsrt, vel_wsrt, yerr=e_vel_wsrt, fmt='.', marker='.', mew=0, label = 'WSRT')
# plt.plot(r, vel, '.', label = 'Noord thesis')
max_coeffs = {}
for ind, photom in enumerate(all_photometry):
if type(photom[5]) == tuple:
disc_v = lambda l: disc_vel(l, photom[7][0](0), photom[5][0], scale, Sigma0_2=photom[7][1](0), h_2=photom[5][1])
else:
disc_v = lambda l: disc_vel(l, photom[7](0), photom[5], scale)
values = map(disc_v, test_points)
disc_max = test_points[values.index(max(values))]
max_coeff = 0.85*spl_gas(disc_max)/disc_v(disc_max)
submax_coeff = 0.6*spl_gas(disc_max)/disc_v(disc_max)
if type(photom[5]) == tuple:
plt.plot(test_points, map(lambda l: disc_vel(l, max_coeff**2 *photom[7][0](0), photom[5][0], scale,
Sigma0_2=max_coeff**2 *photom[7][1](0), h_2=photom[5][1]), test_points), next(linecycler), label=photom[0] + '_MAX')
else:
plt.plot(test_points, map(lambda l: disc_vel(l, max_coeff**2 * photom[7](0), photom[5], scale), test_points), next(linecycler), label=photom[0] + '_MAX')
print '{:15s}: M/L was {:2.2f} and for max it equal {:2.2f}, for submax equal {:2.2f}'.format(photom[0], photom[6], photom[6]*max_coeff**2, photom[6]*submax_coeff**2)
max_coeffs[photom[0]] = [max_coeff**2, submax_coeff**2]
plt.ylim(0, 300)
plt.xlim(0, 300)
plt.legend(loc=(1.01, 0.));
В принципе все смотрится достаточно неплохо, хотя конечно скорее максимум правее. $B$ не проходит по M/L.
In [149]:
from matplotlib.animation import FuncAnimation
fig = plt.gcf()
plt.figure(figsize=(10,6))
ax = plt.gca()
def animate(i):
ax.cla()
if all_photometry[i][7] is not None:
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: max_coeffs[all_photometry[i][0]][0]*tot_dens(all_photometry[i][-1])(l),
star_density_min=lambda l: max_coeffs[all_photometry[i][0]][0]*tot_dens(all_photometry[i][-1])(l),
data_lim=data_lim, color=np.random.rand(3), alpha=0.2, disk_scales=[(all_photometry[i][5], '')])
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: max_coeffs[all_photometry[i][0]][1]*tot_dens(all_photometry[i][-1])(l),
star_density_min=lambda l: max_coeffs[all_photometry[i][0]][1]*tot_dens(all_photometry[i][-1])(l),
data_lim=data_lim, color=np.random.rand(3), alpha=0.2, disk_scales=[(all_photometry[i][5], '')])
ax.axhline(y=1., ls='-', color='grey')
ax.set_title(all_photometry[i][0])
ax.set_ylim(0., 2.5)
return ax
anim = FuncAnimation(plt.gcf(), animate, frames=len(all_photometry), interval=1000)
In [150]:
# anim.save('..\\..\pics\\'+name+'_MAXDISCS.gif', writer='imagemagick', fps=1)
In [151]:
# from IPython.display import HTML
# HTML(anim.to_html5_video())
In [152]:
fig = plt.figure(figsize=[10, 4])
# def y_interp_(r, h):
# return 2.8*np.exp(-r/h)
def y_interp_(r, h):
return 0.20*1e9/(2*math.pi*h**2 * (scale * 1000.)**2) * np.exp(-r/h)
points = np.linspace(0.1, 100., 100.)
plt.plot(r_g_dens, gas_dens, 's-')
plt.plot(points, y_interp_(points, 30.), '--', alpha=0.5)
plt.plot(r_g_dens, map(lambda l: y_interp_(l[0], 30.) + l[1], zip(r_g_dens, gas_dens)), '--', alpha=0.5);
In [153]:
fig = plt.figure(figsize=[10, 6])
total_gas_data_ = zip(r_g_dens, [He_coeff*(y_interp_(l[0], 29.) + l[1]) for l in zip(r_g_dens, gas_dens)])[:7]
ax = plt.gca()
plot_2f_vs_1f(ax=ax, total_gas_data=total_gas_data_[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_V, h=h_disc_V), M_to_L_V, 'V'),
star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_V, h=h_disc_V), M_to_L_V, 'V'),
data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='V maj max/min')
plt.ylim(0., 2.5)
plt.axhline(y=1., ls='-', color='grey')
plot_SF(ax)
plt.grid();
Остальные фотометрии:
In [154]:
from matplotlib.animation import FuncAnimation
fig = plt.gcf()
plt.figure(figsize=(10,6))
ax = plt.gca()
def animate(i):
ax.cla()
if type(all_photometry[i][-1]) == tuple:
total_gas_data_ = zip(r_g_dens, [He_coeff*(y_interp_(l[0], all_photometry[i][-3][0]) + l[1]) for l in zip(r_g_dens, gas_dens)])[:7]
else:
total_gas_data_ = zip(r_g_dens, [He_coeff*(y_interp_(l[0], all_photometry[i][-3]) + l[1]) for l in zip(r_g_dens, gas_dens)])[:7]
plot_2f_vs_1f(ax=ax, total_gas_data=total_gas_data_[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=tot_dens(all_photometry[i][-1]),
star_density_min=tot_dens(all_photometry[i][-1]),
data_lim=data_lim, color=np.random.rand(3), alpha=0.3, disk_scales=disk_scales)
ax.axhline(y=1., ls='-', color='grey')
ax.set_title(all_photometry[i][0])
ax.set_ylim(0., 2.5)
return ax
anim = FuncAnimation(plt.gcf(), animate, frames=len(all_photometry), interval=1000)
In [155]:
# anim.save('..\\..\pics\\'+name+'_molec.gif', writer='imagemagick', fps=1)
In [156]:
# from IPython.display import HTML
# HTML(anim.to_html5_video())
Приподнимается, но совсем чуть-чуть.
И оценки с максимальным диском:
In [157]:
from matplotlib.animation import FuncAnimation
fig = plt.gcf()
plt.figure(figsize=(10,6))
ax = plt.gca()
def animate(i):
ax.cla()
if type(all_photometry[i][5]) == tuple:
total_gas_data_ = zip(r_g_dens, [He_coeff*(y_interp_(l[0], all_photometry[i][-3][0]) + l[1]) for l in zip(r_g_dens, gas_dens)])[:7]
else:
total_gas_data_ = zip(r_g_dens, [He_coeff*(y_interp_(l[0], all_photometry[i][-3]) + l[1]) for l in zip(r_g_dens, gas_dens)])[:7]
plot_2f_vs_1f(ax=ax, total_gas_data=total_gas_data_[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: max_coeffs[all_photometry[i][0]][0]*tot_dens(all_photometry[i][-1])(l),
star_density_min=lambda l: max_coeffs[all_photometry[i][0]][0]*tot_dens(all_photometry[i][-1])(l),
data_lim=data_lim, color=np.random.rand(3), alpha=0.3, disk_scales=[(all_photometry[i][5], '')])
plot_2f_vs_1f(ax=ax, total_gas_data=total_gas_data_[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: max_coeffs[all_photometry[i][0]][1]*tot_dens(all_photometry[i][-1])(l),
star_density_min=lambda l: max_coeffs[all_photometry[i][0]][1]*tot_dens(all_photometry[i][-1])(l),
data_lim=data_lim, color=np.random.rand(3), alpha=0.3, disk_scales=[(all_photometry[i][5], '')])
ax.axhline(y=1., ls='-', color='grey')
ax.set_title(all_photometry[i][0])
ax.set_ylim(0., 1.5)
plot_Q_levels(ax, [1., 1.5, 2., 3.])
return ax
anim = FuncAnimation(plt.gcf(), animate, frames=len(all_photometry), interval=1000)
In [158]:
# anim.save('..\\..\pics\\'+name+'_molec_MAX.gif', writer='imagemagick', fps=1)
In [159]:
# from IPython.display import HTML
# HTML(anim.to_html5_video())
In [123]:
summary_imgs_path = '..\\..\pics\\notebook_summary\\'
def save_model_plot(path):
fig, axes = plt.subplots(1, 5, figsize=[40,7])
fig.tight_layout()
axes[0].imshow(ImagePIL.open('ngc3898_SDSS_labeled.jpeg'), aspect='auto')
axes[0].set_title(name)
try:
axes[1].errorbar(r_sig_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='.', mew=0, color='red', label='$\sigma_{los}^{maj}$')
axes[1].plot(points, map(sig_R_maj_min, points))
axes[1].plot(points, map(sig_R_maj_max, points))
axes[1].plot(points, map(sig_R_maj_maxmaxtrue, points))
except Exception:
pass
try:
axes[1].errorbar(r_sig_mi, sig_mi, yerr=e_sig_mi, fmt='.', marker='.', mew=0, color='b', label='$\sigma_{los}^{min}$')
axes[1].plot(points, map(sig_R_minor_min, points), '--')
axes[1].plot(points, map(sig_R_minor_max, points), '--')
except Exception:
pass
axes[1].set_ylim(0,250)
axes[1].set_xlim(0, 105)
axes[1].grid()
axes[1].legend()
axes[1].set_title('Dispersions')
for photom in all_photometry:
axes[2].plot(r_g_dens, map(tot_dens(photom[-1]), r_g_dens), '-', label='{} [M/L={:2.2f}]'.format(photom[0], photom[-2]))
axes[2].set_xlim(0, 150)
axes[2].set_ylim(0, 300)
axes[2].set_title('Photometry')
axes[2].grid()
axes[2].legend()
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)')
axes[3].set_title('Gas')
axes[3].grid()
axes[3].set_xlim(0, 200)
axes[3].legend()
#change this
@save_model(models_path+'n3898_modelRmax.npy')
def plot_2f_vs_1f_(*args, **kwargs):
plot_2f_vs_1f(*args, **kwargs)
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:10],
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',
ML = 7.80,
CO = lambda l: y_interp_(l, h_disc_R),
verbose=True)
@save_model(models_path+'n3898_modelR2dmax.npy')
def plot_2f_vs_1f_(*args, **kwargs):
plot_2f_vs_1f(*args, **kwargs)
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:10],
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=mu_inner_approx, 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=mu_inner_approx, 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',
ML = 2.93,
CO = lambda l: y_interp_(l, h_approx))
axes[4].set_ylim(0., 2.5)
axes[4].set_xlim(0., 130.)
axes[4].axhline(y=1., ls='-', color='grey')
plot_SF(axes[4])
axes[4].grid()
axes[4].set_title('Instability')
plt.savefig(path+name+'.png', format='png', bbox_inches='tight');
save_model_plot(summary_imgs_path)
Schaye (2004), 'cold gas phase': $$\Sigma_g > 6.1 f_g^{0.3} Z^{-0.3} I^{0.23}$$ или при constant metallicity of 0.1 $Z_{sun}$ and interstellar flux of ionizing photons 10^6 cm−2 s−1: $$\Sigma_g > 6.1 \frac{\Sigma_g}{\Sigma_g + \Sigma_s}$$
In [119]:
plt.plot(zip(*total_gas_data)[0], zip(*total_gas_data)[1], 'o-')
for photom in all_photometry:
if photom[7] is not None:
dens_s04 = [Sigma_crit_S04(l[0], l[1], tot_dens(photom[7])) for l in total_gas_data]
plt.plot(zip(*total_gas_data)[0], dens_s04, '--', label=photom[0])
plt.legend()
plt.ylim(0, 10.);
Видимо везде неустойчиво.
Hunter et al (1998), 'competition with shear' according to Leroy: $$\Sigma_A = \alpha_A\frac{\sigma_g A}{\pi G}$$
In [120]:
fig, (ax1, ax2) = plt.subplots(1,2,figsize=[16, 5])
ax1.plot(test_points, [oort_a(x, gas_approx) for x in test_points], '-', label='poly')
ax1.plot(test_points, [oort_a(x, spl_gas) for x in test_points], '-', label='spline')
ax1.set_xlabel('$R, arcsec$')
ax1.set_ylabel('$d\Omega/dr$', fontsize=15)
ax1.legend()
ax1.set_ylim(0, 100.)
dens_A = [Sigma_crit_A(l, spl_gas, 2., 6.) for l in zip(*total_gas_data)[0]]
ax2.plot(zip(*total_gas_data)[0], dens_A, '--')
ax2.plot(zip(*total_gas_data)[0], zip(*total_gas_data)[1], 'o-')
ax2.set_ylim(0, 10.);
Cудя по всему тоже везде неустойчиво.
Интересный вариант для тех галактик, в которых есть данные по газу. Разница между скоростями вращения звезд и газа вокруг центра галактики называется ассиметричным сдвигом и описывается следующим уравнением (Binney & Tremaine 1987): $$v_{\mathrm{c}}^{2}-\bar{v}_{\varphi}^{2}=\sigma_{R}^{2}\left(\frac{\sigma_{\varphi}^{2}}{\sigma_{R}^{2}}-1-\frac{\partial\ln\Sigma_{\mathrm{s}}}{\partial\ln R}-\frac{\partial\ln\sigma_{R}^{2}}{\partial\ln R}\right)\,$$ Отношение ${\displaystyle \frac{\sigma_{\varphi}^{2}}{\sigma_{R}^{2}}}$ знаем из соответствующего уравнения. Поймем, как в этом выражении вычисляется логарифмическая производная ${\displaystyle \frac{\partial\ln\Sigma_{\mathrm{s}}}{\partial\ln R}}$. Если отношение массы к светимости принять постоянной вдоль радиуса величиной, то в производной ${\displaystyle \frac{\partial\ln\Sigma_{\mathrm{s}}}{\partial\ln R}}$ можно использовать поверхностную яркость звездного диска вместо поверхностной плотности $\Sigma_{\mathrm{s}}$ в тех полосах, которые трассируют старое звездное население. Это означает, что логарифмическая производная должна быть заменена отношением $-{\displaystyle \frac{R}{h_{\text{d}}}}\,,$ где $h_{\text{d}}$ --- экспоненциальный масштаб диска. Вычисление $\frac{\partial\ln\sigma_{R}^{2}}{\partial\ln R}$ из кинематического масштаба равно $-\frac{2R}{h_{kin}}$
In [121]:
def sigR2Evaluation(R, h, h_kin, p_star, p_gas):
'''Вычисление sigmaR^2 в случае, если уже известен кинетический масштаб.'''
return (p_gas(R) ** 2 - p_star(R) ** 2 ) / ( sigPhi_to_sigR_real(R) - 1 + R / h + R / h_kin )
def asymmetricDriftEvaluation(r_pc, h, path, p_star, p_gas, upperLimit):
'''Вычисление ассиметричного сдвига на основе формулы (21) из методички. Логарифмическая производная от радиальной
дисперсии скоростей считается как предложено в статье Silchenko et al. 2011, экспонентой фитируется для R > 1h.
Сами значения считаются только для тех точек, есть данные и по газу и по звездам.'''
eps = 0.1
h_kin = 0
h_kin_next = h
sigR2 = []
upper = upperLimit
r_gt_1h = filter(lambda x: x > h and x <= upper, r_pc)
expfit = poly1d(1)
h_disc = h
print '#!!!!!!!!!!!!# Asymmetric drift evaluation procedure with eps = ' + str(eps) + ' starts.'
while(abs(h_kin - h_kin_next) > eps):
h_kin = h_kin_next
sigR2[:] = []
for R in r_gt_1h:
sigR2.append((p_gas(R) ** 2 - p_star(R) ** 2 ) / ( sigPhi_to_sigR_real(R) - 1 + R / h + R / h_kin ))
sigR2 = map(math.log, sigR2)
expfit = poly1d(polyfit(r_gt_1h, sigR2, deg=1))
h_kin_next = (-1 / expfit.coeffs[0])
print '#!!!!!!!!!!!!# Next approx h_kin =', h_kin_next
h_kin = h_kin_next
sigR2[:] = []
for R in r_pc:
sigR2.append((p_gas(R) ** 2 - p_star(R) ** 2 ) / ( sigPhi_to_sigR_real(R) - 1 + R / h + R / h_kin ))
sigR20 = math.exp(expfit.coeffs[1])
# rexp_sigR2 = evalStartExp(r_pc, sigR2, lambda x: sigR20 * math.exp(-x / h_kin))
return sigR20, h_kin, [sigR2Evaluation(R, h, h_kin, p_star, p_gas) for R in r_pc]
sigR20, h_kin, sigR2 = asymmetricDriftEvaluation(r_sig_ma, 30., '', star_approx, spl_gas, sig_maj_lim)
In [122]:
import scipy.interpolate
fig = plt.figure(figsize=[10, 7])
plt.plot(points, map(sig_R_maj_minmin, points), label = 'maj minmin')
plt.plot(points, map(sig_R_maj_min, points), label = 'maj min')
plt.plot(points, map(sig_R_maj_max, points), label = 'maj max')
plt.plot(points, map(sig_R_maj_maxmax, points), label = 'maj maxmax')
plt.plot(points, map(sig_R_maj_maxmaxtrue, points), label = 'maj maxmaxtrue')
plt.plot(r_sig_ma, np.sqrt(sigR2), 'o')
plt.plot(points, map(lambda l: np.sqrt(sigR20 * math.exp(-l / h_kin)), points), '--', alpha=0.5)
y_interp = scipy.interpolate.interp1d(r_sig_ma, np.sqrt(sigR2))
@flat_end(sig_maj_lim)
def ad_interp_(r):
return y_interp(r)
plt.plot(points[2:], map(ad_interp_, points[2:]), '--', alpha=0.5)
plt.legend()
plt.ylim(50,450)
plt.xlim(0,100);
In [123]:
fig = plt.figure(figsize=[10, 6])
ax = plt.gca()
plot_2f_vs_1f(ax=ax, 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=ad_interp_,
sigma_min=ad_interp_,
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 AD')
plt.ylim(0., 2.5)
plt.axhline(y=1., ls='-', color='grey')
plot_SF(ax)
plt.grid();
Если брать экспонентой:
In [124]:
fig = plt.figure(figsize=[10, 6])
ax = plt.gca()
plot_2f_vs_1f(ax=ax, 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=lambda l: np.sqrt(sigR20 * math.exp(-l / h_kin)),
sigma_min=lambda l: np.sqrt(sigR20 * math.exp(-l / h_kin)),
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 AD')
plt.xlim(0., 300.)
plt.ylim(0., 2.5)
plt.axhline(y=1., ls='-', color='grey')
plot_SF(ax)
plt.grid();
В работах van der Hulst (2016) и Bigiel, Blitz (2012) есть экспоненциальные соотношения для H2+HI (см. заметки).
Можно попробовать использовать это для оценки молекулярной компоненты газа:
In [125]:
r25 = h_disc_B*(25. - mu0d_B)/1.0857
r25, h_disc_B, (25. - mu0d_B)/1.0857
Out[125]:
In [126]:
from scipy.optimize import curve_fit
def func1(x, a):
return a * np.exp(-1.95 * x/r25)
def func2(x, a):
return a * np.exp(-1.65 * x/r25)
popt, pcov = curve_fit(func1, r_g_dens[11:], gas_dens[11:])
points_ = np.linspace(100., max(r_g_dens), 100.)
plt.plot(points_, func1(points_, *popt), '--', alpha=0.3)
popt, pcov = curve_fit(func2, r_g_dens[11:], gas_dens[11:])
points_ = np.linspace(100., max(r_g_dens), 100.)
plt.plot(points_, func2(points_, *popt), '--', alpha=0.3)
for i in range(int(max(r_g_dens)/r25)+1):
if i%2 == 0:
plt.axvspan(i*r25, (i+1)*r25, color='grey', alpha=0.1)
plt.plot(r_g_dens, gas_dens, 's-');
In [127]:
def func(x, a, b):
return a * np.exp(-b * x/r25)
popt, pcov = curve_fit(func, r_g_dens[7:], gas_dens[7:])
print popt[1]
points_ = np.linspace(100., max(r_g_dens), 100.)
plt.plot(points_, func(points_, *popt), '--', alpha=0.3)
for i in range(int(max(r_g_dens)/r25)+1):
if i%2 == 0:
plt.axvspan(i*r25, (i+1)*r25, color='grey', alpha=0.1)
plt.plot(r_g_dens, gas_dens, 's-');
In [128]:
def func(x, a, b):
return a * np.exp(-b * x/r25)
for i in range(7, len(r_g_dens)-5):
popt, pcov = curve_fit(func, r_g_dens[i:], gas_dens[i:])
print popt[1]
points_ = np.linspace(100., max(r_g_dens), 100.)
plt.plot(points_, func(points_, *popt), '--', alpha=0.3)
for i in range(int(max(r_g_dens)/r25)+1):
if i%2 == 0:
plt.axvspan(i*r25, (i+1)*r25, color='grey', alpha=0.1)
plt.plot(r_g_dens, gas_dens, 's-')
plt.ylim(0, 6.);
Тут учитывается толщина диска:
In [246]:
plot_RF13_vs_2F(r_g_dens=r_g_dens, HI_gas_dens=gas_dens, CO_gas_dens=[y_interp_(l, h_disc_R) for l in r_g_dens],
epicycl=epicyclicFreq_real, sound_vel=sound_vel, sigma_R_max=sig_R_maj_max, sigma_R_min=sig_R_maj_min,
star_density=lambda l: surf_density(mu_disc(l, mu0=mu0d_R, h=h_disc_R), 7.80, 'R'),
alpha_max=0.7, alpha_min=0.3, scale=scale, gas_approx=spl_gas, thin=False)
А тут нет:
In [247]:
plot_RF13_vs_2F(r_g_dens=r_g_dens, HI_gas_dens=gas_dens, CO_gas_dens=[y_interp_(l, h_disc_R) for l in r_g_dens],
epicycl=epicyclicFreq_real, sound_vel=sound_vel, sigma_R_max=sig_R_maj_max, sigma_R_min=sig_R_maj_min,
star_density=lambda l: surf_density(mu_disc(l, mu0=mu0d_R, h=h_disc_R), 7.80, 'R'),
alpha_max=0.7, alpha_min=0.3, scale=scale, gas_approx=spl_gas, thin=True)
plt.savefig('..\\..\pics\\RF13\\'+name+'.png', format='png', bbox_inches='tight');
Видно, что согласие достаточно хорошее.
Влияние скорости звука:
In [239]:
%%time
plot_param_depend(ax=plt.gca(), data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='R maxdisc', N = 20,
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:9],
epicycl=epicyclicFreq_real,
gas_approx=spl_gas,
sound_vel=list(np.linspace(4., 15., 20)),
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'))
r25 = h_disc_R*(25. - mu0d_R)/1.0857
plot_2f_vs_1f(ax=plt.gca(), 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=[50*np.exp(-2*l/r25) if l < r25 else 50*np.exp(-2.) for l in r_g_dens[1:7]],
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='m', alpha=0.3, disk_scales=disk_scales, label='R Noord maxdisc')
plt.savefig('..\\..\pics\\cg\\'+name+'.png', format='png', bbox_inches='tight');
Влияние убирания молек. газа:
In [230]:
ax = plt.gca()
plot_2f_vs_1f(ax=ax, 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=ax, total_gas_data=zip(r_g_dens, [He_coeff*(0. + 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='m', alpha=0.1, disk_scales=disk_scales, label='R Noord maxdisc without H2')
plot_2f_vs_1f(ax=ax, total_gas_data=zip(r_g_dens, [He_coeff*(y_interp_(l[0], 2*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 scale=2h')
ax.set_ylim(0., 2.5)
ax.set_xlim(0., 130.)
ax.axhline(y=1., ls='-', color='grey')
plot_SF(ax)
ax.grid()
plt.savefig('..\\..\pics\\He\\'+name+'.png', format='png', bbox_inches='tight');
Влияние изменения M/L:
In [293]:
%%time
plot_param_depend(ax=plt.gca(), data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='R maxdisc', N = 10,
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:9],
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=map(lambda c: lambda l: surf_density(mu_disc(l, mu0=mu0d_R, h=h_disc_R), c, 'R'), np.linspace(1., 12., 10)),
star_density_min=map(lambda c: lambda l: surf_density(mu_disc(l, mu0=mu0d_R, h=h_disc_R), c, 'R'), np.linspace(1., 12., 10)));
Влияние масштаба распределения молекулярного газа (т.е. по сути, какая скорость убывания от $h$):
In [294]:
%%time
plot_param_depend(ax=plt.gca(), data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='R maxdisc', N = 5,
total_gas_data=[zip(r_g_dens, [He_coeff*(y_interp_(l[0], c*h_disc_R) + l[1]) for l in zip(r_g_dens, gas_dens)])[1:9] for c in
[0.25, 0.5, 1., 1.5, 2., 3.]],
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'));
Влияние величины молекулярного газа:
In [295]:
%%time
plot_param_depend(ax=plt.gca(), data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='R maxdisc', N = 7,
total_gas_data=[zip(r_g_dens, [He_coeff*(c*y_interp_(l[0], h_disc_R)/y_interp_(0., h_disc_R) + l[1]) for l in zip(r_g_dens, gas_dens)])[1:9]
for c in [1., 5., 10., 16.5, 20., 50., 100.]],
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'));
Замена spl_gas на gas_approx:
In [296]:
%%time
plot_param_depend(ax=plt.gca(), data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='R maxdisc', N = 2,
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:9],
epicycl=epicyclicFreq_real,
gas_approx=[spl_gas, gas_approx],
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'));
Разные реалистичные дисперсии:
In [297]:
%%time
plot_param_depend(ax=plt.gca(), data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='R maxdisc', N = 10,
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:9],
epicycl=epicyclicFreq_real,
gas_approx=spl_gas,
sound_vel=[[c*np.exp(-2*l/r25) if l < r25 else c*np.exp(-2.) for l in r_g_dens[1:]] for c in list(np.linspace(6., 100., 10))],
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'));
Необходимо узнать, как влияет разброс у гле наклона на итоговый результат. К сожалению кроме как вручную это сложно сделать.
In [258]:
# fig = plt.figure(figsize=[10,8])
# gas_approxes = []
# spl_gases = []
# for i in [36., 40.]:
# incl = i
# sin_i, cos_i = np.sin(incl*np.pi/180.), np.cos(incl*np.pi/180.)
# # Данные по кинематике газа Struve, WSRT (не исправлено за наклон)
# r_wsrt, vel_wsrt, e_vel_wsrt = zip(*np.loadtxt("v_gas_WSRT.dat", float))
# r_wsrt, vel_wsrt, e_vel_wsrt = correct_rotation_curve(r_wsrt, vel_wsrt, e_vel_wsrt, 0.0, v0, incl)
# # Данные по кинематике газа Noordermee 2007, WSRT (не исправлено за наклон?)
# r_noord, vel_noord, e_vel_noord = zip(*np.loadtxt("v_gas_noord.dat", float))
# r_noord, vel_noord, e_vel_noord = correct_rotation_curve(r_noord, vel_noord, e_vel_noord, 0.0, v0, incl)
# plt.plot(r_wsrt, vel_wsrt, '.-', label="gas Struve")
# plt.plot(r_noord, vel_noord, '.-', label="gas Noordermeer 2007")
# gas_approx = poly1d(polyfit(r_ma_n, vel_ma_n, deg=5))
# test_points = np.linspace(min(r_ma_n), max(r_ma_n), 100)
# plt.plot(test_points, gas_approx(test_points), '--', label='approx')
# gas_approxes.append(gas_approx)
# spl_gas = inter.UnivariateSpline(r_noord, vel_noord, k=3, s=10000.)
# plt.plot(test_points, spl_gas(test_points), '-', label='spl')
# spl_gases.append(spl_gas)
# plt.ylim(0, 450)
# plt.legend(loc='lower right');
In [259]:
# fig = plt.figure(figsize=[12, 8])
# for ind, i in enumerate([36., 40.]):
# incl = i
# sin_i, cos_i = np.sin(incl*np.pi/180.), np.cos(incl*np.pi/180.)
# plt.plot(test_points, [epicyclicFreq_real(gas_approxes[ind], x, scale) for x in test_points], '-', label='poly')
# plt.plot(test_points, [epicyclicFreq_real(spl_gases[ind], x, scale) for x in test_points], '-', label='spline')
# print epicyclicFreq_real(spl_gases[ind], 10., scale)
# def epicyclicFreq_real_(spl_gas, x, scale):
# '''продливаем дальше без производной на плато'''
# if x < 60.:
# return epicyclicFreq_real(spl_gas, x, scale)
# else:
# return sqrt(2)*arctanlaw(x, m,c,d)/(x*scale)
# # TODO: check scale multiplication
# # plt.plot(np.linspace(1., 100., 100), [epicyclicFreq_real_(gas_approx, x, scale) for x in np.linspace(1., 100., 100)], '-', label='contin')
# plt.xlabel('$R, arcsec$')
# plt.ylabel('$\kappa,\, km/s/kpc$', fontsize=15)
# plt.ylim(0, 100)
# plt.legend();
In [260]:
# print 145.78/159.43, np.sin(36.*np.pi/180.)/np.sin(40.*np.pi/180.)
In [266]:
plt.errorbar(r_sig_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='.', mew=0, color='red', label='$\sigma_{los}^{maj}$')
for ind, i in enumerate([53., 69.]):
incl = i
sin_i, cos_i = np.sin(incl*np.pi/180.), np.cos(incl*np.pi/180.)
# plt.plot(points, map(sig_R_maj_minmin, points), label = 'minmin')
plt.plot(points, map(sig_R_maj_min, points), label = 'min')
plt.plot(points, map(sig_R_maj_max, points), label = 'max')
# plt.plot(points, map(sig_R_maj_maxmax, points), label = 'maxmax')
# plt.plot(points, map(sig_R_maj_maxmaxtrue, points), label = 'maxmaxtrue')
plt.legend()
plt.ylim(0,400)
plt.xlim(0,100);
In [267]:
# 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')
In [276]:
max_coeffs_incl = []
for ind, i in enumerate([53., 69.]):
incl = i
sin_i, cos_i = np.sin(incl*np.pi/180.), np.cos(incl*np.pi/180.)
max_coeffs = {}
for photom in all_photometry:
if type(photom[5]) == tuple:
disc_v = lambda l: disc_vel(l, photom[7][0](0), photom[5][0], scale, Sigma0_2=photom[7][1](0), h_2=photom[5][1])
else:
disc_v = lambda l: disc_vel(l, photom[7](0), photom[5], scale)
values = map(disc_v, test_points)
disc_max = test_points[values.index(max(values))]
max_coeff = 0.85*spl_gas(disc_max)/disc_v(disc_max)
submax_coeff = 0.6*spl_gas(disc_max)/disc_v(disc_max)
print '{:15s}: M/L was {:2.2f} and for max it equal {:2.2f}, for submax equal {:2.2f}'.format(photom[0], photom[6], photom[6]*max_coeff**2, photom[6]*submax_coeff**2)
max_coeffs[photom[0]] = [photom[6]*max_coeff**2, photom[6]*submax_coeff**2]
max_coeffs_incl.append(max_coeffs)
In [277]:
fig = plt.figure(figsize=[10, 6])
ax = plt.gca()
for ind, i in enumerate([53., 69.]):
incl = i
sin_i, cos_i = np.sin(incl*np.pi/180.), np.cos(incl*np.pi/180.)
plot_2f_vs_1f(ax=ax, 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), max_coeffs_incl[ind]['Noorder R'][0], 'R'),
star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_R, h=h_disc_R), max_coeffs_incl[ind]['Noorder R'][0], 'R'),
data_lim=data_lim, color=cm.rainbow(np.linspace(0, 1, 2))[ind], alpha=0.3, disk_scales=disk_scales, label='R Noord maxdisc inc={}'.format(incl))
print max_coeffs_incl[ind]['Noorder R'][0]
plot_2f_vs_1f(ax=ax, 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=max_coeffs_incl[ind]['Gutierrez R (new)'][0], band='R'),
lambda l2: surf_density(mu=mu_disc(l2, mu0=mu0_out, h=h_out), M_to_L=max_coeffs_incl[ind]['Gutierrez R (new)'][0], 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=max_coeffs_incl[ind]['Gutierrez R (new)'][0], band='R'),
lambda l2: surf_density(mu=mu_disc(l2, mu0=mu0_out, h=h_out), M_to_L=max_coeffs_incl[ind]['Gutierrez R (new)'][0], band='R')))(l),
data_lim=data_lim, color=cm.rainbow(np.linspace(0, 1, 2))[ind], alpha=0.2, disk_scales=disk_scales, label='R Gutierrez 2d maxdisc inc={}'.format(incl))
plt.ylim(0., 2.5)
plt.axhline(y=1., ls='-', color='grey')
plot_SF(ax)
plt.grid()
plt.savefig('..\\..\pics\\incl_summary\\'+name+'.png', format='png', bbox_inches='tight');
In [278]:
incl = 61.
sin_i, cos_i = np.sin(incl*np.pi/180.), np.cos(incl*np.pi/180.)
In [298]:
fig, axes = plt.subplots(1, 5, figsize=[40,7])
fig.tight_layout()
try:
axes[0].errorbar(r_sig_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='.', mew=0, color='red', label='$\sigma_{los}^{maj}$')
axes[0].plot(points, map(sig_R_maj_min, points))
axes[0].plot(points, map(sig_R_maj_max, points))
axes[0].plot(points, map(sig_R_maj_maxmaxtrue, points))
except Exception:
pass
axes[1].set_ylim(0,250)
axes[1].set_xlim(0, 105)
axes[1].grid()
axes[1].legend()
axes[1].set_title('Dispersions')
for photom in all_photometry:
axes[2].plot(r_g_dens, map(tot_dens(photom[-1]), r_g_dens), '-', label='{} [M/L={:2.2f}]'.format(photom[0], photom[-2]))
axes[2].set_xlim(0, 150)
axes[2].set_ylim(0, 300)
axes[2].set_title('Photometry')
axes[2].grid()
axes[2].legend()
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)')
axes[3].set_title('Gas')
axes[3].grid()
axes[3].set_xlim(0, 200)
axes[3].legend()
#change this
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')
axes[4].set_ylim(0., 2.5)
axes[4].set_xlim(0., 130.)
axes[4].axhline(y=1., ls='-', color='grey')
plot_SF(axes[4])
axes[4].grid()
axes[4].set_title('Instability')
Out[298]:
In [ ]:
In [299]:
disk_scales2 = []
# disk_scales2.append(disk_scales[2])
disk_scales2.append(disk_scales[1])
axes[2].plot([0, 1], [-1, -2], 'v-', color='b', label=r'$Q_g^{-1}$')
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_scales2, label='R Noordermeer 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_scales2, label='R Gutierrez 2d maxdisc')
axes[4].set_ylim(0., 2.5)
axes[4].set_xlim(0., 130.)
axes[4].axhline(y=1., ls='-', color='grey')
plot_SF(axes[4])
axes[4].grid()
axes[4].set_title('Instability')
ax = axes[2]
ax.set_ylim(0., 1.0)
ax.set_xlim(0., 100.)
ax.axhline(y=1., ls='-', color='grey')
ax.plot([5., 60.], [0., 0.], '-', lw=7., color='red')
ax.grid()
ax.set_ylabel(r'$Q^{-1}$', fontsize=20)
ax.set_xlabel(r'$R,\,arcsec$', fontsize=20)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
plt.legend(fontsize=20)
plt.show()
In [300]:
fig = plt.figure(figsize=[12, 8])
ax = plt.gca()
ax.errorbar(r_sig_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='.', mew=0, color='red', label='$\sigma_{los}^{maj}$')
ax.plot(points, map(sig_R_maj_min, points), label='min')
ax.plot(points, map(sig_R_maj_max, points), label='max')
ax.legend(fontsize=20)
ax.set_xlabel(r'$R,\,arcsec$', fontsize=20)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
ax.set_ylim(0,350)
ax.set_xlim(0, 100)
ax.set_ylabel(r'$\sigma, km/s$', fontsize=20)
ax.grid()
plt.show()
In [301]:
fig = plt.figure(figsize=[12, 8])
ax = plt.gca()
ax.plot(r_g_dens, gas_dens, 'd-', label=r'$\rm{HI}$', ms=10)
ax.plot(r_g_dens, [y_interp_(l, h_disc_R) for l in r_g_dens], '--', label=r'$\rm{H_2}$')
ax.plot(r_g_dens, [He_coeff*(y_interp_(l[0], h_disc_R) + l[1]) for l in zip(r_g_dens, gas_dens)], 'o-', label=r'$\rm{HI+H_2}$', ms=10)
ax.set_title('Gas')
ax.grid()
ax.set_xlim(0, 200)
ax.legend(fontsize=20)
ax3 = plt.gca()
ax3.set_ylabel(r'$\Sigma,\,M_{sun}/pc^2$', fontsize=20)
ax3.set_xlabel(r'$R,\,arcsec$', fontsize=20)
for tick in ax3.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax3.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax3.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax3.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
plt.show()
In [302]:
fig = plt.figure(figsize=[11, 5])
disk_scales2 = []
disk_scales2.append((36.2 , 'R'))
disk_scales2.append((19.11, 'R2'))
ax = plt.gca()
ax.plot([0, 1], [-1, -2], 'v-', color='b', label=r'$Q_g^{-1}$')
plot_2f_vs_1f(ax=ax, 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_scales2, label='R Noordermeer maxdisc')
plot_2f_vs_1f(ax=ax, 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_scales2, label='R Gutierrez 2d maxdisc')
ax.set_ylim(0., 1.)
ax.set_xlim(0., 100.)
ax.plot([70., 80.], [0., 0.], '-', lw=5., color='red')
ax.plot([190., 210.], [0., 0.], '-', lw=5., color='red')
ax.plot([0., 70.], [0., 0.], '-', lw=3., color='red')
ax.grid()
ax.set_ylabel(r'$Q^{-1}$', fontsize=20)
ax.set_xlabel(r'$R,\,arcsec$', fontsize=20)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
plt.legend(fontsize=16)
plt.show()
И как пример для статьи:
In [303]:
fig = plt.figure(figsize=[16, 16])
ax = plt.subplot2grid((3,2), (0, 0))
ax.imshow(ImagePIL.open('ngc3898_SDSS_labeled.jpeg'), aspect='auto')
ax.set_xticks([])
ax.set_yticks([])
ax = plt.subplot2grid((3,2), (0, 1))
ax.errorbar(r_sig_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='.', mew=0, color='red', label='$\sigma_{los}^{maj}$')
ax.plot(points, map(sig_R_maj_min, points), label='min')
ax.plot(points, map(sig_R_maj_max, points), label='max')
ax.legend(fontsize=20)
ax.set_xlabel(r'$R,\,arcsec$', fontsize=20)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
ax.set_ylim(0,350)
ax.set_xlim(0, 100)
ax.set_ylabel(r'$\sigma, km/s$', fontsize=20)
ax.grid()
ax = plt.subplot2grid((3,2), (1, 0))
ax.plot(r_g_dens, gas_dens, 'd-', label=r'$\rm{HI}$', ms=10)
ax.plot(r_g_dens, [y_interp_(l, h_disc_R) for l in r_g_dens], '--', label=r'$\rm{H_2}$')
ax.plot(r_g_dens, [He_coeff*(y_interp_(l[0], h_disc_R) + l[1]) for l in zip(r_g_dens, gas_dens)], 'o-', label=r'$\rm{HI+H_2}$', ms=10)
ax.set_title('Gas')
ax.grid()
ax.set_xlim(0, 200)
ax.legend(fontsize=20)
ax.set_ylabel(r'$\Sigma,\,M_{sun}/pc^2$', fontsize=20)
ax.set_xlabel(r'$R,\,arcsec$', fontsize=20)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
ax = plt.subplot2grid((3,2), (1, 1))
ax.plot(test_points, spl_gas(test_points), '-', label='spline')
ax.plot(test_points, 0.85*spl_gas(test_points), '--', label='max disc')
ax.errorbar(r_wsrt, vel_wsrt, yerr=e_vel_wsrt, fmt='.', marker='.', mew=0, label = 'WSRT')
# ax.plot(r, vel, '.', label = 'Noord thesis')
max_coeffs = {}
for ind, photom in enumerate(all_photometry):
if type(photom[5]) == tuple:
disc_v = lambda l: disc_vel(l, photom[7][0](0), photom[5][0], scale, Sigma0_2=photom[7][1](0), h_2=photom[5][1])
else:
disc_v = lambda l: disc_vel(l, photom[7](0), photom[5], scale)
values = map(disc_v, test_points)
disc_max = test_points[values.index(max(values))]
max_coeff = 0.85*spl_gas(disc_max)/disc_v(disc_max)
submax_coeff = 0.6*spl_gas(disc_max)/disc_v(disc_max)
if type(photom[5]) == tuple:
ax.plot(test_points, map(lambda l: disc_vel(l, max_coeff**2 *photom[7][0](0), photom[5][0], scale,
Sigma0_2=max_coeff**2 *photom[7][1](0), h_2=photom[5][1]), test_points), next(linecycler), label=photom[0] + '_MAX')
else:
ax.plot(test_points, map(lambda l: disc_vel(l, max_coeff**2 * photom[7](0), photom[5], scale), test_points), next(linecycler), label=photom[0] + '_MAX')
max_coeffs[photom[0]] = [max_coeff**2, submax_coeff**2]
ax.set_ylim(0, 300)
ax.set_xlim(0, 300)
ax.legend(fontsize=10, loc='lower right')
ax.set_ylabel(r'$v,\,km/s$', fontsize=20)
ax.set_xlabel(r'$R,\,arcsec$', fontsize=20)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
ax = plt.subplot2grid((3,2), (2, 0), colspan=2)
disk_scales2 = []
disk_scales2.append((36.2 , 'R'))
disk_scales2.append((19.11, 'R2'))
ax.plot([0, 1], [-1, -2], 'v-', color='b', label=r'$Q_g^{-1}$')
plot_2f_vs_1f(ax=ax, 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_scales2, label='R Noordermeer maxdisc')
plot_2f_vs_1f(ax=ax, 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_scales2, label='R Gutierrez 2d maxdisc')
ax.set_ylim(0., 1.)
ax.set_xlim(0., 100.)
ax.plot([70., 80.], [0., 0.], '-', lw=5., color='red')
ax.plot([190., 210.], [0., 0.], '-', lw=5., color='red')
ax.plot([0., 70.], [0., 0.], '-', lw=3., color='red')
ax.grid()
ax.set_ylabel(r'$Q^{-1}$', fontsize=20)
ax.set_xlabel(r'$R,\,arcsec$', fontsize=20)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
plt.show()
Другая версия:
In [356]:
fig = plt.figure(figsize=[16, 21])
ax = plt.subplot2grid((9,2), (0, 0), rowspan=2)
ax.imshow(ImagePIL.open('ngc3898_SDSS_labeled.jpeg'), aspect='auto')
ax.set_xticks([])
ax.set_yticks([])
ax = plt.subplot2grid((9,2), (0, 1), rowspan=2)
ax.plot(test_points, spl_gas(test_points), '-', label='spline')
ax.plot(test_points, 0.85*spl_gas(test_points), '--', label='max disc')
ax.errorbar(r_wsrt, vel_wsrt, yerr=e_vel_wsrt, fmt='.', marker='.', mew=0, label = 'WSRT')
ax.grid()
# ax.plot(r, vel, '.', label = 'Noord thesis')
max_coeffs = {}
for ind, photom in enumerate(all_photometry):
if photom[0] not in ['Noorder R', 'Gutierrez R (new)']:
continue
if type(photom[5]) == tuple:
disc_v = lambda l: disc_vel(l, photom[7][0](0), photom[5][0], scale, Sigma0_2=photom[7][1](0), h_2=photom[5][1])
else:
disc_v = lambda l: disc_vel(l, photom[7](0), photom[5], scale)
values = map(disc_v, test_points)
disc_max = test_points[values.index(max(values))]
max_coeff = 0.85*spl_gas(disc_max)/disc_v(disc_max)
submax_coeff = 0.6*spl_gas(disc_max)/disc_v(disc_max)
if type(photom[5]) == tuple:
ax.plot(test_points, map(lambda l: disc_vel(l, max_coeff**2 *photom[7][0](0), photom[5][0], scale,
Sigma0_2=max_coeff**2 *photom[7][1](0), h_2=photom[5][1]), test_points), next(linecycler), label=photom[0] + '_MAX')
else:
ax.plot(test_points, map(lambda l: disc_vel(l, max_coeff**2 * photom[7](0), photom[5], scale), test_points), next(linecycler), label=photom[0] + '_MAX')
max_coeffs[photom[0]] = [max_coeff**2, submax_coeff**2]
ax.set_ylim(0, 300)
ax.set_xlim(0, 300)
ax.legend(fontsize=10, loc='lower right')
ax.set_ylabel(r'$v,\,km/s$', fontsize=20)
ax.set_xlabel(r'$R,\,arcsec$', fontsize=20, labelpad=-50.1)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
ax = plt.subplot2grid((9,2), (2, 0), colspan=2, rowspan=2)
ax.errorbar(r_sig_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='.', mew=0, color='red', label='$\sigma_{los}^{maj}$')
ax.plot(points, map(sig_R_maj_min, points), label='min')
ax.plot(points, map(sig_R_maj_max, points), label='max')
ax.legend(fontsize=20)
ax.set_xlabel(r'$R,\,arcsec$', fontsize=20)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
ax.set_ylim(0,350)
ax.set_xlim(0, 100)
ax.set_ylabel(r'$\sigma, km/s$', fontsize=20)
ax.grid()
ax = plt.subplot2grid((9,2), (4, 0), colspan=2, rowspan=2)
ax.plot(r_g_dens, gas_dens, 'd-', label=r'$\rm{HI}$', ms=10)
ax.plot(r_g_dens, [y_interp_(l, h_disc_R) for l in r_g_dens], '--', label=r'$\rm{H_2}$')
ax.plot(r_g_dens, [He_coeff*(y_interp_(l[0], h_disc_R) + l[1]) for l in zip(r_g_dens, gas_dens)], 'o-', label=r'$\rm{HI+H_2}$', ms=10)
# ax.set_title('Gas')
ax.grid()
ax.set_xlim(0, 100)
ax.legend(fontsize=20)
ax.set_ylabel(r'$\Sigma,\,M_{sun}/pc^2$', fontsize=20)
ax.set_xlabel(r'$R,\,arcsec$', fontsize=20)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
ax = plt.subplot2grid((9,2), (6, 0), colspan=2, rowspan=3)
disk_scales2 = []
disk_scales2.append((36.2 , 'R'))
disk_scales2.append((19.11, 'R2'))
ax.plot([0, 1], [-1, -2], 'v-', color='b', label=r'$Q_g^{-1}$')
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.errorbar(rr, (np.array(invQeff_min)+np.array(invQeff_max))/2., yerr=(np.array(invQeff_max)-np.array(invQeff_min))/2., color=color, alpha=0.6, capthick=2)
ax.plot(rr, (np.array(invQeff_min)+np.array(invQeff_max))/2., 'd-', color=color, alpha=0.6, label=label)
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()
plot_2f_vs_1f_(ax=ax, 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_scales2, label='R Noordermeer maxdisc')
plot_2f_vs_1f_(ax=ax, 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_scales2, label='R Gutierrez 2d maxdisc')
for q_ in [1., 1.5, 2., 3.]:
ax.axhline(y=1./q_, lw=10, alpha=0.15, color='grey')
ax.set_ylim(0., 1.)
ax.set_xlim(0., 100.)
ax.plot([70., 80.], [0., 0.], '-', lw=5., color='red')
ax.plot([190., 210.], [0., 0.], '-', lw=5., color='red')
ax.plot([0., 70.], [0., 0.], '-', lw=3., color='red')
ax.grid()
ax.legend(fontsize=20)
ax.set_ylabel(r'$Q^{-1}$', fontsize=20)
ax.set_xlabel(r'$R,\,arcsec$', fontsize=20)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
plt.tight_layout(pad=0.0, w_pad=0.0, h_pad=-1.0)
fig.subplots_adjust(hspace=0)
plt.show()
Картинка с газом и в лог-шкале:
In [144]:
fig = plt.figure(figsize=[16,6])
ax = plt.gca()
disk_scales2 = []
disk_scales2.append((36.2 , 'R'))
disk_scales2.append((19.11, 'R2'))
ax.plot([0, 1], [-1, -2], 'v-', color='b', label=r'$Q_g^{-1}$')
plot_2f_vs_1f(ax=ax, 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_scales2, label='R Noordermeer maxdisc')
plot_2f_vs_1f(ax=ax, 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_scales2, label='R Gutierrez 2d maxdisc')
ax.set_ylim(0.03, 1.1)
ax.set_xlim(0., 100.)
ax.plot([70., 80.], [0., 0.], '-', lw=5., color='red')
ax.plot([190., 210.], [0., 0.], '-', lw=5., color='red')
ax.plot([0., 70.], [0., 0.], '-', lw=3., color='red')
ax.grid()
for q_ in [1., 1.5, 2., 3.]:
ax.axhline(y=1./q_, lw=7, alpha=0.15, color='grey')
ax.set_ylabel(r'$Q^{-1}$', fontsize=20)
ax.set_xlabel(r'$R,\,arcsec$', fontsize=20)
ax.fill_between([81., 130.], [0., 0.], [2., 2.], hatch='.', color='gray', alpha=0.08)
ax2 = ax.twinx()
ax2.fill_between(r_g_dens, [0.]*len(r_g_dens), [He_coeff*(y_interp_(l[0], h_disc_R) + l[1]) for l in zip(r_g_dens, gas_dens)], alpha=0.08, color='y')
ax2.set_ylabel(r'$\Sigma_g$', fontsize=20)
ax2.tick_params('y')
ax2.set_xlim(0., 100.)
ax2.set_ylim(0., 10.)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
ax.set_yscale('log')
# ax2.set_yscale('log')
plt.show()
In [124]:
import matplotlib as mpl
mpl.style.use('classic')
In [127]:
fig = plt.figure(figsize=[16, 12])
ax = plt.subplot2grid((2,2), (0, 0))
# ax.imshow(ImagePIL.open('ngc3898_SDSS_labeled.jpeg'), aspect='auto')
ax.imshow(ImagePIL.open('sdss3898.png'), aspect='auto')
ax.set_xticks([])
ax.set_yticks([])
ax = plt.subplot2grid((2,2), (1, 0))
ax.errorbar(r_sig_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='.', mew=0, color='red', label='$\sigma_{los}^{maj}$', ms=15)
ax.plot(points, map(sig_R_maj_min, points), label=r'$\sigma_R^{min}$')
ax.plot(points, map(sig_R_maj_max, points), label=r'$\sigma_R^{max}$')
ax.legend(fontsize=30)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(23)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(23)
ax.set_ylim(0, 430)
ax.set_xlim(0, 100)
ax.set_ylabel(r'$\sigma, km/s$', fontsize=30)
ax.set_xlabel(r'$R,\,arcsec$', fontsize=30)
ax.grid()
ax = plt.subplot2grid((2,2), (1, 1))
ax.plot(r_g_dens, gas_dens, 'd-', label=r'$\rm{HI}$', ms=8)
ax.plot(r_g_dens, [y_interp_(l, h_disc_R) for l in r_g_dens], '--', label=r'$\rm{H_2}$')
ax.plot(r_g_dens, [He_coeff*(y_interp_(l[0], h_disc_R) + l[1]) for l in zip(r_g_dens, gas_dens)], 'o-', label=r'$\rm{HI+H_2}$', ms=8)
plt.semilogy(test_points, [surf_density(mu_disc(l, mu0=mu0d_R, h=h_disc_R), 7.80, 'R') for l in test_points], '-.', lw=3, color='m')
plt.semilogy(test_points, [tot_dens((lambda l1: surf_density(mu=mu_disc(l1, mu0=mu_inner_approx, 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) for l in test_points], '--', lw=3, color='g')
ax.grid()
ax.set_xlim(0, 100)
ax.set_ylim(0.1, 2590.5)
ax.legend(fontsize=20)
ax.set_ylabel(u'$\Sigma,\,M_\u2609/\mathrm{pc}^2$', fontsize=30, labelpad=-7)
ax.set_xlabel(r'$R,\,arcsec$', fontsize=30)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(23)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(23)
ax.yaxis.set_tick_params(length=10, which='major')
ax.yaxis.set_tick_params(length=7, which='minor')
ax = plt.subplot2grid((2,2), (0, 1))
ax.errorbar(r_wsrt, vel_wsrt, yerr=e_vel_wsrt, fmt='.', marker='.', mew=1, label = 'Noordermeer (2005)', ms=15, color='red')
ax.plot(test_points, spl_gas(test_points), '-', lw=3)
# ax.plot(test_points, 0.85*spl_gas(test_points), '--', label='max disc limit')
# ax.plot(r, vel, '.', label = 'Noord thesis')
max_coeffs = {}
for ind, photom in enumerate(all_photometry):
if photom[0] in ['Gutierrez R (new)', 'Noorder R']:
if type(photom[5]) == tuple:
disc_v = lambda l: disc_vel(l, photom[7][0](0), photom[5][0], scale, Sigma0_2=photom[7][1](0), h_2=photom[5][1])
else:
disc_v = lambda l: disc_vel(l, photom[7](0), photom[5], scale)
values = map(disc_v, test_points)
disc_max = test_points[values.index(max(values))]
max_coeff = 0.85*spl_gas(disc_max)/disc_v(disc_max)
submax_coeff = 0.6*spl_gas(disc_max)/disc_v(disc_max)
if type(photom[5]) == tuple:
ax.plot(test_points, map(lambda l: disc_vel(l, max_coeff**2 *photom[7][0](0), photom[5][0], scale,
Sigma0_2=max_coeff**2 *photom[7][1](0), h_2=photom[5][1]), test_points), '--', label=r'Gutierrez et al. (2011) $R$', lw=3)
else:
ax.plot(test_points, map(lambda l: disc_vel(l, max_coeff**2 * photom[7](0), photom[5], scale), test_points), '-.', label=r'Noordermeer at al. (2007) $R$', lw=3, color='m')
ax.grid()
ax.set_ylim(0, 300)
ax.set_xlim(0, 200)
ax.legend(fontsize=18, loc='lower right')
ax.set_ylabel(r'$v,\,km/s$', fontsize=30, labelpad=-5)
# ax.set_xlabel(r'$R,\,arcsec$', fontsize=20)
ax.set_xticklabels([])
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(23)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(23)
# plt.tight_layout(pad=0.1, w_pad=0.1, h_pad=0.1)
fig.subplots_adjust(wspace=0.2, hspace=0.02)
plt.savefig(paper_imgs_dir+'3898_obs_data.eps', format='eps', bbox_inches='tight')
plt.savefig(paper_imgs_dir+'3898_obs_data.png', format='png', bbox_inches='tight')
plt.savefig(paper_imgs_dir+'3898_obs_data.pdf', format='pdf', dpi=150, bbox_inches='tight')
plt.show()
Еще одно замечание, по поводу того что мне надо было посмотреть - как я и подозревал, у Rafikov написано следующее "For the local analysis we employ the WKB approximation (or tight-winding approximation) which requires that kr ≫ 1 and allows us to neglect terms proportional to 1/r compared to the terms proportional to k." Т.е. применялось ВКБ и похоже мне действительно придется смотреть, какие у меня там длины волн в максимуме (любой адекватный рецензент про это спросит). $$kr \gg 1$$ $$k=\frac{2\pi}{\lambda}$$ $${\displaystyle \bar{k}\equiv\frac{k\,\sigma_{\mathrm{s}}}{\kappa}}$$ $$k\times r = \frac{\bar{k}\varkappa}{\sigma}\times(r\times scale)$$
Проверим для Noordermeer-a:
In [304]:
plot_WKB_dependencies(r_g_dens=r_g_dens[1:7],
gas_dens=[He_coeff*(y_interp_(l[0], h_disc_R) + l[1]) for l in zip(r_g_dens, gas_dens)][1:7],
epicycl=lambda l: epicyclicFreq_real(spl_gas, l, scale),
sound_vel=sound_vel,
star_density=[surf_density(mu_disc(l, mu0=mu0d_R, h=h_disc_R), 7.80, 'R') for l in r_g_dens[1:7]],
sigma=sig_R_maj_max,
scale=scale,
krange=arange(0.01, 1000, 0.01))
plt.xlabel(r'$k\times r$', fontsize=20)
plt.title('WKB check')
plt.xlim(0, 10);
# plt.xlim(0, 5);
Картинка с максимумами для разных длин волн:
In [332]:
print 'Qs\t Qg\t s'
plot_k_dependencies(r_g_dens=r_g_dens[1:7],
gas_dens=[He_coeff*(y_interp_(l[0], h_disc_R) + l[1]) for l in zip(r_g_dens, gas_dens)][1:7],
epicycl=lambda l: epicyclicFreq_real(spl_gas, l, scale),
sound_vel=sound_vel,
star_density=[surf_density(mu_disc(l, mu0=mu0d_R, h=h_disc_R), 7.80, 'R') for l in r_g_dens[1:7]],
sigma=sig_R_maj_max,
krange=arange(0.01, 1000, 0.01),
show=False)
plt.xlabel(r'$\bar{r}$', fontsize=20)
plt.axvline(x=1, alpha=0.1)
plt.xlim(0, 10);
Максимум на одной и той же длине волны у всех, поскольку такая зависимость - два пика соответствуют двум слагаемым и положение первого из них практически не меняется.
Проверим для Gutierrez-a:
In [306]:
plot_WKB_dependencies(r_g_dens=r_g_dens[1:7],
gas_dens=[He_coeff*(y_interp_(l[0], h_approx) + l[1]) for l in zip(r_g_dens, gas_dens)][1:7],
epicycl=lambda l: epicyclicFreq_real(spl_gas, l, scale),
sound_vel=sound_vel,
star_density=[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) for l in r_g_dens[1:7]],
sigma=sig_R_maj_max,
scale=scale,
krange=arange(0.01, 1000, 0.01))
plt.xlabel(r'$k\times r$', fontsize=20)
plt.title('WKB check')
plt.xlim(0, 100);
# plt.xlim(0, 5);
Картинка с максимумами для разных длин волн:
In [307]:
plot_k_dependencies(r_g_dens=r_g_dens[1:7],
gas_dens=[He_coeff*(y_interp_(l[0], h_approx) + l[1]) for l in zip(r_g_dens, gas_dens)][1:7],
epicycl=lambda l: epicyclicFreq_real(spl_gas, l, scale),
sound_vel=sound_vel,
star_density=[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) for l in r_g_dens[1:7]],
sigma=sig_R_maj_max,
krange=arange(0.01, 1000, 0.01))
plt.xlabel(r'$\bar{r}$', fontsize=20)
plt.axvline(x=1, alpha=0.1)
plt.xlim(0, 100);
In [ ]:
In [ ]:
Картинка неустойчивости, как из статей (квадрат частоты $w^2$).
TODO: почему такие странные провалы? (смотри ниже) Это не очень важно, но странно
In [308]:
from scipy.special import i0e, i1e
from scipy.optimize import fsolve
qg = 6.16
qs = 11.06
s = 0.03587
tp = np.arange(0.1, 100., 0.1)
def twofl_w2(l):
global qg, qs, s, dimlK
return 2. / dimlK / qs * (1 - i0e(dimlK ** 2)) * (1. / (1 - l)) + 2*s*dimlK / qg / (1 + dimlK**2 * s**2 - l) - 1.
sol = []
for d in tp:
dimlK = d
initial_guess = 0.
solution = fsolve(twofl_w2, initial_guess)
sol.append(solution)
def star_w2(qs, dimlK):
return 1. - 2. / dimlK / qs * (1 - i0e(dimlK ** 2))
cold_gas = lambda l: l**2 * s**2 + 1. - 2./qg * s * l
fig = plt.figure(figsize=[8, 6])
plt.axhline(y=0, ls='-.', alpha=0.5)
plt.plot(tp, map(cold_gas, tp), '--', label='G')
plt.plot(tp, [star_w2(1.5, x) for x in tp], '--', label='S')
plt.plot(tp, sol, '--', label='S+G')
plt.xlim(0, 100)
# plt.ylim(-0.2, 1.5)
plt.legend()
# plt.xticks(np.arange(0, 10, 0.5))
plt.ylim(0, 4);
Видимо просто неустойчивости, их можно игнорировать в разумных пределах (если брать другое начальное приближение, то все нормально решается, но автоматически этого не понять).
Функция, предсказывающая минимальное $Q_s$ для неустойчивости:
In [309]:
qg = 6.16
qs = 11.06
s = 0.03587
def predict_hydro_unstable_Qs(Qg, s):
'''Рассчитывает величину Qs, при которой достигается минимальная неустойчивость.
Для этого выражения для двухжидкостной неуст в кинематич. приближении приравнивается 1
и максимум выражения ниже и есть искомое.'''
fuu = lambda l: 2*k*Qg*(1+k**2*s**2)/(1+k**2)/(Qg+Qg*k**2*s**2 - 2*s*k)
candidate_Qs = max([fuu(k) for k in np.arange(0.01, 60000., 1.)])
Qeff = findInvKinemQeffBrute(candidate_Qs, Qg, s, np.arange(0.01, 60000., 1.))[0]
assert abs(Qeff-1) < 0.2
return candidate_Qs
In [310]:
predict_hydro_unstable_Qs(qg, s)
Out[310]:
И проверка этой функции для ряда параметров:
In [311]:
Qgs = []
Qss = []
invQeff = []
for r, gd in zip(r_g_dens, gas_dens)[1:7]: #только до 90" берем, дальне нет данных по дисперсии
Qgs.append(Qg(epicycl=epicyclicFreq_real(gas_approx, r, scale), sound_vel=sound_vel, gas_density=gd*He_coeff))
Qss.append(predict_hydro_unstable_Qs(Qgs[-1], sound_vel/sig_R_maj_max(r)))
qeff = findInvHydroQeffBrentq(Qss[-1], Qgs[-1], sound_vel/sig_R_maj_max(r), np.arange(0.01, 60000., 1.))
print 'Qs = {:2.2f}; Qg = {:2.2f}; Qeff = {:2.2f}'.format(Qss[-1], Qgs[-1], 1./qeff[1])
invQeff.append(qeff)
In [ ]:
In [ ]:
Функция, предсказывающая минимальное $Q_g$ для неустойчивости:
In [398]:
qg = 6.16
qs = 11.06
s = 0.03587
def predict_kinem_unstable_Qg(Qs, s, alpha):
'''Рассчитывает величину Qs, при которой достигается минимальная неустойчивость.
Для этого выражения для двухжидкостной неуст в кинематич. приближении приравнивается alpha
и максимум выражения ниже и есть искомое.'''
fuu = lambda k: 2*s*k/((1+k**2*s**2)*(alpha-2. / k / Qs * (1 - i0e(k ** 2))))
candidate_Qg = max([fuu(k) for k in np.arange(0.01, 60000., 1.)])
Qeff = findInvKinemQeffBrentq(Qs, candidate_Qg, s, np.arange(0.01, 60000., 1.))[1]
assert abs(Qeff-alpha) < 0.2
return candidate_Qg
In [399]:
predict_kinem_unstable_Qg(qg, s, 1.)
Out[399]:
In [ ]:
Посмотрим, сколько нужно $\Sigma_{H_2}$, чтобы была одножидкостная неустойчивость для разных $\alpha$:
In [416]:
colors = cm.rainbow(np.linspace(0, 1, 4))
for index, alpha in enumerate([0.33, 0.5, 0.75, 1.]):
CO_tmp = []
for ind, r in enumerate(r_g_dens[1:7]):
Sigma_CO = (alpha - 1./Qg(epicycl=epicyclicFreq_real(spl_gas, r, scale), sound_vel=sound_vel,
gas_density=gas_dens[1:7][ind]))/(He_coeff*1./Qg(epicycl=epicyclicFreq_real(spl_gas, r, scale), sound_vel=sound_vel, gas_density=1.))
CO_tmp.append(Sigma_CO)
plt.plot(r_g_dens[1:7], CO_tmp, '-', color=colors[index], label=alpha)
plt.legend()
plt.xlabel(r'$R$', fontsize=20)
plt.ylabel(r'$\Sigma_{H_2}$', fontsize=20);
In [ ]:
In [ ]:
In [ ]:
In [312]:
plot_2f_vs_1f(ax=ax, 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_scales2, label='R Noordermeer maxdisc')
In [313]:
%%time
res = []
r = r_g_dens[5]
gd = [He_coeff*(y_interp_(l[0], h_disc_R) + l[1]) for l in zip(r_g_dens, gas_dens)][4]
sd = surf_density(mu_disc(r, mu0=mu0d_R, h=h_disc_R), 7.80, 'R')
print r, gd, sd
for s_v in np.arange(0.1, 10., 1.):
for sig_ in np.arange(50., 350., 25.):
Qgs = Qg(epicycl=epicyclicFreq_real(spl_gas, r, scale), sound_vel=s_v, gas_density=gd)
Qss = Qs(epicycl=epicyclicFreq_real(spl_gas, r, scale), sigma=sig_, star_density=sd)
qeff = findInvKinemQeffBrentq(Qss, Qgs, s_v/sig_, np.arange(0.01, 60000., 1.))
res.append(qeff[1])
In [314]:
z = np.array(res).reshape(len(np.arange(0.1, 10., 1.)), len(np.arange(50., 350., 25.)))
plt.pcolormesh(z)
plt.colorbar()
curves = 10
m = max([max(row) for row in z])
levels = np.linspace(0, m, 10)
# plt.contour(z, colors="white", levels=levels)
plt.show()
In [315]:
def fb(x, y):
return x*y
xmax = ymax = 100
z = numpy.array([[fb(x, y) for x in range(xmax)] for y in range(ymax)])
plt.pcolormesh(z)
plt.colorbar()
curves = 10
m = max([max(row) for row in z])
levels = numpy.arange(0, m, (1 / float(curves)) * m)
plt.contour(z, colors="white", levels=levels)
plt.show()
In [ ]:
In [ ]:
In [316]:
def plot_WKB_dependency(Qs=None, Qg=None, s=None, krange=None, ax=None, label=None, color=None, r=None, scale=None, epicycl=None, sound_vel=None):
'''рисуется зависимость между (k x r) и двухжидкостной неустойчивостью, показан максимум (см. уравнение выше), чтобы проверить справедливость WKB'''
TFcriteria = []
sigma = sound_vel/s
_tmp = [inverse_kinem_Qeff_from_k(dimlK, Qg=Qg, Qs=Qs, s=s) for dimlK in krange]
root_for_max, max_val = findInvKinemQeffBrentq(Qs, Qg, s, krange)
factor = epicycl*r*scale/sigma
ax.plot(np.array(krange)*factor, _tmp, '-', label=label, color=color)
ax.plot(root_for_max*factor, max_val, 'o', color=color)
def plot_WKB_dependencies(r_g_dens=None, gas_dens=None, epicycl=None,
sound_vel=None, star_density=None, sigma=None, krange=None, scale=None):
'''рисуем много зависимостей сразу для проверки WKB'''
Qgs, Qss, s_params = [], [], []
maxk = 30.
fig = plt.figure(figsize=[16,8])
ax = plt.gca()
colors = cm.rainbow(np.linspace(0, 1, len(r_g_dens)))
for r, gd, sd, color in zip(r_g_dens, gas_dens, star_density, colors):
Qgs.append(Qg(epicycl=epicycl(r), sound_vel=sound_vel, gas_density=gd))
Qss.append(Qs(epicycl=epicycl(r), sigma=sigma(r), star_density=sd))
s_params.append(sound_vel/sigma(r))
plot_WKB_dependency(Qs=Qss[-1], Qg=Qgs[-1], s=s_params[-1], krange=krange, ax=ax, label=str(r),
color=color, r=r, scale=scale, epicycl=epicycl(r), sound_vel=sound_vel)
maxk = max(maxk, findInvKinemQeffBrentq(Qss[-1], Qgs[-1], s_params[-1], krange)[0]) #not optimal
plt.legend()
plt.xlim(0, maxk+100.)
def inverse_hydro_Qeff_from_k(dimlK, Qg=None, Qs=None, s=None):
return 2.*dimlK / Qs / (1 + dimlK**2) + 2*s*dimlK / Qg / (1 + dimlK**2 * s**2)
def inverse_kinem_Qeff_from_k(dimlK, Qg=None, Qs=None, s=None):
return 2. / dimlK / Qs * (1 - i0e(dimlK ** 2)) + 2*s*dimlK / Qg / (1 + dimlK**2 * s**2)
In [317]:
def plot_k_dependency(Qs=None, Qg=None, s=None, krange=None, ax=None, label=None, color=None):
'''рисуется зависимость между волновыми числами и двухжидкостной неустойчивостью, показан максимум'''
print Qs, Qg, s
TFcriteria = []
_tmp = [inverse_kinem_Qeff_from_k(dimlK, Qg=Qg, Qs=Qs, s=s) for dimlK in krange]
root_for_max, max_val = findInvKinemQeffBrentq(Qs, Qg, s, krange)
ax.plot(krange, _tmp, '-', label=label, color=color)
ax.plot(root_for_max, max_val, 'o', color=color)
def plot_k_dependencies(r_g_dens=None, gas_dens=None, epicycl=None,
sound_vel=None, star_density=None, sigma=None, krange=None):
'''рисуем много зависимостей сразу'''
Qgs, Qss, s_params = [], [], []
maxk = 30.
fig = plt.figure(figsize=[16,8])
ax = plt.gca()
colors = cm.rainbow(np.linspace(0, 1, len(r_g_dens)))
for r, gd, sd, color in zip(r_g_dens, gas_dens, star_density, colors):
Qgs.append(Qg(epicycl=epicycl(r), sound_vel=sound_vel, gas_density=gd))
Qss.append(Qs(epicycl=epicycl(r), sigma=sigma(r), star_density=sd))
s_params.append(sound_vel/sigma(r))
plot_k_dependency(Qs=Qss[-1], Qg=Qgs[-1], s=s_params[-1], krange=krange, ax=ax, label=str(r), color=color)
maxk = max(maxk, findInvKinemQeffBrentq(Qss[-1], Qgs[-1], s_params[-1], krange)[0]) #not optimal
plt.legend()
plt.xlim(0, maxk+100.)
In [318]:
%%time
fig = plt.figure(figsize=[16, 8])
ax = plt.gca()
count = 0
colors = cm.rainbow(np.linspace(0, 1, 27))
for Qs_ in [0.1, 1., 10.]:
for Qg_ in [0.1, 1., 10.]:
for s_ in [0.1, 1., 10.]:
print Qs_, Qg_, s_
plot_k_dependency(Qs=Qs_, Qg=Qg_, s=s_, krange=arange(0.001, 10, 0.001), ax=ax, color=colors[count])
count+=1
plt.xlabel(r'$\bar{r}$', fontsize=20)
plt.axvline(x=1, alpha=0.1)
plt.xlim(0, 5);
In [320]:
# r 15-90
# Qs Qg s
# 5.59185852878 29.5792958425 0.0294726918078
# 3.92966549184 13.6304146876 0.0269366023417
# 2.79468647063 7.05537280878 0.0299828699262
# 3.2978809954 6.07349693121 0.0341363619448
# 4.42732739473 5.52623632628 0.0323983342832
# 5.1828928061 5.1313317673 0.0323983342832
In [321]:
326*6/5./4.32/3.14
Out[321]:
In [322]:
epicyclicFreq_real(spl_gas, 15., scale)
Out[322]:
In [323]:
surf_density(mu_disc(15., mu0=mu0d_R, h=h_disc_R), 7.80, 'R')
Out[323]:
In [324]:
326.*200./818./4.32/3.14
Out[324]:
In [ ]:
In [325]:
def plot_k_dependency(Qs=None, Qg=None, s=None, krange=None, ax=None, label=None, color=None):
'''рисуется зависимость между волновыми числами и двухжидкостной неустойчивостью, показан максимум'''
print '{:7.3f} {:7.3f} {:7.3f} {:7.3f} {:7.3f} {:9.3f} {:7.3f} {:7.3f} {:7.6f} {:7.6f} {:7.6f} {:7.6f}'.format(Qs, Qg, s, Qs/Qg, Qs/s, Qg/s, Qs*s,
Qg*s/Qs, s/Qg/(1+s**2)/Qs, Qg*(s**4) + Qs*(s**2), 2*Qs*s-Qg,
Qg*s**4 + Qs*s - 2*Qg*s**2 - 2*Qs*s**3)
TFcriteria = []
_tmp = [inverse_kinem_Qeff_from_k(dimlK, Qg=Qg, Qs=Qs, s=s) for dimlK in krange]
root_for_max, max_val = findInvKinemQeffBrentq(Qs, Qg, s, krange)
ax.plot(krange, _tmp, '-', label=label, color=color)
ax.plot(root_for_max, max_val, 'o', color=color)
def plot_k_dependencies(r_g_dens=None, gas_dens=None, epicycl=None,
sound_vel=None, star_density=None, sigma=None, krange=None):
'''рисуем много зависимостей сразу'''
Qgs, Qss, s_params = [], [], []
maxk = 30.
fig = plt.figure(figsize=[16,8])
ax = plt.gca()
colors = cm.rainbow(np.linspace(0, 1, len(r_g_dens)))
for r, gd, sd, color in zip(r_g_dens, gas_dens, star_density, colors):
Qgs.append(Qg(epicycl=epicycl(r), sound_vel=sound_vel, gas_density=gd))
Qss.append(Qs(epicycl=epicycl(r), sigma=sigma(r), star_density=sd))
s_params.append(sound_vel/sigma(r))
plot_k_dependency(Qs=Qss[-1], Qg=Qgs[-1], s=s_params[-1], krange=krange, ax=ax, label=str(r), color=color)
maxk = max(maxk, findInvKinemQeffBrentq(Qss[-1], Qgs[-1], s_params[-1], krange)[0]) #not optimal
plt.legend()
plt.xlim(0, maxk+100.)
In [326]:
print 'Qs\t Qg\t s'
plot_k_dependencies(r_g_dens=r_g_dens[1:7],
gas_dens=[He_coeff*(y_interp_(l[0], h_disc_R) + l[1]) for l in zip(r_g_dens, gas_dens)][1:7],
epicycl=lambda l: epicyclicFreq_real(spl_gas, l, scale),
sound_vel=sound_vel,
star_density=[surf_density(mu_disc(l, mu0=mu0d_R, h=h_disc_R), 7.80, 'R') for l in r_g_dens[1:7]],
sigma=sig_R_maj_max,
krange=arange(0.01, 1000, 0.01))
plt.xlabel(r'$\bar{r}$', fontsize=20)
plt.axvline(x=1, alpha=0.1)
plt.xlim(0, 10);
Замена sig_R_maj_max на sig_R_maj_min картины не меняет, только приподнимает как целое.
In [327]:
print 'Qs\t Qg\t s'
plot_k_dependencies(r_g_dens=r_g_dens[1:7],
gas_dens=[He_coeff*(y_interp_(l[0], h_disc_R) + l[1]) for l in zip(r_g_dens, gas_dens)][1:7],
epicycl=lambda l: epicyclicFreq_real(spl_gas, l, scale),
sound_vel=sound_vel,
star_density=[surf_density(mu_disc(l, mu0=mu0d_R, h=h_disc_R), 7.80, 'R') for l in r_g_dens[1:7]],
sigma=sig_R_maj_min,
krange=arange(0.01, 1000, 0.01))
plt.xlabel(r'$\bar{r}$', fontsize=20)
plt.axvline(x=1, alpha=0.1)
plt.xlim(0, 10);
Ополовинивание молгаза
In [328]:
print 'Qs\t Qg\t s'
plot_k_dependencies(r_g_dens=r_g_dens[1:7],
gas_dens=[He_coeff*(y_interp_(l[0], h_disc_R/2.) + l[1]) for l in zip(r_g_dens, gas_dens)][1:7],
epicycl=lambda l: epicyclicFreq_real(spl_gas, l, scale),
sound_vel=sound_vel,
star_density=[surf_density(mu_disc(l, mu0=mu0d_R, h=h_disc_R), 7.80, 'R') for l in r_g_dens[1:7]],
sigma=sig_R_maj_min,
krange=arange(0.01, 1000, 0.01))
plt.xlabel(r'$\bar{r}$', fontsize=20)
plt.axvline(x=1, alpha=0.1)
plt.xlim(0, 10);
Убрать молек. газ
In [329]:
print 'Qs\t Qg\t s'
plot_k_dependencies(r_g_dens=r_g_dens[1:7],
gas_dens=[He_coeff*(0.+ l[1]) for l in zip(r_g_dens, gas_dens)][1:7],
epicycl=lambda l: epicyclicFreq_real(spl_gas, l, scale),
sound_vel=sound_vel,
star_density=[surf_density(mu_disc(l, mu0=mu0d_R, h=h_disc_R), 7.80, 'R') for l in r_g_dens[1:7]],
sigma=sig_R_maj_min,
krange=arange(0.01, 1000, 0.01))
plt.xlabel(r'$\bar{r}$', fontsize=20)
plt.axvline(x=1, alpha=0.1)
plt.xlim(0, 10);
А вот взять не макс. фотометрию - меняет
In [330]:
print 'Qs\t Qg\t s'
plot_k_dependencies(r_g_dens=r_g_dens[1:7],
gas_dens=[He_coeff*(y_interp_(l[0], h_disc_R) + l[1]) for l in zip(r_g_dens, gas_dens)][1:7],
epicycl=lambda l: epicyclicFreq_real(spl_gas, l, scale),
sound_vel=sound_vel,
star_density=[surf_density(mu_disc(l, mu0=mu0d_R, h=h_disc_R), 4.0, 'R') for l in r_g_dens[1:7]],
sigma=sig_R_maj_min,
krange=arange(0.01, 1000, 0.01))
plt.xlabel(r'$\bar{r}$', fontsize=20)
plt.axvline(x=1, alpha=0.1)
plt.xlim(0, 10);
In [ ]: