Галактика найдена из пересечения HYPERLEDA и https://arxiv.org/pdf/1608.06735v1.pdf.
In [1]:
from IPython.display import HTML
from IPython.display import Image
import os
%pylab
%matplotlib inline
%run ../../utils/load_notebook.py
In [2]:
from photometry import *
In [3]:
from instabilities import *
In [4]:
from utils import *
In [5]:
name = 'N4725'
gtype = 'SABa' #LEDA, 'SBbc' from Heraudeau98
incl = 50. #mean value
distance = 18.2 # Noordermeer
scale = 0.088 #kpc/arcsec according to Noordermeer
data_path = '../../data/n4725_u7989'
sin_i, cos_i = np.sin(incl*np.pi/180.), np.cos(incl*np.pi/180.)
In [6]:
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')
In [7]:
os.chdir(data_path)
# Данные из NED
HTML('<iframe src=http://ned.ipac.caltech.edu/cgi-bin/objsearch?objname=ngc+4725&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[7]:
In [8]:
# Данные из HYPERLEDA
HTML('<iframe src=http://leda.univ-lyon1.fr/ledacat.cgi?o=ngc4725 width=1000 height=350></iframe>')
Out[8]:
In [9]:
#SDSS
Image('n4725_SDSS.jpg', width=500)
Out[9]:
Проблема в том, что эта галактика есть в DR7, но в более поздних релихах я не смог ее найти, соответственно не знаю масштаба корректного.
Однако я понял, что сама картинка из выборки http://cosmo.nyu.edu/hogg/rc3/ и там есть с маштабом изображение:
In [10]:
Image('n4725_SDSS_labeled.jpg', width=500)
Out[10]:
In [11]:
#JHK
Image('n4725_2MASS.jpg', width=300)
Out[11]:
In [12]:
Image('noord_p113_cite.png')
Out[12]:
Дисперсии скоростей и кривая вращения - есть в Heraudeau 1999 http://adsabs.harvard.edu/cgi-bin/bib_query?1999A%26AS..136..509H до ~50'' (1 разрез), PA=$35^{\circ}$
In [13]:
# Данные по звездной кинематике Heraudeau+1999 вдоль большой полуоси (не исправленные за наклон?) - из HYPERLEDA
r_ma, vel_ma, e_vel_ma, sig_ma, e_sig_ma = zip(*np.loadtxt("her99_kinem.dat", float))
fig = plt.figure(figsize=[8,5])
plt.errorbar(r_ma, vel_ma, e_vel_ma, fmt='.', marker='.', mew=0, label="Heraudeau 1998 stars maj")
plt.legend()
plt.ylim(-350., 350.)
plt.show()
In [14]:
Image('her99_rot.png') #оригинал
Out[14]:
Приближение:
In [15]:
r_ma_b, vel_ma_b, e_vel_b = zip(*sorted(zip(np.abs(r_ma), np.abs(vel_ma), e_vel_ma)))
In [16]:
fig = plt.figure(figsize=[8,4])
plt.errorbar(r_ma_b, vel_ma_b, yerr=e_vel_b, fmt='.', marker='.', mew=0, color='blue', label = 'Her98 star maj')
test_points = np.linspace(0.0, max(r_ma_b), 100)
poly_star = poly1d(polyfit(r_ma_b, vel_ma_b, deg=7))
plt.plot(test_points, poly_star(test_points), '-', label='poly deg=7')
def w(arr):
return map(lambda l: 1/(1. + l**2), arr)
import scipy.interpolate as inter
spl = inter.UnivariateSpline(r_ma_b, vel_ma_b, k=3, s=10000., w=w(e_vel_b))
plt.plot(test_points, spl(test_points), '-', label='spl s=10000 w^2')
spl = inter.UnivariateSpline(r_ma_b, vel_ma_b, k=3, s=10000.)
plt.plot(test_points, spl(test_points), '-', label='spl s=10000')
plt.legend(loc='upper left')
plt.ylim(0, 170)
plt.show()
C весами плохо получается, полином и обычный почти совпадают - берем их:
In [17]:
star_approx = spl
In [18]:
r_sig_ma = r_ma #Heraudeau+1999
fig = plt.figure(figsize=[6., 4.])
plt.errorbar(r_sig_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='.', mew=0, color='blue', label=r'$maj\, Heraudeau $')
plt.xlabel('$R$')
plt.ylabel('$\sigma$')
plt.ylim(0, 350)
plt.legend();
In [19]:
Image('her99_disp.png') #из статьи
Out[19]:
In [20]:
fig = plt.figure(figsize=[6., 4.])
plt.errorbar(map(abs, r_sig_ma), sig_ma, yerr=e_sig_ma, fmt='.', marker='.', mew=0, color='blue', label=r'$maj\, Heraudeau $')
plt.xlabel('$R$')
plt.ylabel('$\sigma$')
plt.ylim(0, 250)
plt.legend();
Для большой оси: $\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 [21]:
spl_maj = inter.UnivariateSpline(r_sig_ma, sig_ma, k=3, s=10000.)
sig_maj_lim = max(r_sig_ma)
points = np.linspace(0.1, max(r_ma)+15., 100)
In [22]:
# TODO: move to external file
def flat_end(argument):
'''декоратор для того, чтобы продолжать функцию на уровне последнего значения'''
def real_decorator(function):
def wrapper(*args, **kwargs):
if args[0] < argument:
return function(*args, **kwargs)
else:
return function(argument, *args[1:], **kwargs)
return wrapper
return real_decorator
@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))
Используем соотношение $\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 [23]:
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 [24]:
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,80);
Видно, что настоящее еще и больше.
In [25]:
Image('HI_rot.png') #HI WHISP from van Eymeren 2011, PA=32.45, R25=25.91kpc, разные цвета - это две стороны, черная - усредненная
Out[25]:
Достаточно странно расположились черные точки - не между двумя другими, выглядит неверным
In [26]:
R25 = 25.91
# Данные по кинематике газа van Eymeren 2011 в HI
r_hi, vel_hi = zip(*np.loadtxt("HI_rot.dat", float, delimiter=','))
plt.plot(r_hi, vel_hi, 'd', label='HI Eymeren+2011')
plt.ylim(0, 300)
plt.xlim(-0.5, 2.)
Out[26]:
Из Ноордермеера:
In [27]:
Image('noord_rot.png')
Out[27]:
In [28]:
r_n, vel_n = zip(*np.loadtxt("noord_rot.dat", float, delimiter=','))
plt.plot(r_n, vel_n, 'd', label='HI Noord+2005')
plt.ylim(850, 1550)
plt.xlim(-10., 11.)
Out[28]:
In [29]:
vel_n = map(lambda l: l-1208., vel_n)
r_n, vel_n = zip(*sorted(zip(np.abs(r_n), np.abs(vel_n))))
r_n = [l*60 for l in r_n]
Посмотрим на согласие между измерениями:
In [30]:
plt.plot(r_n, vel_n, 'd', label='HI Noord+2005')
plt.plot([l*R25/calc_scale(18.2) for l in r_hi], vel_hi, 'd', label='HI Eymeren+2011')
plt.legend();
Похоже, что у Ноордермеера не нормировано на угол. Проверим:
In [31]:
plt.plot(r_n, [l/sin_i for l in vel_n], 'd', label='HI Noord+2005')
plt.plot([l*R25/scale for l in r_hi], vel_hi, 'd', label='HI Eymeren+2011')
plt.legend(loc='lower right');
Да, так и есть - данные хорошо совпали для исправленных за угол измерений Ноордермеера (даже для такой низкой точности снятия данных с рисунка). Странно, что у более новых данных протяженность значительно меньше. Для построения кривой возьмем более оба набора до $200^{''}$:
In [32]:
r_hi = [l*R25/calc_scale(18.2) for l in r_hi]
vel_n = [l/sin_i for l in vel_n]
In [33]:
fig = plt.figure(figsize=[10,6])
_1,_2, = [0.0,],[0.0,]
_1.extend(r_hi)
_2.extend(vel_hi)
_1,_2 = zip(*sorted(zip(_1,_2)))
gas_approx = poly1d(polyfit(_1, _2, deg=9))
test_points = np.linspace(0, max(r_hi), 100)
plt.plot(test_points, gas_approx(test_points), '--', label='poly approx')
spl_gas = inter.UnivariateSpline(_1, _2, k=3, s=2000.)
plt.plot(test_points, spl_gas(test_points), '-', label='spline')
plt.plot(r_hi, vel_hi, 'd', label='HI Eymeren+2011')
gas_approx = poly1d(polyfit(r_n[:-20], vel_n[:-20], deg=9))
plt.plot(test_points, gas_approx(test_points), '--', label='poly approx N')
spl_gas = inter.UnivariateSpline(r_n[:-20], vel_n[:-20], k=3, s=10000.)
plt.plot(test_points, spl_gas(test_points), '-', label='spline N')
plt.plot(r_n[:-20], vel_n[:-20], 'v', label='HI Noord+2005')
plt.ylim(0, 300)
plt.legend(loc='lower right');
Похоже для данных Eymeren вторая точка все портит и вообще точность в нужной нам области хромает.
In [34]:
fig = plt.figure(figsize=[10,6])
_1,_2, = [0.0,],[0.0,]
_1.extend(r_hi[1:])
_2.extend(vel_hi[1:])
_1,_2 = zip(*sorted(zip(_1,_2)))
gas_approx = poly1d(polyfit(_1, _2, deg=9))
test_points = np.linspace(0, max(r_hi), 100)
plt.plot(test_points, gas_approx(test_points), '--', label='poly approx')
spl_gas = inter.UnivariateSpline(_1, _2, k=3, s=2000.)
plt.plot(test_points, spl_gas(test_points), '-', label='spline')
plt.plot(_1, _2, 'd', label='HI Eymeren+2011')
gas_approx = poly1d(polyfit(r_n[:-20], vel_n[:-20], deg=9))
plt.plot(test_points, gas_approx(test_points), '--', label='poly approx N')
spl_gas_N = inter.UnivariateSpline(r_n[:-20], vel_n[:-20], k=3, s=10000.)
plt.plot(test_points, spl_gas_N(test_points), '-', label='spline N')
plt.plot(r_n[:-20], vel_n[:-20], 'v', label='HI Noord+2005')
plt.ylim(0, 300)
plt.legend(loc='lower right');
In [35]:
fig = plt.figure(figsize=[10,6])
_1,_2, = [0.0,],[0.0,]
_1.extend(r_hi[1:])
_1.extend(r_n[1::2])
_2.extend(vel_hi[1:])
_2.extend(vel_n[1::2])
_1,_2 = zip(*sorted(zip(_1,_2)))
gas_approx = poly1d(polyfit(_1, _2, deg=9))
test_points = np.linspace(0, max(r_hi), 100)
plt.plot(test_points, gas_approx(test_points), '--', label='poly approx')
spl_gas = inter.UnivariateSpline(_1, _2, k=3, s=5000.)
plt.plot(test_points, spl_gas(test_points), '-', label='spline')
plt.plot(_1, _2, 'd', label='HI Eymeren+2011')
plt.plot(r_n, vel_n, 'v', label='HI Noord+2005')
plt.ylim(0, 300)
plt.xlim(0, 300)
plt.legend(loc='lower right');
Если убрать эту точку - то сплайны почти совпадают, но все равно точность недостаточная.
TODO: поискать более точные данные - например CO
Для случая бесконечного тонкого диска: $$\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 [36]:
test_points = np.linspace(0, 350, 1000)
fig = plt.figure(figsize=[8, 4])
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.plot(test_points, [epicyclicFreq_real(spl_gas_N, x, scale) for x in test_points], '-', label='spline N')
plt.xlabel('$R, arcsec$')
plt.ylabel('$\kappa,\, km/s/kpc$', fontsize=15)
plt.ylim(0, 300)
plt.legend();
Достаточно сложно, учитывая что нас интересуют первые 50 секунд. Будем считать, что нас интересуют сплайны.
(необходимо иметь в виду, что в первой работе общий профиль исправлен за гелий следующим образом $\Sigma_g = 1.36\times(\Sigma_{H_2} + \Sigma_{HI})$ и переход CO-to-H2 тоже конкретный)
In [37]:
# у них расстояние в работе 26.8
gas_scale = calc_scale(26.8)
print gas_scale
In [38]:
r_g_dens1, gas_dens1 = zip(*np.loadtxt("gas_dens_err.txt", float, delimiter=';'))
hi_yerr = np.abs((np.column_stack([gas_dens1[58::3], gas_dens1[59::3]])-np.column_stack([gas_dens1[57::3], gas_dens1[57::3]]))).T
h2_yerr = np.abs((np.column_stack([gas_dens1[1:57:3], gas_dens1[2:57:3]])-np.column_stack([gas_dens1[0:57:3], gas_dens1[0:57:3]]))).T
In [39]:
fig = plt.figure(figsize=[10, 7.5])
plt.errorbar([l/gas_scale for l in r_g_dens1[57::3]], gas_dens1[57::3], hi_yerr, fmt='bo', capsize=5)
plt.errorbar([l/gas_scale for l in r_g_dens1[0:57:3]], gas_dens1[0:57:3], h2_yerr, fmt='ro', capsize=5)
ax = plt.gca()
ax.set_yscale("log", nonposy='clip')
plt.ylim(0.01, 1000)
plt.xlim(0);
In [40]:
Image('u7989_gas_dens.png')
Out[40]:
Хорошо совпадает. Дальше будем работать с теми же что и раньше, а эти прибережем для финальной картинки. Единственное что все (и ошибки) надо исправить за другое расстояние и $X_{CO}$
In [41]:
gas_dens1 = [l*(26.8/distance)**2 for l in gas_dens1]
hi_yerr = hi_yerr*(26.8/distance)**2
h2_yerr = h2_yerr*(26.8/distance)**2
r_g_dens1 = [l/gas_scale for l in r_g_dens1]
In [42]:
r_g_dens, gas_dens = zip(*np.loadtxt("gas_dens.dat", float, delimiter=','))
plt.semilogy([l/gas_scale for l in r_g_dens[:58]], gas_dens[:58], 'd', color='blue')
plt.semilogy([l/gas_scale for l in r_g_dens[58:77]], gas_dens[58:77], 'd', color='red')
plt.semilogy([l/gas_scale for l in r_g_dens[77:]], gas_dens[77:], 'd', color='black')
plt.ylim(0.01, 1000);
In [43]:
plt.semilogy([l/gas_scale for l in r_g_dens[77:]], gas_dens[77:], 'o');
In [44]:
plt.plot([l/gas_scale for l in r_g_dens[58:77]], gas_dens[58:77], '-o', color='r')
plt.plot([l/gas_scale for l in r_g_dens[77:]], gas_dens[77:], '-o') #только сумма
plt.xlim(0, 200);
Ноордермеер:
In [45]:
Image('noord_gdens.png')
Out[45]:
In [46]:
r_g_n, gas_dens_n = zip(*np.loadtxt("HI_dens.dat", float, delimiter=','))
In [47]:
plt.plot(r_g_n, gas_dens_n, 'o')
plt.xlim(0, 450)
plt.ylim(0, 5)
Out[47]:
Сравнение:
In [48]:
plt.semilogy([l/gas_scale for l in r_g_dens[:58]], gas_dens[:58], 'o', color='blue', label='Hulst 2016')
plt.semilogy(r_g_n, gas_dens_n, 's', label='Noord HI', color='red')
plt.legend()
Out[48]:
Как и было изначально видно - пики заметно смещены (теперь нет, т.к. я это учел)
Разгадка кроется, кмк, в том что в работе 2016 года расстояние взято равным 26.8 Мпк, а в NED это соответствует (Virgo + GA + Shapley) и масштаб 0.130 kpc/arcseс.
Также это подтерждается сравнением с верхней шкалой на рисунке, где $R/R_{25}$, $R_{25}=321.^{"}$ и примерно по соотношенияю можно оценить
$$20.5\,kpc/(0.5*321)\,arcsec = 0.1308\, kpc/arcsec$$
In [49]:
plt.semilogy([l/0.130 for l in r_g_dens[:58]], gas_dens[:58], 'o', color='blue', label='Hulst 2016')
plt.semilogy(r_g_n, gas_dens_n, 's', label='Noord HI', color='red')
plt.legend();
И теперь с учетом спирали они стали похожи друг на друга.
Используем более современный газ $\rm{HI}$ + $\rm{HII}$:
In [50]:
r_HI_dens, HI_dens = [l/gas_scale for l in r_g_dens[:58]], gas_dens[:58]
r_mol_dens, mol_dens = [l/gas_scale for l in r_g_dens[58:77]], gas_dens[58:77]
r_g_dens, gas_dens = [l/gas_scale for l in r_g_dens[77:]], gas_dens[77:] #используем только полный газ
Последнее - проверим, что действительно исправлялось как сказанно в работе - домножением на 1.36 суммы газов (а заодно сравним, что если мы будем домножать только атомарный водород):
In [51]:
from scipy.interpolate import interp1d
tmp1_ = interp1d(r_HI_dens, HI_dens)
for r_, d_ in zip(r_mol_dens, mol_dens):
plt.scatter(r_, 1.36*(d_ + tmp1_(r_)), color='b')
plt.scatter(r_, 1.36*d_ + tmp1_(r_), color='m')
plt.plot(r_g_dens, gas_dens, '-')
plt.ylim(0, 20)
plt.xlim(0, 200);
Да, действительно домножается сумма и в принципе расхождение не такое большое.
Надо исправить за другое расстояние (26.8) и $X_{CO}=1.9$. Сейчас исправим только за расстояние, т.к. фактор у нас пока такой же:
In [52]:
HI_dens = [l*(26.8/distance)**2 for l in HI_dens]
mol_dens = [l*(26.8/distance)**2 * (X_CO/1.9) for l in mol_dens]
gas_dens = [l*(26.8/distance)**2 for l in gas_dens]
In [53]:
import scipy.interpolate
y_interp = scipy.interpolate.interp1d(list(r_mol_dens), list(mol_dens))
def y_interp_(r):
if r <= min(r_mol_dens):
return y_interp(min(r_mol_dens))
elif r < max(r_mol_dens):
return y_interp(r)
else:
return 0.
In [54]:
plt.plot(r_HI_dens, HI_dens)
plt.plot(r_mol_dens, mol_dens)
plt.plot(r_g_dens, gas_dens);
In [55]:
plt.plot(map(lambda l: l*scale, r_HI_dens), HI_dens, '.-')
plt.plot(map(lambda l: l*scale, r_g_n), gas_dens_n, 'o-');
И проинтегрировать массы:
In [56]:
import scipy.integrate as integrate
import scipy.interpolate
tmp_ = scipy.interpolate.interp1d(r_HI_dens, HI_dens)
result = integrate.quad(lambda l: 2*np.pi*l*tmp_(l), r_HI_dens[0], r_HI_dens[-1])
print (scale * 1000.)**2 * result[0]/1e9
In [57]:
import scipy.integrate as integrate
import scipy.interpolate
tmp_1 = filter(lambda l: l[0] < 360., zip(r_HI_dens, HI_dens))
tmp_ = scipy.interpolate.interp1d(zip(*tmp_1)[0], zip(*tmp_1)[1])
result = integrate.quad(lambda l: 2*np.pi*l*tmp_(l), zip(*tmp_1)[0][0], zip(*tmp_1)[0][-1])
print (scale * 1000.)**2 * result[0]/1e9
In [58]:
tmp_ = scipy.interpolate.interp1d(r_mol_dens, mol_dens)
result = integrate.quad(lambda l: 2*np.pi*l*tmp_(l), r_mol_dens[0], r_mol_dens[-1])
print (scale * 1000.)**2 * result[0]/1e9
In [59]:
tmp_ = scipy.interpolate.interp1d(r_g_n, gas_dens_n)
result = integrate.quad(lambda l: 2*np.pi*l*tmp_(l), r_g_n[0], r_g_n[-1])
print (scale * 1000.)**2 * result[0]/1e9
In [60]:
all_photometry = []
S4G данные из GALFIT (есть бар в модели!):
In [61]:
r_eff_s4g = 10.14
# mu_eff_s4g = ...
n_s4g = 2.212
mu0d_s4g = 20.337
h_disc_s4g = 73.20
Тут нужно учитывать, что эти параметры в AB-mag и нуждаются в доп. исправлении.
In [62]:
M_to_L_s4g = s4g_mass_to_light(-21.762, -21.267)
M_to_L_s4g
Out[62]:
In [63]:
p_ = np.arange(0.1, 200., 0.1)
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.legend();
Достаточно маленькие.
In [64]:
all_photometry.append(('S4G 3.6', r_eff_s4g, None, n_s4g, mu0d_s4g, h_disc_s4g, M_to_L_s4g,
lambda l: s4g_surf_density(mu_disc(l, mu0=mu0d_s4g, h=h_disc_s4g), M_to_L_s4g)))
In [65]:
mu0d_J = 17.78
h_disc_J = 49.99
mu0d_H = 17.11
h_disc_H = 50.28
mu0d_K = 17.01
h_disc_K = 54.66
In [66]:
b_v_color = 0.012 #TODO: не знаем на самом деле какой цвет (вот тут https://arxiv.org/pdf/1102.1724v1.pdf указано 0.012)
# тут https://arxiv.org/pdf/astro-ph/0610688v2.pdf есть B и V в Янских
M_to_L_J = bell_mass_to_light(b_v_color, 'J', 'B-V')
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))
M_to_L_H = bell_mass_to_light(b_v_color, 'H', 'B-V')
surf_H = [surf_density(mu=mu_disc(l, mu0=mu0d_H, h=h_disc_H), 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_v_color, 'K', 'B-V')
surf_K = [surf_density(mu=mu_disc(l, mu0=mu0d_K, h=h_disc_K), 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()
Out[66]:
In [67]:
p_ = np.arange(0.1, 200., 0.1)
b_r_color = 0.55 #вот отсюда древнее B-R https://ui.adsabs.harvard.edu/#abs/1995AJ....109..543B/abstract
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))
M_to_L_H = bell_mass_to_light(b_r_color, 'H', 'B-R')
surf_H = [surf_density(mu=mu_disc(l, mu0=mu0d_H, h=h_disc_H), 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=mu0d_K, h=h_disc_K), 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()
Out[67]:
Разница, как мы видим, не столь существенная.
TODO: разобраться с цветом
In [68]:
all_photometry.append(('Heidt J', None, None, None, 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')))
all_photometry.append(('Heidt H', None, None, None, mu0d_H, h_disc_H, M_to_L_H,
lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_H, h=h_disc_H), M_to_L=M_to_L_H, band='H')))
all_photometry.append(('Heidt K', None, None, None, mu0d_K, h_disc_K, M_to_L_K,
lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_K, h=h_disc_K), M_to_L=M_to_L_K, band='K')))
In [69]:
# for 3.6 band
dist_36 = 13.24 #Mpc
r_eff_36 = np.power(10., 3.)/scale/1000./(dist_36/20.4) #pc
mu_eff_36 = 17.49
n_36 = 3.61
mu0d_36 = 19.65
h_disc_36 = np.power(10., 3.66)/scale/1000./(dist_36/20.4) #because original is Log(), pc
Опять же, надо исправлять.
In [70]:
surf_36 = [s4g_surf_density(mu=m, M_to_L=M_to_L_s4g) for m in [mu_disc(l, mu0=mu0d_36, h=h_disc_36) for l in p_]]
plt.plot(p_, surf_36, '-', label='3.6 [M/L={:2.2f}]'.format(0.6))
plt.legend()
Out[70]:
Вполне неплохо. Можем сравнить еще с калибровками McGaugh:
In [71]:
bv_color = 0.9
mcgaugh_mass_to_light(bv_color, 'mu36')
Out[71]:
Ну да, тут все похоже честно. Смущает конечно разница в центральной яркости с S4G, причем я не смог ничего найти насчет депроецировано или нет. Если депроецировать - становятся похожи:
In [72]:
mu_face_on(mu0d_36, cos_i)
Out[72]:
Добавим и депроецированные и нет:
In [73]:
all_photometry.append(('infra 3.6', r_eff_36, mu_eff_36, n_36, mu0d_36, h_disc_36, M_to_L_s4g,
lambda l: s4g_surf_density(mu=mu_disc(l, mu0=mu0d_36, h=h_disc_36), M_to_L=M_to_L_s4g)))
all_photometry.append(('infra 3.6 face-on', r_eff_36, mu_eff_36, n_36, mu_face_on(mu0d_36, cos_i), h_disc_36, M_to_L_s4g,
lambda l: s4g_surf_density(mu=mu_disc(l, mu0=mu_face_on(mu0d_36, cos_i), h=h_disc_36), M_to_L=M_to_L_s4g)))
Суммарная картинка:
In [74]:
for photom in all_photometry:
plt.plot(p_, map(photom[-1], p_), '-', label=photom[0])
plt.ylabel('$M_{sun}/{pc}^2$', fontsize=15.)
plt.legend(loc='best')
plt.legend()
Out[74]:
Неплохо.
In [75]:
show_all_photometry_table(all_photometry, scale)
Можно провести тест-сравнение с кривой вращения тонкого диска при заданной фотометрии, если она слишком массивная - то не брать ее (это ограничение сверху).
In [76]:
fig = plt.figure(figsize=[10,6])
plt.plot(_1, _2, 'd', label='HI Eymeren+2011')
plt.plot(test_points, spl_gas_N(test_points), '-', label='spline N')
plt.plot(r_n[:-20], vel_n[:-20], 'v', label='HI Noord+2005')
for photom in all_photometry:
plt.plot(test_points, map(lambda l: disc_vel(l, photom[7](0), photom[5], scale), test_points), '--', label=photom[0])
plt.ylim(0, 300)
plt.xlim(0, 300)
plt.legend(loc='best');
Все достаточно хорошо.
$H_{\alpha}$ , $UV$
есть $H_{\alpha}$ в Hameed 2005 http://iopscience.iop.org/article/10.1086/430211/pdf
ВАЖНО: в этой работе расстояние вообще 12.4 Мпк, соответственно неправильно пересчитываются масштабы
TODO: разобраться с расстояниями-масштабами, проверить в других галактиках
In [77]:
Image('n4725_halpha.png')
Out[77]:
In [78]:
Image('n4725_halpha_dist.png')
Out[78]:
Совмещенная картинка $H_{\alpha}$ + SDSS: (на самом деле совместил не очень качествено - видно например по трем звездам в центре)
In [79]:
Image('n4725_halpha_plus_sdss.png')
Out[79]:
В обратную сторону совмещенная:
In [80]:
Image('n4725_SDSS_labeled_plus_halpha.jpg')
Out[80]:
In [81]:
#SDSS
print 251*300/238. #arcsec
print 251*300/238.*0.130 #kpc, for distance 26.8 Mpc
#Hameed
print 161./9. #kpc
#distance ratio
print 26.8/12.4
print 41./17.9
Т.е. можно считать, что для исправления надо примерно *2.2 (чтобы привести к такому же масштабу, как SDSS), но тогда не очень точно получается спираль внешняя. Проще всего еще раз перемерить для SDSS картинки: (также здесь отмечены звезды, для которых я мерял выше расстояния)
In [82]:
Image('n4725_SDSS_labeled_sizes.jpg', width=600)
Out[82]:
In [83]:
def plot_SF(ax):
# как было
# ax.plot([0., 35.5/2./9./scale], [0., 0.], '-', lw=7., color='red')
# ax.plot([107./2./9./scale, 177./2./9./scale], [0., 0.], '-', lw=7., color='red')
# ax.plot([146./9./scale, 169./9./scale], [0., 0.], '-', lw=7., color='b') #внешняя спираль
ax.plot([0., 300./238 * 161./4.], [0., 0.], '-', lw=7., color='red') #очень примерно, потому что в SDSS не видно, но оно там есть
ax.plot([300./238 * 161./4., 300./238 *268./2.], [0., 0.], '--', lw=6., color='red', alpha=0.5) #в спитцере видно слабое
ax.plot([300./238 * 161./2., 300./238 *268./2.], [0., 0.], '-', lw=7., color='red')
ax.plot([290./238 *220., 300./238 *240], [0., 0.], '-', lw=7., color='b') #внешняя спираль, ширина условна
plot_SF(plt.gca())
plt.xlim(0, 350)
plt.ylim(0, 200)
Out[83]:
Видно, что теперь внешняя спираль находится на втором горбе данных, как и должно быть.
Есть еще ИК изображение от SPITZER http://www.spitzer.caltech.edu/images/2355-sig05-011-NGC-4725
Red represents warm dust clouds illuminated by newborn stars, while blue indicates older, cooler stellar populations.
In [84]:
Image('n4725_SPITZER.jpg', width=500)
Out[84]: