In [1]:
from IPython.display import HTML
from IPython.display import Image
import os
# import scipy.interpolate as inter
%pylab
%matplotlib inline
%run ../../utils/load_notebook.py
In [2]:
from photometry import *
In [3]:
from instabilities import *
In [4]:
from utils import *
In [5]:
name = 'N2985'
gtype = '(R)SA(r)ab' #TODO: откуда
incl = 36. #adopted by Noordermeer+08, LEDA - 37.9, Epinat+08 - 36.
distance = 21.1 # Noordermeer
scale = 0.102 #kpc/arcsec according Noordermeer
data_path = '../../data/ngc2985'
sin_i, cos_i = np.sin(incl*np.pi/180.), np.cos(incl*np.pi/180.)
In [6]:
# plt.rcParams['image.cmap'] = 'hsv'
# plt.rcParams['axes.prop_cycle'] = plt.cycler('color', plt.get_cmap('hsv')(np.linspace(0, 1.0, 12)))
In [7]:
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')
TODO: add
In [8]:
# Данные из NED
HTML('<iframe src=http://ned.ipac.caltech.edu/cgi-bin/objsearch?objname=ngc+2985&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[8]:
In [9]:
# Данные из HYPERLEDA
HTML('<iframe src=http://leda.univ-lyon1.fr/ledacat.cgi?o=ngc2985 width=1000 height=350></iframe>')
Out[9]:
In [10]:
os.chdir(data_path)
Image('noordermeer_data/n2985_cite_p39.png')
Out[10]:
In [11]:
Image('noordermeer_data/n2985_cite_p110.png')
Out[11]:
In [12]:
Image('noordermeer_data/n2985_cite_pp187_188.png')
Out[12]:
In [13]:
#2MASS
Image('ngc2985_JHK.jpg', width=300)
Out[13]:
DSS картинка 10х10 arcmin:
In [14]:
Image('n2985_DSS_10x10min.jpg')
Out[14]:
Возможно картинка из Хаббла:
да, это действительно она в обработке - на Hubble Legacy Archive https://hla.stsci.edu/ вот в этом снимке видно
In [15]:
Image('n2985_hubble.jpg', width=300)
Out[15]:
Извлеченная картинка самостоятельно:
In [16]:
Image('HST_color.png')
Out[16]:
Совмещенное с DSS изображение:
In [17]:
Image('hst_dss_combined.jpg')
Out[17]:
In [18]:
Image('noordermeer_data/n2985_rc.png')
Out[18]:
In [19]:
r, vel = zip(*np.loadtxt("noordermeer_data/n2985_rc_noorderm.dat", float, delimiter=','))
fig = plt.figure(figsize=[10,8])
plt.plot(r, vel, 's-')
plt.legend();
In [20]:
Image('noordermeer_data/n2985_photom.png')
Out[20]:
TODO: add links and check exact values
TODO: extract $\sigma_g$ values
In [21]:
Image('Dumas_gas_disp.png')
Out[21]:
In [22]:
# Данные по звездной кинематике Noordermeer вдоль большей полуоси (не исправленные за наклон?)
r_ma, vel_ma, e_vel_ma, sig_ma, e_sig_ma = zip(*np.loadtxt("v_stars_ma.dat", float))
# Данные по звездной кинематике HSM99 вдоль большой полуоси (не исправленные за наклон)
r_ma_hsm, vel_ma_hsm, e_vel_ma_hsm, sig_ma_hsm, e_sig_ma_hsm = zip(*np.loadtxt("v_stars_her.dat", float))
# Данные по звездной кинематике Noord вдоль большой полуоси (исправленные за наклон)
r_ma_n, vel_ma_n, e_vel_ma_n, sig_ma_n, e_sig_ma_n = zip(*np.loadtxt("v_stars_noord_1.dat", float))
# Данные по звездной кинематике Noordermeer вдоль малой полуоси (исправленные за наклон?)
r_mi, vel_mi, e_vel_mi, sig_mi, e_sig_mi = zip(*np.loadtxt("v_stars_mi.dat", float))
fig = plt.figure(figsize=[10,8])
plt.errorbar(r_ma, vel_ma, e_vel_ma, fmt='-.', marker='.', mew=0, label="Noord stars maj")
plt.errorbar(r_mi, vel_mi, e_vel_mi, fmt='-.', marker='.', mew=0, label="Noord stars min")
plt.errorbar(r_ma_hsm, vel_ma_hsm, e_vel_ma_hsm, fmt='-.', marker='.', mew=0, label="HSM99 stars maj")
plt.errorbar(r_ma_n, vel_ma_n, e_vel_ma_n, fmt='-.', marker='.', mew=0, label="Noord corr stars maj")
plt.legend();
Малая ось точно не валидна, другие данные похожи. Перегнем:
In [23]:
r_ma, vel_ma, e_vel_ma, sig_ma, e_sig_ma = zip(*sorted(zip(np.abs(r_ma), np.abs(vel_ma), e_vel_ma, sig_ma, e_sig_ma)))
r_ma_hsm, vel_ma_hsm, e_vel_ma_hsm, sig_ma_hsm, e_sig_ma_hsm = zip(*sorted(zip(np.abs(r_ma_hsm), np.abs(vel_ma_hsm), e_vel_ma_hsm, sig_ma_hsm, e_sig_ma_hsm)))
In [24]:
# Данные по звездной кинематике Gerssen вдоль большой полуоси (не исправленные за наклон?)
r_g, vel_g = zip(*np.loadtxt("gerssen_2000_vel.dat", float, delimiter=','))
fig = plt.figure(figsize=[10,8])
plt.errorbar(r_ma, vel_ma, e_vel_ma, fmt='-.', marker='.', mew=0, label="Noord stars maj")
plt.errorbar(r_ma_hsm, vel_ma_hsm, e_vel_ma_hsm, fmt='-.', marker='.', mew=0, label="HSM99 stars maj")
plt.errorbar(r_ma_hsm, map(lambda l: l/sin_i, vel_ma_hsm), e_vel_ma_hsm, fmt='-.', marker='.', mew=0, label="HSM99 corr stars maj")
plt.errorbar(r_ma_n, vel_ma_n, e_vel_ma_n, fmt='-.', marker='.', mew=0, label="Noord corr stars maj")
plt.plot(r_g, vel_g, 'v', label='Gerssen')
plt.legend()
plt.ylim(0, 300);
Данные Герссена тоже хорошо ложатся.
TODO: проверить, что данные валидны (в работе Ноордермеера 2008 видно, что надо использовать верхнюю)
In [25]:
fig = plt.figure(figsize=[8,4])
plt.errorbar(r_ma_n, vel_ma_n, yerr=e_vel_ma_n, fmt='.', marker='.', mew=0, color='blue', label = 'Her98 star maj')
test_points = np.linspace(0.0, max(r_ma_n), 100)
poly_star = poly1d(polyfit(r_ma_n, vel_ma_n, deg=4))
plt.plot(test_points, poly_star(test_points), '-', label='poly deg=4')
def w(arr):
return map(lambda l: 1/(l**2), arr)
import scipy.interpolate as inter
spl = inter.UnivariateSpline(r_ma_n, vel_ma_n, k=3, s=10000., w=w(e_vel_ma_n))
plt.plot(test_points, spl(test_points), '-', label='spl s=10000 w^2')
spl = inter.UnivariateSpline(r_ma_n, vel_ma_n, k=3, s=2000.)
plt.plot(test_points, spl(test_points), '-', label='spl s=2000')
plt.legend(loc='lower right')
plt.ylim(0, 300);
In [26]:
star_approx = spl
Прежде, чем оценивать, проведем несколько технических проверок.
Проверка 1: верно ли, что можно использовать дисперсии вдоль большой оси из v_stars_ma.dat вместо s_stars_maN.dat?
Можно, но там точность хуже - лучше использовать вторые данные, ибо они в явном виде даны в статье (таблица 4) и снятые мной точки с графика лучше ложатся на данные s_stars_maN.dat.
Проверка 2: надо ли исправлять малую ось?
Да, надо - см. рисунок 5, там явно растянута малая полуось.
In [27]:
Image('noordermeer_data/n08_disp_ratio.png', width=400)
Out[27]:
Также похоже, что в работе использовано немного другое значение масштаба - около 0.102 (выяснил эмпирически).
Для большой оси: $\sigma^2_{maj} = \sigma^2_{\varphi}\sin^2 i + \sigma^2_{z}\cos^2 i$, следовательно примерные ограничения $$\sigma_{maj} < \frac{\sigma_{maj}}{\sqrt{\sin^2 i + 0.49\cos^2 i}}< \sigma_R = \frac{\sigma_{maj}}{\sqrt{f\sin^2 i + \alpha^2\cos^2 i}} ~< \frac{\sigma_{maj}}{\sqrt{0.5\sin^2 i + 0.09\cos^2 i}} < \frac{\sqrt{2}\sigma_{maj}}{\sin i} (или \frac{\sigma_{maj}}{\sqrt{f}\sin i}),$$ или можно более точную оценку дать, если построить $f$ (сейчас $0.5 < f < 1$).
Для малой оси: $\sigma^2_{min} = \sigma^2_{R}\sin^2 i + \sigma^2_{z}\cos^2 i$ и ограничения $$\sigma_{min} < \frac{\sigma_{min}}{\sqrt{\sin^2 i + 0.49\cos^2 i}} < \sigma_R = \frac{\sigma_{min}}{\sqrt{\sin^2 i + \alpha^2\cos^2 i}} ~< \frac{\sigma_{min}}{\sqrt{\sin^2 i + 0.09\cos^2 i}} < \frac{\sigma_{min}}{\sin i}$$
Соответственно имеем 5 оценок из maj и 4 оценки из min.
In [28]:
# Исправляем значения вдоль малой оси на синус угла:
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))
#Данные Герссена по большой и малой оси (наверное уже раздвинутые)
r_sig_mi_g, sig_mi_g = zip(*np.loadtxt("gerssen_2000_sig_min.dat", float, delimiter=','))
r_sig_ma_g, sig_ma_g = zip(*np.loadtxt("gerssen_2000_sig_maj.dat", float, delimiter=','))
fig = plt.figure(figsize=[11, 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_ma_hsm), sig_ma_hsm, yerr=e_sig_ma_hsm, fmt='.', marker='.', mew=0, color='red', label=r'$HSM99\, maj$')
# plt.errorbar(r_sig_mi, sig_mi, yerr=e_sig_mi, fmt='.', marker='.', mew=0, color='black', label='$\sigma_{los}^{min} Noord$')
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.plot(r_sig_mi_g, sig_mi_g, 'v', color='m', label='$\sigma_{los}^{min} Gerssen$')
plt.plot(r_sig_ma_g, sig_ma_g, 'v', color='g', label='$\sigma_{los}^{maj} Gerssen$')
plt.xlabel('$R$')
plt.ylabel('$\sigma$')
plt.legend();
Видно, что в пределах первых 20 секунд HSM хорошо следует Ноордермееру (дальше тоже, но ошибки больше). Данные Герссена плохо ложатся на малую полуось.
TODO: поискать еще данных
TODO: почему в дипломе наблюдательные точки дисперсий такие маленькие?
In [29]:
spl_min = inter.UnivariateSpline(r_sig_mi, sig_mi, k=3, s=100.)
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 [30]:
@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 [31]:
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 [32]:
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 [33]:
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);
Для настоящей maxmaxtrue почти не отличается от maxmax.
Добавим оценку из Ноордермеера и Герссена, а также сравним major vs minor оценки:
In [34]:
fig = plt.figure(figsize=[10, 7])
r_sig_n, sig_R_n = zip(*np.loadtxt("noordermeer_data/n08_sigR.dat", float, delimiter=','))
plt.plot(map(lambda l: l/0.102, r_sig_n), sig_R_n, 'o-', label='Noord08 sigR, alpha=0.82')
plt.plot(points, map(lambda l: 149.*np.exp(-l/73.), points), '.-', label='Gerssen2000 sigR, alpha=0.85')
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(0,400)
plt.xlim(0,100);
Видно, что оценка Ноордрмеера лежит около minor minmin, тогда как оценка Герссена куда хуже и вообще плоха. Также видно, что почти везде оценки из минимальной дисперсии точнее, чем из максимальной, однако там короче сами данные (почти в два раза). Можно попытаться в центральной части использовать более точные оценки.
TODO: посмотреть дисперсиии тут http://adsabs.harvard.edu/cgi-bin/bib_query?2008MNRAS.390.1089P
TODO: поискать еще данные и добавить статьи
TODO: добавить
In [35]:
fig = plt.figure(figsize=[10,8])
# TODO: проверить как сняты данные
# Noordermeer+2007 ionized gas + HI
r_wsrt, vel_wsrt, e_vel_wsrt = zip(*np.loadtxt("v_gas_WSRT.dat", float))
# Данные по кинематике газа Epinat+2008 в Halpha
r_ha, dr_ha,_,_, vel_ha, e_vel_ha, _,_ = zip(*np.loadtxt("v_gasHa.dat", str))
r_ha, dr_ha, vel_ha, e_vel_ha = np.array(r_ha, dtype='float'), np.array(dr_ha, dtype='float'), np.array(vel_ha, dtype='float'), np.array(e_vel_ha, dtype='float')
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, xerr=dr_ha, fmt='.', marker='d', mew=0, label = 'Halpha')
plt.plot(r, vel, '.-', label = 'Noord thesis')
plt.ylim(0, 300)
plt.legend(loc='lower right');
In [36]:
fig = plt.figure(figsize=[10,8])
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, xerr=dr_ha, fmt='.', marker='d', mew=0, label = r'$H_{\alpha}$')
plt.plot(r, vel, '.-', label = 'Noord thesis')
plt.ylim(0, 300)
plt.xlim(0, 150.)
plt.legend(loc='lower right');
Видно, что я неплохо снял данные из диссертации и что в целом все данные следуют друг другу. Опять, как и в 3898 $H_{\alpha}$ ниже HI.
Приблизим:
In [37]:
fig = plt.figure(figsize=[10,6])
gas_approx = poly1d(polyfit(r_wsrt, vel_wsrt, deg=17))
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.plot(r, vel, '.', label = 'Noord thesis')
plt.ylim(0, 300)
plt.legend(loc='lower right');
In [38]:
fig = plt.figure(figsize=[10,6])
spl_gas = inter.UnivariateSpline([0.]+list(r_wsrt[1:]), [0.]+list(vel_wsrt[1:]), k=3, s=100.)
plt.errorbar(r_wsrt, vel_wsrt, yerr=e_vel_wsrt, fmt='.', marker='.', mew=0, label = 'WSRT')
plt.plot(test_points, spl_gas(test_points), '-', label='spline')
plt.ylim(0, 300)
plt.legend(loc='lower right');
Непонятно конечно, что делать с таким обилием точек как раз в нужной области - по ним выходит, что все не так уж гладко там.
TODO: разобраться
Для случая бесконечного тонкого диска: $$\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 [39]:
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();
Плотность HI:
In [40]:
r_g_dens, gas_dens = zip(*np.loadtxt("gas_density.dat", float))
plt.plot(r_g_dens, gas_dens, 's-');
По $\rm{CO}$ есть профиль интенсивности в работе THE FCRAO EXTRAGALACTIC CO SURVEY. I. THE DATA (1995) https://ui.adsabs.harvard.edu/#abs/1995ApJS...98..219Y/abstract
Ее можно пересчитать в профиль плотности по формуле из van der Hulst 2016 $$\Sigma_{H_2}[M_{\circ}\,{pc}^{−2}] = 3.2\times I_{CO} [K\, km\, s^{−1}]$$
Только надо исправить, т.к. эта формула для $X_{CO} = 1.9$, а еще расстояние там 28.6. Раз есть пк, то значит делится на квадрат расстояния, а значит как обратный квадрат должно зависеть.
In [41]:
Image('2985_CO_data.png')
Out[41]:
Я зря снимал данные с графика, т.к. там есть табличные данные и они не совпадают с картинкой (точки парами должны быть по идее друг над другом). Возьмем среднее между точками:
In [42]:
r_H2_dens_, H2_dens_ = zip(*[(1.134957189032626, 5.274475407021447),
(43.71959002993408, 1.9954912250580392),
(46.072588290621894, 1.421413695030643),
(88.78352654877163, 0.5450392278349678),
(91.457215159636, 0.737871470054718)])
plt.plot(r_H2_dens_, H2_dens_, 's-')
r_H2_dens, H2_dens = zip(*[(0.0, 5.25),
(45., 1.67),
(90., 0.57)])
plt.plot(r_H2_dens, H2_dens, 's-')
plt.xlim(0, 120.);
In [43]:
H2_dens = [3.2*(X_CO/1.9)*(28.6/distance)**2 * l for l in H2_dens]
In [44]:
plt.plot(r_H2_dens, H2_dens, 'o-')
plt.plot(r_g_dens, gas_dens, 's-');
Вполне похоже на правду.
In [45]:
import scipy.integrate as integrate
import scipy.interpolate
tmp_ = scipy.interpolate.interp1d(r_H2_dens, H2_dens)
result = integrate.quad(lambda l: 2*np.pi*l*tmp_(l), r_H2_dens[0], r_H2_dens[-1])
print (scale * 1000.)**2 * result[0]/1e9
Диплом: B, R - маленькие, около 60, макс. диск ~ 450 (M/L=6), другие R и J - больше 1000, похоже два диска
In [46]:
all_photometry = []
In [47]:
from wand.image import Image as WImage
img = WImage(filename='ngc2985.pdf', resolution=200)
img[:, 150:1800]
Out[47]:
In [48]:
img = WImage(filename='ngc2985.pdf[1]', resolution=200)
img[:, 150:1300]
Out[48]:
In [49]:
all_photometry = []
Фотометрия Ноордермеера:
In [50]:
Image('noordermeer_data/n2985_photom.png')
Out[50]:
Снятые R-данные:
In [51]:
r_phot2, r_band = zip(*np.loadtxt('Rband_ngc2985.dat', float))
plt.plot(r_phot2, r_band, 's')
plt.ylim(30, 15)
Out[51]:
Снятые вместе данные:
In [52]:
r_phot, mu_phot = zip(*np.loadtxt('noordermeer_data/n2985_noord_photoRB.dat', float, delimiter=','))
plt.plot(r_phot, mu_phot, 's')
plt.plot(r_phot2, r_band, '.')
plt.ylim(30, 15);
Согласуются.
In [53]:
M_R = -20.89 #10.80 - это правильно? надо брать абсолютные? в дипломе были относительные, тут разница уже существенная
M_B = -20.31 #11.43
In [54]:
print 'Abs B : {:2.2f}; R: {:2.2f}.'.format(bell_mass_to_light(M_B-M_R, 'B', 'B-R'), bell_mass_to_light(M_B-M_R, 'R', 'B-R'))
print 'Rel B : {:2.2f}; R: {:2.2f}.'.format(bell_mass_to_light(11.43-10.80, 'B', 'B-R'), bell_mass_to_light(11.43-10.80, 'R', 'B-R'))
Тут разница не очень большая, что именно брать - и то и то маленькое.
In [55]:
# R-band
r_eff_R = 25.1
mu_eff_R = 20.48 # уточнить это ли число
n_R = 3.9
mu0d_R = 21.16 # и тут тоже
h_disc_R = 52.2
mu_eff_Rc = 20.41 # уточнить это ли число
mu0d_Rc = 21.32 # и тут тоже
In [56]:
# B-band
r_eff_B = 25.1
mu_eff_B = 21.98 # уточнить это ли число
n_B = 3.9
mu0d_B = 21.95 # и тут тоже
h_disc_B = 57.6
mu_eff_Bc = 21.86 # уточнить это ли число
mu0d_Bc = 22.06 # и тут тоже
In [57]:
p_ = np.arange(0., 300., 0.1)
fig = plt.figure(figsize=[8, 5])
plt.plot(r_phot, mu_phot, 's')
plt.plot(p_, total_mu_profile([mu_bulge(l, mu_eff=mu_eff_R, r_eff=r_eff_R, n=n_R) for l in p_],
[mu_disc(l, mu0=mu0d_R, h=h_disc_R) for l in p_]), '-', label='sum R', color='#007700')
plt.plot(p_, total_mu_profile([mu_bulge(l, mu_eff=mu_eff_B, r_eff=r_eff_B, n=n_B) for l in p_],
[mu_disc(l, mu0=mu0d_B, h=h_disc_B) for l in p_]), '-', label='sum B', color='#000077')
plt.axvline(x = h_disc_B*(26.5 - mu0d_B)/1.0857, ls='--', alpha=0.5) # как указать расстояние, соответствующее 26.5m
plt.axhline(y=26.5, ls='--', alpha=0.5)
plt.xlim(-10, 300)
plt.ylim(30, 15)
plt.legend();
Отлично, похоже на правду. Массовые модели:
In [58]:
b_r_color = M_B-M_R
M_to_L_R = bell_mass_to_light(b_r_color, 'R', 'B-R')
surf_R = [surf_density(mu=mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L=M_to_L_R, band='R') for l in p_]
plt.plot(p_, surf_R, '-', label='R [M/L={:2.2f}]'.format(M_to_L_R))
M_to_L_B = bell_mass_to_light(b_r_color, 'B', 'B-R')
surf_B = [surf_density(mu=mu_disc(l, mu0=mu0d_Bc, h=h_disc_B), M_to_L=M_to_L_B, band='B') for l in p_]
plt.plot(p_, surf_B, '-', label='B [M/L={:2.2f}]'.format(M_to_L_B))
diploma_color = 11.43-10.80
M_to_L_Rc = bell_mass_to_light(diploma_color, 'R', 'B-R')
surf_R = [surf_density(mu=mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L=M_to_L_Rc, band='R') for l in p_]
plt.plot(p_, surf_R, '-', label='R [M/L={:2.2f}]'.format(M_to_L_R))
M_to_L_Bc = bell_mass_to_light(diploma_color, 'B', 'B-R')
surf_B = [surf_density(mu=mu_disc(l, mu0=mu0d_Bc, h=h_disc_B), M_to_L=M_to_L_Bc, band='B') for l in p_]
plt.plot(p_, surf_B, '-', label='B [M/L={:2.2f}]'.format(M_to_L_B))
print b_r_color, diploma_color
plt.legend();
Тут все как в дипломе, очень маленькие значения получились. Разница в цвете никакая, можно не учитывать.
In [59]:
all_photometry.append((BAD_MODEL_PREFIX+'Noorder R', r_eff_R, mu_eff_Rc, n_R, mu0d_Rc, h_disc_R, M_to_L_Rc,
lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L=M_to_L_Rc, band='R')))
all_photometry.append((BAD_MODEL_PREFIX+'Noorder B', r_eff_B, mu_eff_Bc, n_B, mu0d_Bc, h_disc_B, M_to_L_Bc,
lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_Bc, h=h_disc_B), M_to_L=M_to_L_Bc, band='B')))
Méndez-Abreu фотометрия в $J$:
In [60]:
# J-band
r_eff_J = 13.2
mu_eff_J = 17.94
n_J = 2.92
mu0d_J = 18.22
h_disc_J = 25.8
In [61]:
M_to_L_J = bell_mass_to_light(b_r_color, 'J', 'B-R')
surf_J = [surf_density(mu=mu_disc(l, mu0=mu0d_J, h=h_disc_J), M_to_L=M_to_L_J, band='J') for l in p_]
plt.plot(p_, surf_J, '-', label='J [M/L={:2.2f}]'.format(M_to_L_J))
plt.legend();
Получилось не то, что в дипломе - надо понять, почему.
TODO: понять почему отличается (M/L в дипломе другое)
In [62]:
all_photometry.append((BAD_MODEL_PREFIX+'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$ (в пдф-ке неверно записано J), модель достаточно странная - ибо нет балджа.
In [63]:
h_in=18.1
h_out=81.0
h_brk=69.
mu_in=18.57
mu_out=21.76
mu_brk=21.8
In [64]:
Image('gutierrez_Rphotom.png', width=500)
Out[64]:
In [65]:
fig = plt.figure(figsize=[5, 5])
_, _, r_phot2, mu_phot2 = zip(*np.loadtxt("gutierrez_Rphotom.dat", str)) #данные из таблицы в Vizier
plt.plot(r_phot2, mu_phot2, '.')
plt.plot(p_, [mu_disc(l, mu0=mu_in, h=h_in) for l in p_], '--', label='in', color='#007700')
plt.plot(p_, [mu_disc(l, mu0=mu_out, h=h_out) for l in p_], '--', label='out', color='#000077')
plt.plot(p_, total_mu_profile([mu_disc(l, mu0=mu_in, h=h_in) for l in p_],
[mu_disc(l, mu0=mu_out, h=h_out) for l in p_]), '-', label='sum', color='#770000')
plt.xlim(0., 320)
plt.ylim(26.5, 16)
plt.legend()
plt.axvline(x = h_brk, linestyle='-.');
Похоже. Возьмем оба диска:
In [66]:
surf_R2 = [surf_density(mu=m, M_to_L=M_to_L_R, band='R') for m in
total_mu_profile([mu_disc(l, mu0=mu_in, h=h_in) for l in p_],
[mu_disc(l, mu0=mu_out, h=h_out) for l in p_])]
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=mu_in, h=h_in), 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)) #только внутренний диск
plt.plot(p_, np.array([surf_density(mu=mu_disc(l, mu0=mu_in, h=h_in), M_to_L=M_to_L_R, band='R') for l in p_])+
np.array([surf_density(mu=mu_disc(l, mu0=mu_out, h=h_out), M_to_L=M_to_L_R, band='R') for l in p_]), '--') #проверяю, что одно и то же получается
plt.legend();
Опять отличается, не ясно почему (M/L на этот раз похоже), хотя все вместе вполне нормальное.
TODO: понять почему
In [67]:
all_photometry.append(('Gutierrez R 2d', None, None, None, (mu_in, mu_out), (h_in, h_out), M_to_L_R,
(lambda l: surf_density(mu=mu_disc(l, mu0=mu_in, h=h_in), M_to_L=M_to_L_R, band='R'),
lambda l: surf_density(mu=mu_disc(l, mu0=mu_out, h=h_out), M_to_L=M_to_L_R, band='R'))))
Декомпозиция в Sofue 2016 из кривой вращения:
In [68]:
M_bulge = 0.39 #± 0.07
a_b = 0.56 #± 0.14 kpc
M_disk = 1.1 #± 0.1
a_d = 1.1 #± 0.2 kpc
Для диска верно: $M_d = 2\pi a_d^2 \Sigma_0$
In [69]:
M_disk * 10.**10 / (2*np.pi) / (a_d * 1000.)**2
Out[69]:
In [70]:
Sigma_0_sofue = M_disk * 10.**10 / (2*np.pi) / (a_d * 1000.)**2
surf_sofue = [Sigma_0_sofue*np.exp(-l*scale/a_d) for l in p_]
plt.plot(p_, surf_sofue, '-', label='Sofue')
plt.legend();
Наконец, JHK из работы Heidt 2001 http://www.aanda.org/articles/aa/pdf/2001/10/aa10227.pdf:
In [71]:
mudJ, mudH, mudK = 18.05, 17.27, 17.32
hJ, hH, hK = 30.63, 26.10, 31.08
mueJ, mueH, mueK = 17.73, 17.03, 17.23
reJ, reH, reK = 12.08, 12.81, 14.06
nJ, nH, nK = 1./0.35, 1./0.35, 1./0.33
In [72]:
M_to_L_J = bell_mass_to_light(b_r_color, 'J', 'B-R')
surf_J2 = [surf_density(mu=mu_disc(l, mu0=mudJ, h=hJ), M_to_L=M_to_L_J, band='J') for l in p_]
plt.plot(p_, surf_J2, '-', label='J2 [M/L={:2.2f}]'.format(M_to_L_J))
M_to_L_H = bell_mass_to_light(b_r_color, 'H', 'B-R')
surf_H = [surf_density(mu=mu_disc(l, mu0=mudH, h=hH), M_to_L=M_to_L_H, band='H') for l in p_]
plt.plot(p_, surf_H, '-', label='H [M/L={:2.2f}]'.format(M_to_L_H))
M_to_L_K = bell_mass_to_light(b_r_color, 'K', 'B-R')
surf_K = [surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=M_to_L_K, band='K') for l in p_]
plt.plot(p_, surf_K, '-', label='K [M/L={:2.2f}]'.format(M_to_L_K))
plt.legend();
На удивление согласуется друг с другом.
In [73]:
all_photometry.append(('Heidt J', reJ, mueJ, nJ, mudJ, hJ, M_to_L_J,
lambda l: surf_density(mu=mu_disc(l, mu0=mudJ, h=hJ), M_to_L=M_to_L_J, band='J')))
all_photometry.append(('Heidt K', reK, mueK, nK, mudK, hK, M_to_L_K,
lambda l: surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=M_to_L_K, band='K')))
all_photometry.append(('Heidt H', reH, mueH, nH, mudH, hH, M_to_L_H,
lambda l: surf_density(mu=mu_disc(l, mu0=mudH, h=hH), M_to_L=M_to_L_H, band='H')))
В работе Герссена 2000го года есть фотометрия в $I$ и декомпозиция на диск-балдж, но там как-то произвольно выбраны величины на оси:
In [74]:
Image('gerssen_I_photom.png')
Out[74]:
Единственное что понятно - что $h_{disc} = 30 \pm 4$.
S4G данные из GALFIT (дисков два, непонятно в чем отличие):
In [75]:
r_eff_s4g = 6.32
# mu_eff_s4g = ...
n_s4g = 2.822
mu0d_s4g = 18.552
h_disc_s4g = 12.78
mu0d_s4g_2 = 20.845
h_disc_s4g_2 = 48.89
Тут нужно учитывать, что эти параметры в AB-mag и нуждаются в доп. исправлении.
In [76]:
M_to_L_s4g = s4g_mass_to_light(-21.746, -21.276)
M_to_L_s4g
Out[76]:
In [77]:
surf_s4g = [s4g_surf_density(mu_disc(l, mu0=mu0d_s4g, h=h_disc_s4g), M_to_L_s4g) for l in p_]
plt.plot(p_, surf_s4g, '-', label='S4G [M/L={:2.2f}]'.format(M_to_L_s4g))
plt.plot(p_, [s4g_surf_density(mu_disc(l, mu0=mu0d_s4g_2, h=h_disc_s4g_2), M_to_L_s4g) for l in p_], '-')
plt.legend();
Достаточно массивный.
In [78]:
all_photometry.append(('S4G 2d', r_eff_s4g, None, n_s4g, (mu0d_s4g, mu0d_s4g_2), (h_disc_s4g, h_disc_s4g_2), M_to_L_s4g,
(lambda l: s4g_surf_density(mu_disc(l, mu0=mu0d_s4g, h=h_disc_s4g), M_to_L_s4g),
lambda l: s4g_surf_density(mu_disc(l, mu0=mu0d_s4g_2, h=h_disc_s4g_2), M_to_L_s4g))))
Cаша Мосенков прогнал GALFIT декомпозиию для двух моделей - с двумя дисками и с одним. Результаты (пиксели в секунды переводятся домножением на 0.75, центр. поверх яркость по m0d = mag + 2.5log10(2.math.pi*h^2)):
In [79]:
h_AM = 31.4388 * 0.75
mu0d_AM = 12.0660 + 2.5*np.log10(2.*math.pi*h_AM**2)
h_AM, mu0d_AM
Out[79]:
In [80]:
surf_s4g_AM = [s4g_surf_density(mu_disc(l, mu0=mu0d_AM, h=h_AM), M_to_L_s4g) for l in p_]
plt.plot(p_, surf_s4g_AM, '-', label='S4G AM [M/L={:2.2f}]'.format(0.66))
plt.legend();
Слишком маленькая.
Двухдисковая:
In [81]:
h_AM_in = 15.6155 * 0.75
mu0d_AM_in = 10.9384 + 2.5*np.log10(2.*math.pi*h_AM_in**2)
h_AM_out = 62.9739 * 0.75
mu0d_AM_out = 10.3740 + 2.5*np.log10(2.*math.pi*h_AM_out**2)
h_AM_in, mu0d_AM_in, h_AM_out, mu0d_AM_out
Out[81]:
In [82]:
surf_s4g_AM_in = [s4g_surf_density(mu=mu_disc(l, mu0=mu0d_AM_in, h=h_AM_in), M_to_L=M_to_L_s4g) for l in p_]
plt.plot(p_, surf_s4g_AM_in, '-', label='S4G AM [M/L={:2.2f}]'.format(0.66))
surf_s4g_AM_out = [s4g_surf_density(mu=mu_disc(l, mu0=mu0d_AM_out, h=h_AM_out), M_to_L=M_to_L_s4g) for l in p_]
plt.plot(p_, surf_s4g_AM_out, '-', label='S4G AM [M/L={:2.2f}]'.format(0.66))
plt.legend();
Как видим - поменялось не так уж сильно, раза в полтора, но итог очень тяжеловесный.
In [83]:
all_photometry.append(('S4G_AM 2d', None, None, None, (mu0d_AM_in, mu0d_AM_out), (h_AM_in, h_AM_out), M_to_L_s4g,
(lambda l: s4g_surf_density(mu_disc(l, mu0=mu0d_AM_in, h=h_AM_in), M_to_L_s4g),
lambda l: s4g_surf_density(mu_disc(l, mu0=mu0d_AM_out, h=h_AM_out), M_to_L_s4g))))
Финальная сводная картинка:
In [84]:
fig = plt.figure(figsize=[8, 5])
plt.plot(p_, surf_R, '-', label='R [M/L={:2.2f}]'.format(M_to_L_R))
plt.plot(p_, surf_B, '-', label='B [M/L={:2.2f}]'.format(M_to_L_B))
plt.plot(p_, surf_R2, '-', label='R2 [M/L={:2.2f}]'.format(M_to_L_R))
plt.plot(p_, surf_J, '-', label='J [M/L={:2.2f}]'.format(M_to_L_J))
plt.plot(p_, surf_sofue, '-', label='Sofue')
plt.plot(p_, surf_J2, '-', label='J2 [M/L={:2.2f}]'.format(M_to_L_J))
plt.plot(p_, surf_H, '-', label='H [M/L={:2.2f}]'.format(M_to_L_H))
plt.plot(p_, surf_K, '-', label='K [M/L={:2.2f}]'.format(M_to_L_K))
plt.plot(p_, surf_s4g, '-', label='3.6 [M/L={:2.2f}]'.format(M_to_L_s4g))
plt.plot(p_, surf_s4g_AM_in, '-', label='3.6 AM [M/L={:2.2f}]'.format(M_to_L_s4g))
#maximal disk
plt.plot(p_, [surf_density(mu=mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L=6., band='R') for l in p_], '-', label='R [M/L={:2.2f}]'.format(6.))
plt.legend()
plt.xlim(0, 150.);
In [85]:
all_photometry.append(('Noorder R_max', r_eff_R, mu_eff_Rc, n_R, mu0d_Rc, h_disc_R, 6.,
lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L=6., band='R')))
Видно, что фотометрия Ноордрмеера совсем какая-то левая, пять остальных примерно согласуются со значениями в 500-800 и есть две очень большие, одна из которых S4G (а одна совсем большая).
На самом деле во всем этом еще важна скорость падения! Видно, что в интересующем нас диапазоне эти фотметрии по сути разбиваются на две группы.
In [86]:
show_all_photometry_table(all_photometry, scale)
Можно провести тест-сравнение с кривой вращения тонкого диска при заданной фотометрии, если она слишком массивная - то не брать ее (это ограничение сверху).
In [87]:
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')
plt.plot(r, vel, '.', label = 'Noord thesis')
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');
Видно, что все двухдисковые модели похожи друг на друга, но тут опять же сложно выявить максимум.
Hameed & Devereux 2005 http://iopscience.iop.org/article/10.1086/430211/pdf $H_{\alpha}$ (расстояние у него 22.4 Мпк, тогда как у нас scale для ~19 Мпк):
TODO: поискать еще данных и измерить корректно размер
In [88]:
Image('ngc2985_Halpha.jpg')
Out[88]:
Вручную примерно так:
In [89]:
Image('ngc2985_Halpha_dist.jpg')
Out[89]:
XDSS и $H_{alpha}$ из http://adsabs.harvard.edu/cgi-bin/bib_query?2008MNRAS.390..466E (расстояние 21.1 Мпк)
In [90]:
Image('halpha_xdss_u5253.png')
Out[90]:
На кривой вращения видно, что светятся области 15-30 и 10-45+(до 70).
И кривая вращения из той же работы что и выше (обзор GHASP):
In [91]:
Image('ghasp_2008_halpha_vel.png')
Out[91]:
In [92]:
def plot_SF(ax):
ax.plot([10., 70.], [0., 0.], '-', lw=7., color='b')
# ax.plot([10., 7.2/scale], [0., 0.], '-', lw=7., color='red')
ax.plot([10., 7.2/(scale*22.4/distance)], [0., 0.], '-', lw=7., color='r') #TODO: исправить менее грубо
plot_SF(plt.gca())
plt.xlim(0, 350)
plt.ylim(0, 200)
Out[92]:
Согласуются примерно данные (хотя по яркости кривой вращения в Halpha так и не скажешь).
Изображение из AINUR: Atlas of Images of NUclear Rings (DSS изображение из http://adsabs.harvard.edu/cgi-bin/bib_query?2010MNRAS.402.2462C)- но это самый центр, не нужно
In [93]:
Image('AINUR_image.png')
Out[93]:
In [94]:
Image('diplom_results.png')
Out[94]:
Устойчиво, когда > 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 [95]:
sound_vel = 6 #скорость звука в газе, км/с
data_lim = min(max(r_sig_ma), max(r_wsrt)) #где заканчиваются данные
In [96]:
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(spl_gas, 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=4., gas_density=l[1]) for l in
zip([epicyclicFreq_real(spl_gas, x, scale) for x in r_g_dens], gas_dens)], 's-', label='$Q_{gas}$ & c=4')
plt.plot(r_g_dens, [1./Qg(epicycl=l[0], sound_vel=sound_vel, gas_density=He_coeff*l[1]) for l in
zip([epicyclicFreq_real(spl_gas, x, scale) for x in r_g_dens], gas_dens)], 's-', label='$Q_{gas}$ & He_coeff')
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(l_, M_to_L_K, 'K') for l_ in [mu_disc(ll, mu0=mudK, h=hK) for ll in r_g_dens]])], 's-', label='$Q_{star}^{min}$')
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(l_, M_to_L_K, 'K') for l_ in [mu_disc(ll, mu0=mudK, h=hK) for ll in r_g_dens]])], 's-', label='$Q_{star}^{max}$')
plt.axhline(y=1, ls='--')
plt.legend()
plot_SF(plt.gca())
plt.ylabel('$Q^{-1}$', fontsize=15);
Значения одножидкостного похожи на диплом (без 1.6).
НЕ ИСПРАВЛЕНО ЗА 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 [97]:
total_gas_data = zip(r_g_dens, map(lambda l: He_coeff*l, gas_dens))[:5]
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_Rc, h=h_disc_R), 6., 'R'),
star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu_in, h=h_in), M_to_L_R, 'R'),
data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='R/R_max maj max/min')
plt.ylim(0., 2.5)
plt.axhline(y=1., ls='-', color='grey')
plot_SF(ax)
plt.grid();
In [98]:
from matplotlib.animation import FuncAnimation
fig = plt.gcf()
plt.figure(figsize=(10,6))
ax = plt.gca()
def animate(i):
ax.cla()
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 [99]:
# anim.save('..\\..\pics\\'+name+'.gif', writer='imagemagick', fps=1)
In [100]:
# from IPython.display import HTML
# HTML(anim.to_html5_video())
Существует ограничение на максимальный диск в ~0.85 (изотермическое гало) и на субмаксимальный в 0.55-0.6 (NFW гало). Попробуем дотянуть фотметрию до максимальных дисков и посмотрим, как изменятся M/L (скорость зависит как корень из M/L):
In [101]:
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.));
Для R многовато, остальные более-менее. Видно, что по сути разбивается на две группы - Ноордермеер vs все остальные.
In [102]:
from matplotlib.animation import FuncAnimation
fig = plt.gcf()
plt.figure(figsize=(10,6))
ax = plt.gca()
def animate(i):
ax.cla()
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)
plot_SF(ax)
return ax
anim = FuncAnimation(plt.gcf(), animate, frames=len(all_photometry), interval=1000)
In [103]:
# anim.save('..\\..\pics\\'+name+'_MAXDISCS.gif', writer='imagemagick', fps=1)
In [104]:
# from IPython.display import HTML
# HTML(anim.to_html5_video())
In [105]:
import scipy.interpolate
fig = plt.figure(figsize=[10, 4])
y_interp = scipy.interpolate.interp1d(list(r_H2_dens)[:-1] + [91.], list(H2_dens)[:-1] + [0.0])
def y_interp_(r):
if r < min(r_H2_dens):
return y_interp(min(r_H2_dens))
elif r < 91:
return y_interp(r)
else:
return 0.
points = np.linspace(1.2, 200., 100.)
plt.plot(r_H2_dens, H2_dens, 'o-')
plt.plot(r_g_dens, gas_dens, 's-')
plt.plot(points, map(y_interp_, points), '--', alpha=0.5)
plt.plot(r_g_dens, map(lambda l: y_interp_(l[0]) + l[1], zip(r_g_dens, gas_dens)), '--', alpha=0.5)
plt.ylim(0);
In [106]:
fig = plt.figure(figsize=[10, 6])
total_gas_data_ = zip(r_g_dens, [He_coeff*(y_interp_(l[0]) + l[1]) for l in zip(r_g_dens, gas_dens)])[1: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_Rc, h=h_disc_R), 6., 'R'),
star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), 6., 'R'),
data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='R_max AD disp')
plt.ylim(0., 2.5)
plt.axhline(y=1., ls='-', color='grey')
plot_SF(ax)
plt.grid();
Стало более полого, но все еще недостаточно.
Для максимальных дисков:
In [107]:
from matplotlib.animation import FuncAnimation
fig = plt.gcf()
plt.figure(figsize=(10,6))
ax = plt.gca()
def animate(i):
ax.cla()
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)
plot_SF(ax)
return ax
anim = FuncAnimation(plt.gcf(), animate, frames=len(all_photometry), interval=1000)
In [108]:
# anim.save('..\\..\pics\\'+name+'_MAXDISCS_all_gas.gif', writer='imagemagick', fps=1)
In [109]:
# from IPython.display import HTML
# HTML(anim.to_html5_video())
Видно, что в смысле звездообразования лучше всего фотометрия Heidt подходит.
In [110]:
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('HST_color_dist.jpg'), aspect='auto', cmap='Greys')
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), next(linecycler), 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]) + l[1]) for l in zip(r_g_dens, gas_dens)], '*-')
axes[3].plot(r_g_dens, [y_interp_(l) for l in r_g_dens], '--')
axes[3].set_title('Gas')
axes[3].grid()
axes[3].set_xlim(0, 200)
axes[3].legend()
#change this
@save_model(models_path+'n2985_modelKmax.npy')
def plot_2f_vs_1f_(*args, **kwargs):
plot_2f_vs_1f(*args, **kwargs)
plot_2f_vs_1f_(ax=axes[4], total_gas_data=total_gas_data_,
epicycl=epicyclicFreq_real,
gas_approx=spl_gas,
sound_vel=sound_vel,
scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.27, band='K'),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.27, band='K'),
data_lim=data_lim, color='g', alpha=0.2, disk_scales=disk_scales, label='K Heidt maxdisc',
Ml = 1.27,
CO = lambda l: y_interp_(l))
@save_model(models_path+'n2985_model36max.npy')
def plot_2f_vs_1f_(*args, **kwargs):
plot_2f_vs_1f(*args, **kwargs)
plot_2f_vs_1f_(ax=axes[4], total_gas_data=total_gas_data_,
epicycl=epicyclicFreq_real,
gas_approx=spl_gas,
sound_vel=sound_vel,
scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu0d_s4g, h=h_disc_s4g), 1.23),
lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu0d_s4g_2, h=h_disc_s4g_2), 1.23)))(l),
star_density_min=lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu0d_s4g, h=h_disc_s4g), 1.23),
lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu0d_s4g_2, h=h_disc_s4g_2), 1.23)))(l),
data_lim=data_lim, color='y', alpha=0.2, disk_scales=disk_scales, label='S4G 2d maxdisc',
ML = 1.21,
CO = lambda l: y_interp_(l))
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 [109]:
plt.plot(zip(*total_gas_data)[0], zip(*total_gas_data)[1], 'o-')
for photom in all_photometry:
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, 50.);
Видимо везде неустойчиво.
Hunter et al (1998), 'competition with shear' according to Leroy: $$\Sigma_A = \alpha_A\frac{\sigma_g A}{\pi G}$$
In [110]:
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, 50.)
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, 50.);
Похоже неустойчиво почти всюду.
Интересный вариант для тех галактик, в которых есть данные по газу. Разница между скоростями вращения звезд и газа вокруг центра галактики называется ассиметричным сдвигом и описывается следующим уравнением (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 [111]:
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 [112]:
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)
ad_interp = scipy.interpolate.interp1d(r_sig_ma, np.sqrt(sigR2))
@flat_end(sig_maj_lim)
def ad_interp_(r):
return ad_interp(r)
plt.plot(points[2:], map(ad_interp_, points[2:]), '--', alpha=0.5)
plt.legend()
plt.ylim(0,450)
plt.xlim(0,100);
In [113]:
fig = plt.figure(figsize=[10, 6])
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=ad_interp_,
sigma_min=ad_interp_,
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.42, band='K'),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.42, band='K'),
data_lim=data_lim, color='g', alpha=0.2, disk_scales=disk_scales, label='K Heidt maxdisc AD')
plt.ylim(0., 2.5)
plt.axhline(y=1., ls='-', color='grey')
plot_SF(ax)
plt.grid();
Не сильно отличается, что вообще-то ожидаемо.
Если брать экспонентой:
In [114]:
fig = plt.figure(figsize=[10, 6])
ax = plt.gca()
plot_2f_vs_1f(ax=ax, total_gas_data=zip(r_g_dens[1:], [He_coeff*(y_interp_(l[0]) + 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=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.42, band='K'),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.42, band='K'),
data_lim=data_lim, color='g', alpha=0.2, disk_scales=disk_scales, label='K Heidt maxdisc AD')
plt.xlim(0., 500.)
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 [115]:
r25 = h_disc_B*(25. - mu0d_B)/1.0857
r25, h_disc_B, (25. - mu0d_B)/1.0857
Out[115]:
In [116]:
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[7:], gas_dens[7:])
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[7:], gas_dens[7:])
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 [117]:
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 [118]:
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-');
Но тут надо иметь в виду, что возможно у Ноордермеера неадекватные оценки диска и $r_{25}$ вполне может быть раза в два меньше.
Из https://ui.adsabs.harvard.edu/#abs/2016MNRAS.460.1106W/abstract: два возможных вида связи между молекулярным и атомарным газом $R_{mol} = \Sigma_{H_2}/\Sigma_{HI}$:
$$R_{mol} = \Sigma_{star}/81$$или $$R_{mol} = \left(\frac{P_h}{1.7 \times 10^4 cm^{-3}K k_B } \right)^{0.8},\, P_h = \frac{\pi}{2}G\Sigma_g(\Sigma_g + \frac{\sigma_g}{\sigma_z}\Sigma_{star})$$
In [119]:
def R1(Sigma_star):
return Sigma_star/81.
def h2_gas(r, h_gas_dens):
return R1(star_density(r))*h_gas_dens
star_density=lambda l: surf_density(mu_disc(l, mu0=mudK, h=hK), M_to_L_K, 'K')
plt.plot(r_g_dens, gas_dens, 's-', color='b')
plt.plot(r_g_dens, [h2_gas(_[0], _[1]) for _ in zip(r_g_dens, gas_dens)], 's-', color='r')
plt.plot(r_H2_dens, H2_dens, 'o-', alpha=0.6);
Вторая оценка:
In [120]:
spl_min = inter.UnivariateSpline(r_sig_mi, sig_mi, k=3, s=100.)
def sig_z(r, alpha):
if r < max(r_sig_mi):
return spl_min(r)/sqrt(1/alpha**2 * sin_i**2 + cos_i**2)
else:
return spl_min(max(r_sig_mi))/sqrt(1/alpha**2 * sin_i**2 + cos_i**2)
plt.errorbar(r_sig_mi, sig_mi, yerr=e_sig_mi, fmt='.', marker='.', mew=0, color='red', label='$\sigma_{los}^{min}$')
plt.plot(points, spl_min(points), label = '$\sigma_{los}^{min}\, splinefit$', color='red')
plt.plot(points, map(lambda l: sig_z(l, 0.3), points), label = '0.3')
plt.plot(points, map(lambda l: sig_z(l, 0.5), points), label = '0.5')
plt.plot(points, map(lambda l: sig_z(l, 0.7), points), label = '0.7')
plt.legend()
plt.ylim(0,400)
plt.xlim(0,100);
In [121]:
def R2(r, h_gas_dens, alpha, sound_vel):
G = 6.67408
kB = 3.7529917
Ph = np.pi/2. * G * h_gas_dens * (h_gas_dens + sound_vel/sig_z(r, alpha) * star_density(r))
return np.power(4.363474*Ph/(1.7 * 10000. * kB) , 0.8)
def h2_gas2(r, h_gas_dens, alpha, sound_vel):
return R2(r, h_gas_dens, alpha, sound_vel)*h_gas_dens
plt.plot(r_g_dens[:4], [R2(_[0], _[1], 0.3, 6.) for _ in zip(r_g_dens, gas_dens)[:4]], 's-')
plt.plot(r_g_dens[:4], [R2(_[0], _[1], 0.5, 6.) for _ in zip(r_g_dens, gas_dens)[:4]], 's-')
plt.plot(r_g_dens[:4], [R2(_[0], _[1], 0.7, 6.) for _ in zip(r_g_dens, gas_dens)[:4]], 's-');
In [122]:
plt.plot(r_g_dens, gas_dens, 's-', color='b')
plt.plot(r_g_dens, [h2_gas2(_[0], _[1], 0.3, 6.) for _ in zip(r_g_dens, gas_dens)], 's-')
plt.plot(r_g_dens, [h2_gas2(_[0], _[1], 0.5, 6.) for _ in zip(r_g_dens, gas_dens)], 's-')
plt.plot(r_g_dens, [h2_gas2(_[0], _[1], 0.7, 6.) for _ in zip(r_g_dens, gas_dens)], 's-')
plt.plot(r_g_dens, map(lambda l: 0.44*l, gas_dens), '-', color='m');
И теперь сравнение с настоящим значением:
In [123]:
fig = plt.figure(figsize=[10, 10])
# plt.plot(r_g_dens, gas_dens, 's-', color='b')
plt.semilogy(r_g_dens, [h2_gas2(_[0], _[1], 0.3, 6.) for _ in zip(r_g_dens, gas_dens)], '-', label='0.3')
plt.semilogy(r_g_dens, [h2_gas2(_[0], _[1], 0.5, 6.) for _ in zip(r_g_dens, gas_dens)], '-', label='0.5')
plt.semilogy(r_g_dens, [h2_gas2(_[0], _[1], 0.7, 6.) for _ in zip(r_g_dens, gas_dens)], '-', label='0.7')
# plt.semilogy(r_g_dens, map(lambda l: 0.44*l, gas_dens), '-', color='m', label='x0.44')
plt.semilogy(r_H2_dens, H2_dens, 'o-', color='r')
plt.semilogy(r_g_dens, [h2_gas(_[0], _[1]) for _ in zip(r_g_dens, gas_dens)], '-', label='/81')
plt.legend()
plt.ylim(0.1, 100.)
plt.xlim(0, 100);
На этот раз первая модель лучше всего, остальные слишком маленькие.
Тут учитывается толщина диска:
In [354]:
plot_RF13_vs_2F(r_g_dens=r_g_dens, HI_gas_dens=gas_dens, CO_gas_dens=[y_interp_(l) 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=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.42, band='K'),
alpha_max=0.7, alpha_min=0.3, scale=scale, gas_approx=spl_gas, thin=False)
А тут нет:
In [110]:
plot_RF13_vs_2F(r_g_dens=r_g_dens, HI_gas_dens=gas_dens, CO_gas_dens=[y_interp_(l) 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=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.27, band='K'),
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');
Видно, что согласие достаточно хорошее.
Но при этом выясняется интересная и достаточно занимательная вещь - так как теперь две версии, с которыми можно сравнивать - RW2 и RW3 (двух- и трех-компонентные). Так вот изменим совсем немного COdensity в одной можно получить на 30% меньше значение $Q{RW}$:
In [126]:
print romeo_Qinv2(r=60., epicycl=56.367905714, sound_vel=6., sigma_R=77.1464789152, star_density=119.826043553,
HI_density=4., CO_density=6.6, alpha=0.3, verbose=True, thin=True, He_corr=False)
print('======')
print romeo_Qinv2(r=60., epicycl=56.367905714, sound_vel=6., sigma_R=77.1464789152, star_density=119.826043553,
HI_density=4., CO_density=7.0, alpha=0.3, verbose=True, thin=True, He_corr=False)
print('======')
print romeo_Qinv(r=60., epicycl=56.367905714, sound_vel=6., sigma_R=77.1464789152, star_density=119.826043553,
HI_density=4., CO_density=6.6, alpha=0.3, verbose=True, thin=True, He_corr=False)
print('======')
# change CO_density a little
print '!!!=>',romeo_Qinv(r=60., epicycl=56.367905714, sound_vel=6., sigma_R=77.1464789152, star_density=119.826043553,
HI_density=4., CO_density=7.0, alpha=0.3, verbose=True, thin=True, He_corr=False)
print('======')
print romeo_Qinv(r=60., epicycl=56.367905714, sound_vel_CO=6., sound_vel_HI=11., sigma_R=77.1464789152, star_density=119.826043553,
HI_density=4., CO_density=6.6, alpha=0.3, verbose=True, thin=True, He_corr=False)
print('======')
print romeo_Qinv(r=60., epicycl=56.367905714, sound_vel_CO=6., sound_vel_HI=11., sigma_R=77.1464789152, star_density=119.826043553,
HI_density=4., CO_density=7.0, alpha=0.3, verbose=True, thin=True, He_corr=False)
In [123]:
fig = plt.figure(figsize=[9,9])
ax = plt.gca()
plt.plot([6./77.1464789152], [1.7286/2.6740], 's', color='g')
plt.plot([6./77.1464789152], [1.665/2.6740], 's', color='g')
plt.plot([6./77.1464789152, 6./77.1464789152], [4.580/2.6740, 2.776/2.6740], 'o-', color='m')
plt.plot([11./77.1464789152, 6./77.1464789152], [4.580/2.6740, 2.617/2.6740], 'o-', color='m')
plt.plot([0.013, 0.21], [10., 0.88], '-')
plt.plot([0.017, 0.21], [0.1, 0.88], '-')
plt.ylim(0.1, 10.)
plt.xlim(0.01, 1.)
plt.xlabel('s')
plt.ylabel('q')
ax.set_xscale("log", nonposx='clip')
ax.set_yscale("log", nonposy='clip');
TODO: действительно ли приближение фейлит на транзите через границу?
Влияние скорости звука:
In [356]:
%%time
r25_ = hK*(25. - mudK)/1.0857
plot_2f_vs_1f(ax=plt.gca(), total_gas_data=total_gas_data_,
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:6]],
scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.27, band='K'),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.27, band='K'),
data_lim=data_lim, color='g', alpha=0.2, disk_scales=disk_scales, label='K Heidt maxdisc')
plot_param_depend(ax=plt.gca(), data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='K maxdisc', N = 20,
total_gas_data=total_gas_data_,
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=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.27, band='K'),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.27, band='K'))
plt.savefig('..\\..\pics\\cg\\'+name+'.png', format='png', bbox_inches='tight');
Влияние убирания молек. газа:
In [357]:
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=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.27, band='K'),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.27, band='K'),
data_lim=data_lim, color='g', alpha=0.2, disk_scales=disk_scales, label='K Heidt 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)])[:6],
epicycl=epicyclicFreq_real,
gas_approx=spl_gas,
sound_vel=sound_vel,
scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.27, band='K'),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.27, band='K'),
data_lim=data_lim, color='m', alpha=0.1, disk_scales=disk_scales, label='K Heidt maxdisc without H2')
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 [137]:
%%time
plot_param_depend(ax=plt.gca(), data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='K maxdisc', N = 10,
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=map(lambda c: lambda l: surf_density(mu_disc(l, mu0=mudK, h=hK), c, 'K'), np.linspace(1., 12., 10)),
star_density_min=map(lambda c: lambda l: surf_density(mu_disc(l, mu0=mudK, h=hK), c, 'K'), np.linspace(1., 12., 10)));
Замена spl_gas на gas_approx:
In [138]:
%%time
plot_param_depend(ax=plt.gca(), data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='K maxdisc', N = 2,
total_gas_data=total_gas_data_,
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=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.42, band='K'),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.42, band='K'));
Разные реалистичные дисперсии:
In [139]:
%%time
plot_param_depend(ax=plt.gca(), data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='K maxdisc', N = 2,
total_gas_data=total_gas_data_,
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=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.42, band='K'),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.42, band='K'));
Необходимо узнать, как влияет разброс у гле наклона на итоговый результат. К сожалению кроме как вручную это сложно сделать.
In [112]:
# Вроде бы кривые вращения не надо исправлять
In [113]:
# fig = plt.figure(figsize=[12, 8])
# for ind, i in enumerate([34., 38.]):
# 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 [114]:
# print 145.78/159.43, np.sin(36.*np.pi/180.)/np.sin(40.*np.pi/180.)
In [358]:
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([34., 38.]):
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 [116]:
# plot_2f_vs_1f(ax=axes[4], total_gas_data=total_gas_data_,
# epicycl=epicyclicFreq_real,
# gas_approx=spl_gas,
# sound_vel=sound_vel,
# scale=scale,
# sigma_max=sig_R_maj_max,
# sigma_min=sig_R_maj_min,
# star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.42, band='K'),
# star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.42, band='K'),
# data_lim=data_lim, color='g', alpha=0.2, disk_scales=disk_scales, label='K Heidt maxdisc')
# plot_2f_vs_1f(ax=axes[4], total_gas_data=total_gas_data_,
# epicycl=epicyclicFreq_real,
# gas_approx=spl_gas,
# sound_vel=sound_vel,
# scale=scale,
# sigma_max=sig_R_maj_max,
# sigma_min=sig_R_maj_min,
# star_density_max=lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu0d_s4g, h=h_disc_s4g), 1.35),
# lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu0d_s4g_2, h=h_disc_s4g_2), 1.35)))(l),
# star_density_min=lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu0d_s4g, h=h_disc_s4g), 1.35),
# lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu0d_s4g_2, h=h_disc_s4g_2), 1.35)))(l),
# data_lim=data_lim, color='y', alpha=0.2, disk_scales=disk_scales, label='S4G 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')
In [361]:
max_coeffs_incl = []
for ind, i in enumerate([34., 38.]):
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 [362]:
fig = plt.figure(figsize=[10, 6])
ax = plt.gca()
for ind, i in enumerate([34., 38.]):
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=total_gas_data_,
epicycl=epicyclicFreq_real,
gas_approx=spl_gas,
sound_vel=sound_vel,
scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=max_coeffs_incl[ind]['Heidt K'][0], band='K'),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=max_coeffs_incl[ind]['Heidt K'][0], band='K'),
data_lim=data_lim, color=cm.rainbow(np.linspace(0, 1, 4))[ind], alpha=0.2, disk_scales=disk_scales, label='K Heidt maxdisc inc={}'.format(incl))
plot_2f_vs_1f(ax=ax, total_gas_data=total_gas_data_,
epicycl=epicyclicFreq_real,
gas_approx=spl_gas,
sound_vel=sound_vel,
scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu0d_s4g, h=h_disc_s4g), max_coeffs_incl[ind]['S4G 2d'][0]),
lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu0d_s4g_2, h=h_disc_s4g_2), max_coeffs_incl[ind]['S4G 2d'][0])))(l),
star_density_min=lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu0d_s4g, h=h_disc_s4g), max_coeffs_incl[ind]['S4G 2d'][0]),
lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu0d_s4g_2, h=h_disc_s4g_2), max_coeffs_incl[ind]['S4G 2d'][0])))(l),
data_lim=data_lim, color=cm.rainbow(np.linspace(0, 1, 4))[ind+2], alpha=0.2, disk_scales=disk_scales, label='S4G 2d maxdisc inc={}'.format(incl))
plt.ylim(0., 2.5)
plt.xlim(0., 130.)
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 [363]:
incl = 36.
sin_i, cos_i = np.sin(incl*np.pi/180.), np.cos(incl*np.pi/180.)
In [111]:
import matplotlib as mpl
mpl.style.use('classic')
In [124]:
from matplotlib import container
fig = plt.figure(figsize=[16, 12])
ax = plt.subplot2grid((2,2), (0, 0))
# ax.imshow(ImagePIL.open('n2985_DSS_10x10min.jpg'), aspect='auto', cmap='gray')
ax.imshow(ImagePIL.open('dss2985.png'), aspect='auto', cmap='gray')
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,350)
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.errorbar(r_H2_dens, H2_dens, [0.75*3.2, 0.25*3.2, 0.2*3.2], fmt='s-', label=r'$\rm{H_2}$', ms=8, capsize=5)
# ax.plot(r_H2_dens, H2_dens, 's-', label=r'$\rm{H_2}$', ms=8)
# ax.plot(r_g_dens, [y_interp_(l) for l in r_g_dens], 's-', label=r'$\rm{H_2}$')
sum_dens = [He_coeff*(y_interp_(l[0]) + l[1]) for l in zip(r_g_dens, gas_dens)]
sum_dens[3] = He_coeff*(gas_dens[3] + H2_dens[2])
ax.errorbar(r_g_dens, sum_dens, fmt='o-', label=r'$\rm{HI+H_2}$', ms=8)
plt.semilogy(test_points, [surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.27, band='K') for l in test_points], '--', lw=3, color='g')
plt.semilogy(test_points, [tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu0d_s4g, h=h_disc_s4g), 1.23),
lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu0d_s4g_2, h=h_disc_s4g_2), 1.23)))(l) for l in test_points], '-.', lw=3, color='m')
ax.grid()
ax.set_xlim(0, 100)
ax.set_ylim(1.0, 1550.5)
ax.legend(fontsize=20)
# trick to make H2 second not third in legend
handles, labels = ax.get_legend_handles_labels()
handles2 = [h[0] if isinstance(h, container.ErrorbarContainer) else h for h in handles]
handles2[1] = handles[1]
ax.legend(handles2, labels, fontsize=20)
ax.set_ylabel(u'$\Sigma,\,M_\u2609/\mathrm{pc}^2$', fontsize=30, labelpad=-5)
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 = 'data', ms=15, color='red')
plt.errorbar(r_wsrt, vel_wsrt, yerr=e_vel_wsrt, fmt='.', marker='.', mew=1, label = r'Noordermeer et al. (2007) $\rm{HI}$', color='r', ms=15)
ax.plot(test_points, spl_gas(test_points), '-', lw=3)
# plt.plot(r_wsrt, vel_wsrt, '.', label="gas Struve")
# plt.plot(r_ma_n, vel_ma_n, 's', label='HI from fig')
# 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 ['Heidt K', 'S4G 2d']:
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'$S^4G$ Salo at al. (2015) $3.6$', lw=3, color='m')
else:
ax.plot(test_points, map(lambda l: disc_vel(l, max_coeff**2 * photom[7](0), photom[5], scale), test_points), '--', label=r'Mollenhof & Heidt (2001) $K$ max', lw=3, color='g')
ax.grid(linewidth=0.5)
ax.set_ylim(0, 270)
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+'2985_obs_data.eps', format='eps', bbox_inches='tight')
plt.savefig(paper_imgs_dir+'2985_obs_data.png', format='png', bbox_inches='tight')
plt.savefig(paper_imgs_dir+'2985_obs_data.pdf', format='pdf', dpi=150, bbox_inches='tight')
plt.show()
In [110]:
test_points = np.linspace(0.1, 400., 100)
fig, axes = plt.subplots(figsize=[18, 22], ncols=2, nrows=3)
ax = axes[0][0]
ax.plot(test_points, spl_gas(test_points), '-')
ax.set_title('v_c')
ax.set_ylim(0)
ax = axes[0][1]
ax.plot(test_points, spl_gas.derivative()(test_points), '-')
ax.set_title('dv_c/dR')
ax = axes[1][0]
ax.plot(test_points, [sqrt(2)*(spl_gas(l)/(l*scale))*sqrt(1+spl_gas.derivative()(l)*l/spl_gas(l)) for l in test_points], '-')
ax.set_title('kappa')
# ax.set_ylim(0, 600)
ax = axes[1][1]
ax.plot(test_points, [spl_gas(l)/(l*scale) for l in test_points], '-')
ax.set_title('Omega')
ax = axes[2][0]
ax.plot(test_points, [0.5*(1+spl_gas.derivative()(l)*l/spl_gas(l)) for l in test_points], '-')
ax.set_title('sigPhi2/sigR2')
ax = axes[2][1]
ax.plot(test_points, [sqrt(0.5*(1+spl_gas.derivative()(l)*l/spl_gas(l))) for l in test_points], '-')
ax.plot(test_points, [scale*l*epicyclicFreq_real(spl_gas, l, scale)/(2*spl_gas(l)) for l in test_points], '-')
ax.set_title('kappa/2Omega, sigPhi/sigR')
plt.show()
А вот почему с дипломом не сходится - надо было взять фотометрию в $R$ с M/L=6 и сделать дисперсии похожими на диплом:
In [140]:
def sig_R_tmp(r):
if r < sig_maj_lim:
return spl_maj(r) + 15.
else:
return spl_maj(sig_maj_lim) + 15.
plt.errorbar(r_sig_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='.', mew=0, color='red', label='$\sigma_{los}^{maj}$')
# plt.plot(points, spl_maj(points), label = '$\sigma_{los}^{maj}\, splinefit$', color='red')
plt.plot(points, map(sig_R_tmp, points), label = '$\sigma_R^{max}$', color='blue')
plt.legend()
plt.ylim(0,95)
plt.xlim(0,600);
In [141]:
gas_data = zip(r_g_dens, map(lambda l: l*He_coeff, gas_dens))[:5]
invQg, invQs, invQeff = zip(*get_invQeff_from_data(gas_data=gas_data,
epicycl=epicyclicFreq_real,
gas_approx=spl_gas,
sound_vel=sound_vel,
scale=scale,
sigma=sig_R_maj_maxmax,
star_density=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L=6., band='R')))
invQg, invQs, invQeff_ma = zip(*get_invQeff_from_data(gas_data=gas_data,
epicycl=epicyclicFreq_real,
gas_approx=spl_gas,
sound_vel=sound_vel,
scale=scale,
sigma=sig_R_maj_minmin,
star_density=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), 6., 'R')))
invQg = map(lambda l: l*1.6, invQg)
invQeff = map(lambda l: l*1.6, invQeff)
invQeff_ma = map(lambda l: l*1.6, invQeff_ma)
fig = plt.figure(figsize=[10, 6])
plt.fill_between(r_g_dens[:5], invQeff, invQeff_ma, color='g', alpha=0.3)
plt.plot(r_g_dens[:5], invQeff, 'v-', label='Qeff', color='g', alpha=0.6)
plt.plot(r_g_dens[:5], invQeff_ma, 'v-', color='g', alpha=0.6)
plt.plot(r_g_dens[:5], invQg, 'v-', label='Qg', color='b')
r_dip, Q_dip = zip(*np.loadtxt("diplom_results_data.dat", float, delimiter=','))
plt.plot(r_dip[:3], Q_dip[:3], 's-', color='g', alpha=0.9)
plt.plot(r_dip[3:], Q_dip[3:], 's-', color='b', alpha=0.9)
invQg, invQs, invQeff_ = zip(*get_invQeff_from_data(gas_data=gas_data,
epicycl=epicyclicFreq_real,
gas_approx=spl_gas,
sound_vel=sound_vel,
scale=scale,
sigma=sig_R_tmp,
star_density=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), 6., 'R')))
invQeff_ = map(lambda l: l*1.6, invQeff_)
plt.plot(r_g_dens[:5], invQeff_, 'o-', color='g', alpha=0.6)
plt.ylim(0., 1.)
plt.xlim(0., 150.)
plot_SF(plt.gca())
plot_data_lim(plt.gca(), data_lim)
plot_disc_scale(hK, plt.gca())
plt.legend()
plt.grid();
Plotly:
In [143]:
from plotly import tools
from plotly.offline import init_notebook_mode, iplot
import plotly.plotly as py
import plotly.graph_objs as go
init_notebook_mode(connected=True)
acc_dot = go.Scatter(
x = r_g_dens[:5],
y = invQeff,
mode='lines+markers',
name = 'Qeff',
marker=dict(
size='16',
color = ['#770000']*len(invQeff))
)
qgdat = go.Scatter(
x = r_g_dens[:5],
y = invQg,
mode='lines+markers',
name = 'Qg',
marker=dict(
size='16',
color = ['#007700']*len(invQeff))
)
data = [acc_dot, qgdat]
iplot(data, filename='scatter-plot-with-colorscale')
Проверки про дисперсии разные:
In [150]:
# Исправляем значения вдоль малой оси на синус угла:
def correct_min(R):
return R / cos_i
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)
r_sig_mi, sig_mi, e_sig_mi = zip(*np.loadtxt("s_stars_miN.dat", float))
# r_sig_ma, sig_ma, e_sig_ma = zip(*np.loadtxt("s_stars_maN.dat", float))
r_, s_ = zip(*np.loadtxt("noordermeer_data/n08_disp.dat", float, delimiter=','))
rr_, sr_ = zip(*np.loadtxt("noordermeer_data/n08_disp_ratio.dat", float, delimiter=','))
fig = plt.figure(figsize=[11, 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(r_sig_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='.', mew=0, color='blue', label=r'$\sigma_{los}^{maj} Noord$')
r_mi_extend = map(correct_min, r_sig_mi)
plt.errorbar(r_mi_extend, sig_mi, yerr=e_sig_mi, fmt='.', marker='.', mew=0, color='black', label='$\sigma_{los}^{min} Noord$')
spl_mi1 = inter.UnivariateSpline(r_mi_extend, sig_mi, k=3, s=100.)
plt.plot(points, spl_mi1(points), '--', alpha=0.3)
plt.errorbar(r_sig_mi, sig_mi, yerr=e_sig_mi, fmt='.', marker='.', mew=0, color='m', label='$\sigma_{los}^{min} Noord$')
sc = 0.102
plt.plot(map(lambda l: l/sc, r_), s_, '.', color='r')
plt.plot(map(lambda l: l/sc, rr_), map(lambda l: l[1]*spl_maj(l[0]/sc), zip(rr_, sr_)), '.', color='g')
# plt.plot(points, spl_maj(points), '--')
plt.xticks([l for l in np.linspace(0, 110., 20.)], ['%2.1f'%(l*sc) for l in np.linspace(0, 110., 20.)])
plt.xlabel('$R$')
plt.ylabel('$\sigma$')
plt.grid()
plt.ylim(0)
plt.legend();
In [151]:
Image('noordermeer_data/n08_disp.png', width=300)
Out[151]:
Проверим применимость WKB приближения, т.е. $k\times r \gg 1$:
Для $K$:
In [130]:
plot_WKB_dependencies(r_g_dens=zip(*total_gas_data_)[0],
gas_dens=zip(*total_gas_data_)[1],
epicycl=lambda l: epicyclicFreq_real(spl_gas, l, scale),
sound_vel=sound_vel,
star_density=[surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.42, band='K') for l in zip(*total_gas_data_)[0]],
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.axvline(x=1, alpha=0.1)
plt.xlim(0, 100);
Исходная зависимость:
In [131]:
plot_k_dependencies(r_g_dens=zip(*total_gas_data_)[0],
gas_dens=zip(*total_gas_data_)[1],
epicycl=lambda l: epicyclicFreq_real(spl_gas, l, scale),
sound_vel=sound_vel,
star_density=[surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.42, band='K') for l in zip(*total_gas_data_)[0]],
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);
Для $S^4G$:
In [132]:
plot_WKB_dependencies(r_g_dens=zip(*total_gas_data_)[0],
gas_dens=zip(*total_gas_data_)[1],
epicycl=lambda l: epicyclicFreq_real(spl_gas, l, scale),
sound_vel=sound_vel,
star_density=[tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu0d_s4g, h=h_disc_s4g), 1.35),
lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu0d_s4g_2, h=h_disc_s4g_2), 1.35)))(l) for l in zip(*total_gas_data_)[0]],
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.axvline(x=1, alpha=0.1)
plt.xlim(0, 100);
Исходная зависимость:
In [133]:
plot_k_dependencies(r_g_dens=zip(*total_gas_data_)[0],
gas_dens=zip(*total_gas_data_)[1],
epicycl=lambda l: epicyclicFreq_real(spl_gas, l, scale),
sound_vel=sound_vel,
star_density=[tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu0d_s4g, h=h_disc_s4g), 1.35),
lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu0d_s4g_2, h=h_disc_s4g_2), 1.35)))(l) for l in zip(*total_gas_data_)[0]],
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 [152]:
Image('HST_color_dist_estimation.jpg')
Out[152]:
Для угловых координат (по идее достаточно далеко, чтобы просто евклидово расстояние использовать):
In [153]:
from scipy.spatial.distance import euclidean
point_a = (16.109, 74.77)
point_b = (16.217, 43.07)
point_c = (26.580, 57.24)
print euclidean(point_a, point_b), euclidean(point_a, point_b)/319.
print euclidean(point_b, point_c), euclidean(point_b, point_c)/491.
print euclidean(point_a, point_c), euclidean(point_a, point_c)/507.
Странно, первое расстояние получается больше (а должно быть меньше) и масштаб не тот что нужен.
При этом для координат попиксельных все шикарно сходится:
In [154]:
point_a = (2731., 2874.)
point_b = (2721., 2240.)
point_c = (1775., 2524.)
print euclidean(point_a, point_b), euclidean(point_a, point_b)/319.
print euclidean(point_b, point_c), euclidean(point_b, point_c)/491.
print euclidean(point_a, point_c), euclidean(point_a, point_c)/507.
Ну ок, я не буду разбираться, возьму последние два масштаба, усредню их и рассчитаю полминуты:
In [155]:
30./((0.0402745085126 + 0.0357537175275)/2.) #px
Out[155]:
Похоже, если сравнивать визуально с DSS:
In [156]:
Image('n2985_DSS_dist_estimation.jpg')
Out[156]:
Тогда 15 секунд:
In [157]:
15./((0.0402745085126 + 0.0357537175275)/2.) #px
Out[157]:
их и нанесем:
In [158]:
Image('HST_color_dist.jpg', width=300)
Out[158]:
In [113]:
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,350)
plt.xlim(0,100);
In [114]:
fig, axes = plt.subplots(1, 5, figsize=[40,7])
fig.tight_layout()
axes[0].imshow(ImagePIL.open('HST_color_dist.jpg'), aspect='auto', cmap='Greys')
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), next(linecycler), 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]) + l[1]) for l in zip(r_g_dens, gas_dens)], '*-')
axes[3].plot(r_g_dens, [y_interp_(l) for l in r_g_dens], '--')
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=total_gas_data_,
epicycl=epicyclicFreq_real,
gas_approx=spl_gas,
sound_vel=sound_vel,
scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.42, band='K'),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.42, band='K'),
data_lim=data_lim, color='g', alpha=0.2, disk_scales=disk_scales, label='K Heidt maxdisc')
plot_2f_vs_1f(ax=axes[4], total_gas_data=total_gas_data_,
epicycl=epicyclicFreq_real,
gas_approx=spl_gas,
sound_vel=sound_vel,
scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu0d_s4g, h=h_disc_s4g), 1.35),
lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu0d_s4g_2, h=h_disc_s4g_2), 1.35)))(l),
star_density_min=lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu0d_s4g, h=h_disc_s4g), 1.35),
lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu0d_s4g_2, h=h_disc_s4g_2), 1.35)))(l),
data_lim=data_lim, color='y', alpha=0.2, disk_scales=disk_scales, label='S4G 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[114]:
In [131]:
fig = plt.figure(figsize=[16, 8])
plt.plot(r_g_dens, gas_dens, 'd-', label=r'$\rm{HI}$', ms=10)
plt.plot(r_g_dens, [y_interp_(l) for l in r_g_dens], 's-', label=r'$\rm{H_2}$', ms=10)
plt.plot(r_g_dens, [He_coeff*(y_interp_(l[0]) + l[1]) for l in zip(r_g_dens, gas_dens)], 'o-', label=r'$\rm{HI+H_2}$', ms=10)
plt.grid()
plt.xlim(0, 200)
plt.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 [138]:
# from matplotlib.ticker import MultipleLocator, FormatStrFormatter
# majorLocator = MultipleLocator(20)
# majorFormatter = FormatStrFormatter('%d')
# minorLocator = MultipleLocator(5)
# t = np.arange(0.0, 100.0, 0.1)
# s = np.sin(0.1*np.pi*t)*np.exp(-t*0.01)
# fig, ax = plt.subplots()
# plt.plot(t, s)
# ax.xaxis.set_major_locator(majorLocator)
# ax.xaxis.set_major_formatter(majorFormatter)
# # for the minor ticks, use no labels; default NullFormatter
# ax.xaxis.set_minor_locator(minorLocator)
fig = plt.figure(figsize=[16, 8])
disk_scales2 = []
disk_scales2.append((31.08, 'K'))
disk_scales2.append((12.78, 'S4G'))
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=total_gas_data_,
epicycl=epicyclicFreq_real,
gas_approx=spl_gas,
sound_vel=sound_vel,
scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.42, band='K'),
star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.42, band='K'),
data_lim=data_lim, color='g', alpha=0.2, disk_scales=disk_scales2, label='K Heidt maxdisc')
plot_2f_vs_1f(ax=ax, total_gas_data=total_gas_data_,
epicycl=epicyclicFreq_real,
gas_approx=spl_gas,
sound_vel=sound_vel,
scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu0d_s4g, h=h_disc_s4g), 1.35),
lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu0d_s4g_2, h=h_disc_s4g_2), 1.35)))(l),
star_density_min=lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu0d_s4g, h=h_disc_s4g), 1.35),
lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu0d_s4g_2, h=h_disc_s4g_2), 1.35)))(l),
data_lim=data_lim, color='y', alpha=0.2, disk_scales=disk_scales2, label='S4G 2d maxdisc')
ax.set_ylim(0., 1.0)
ax.set_xlim(0., 130.)
ax.axhline(y=1., ls='-', color='grey')
ax.plot([10., 7.2/(scale*22.4/19.)], [0., 0.], '-', lw=7., color='r')
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 [159]:
# %install_ext http://raw.github.com/jrjohansson/version_information/master/version_information.py
%load_ext version_information
%reload_ext version_information
%version_information numpy, scipy, matplotlib
Out[159]: