NGC 4725 (UGC 7989)

Галактика найдена из пересечения 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


Using matplotlib backend: Qt5Agg
Populating the interactive namespace from numpy and matplotlib

In [2]:
from photometry import *


importing Jupyter notebook from photometry.ipynb
Using matplotlib backend: Qt5Agg
Populating the interactive namespace from numpy and matplotlib

In [3]:
from instabilities import *


importing Jupyter notebook from instabilities.ipynb
Using matplotlib backend: Qt5Agg
Populating the interactive namespace from numpy and matplotlib

In [4]:
from utils import *


importing Jupyter notebook from utils.ipynb

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')


Оглавление

Статьи

TODO: add arcticles

Разное


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}$

TODO: понять, исправлено ли за наклон

Кривая вращения


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]:
(-0.5, 2.0)

Из Ноордермеера:


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]:
(-10.0, 11.0)

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


0.129930204938

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]:
(0, 5)

Сравнение:


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]:
<matplotlib.legend.Legend at 0xcbe5ba8>

Как и было изначально видно - пики заметно смещены (теперь нет, т.к. я это учел)

Разгадка кроется, кмк, в том что в работе 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


27.4434876449
C:\Anaconda\lib\site-packages\scipy\integrate\quadpack.py:364: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
  If increasing the limit yields no improvement it is advised to analyze 
  the integrand in order to determine the difficulties.  If the position of a 
  local difficulty can be determined (singularity, discontinuity) one will 
  probably gain from splitting up the interval and calling the integrator 
  on the subranges.  Perhaps a special-purpose integrator should be used.
  warnings.warn(msg, IntegrationWarning)

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


20.162389469

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


2.5017812108

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


8.15586889845

Данные по фотометрии


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]:
0.67889873586113669

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]:
<matplotlib.legend.Legend at 0xcc50518>

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]:
<matplotlib.legend.Legend at 0xf04f390>

Разница, как мы видим, не столь существенная.

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]:
<matplotlib.legend.Legend at 0x105f3198>

Вполне неплохо. Можем сравнить еще с калибровками McGaugh:


In [71]:
bv_color = 0.9
mcgaugh_mass_to_light(bv_color, 'mu36')


Out[71]:
array([ 0.47,  0.67,  1.03,  0.8 ])

Ну да, тут все похоже честно. Смущает конечно разница в центральной яркости с S4G, причем я не смог ничего найти насчет депроецировано или нет. Если депроецировать - становятся похожи:


In [72]:
mu_face_on(mu0d_36, cos_i)


Out[72]:
20.129831258118912

Добавим и депроецированные и нет:


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]:
<matplotlib.legend.Legend at 0xee0cef0>

Неплохо.


In [75]:
show_all_photometry_table(all_photometry, scale)


+------+-------------------+---------+----------+--------+---------+----------+-------+-------------+-----------+
|      | Name              |   r_eff |   mu_eff |      n |   mu0_d |   h_disc |   M/L | M_d/M_sun   |   Sigma_0 |
|------+-------------------+---------+----------+--------+---------+----------+-------+-------------+-----------|
| 0.00 | S4G 3.6           |   10.14 |   nan    |   2.21 |   20.34 |    73.20 |  0.68 | 7.47E+10.   |       286 |
| 1.00 | Heidt J           |  nan    |   nan    | nan    |   17.78 |    49.99 |  0.75 | 8.55E+10.   |       703 |
| 2.00 | Heidt H           |  nan    |   nan    | nan    |   17.11 |    50.28 |  0.69 | 1.10E+11.   |       891 |
| 3.00 | Heidt K           |  nan    |   nan    | nan    |   17.01 |    54.66 |  0.65 | 1.29E+11.   |       889 |
| 4.00 | infra 3.6         |   17.51 |    17.49 |   3.61 |   19.65 |    80.03 |  0.68 | 1.68E+11.   |       539 |
| 5.00 | infra 3.6 face-on |   17.51 |    17.49 |   3.61 |   20.13 |    80.03 |  0.68 | 1.08E+11.   |       347 |
+------+-------------------+---------+----------+--------+---------+----------+-------+-------------+-----------+

Сравнение с Кривой вращения тонкого диска

Можно провести тест-сравнение с кривой вращения тонкого диска при заданной фотометрии, если она слишком массивная - то не брать ее (это ограничение сверху).


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


316.386554622
41.1302521008
17.8888888889
2.16129032258
2.2905027933

Т.е. можно считать, что для исправления надо примерно *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]:
(0, 200)

Видно, что теперь внешняя спираль находится на втором горбе данных, как и должно быть.

Есть еще ИК изображение от 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]: