Галактика из диплома.
In [1]:
from IPython.display import HTML
from IPython.display import Image
import os
from PIL import Image as ImagePIL
%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 = 'N1167'
gtype = 'S0(A)' #CALIFA
incl = 36.0 #Zasov, in LEDA 49deg, CALIFA 42deg TODO: check, хотя разброс и не влияет на результат сильно
scale = 0.321 #kpc/arcsec from NED
data_path = '../../../data/ngc1167'
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)
In [8]:
# Данные из NED
HTML('<iframe src=http://ned.ipac.caltech.edu/cgi-bin/objsearch?objname=ngc+1167&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=ngc1167 width=1000 height=350></iframe>')
Out[9]:
In [10]:
fig, subplot = plt.subplots(1, 2, figsize=[12, 12])
subplot[0].imshow(np.asarray(ImagePIL.open("ngc1167_JHK.jpg")))
subplot[0].set_title("2MASS image JHK")
subplot[0].set_xticks([]);subplot[0].set_yticks([])
subplot[1].imshow(np.asarray(ImagePIL.open("ngc1167_SDSS.jpeg")))
subplot[1].set_title("SDSS DR9 whole image")
subplot[1].set_xticks([]);subplot[1].set_yticks([]);
Фотометрия в B и R дают согласованные значения центральной поверхностной плотности для диска. Галактика раннего типа. В диске наблюдаются туго закрученные тонкие ветви (SDSS). Излучение Hα только в центре. На профиле яркости в полосе B никаких особенностей. У галактики много газа. Интегрально MHI = 1.65·1010M⊙ . Но газ распределен широко, и значения поверхностной плотности низкие, даже в области газового кольца примерно на 60 arсseс. Газовая (WSRT) и звездная кинематика плохо согласуются (точки на звездной кривой вращения лежат выше точек на газовой кривой вращения). Профиль дисперсии скоростей σR восстанавливать через σlos .
p.37:
UGC 2487 (NGC 1167) is a nice example of an S0 galaxy with an extended and regularly rotating HI disk. The gas seems to be in circular motion at radii up till 80 kpc. The central hole does not reflect a true absence of gas, but is rather the result of absorption against a central continuum source.
p.110:
UGC 2487 (NGC 1167) is the most luminous galaxy in our sample ($M_{lim,R}$ = −23.24). It is an almost perfect superposition of an exponential disk and a Sersic bulge; the residuals with respect to the fitted bulge and disk profiles are less than 0.05 mag $arcsec^{−2}$ everywhere. Some very faint spiral structure is visible in the disk.
Сейфертовская галактика? p.165:
Category II contains galaxies with, for example, mild asymmetries, bar-induced streaming motions or signs of interactions or tidal distortions, as well as galaxies for which the orientation angles cannot be determined with sufficient accuracy. The following galaxies were classed as category II: UGC 2487 (Seyfert nucleus) ...
TODO: снять с большого количества кривых вращения кривую, это #2
да, Сейфертовская, p.183:
UGC 2487 (NGC 1167) is a giant S0 galaxy (MB = −21.88, D B 25 = 54 kpc; see table 3.4) with an extended, highly regular gas disk. We can trace the HI rotation curve out to radii of 80 kpc (10 R-band disk scale lengths) and although there is a small decline in the rotation velocities, they remain well above 300 km/s till the outermost point. The total mass enclosed within the last measured point is Menc = 2.1 · 1012 M, which makes UGC 2487 the most massive galaxy in our sample. The total enclosed mass is larger even than those in the giant Sc galaxies NGC 2916 and UGC 2885 (Rubin et al. 1979; Roelfsema & Allen 1985, note that in both papers a Hubble constant of H0 = 50 km s −1 Mpc−1 is assumed; their derived masses have to be divided by 1.5 when using our value of 75 km s −1 Mpc−1 ); to our knowledge, it is the largest mass ever derived from a rotation curve. Saglia & Sancisi (1988) list a number of other large disk galaxies with extremely high rotation velocities; some of those galaxies may be even more massive than UGC 2487, but since no spatially resolved rotation curves are available for these systems, no accurate values for the total masses can be derived. In any case, UGC 2487 seems member of a class of extremely massive disk galaxies (see also Giovanelli et al. 1986; Carignan et al. 1997), with masses that rival those of the most massive elliptical galaxies (e.g. Bertin et al. 1988; Minniti et al. 1998).
UGC 2487 is also classified as a Seyfert galaxy, explaining the broad emission lines in the nucleus (cf. Filippenko & Sargent 1985). It has a central compact steep spectrum (CSS) radio source (e.g. Sanghera et al. 1995; Giovannini et al. 2001), which is responsible for the H absorption in the center. Away from the bright nucleus, we detect some very faint emission in the optical spectrum. Although this emission seems to follow the general sense of rotation of the galaxy, the emission profiles are broad and do not have well-defined peaks. From these data alone, it is difficult to determine whether this faint emission traces regular rotation in the circumnuclear regions, or whether it is related to outflows from the active nucleus. Thus, this emission gives no useful information on the shape of the potential in the inner regions and we have decided not to use it in the derivation of the rotation curve. A small HII region is detected 3000 away from the center on the approaching side; the emission from this region has regular line profiles and its velocity is consistent with the rotation velocities of the H at the corresponding location.
In [11]:
Image('noord_pics/instab.png')
Out[11]:
In [12]:
r_q_n, q_n = zip(*np.loadtxt("Q1F_from_webplotdigitizer.dat", float, delimiter=','))
fig = plt.figure(figsize=[8,2])
plt.plot(r_q_n, q_n, 'o-')
plt.xlim(0, 80)
plt.ylim(0, 1.2);
TODO: добавить остальные данные Ноордермеера и подчистить их
In [13]:
Image('noord_pics/rot_vel.png', width=400)
Out[13]:
In [14]:
r_ma_n, vel_ma_n = zip(*np.loadtxt("v_gas_from_webplotdigitizer.dat", float, delimiter=','))
plt.plot(r_ma_n, vel_ma_n, 's-')
plt.xlim(0, 300)
plt.ylim(0, 450);
TODO: статья Струве, просмотреть
Данные Засов (данные Катков или Моисеев):
TODO: добавить данные Засова в работу
In [15]:
Image('zasov_pics/rot_vel.png')
Out[15]:
In [16]:
r_ma_z, vel_ma_z = zip(*np.loadtxt("v_stars_from_webplotdigitizer_zasov.dat", float, delimiter=','))
plt.plot(r_ma_z, vel_ma_z, 's')
plt.xlim(0, 120)
plt.ylim(0, 420);
In [17]:
# Данные по звездной кинематике Засова 2012 вдоль большей полуоси, не исправленные за наклон
r_ma, vel_ma, e_vel_ma, sig_ma, e_sig_ma = zip(*np.loadtxt("v_stars_maZ.dat", float))
# Данные по звездной кинематике Засова 2012 вдоль малой полуоси, не исправленные за наклон
r_mi, vel_mi, e_vel_mi, sig_mi, e_sig_mi = zip(*np.loadtxt("v_stars_miZ.dat", float))
plt.plot(r_ma, vel_ma, '.-', label="Zasov 2008, maj")
plt.plot(r_mi, vel_mi, '.-', label="Zasov 2008, min")
plt.legend();
In [20]:
%run ../../utils/rotvelutils #import util functions for rotcurv correction
In [21]:
v0 = 4959.3 #TODO: вспомнить, как нормально называется эта величина
r_ma_b, vel_ma_b, e_vel_b = correct_rotation_curve(r_ma, vel_ma, e_vel_ma, 0.0, v0, incl)
r_mi_b, vel_mi_b, e_vel_mi_b = correct_rotation_curve(r_mi, vel_mi, e_vel_mi, 0.0, v0, incl)
fig = plt.figure(figsize=[10,8])
plt.errorbar(r_ma_b, vel_ma_b, yerr=e_vel_b, fmt='.', marker='.', mew=0, color='blue', label = 'Zasov star maj')
plt.errorbar(r_mi_b, vel_mi_b, yerr=e_vel_mi_b, fmt='.', marker='.', mew=0, color='green', label = 'Zasov star min')
plt.plot(r_ma_z, vel_ma_z, 'o', color='m', label='Zasov from fig')
plt.legend()
plt.ylim(0, 450);
Похоже, но не слишком. Ладно, большой роли в анализе это играть не должно.
Теперь приближения:
In [ ]:
fig = plt.figure(figsize=[10,8])
plt.errorbar(r_ma_b, vel_ma_b, yerr=e_vel_b, fmt='.', marker='.', mew=0, color='blue', label = 'Zasov star maj')
plt.plot(r_ma_z, vel_ma_z, 'o', color='m', label='Zasov from fig')
test_points = np.linspace(0.0, max(r_ma_b), 100)
poly_star = poly1d(polyfit(r_ma_z, vel_ma_z, deg=3))
plt.plot(test_points, poly_star(test_points), '-', label='poly deg=3*')
poly_star = poly1d(polyfit(r_ma_b, vel_ma_b, deg=3))
plt.plot(test_points, poly_star(test_points), '-', label='poly deg=3')
def w(arr):
return map(lambda l: 1/(1. + l**2), arr)
import scipy.interpolate as inter
spl = inter.UnivariateSpline(r_ma_b, vel_ma_b, k=3, s=100., w=w(e_vel_b))
plt.plot(test_points, spl(test_points), '-', label='spl s=100 w^2')
plt.legend(loc='lower right')
plt.ylim(0, 450);
TODO: посмотреть аппроксимации, выбрать наиболее нормальную, границы опять же
In [ ]:
star_approx = poly1d(polyfit(r_ma_b, vel_ma_b, deg=3))
In [ ]:
#TODO: свидетельство того, что надо раздвигать
# Исправляем значения вдоль малой оси на синус угла:
def correct_min(R):
return R / cos_i
plt.errorbar(r_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='.', mew=0, color='blue', label='$\sigma_{los}^{maj}$')
r_mi_extend = map(correct_min, r_mi)
plt.errorbar(r_mi_extend, sig_mi, yerr=e_sig_mi, fmt='.', marker='.', mew=0, color='black', label='$\sigma_{los}^{min}$')
plt.xlabel('$R$')
plt.ylabel('$\sigma$')
plt.legend();
In [ ]:
bind_curve = lambda p: (abs(p[0]), abs(p[1]), p[2])
sig_maj_data = sorted(zip(r_ma, sig_ma, e_sig_ma))
sig_maj_data = map(bind_curve, sig_maj_data)
sig_maj_data.sort()
radii_maj, sig_maj_p, e_sig_maj_p = zip(*sig_maj_data)
sig_min_data = sorted(zip(r_mi_extend, sig_mi, e_sig_mi))
sig_min_data = map(bind_curve, sig_min_data)
sig_min_data.sort()
radii_min, sig_min_p, e_sig_min_p = zip(*sig_min_data)
points = np.arange(0, max(radii_min), 0.1)
plt.errorbar(radii_maj, sig_maj_p, yerr=e_sig_maj_p, fmt='.', marker='.', mew=0, color='blue', label='$\sigma_{los}^{maj}$')
plt.errorbar(radii_min, sig_min_p, yerr=e_sig_min_p, fmt='.', marker='.', mew=0, color='red', label='$\sigma_{los}^{min}$')
plt.legend()
plt.ylim(0,250)
plt.xlim(0,70);
In [ ]:
r_sig_ma, sig_ma, e_sig_ma = radii_maj, sig_maj_p, e_sig_maj_p
r_sig_mi, sig_mi, e_sig_mi = radii_min, sig_min_p, e_sig_min_p
spl_min = inter.UnivariateSpline(r_sig_mi, sig_mi, k=3, s=7000.)
sig_min_lim = max(r_sig_mi)
# hack for smooth at the end
spl_maj = inter.UnivariateSpline(r_sig_ma[:-2] + tuple(np.linspace(53, 60, 10)), sig_ma[:-2] + tuple(np.linspace(90, 80, 10)), k=3, s=10000.)
# 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_sig_ma)+15., 100)
In [ ]:
@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 [ ]:
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 [ ]:
def sigPhi_to_sigR_real(R):
try:
return 0.5 * (1 + R*star_approx.derivative()(R) / star_approx(R))
except:
return 0.5 * (1 + R*star_approx.deriv()(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 [ ]:
plt.errorbar(r_sig_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='.', mew=0, color='red', label='$\sigma_{los}^{maj}$')
plt.plot(points, map(sig_R_maj_minmin, points), label = 'minmin')
plt.plot(points, map(sig_R_maj_min, points), label = 'min')
plt.plot(points, map(sig_R_maj_max, points), label = 'max')
plt.plot(points, map(sig_R_maj_maxmax, points), label = 'maxmax')
plt.plot(points, map(sig_R_maj_maxmaxtrue, points), label = 'maxmaxtrue')
plt.legend()
plt.ylim(0,400)
plt.xlim(0,100);
Сравним major vs minor оценки:
In [ ]:
fig = plt.figure(figsize=[10, 7])
plt.plot(points, map(sig_R_maj_minmin, points), label = 'maj minmin')
plt.plot(points, map(sig_R_maj_min, points), label = 'maj min')
plt.plot(points, map(sig_R_maj_max, points), label = 'maj max')
plt.plot(points, map(sig_R_maj_maxmax, points), label = 'maj maxmax')
plt.plot(points, map(sig_R_maj_maxmaxtrue, points), label = 'maj maxmaxtrue')
plt.plot(points, map(sig_R_minor_minmin, points), '--', label = 'minor minmin')
plt.plot(points, map(sig_R_minor_min, points), '--', label = 'minor min')
plt.plot(points, map(sig_R_minor_max, points), '--', label = 'minor max')
plt.plot(points, map(sig_R_minor_maxmax, points), '--', label = 'minor maxmax')
plt.legend()
plt.ylim(50,400)
plt.xlim(0,100);
Достаточно кривое получилось, ну да не страшно.
Нанесем также из пред работы (правда в предыдущей работе были получены величины $\sigma_R$ на двух участках, которые не особо согласуются друг с другом, что плохо):
In [ ]:
am1 = 0.72
dam1 = 0.09
sm1 = 281.
dsm1 = 28.
hkin1 = 52.6
am2 = 0.3
dam2 = 0.08
sm2 = 222.
dsm2 = 35.
hkin2 = 1343.
def sigR_ger_exp(R, sig_R_0, h_kin):
return sig_R_0*exp(-R/h_kin)
fig = plt.figure(figsize=[10, 7])
points_ = np.arange(6.7, 30., 0.1)
# alpha = am1
# sig_R_0 = sm1
# h_kin = hkin1
# plt.plot(points, [sigR_ger_exp(l, sm1, hkin1) for l in points], ls='-', color='red')
# alpha = am1-dam1
# sig_R_0 = sm1-dsm1
# plt.plot(points, [sigR_ger_exp(l) for l in points], ls='--', color='red')
# alpha = am1+dam1
# sig_R_0 = sm1+dsm1
# plt.plot(points, [sigR_ger_exp(l) for l in points], ls='--', color='red')
plt.fill_between(points_, [sigR_ger_exp(l, sm1-dsm1, hkin1) for l in points_], [sigR_ger_exp(l, sm1+dsm1, hkin1) for l in points_], color='red', alpha=0.3)
points_ = np.arange(30., 60., 0.1)
# alpha = am2
# sig_R_0 = sm2
# h_kin = hkin2
# plt.plot(points, [sigR_ger_exp(l, sm2, hkin2) for l in points], ls='-', color='red')
# alpha = am2-dam2
# sig_R_0 = sm2-dsm2
# plt.plot(points, [sigR_ger_exp(l) for l in points], ls='--', color='red')
# alpha = am2+dam2
# sig_R_0 = sm2+dsm2
# plt.plot(points, [sigR_ger_exp(l) for l in points], ls='--', color='red')
plt.fill_between(points_, [sigR_ger_exp(l, sm2-dsm2, hkin2) for l in points_], [sigR_ger_exp(l, sm2+dsm2, hkin2) for l in points_], color='red', alpha=0.3)
plt.plot(points, map(sig_R_maj_minmin, points), label = 'maj minmin')
plt.plot(points, map(sig_R_maj_min, points), label = 'maj min')
plt.plot(points, map(sig_R_maj_max, points), label = 'maj max')
plt.plot(points, map(sig_R_maj_maxmax, points), label = 'maj maxmax')
plt.plot(points, map(sig_R_maj_maxmaxtrue, points), label = 'maj maxmaxtrue')
plt.plot(points, map(sig_R_minor_minmin, points), '--', label = 'minor minmin')
plt.plot(points, map(sig_R_minor_min, points), '--', label = 'minor min')
plt.plot(points, map(sig_R_minor_max, points), '--', label = 'minor max')
plt.plot(points, map(sig_R_minor_maxmax, points), '--', label = 'minor maxmax')
plt.legend()
plt.ylim(50,400)
plt.xlim(0,100);
Не смотря на то, что расходится - не такое уж и плохое согласие.
In [ ]:
# Данные по кинематике газа Struve, WSRT (не исправлено за наклон)
r_wsrt, vel_wsrt, e_vel_wsrt = zip(*np.loadtxt("v_gas_WSRT.dat", float))
r_wsrt, vel_wsrt, e_vel_wsrt = correct_rotation_curve(r_wsrt, vel_wsrt, e_vel_wsrt, 0.0, v0, incl)
# Данные по кинематике газа Noordermee 2007, WSRT (не исправлено за наклон?)
r_noord, vel_noord, e_vel_noord = zip(*np.loadtxt("v_gas_noord.dat", float))
r_noord, vel_noord, e_vel_noord = correct_rotation_curve(r_noord, vel_noord, e_vel_noord, 0.0, v0, incl)
gas_approx = poly1d(polyfit(r_ma_n, vel_ma_n, deg=5))
test_points = np.linspace(min(r_ma_n), max(r_ma_n), 100)
plt.plot(test_points, gas_approx(test_points), '--', label='approx')
spl_gas = inter.UnivariateSpline(r_noord, vel_noord, k=3, s=10000.)
plt.plot(test_points, spl_gas(test_points), '-', label='spl')
plt.plot(r_wsrt, vel_wsrt, '.-', label="gas Struve")
plt.plot(r_noord, vel_noord, '.-', label="gas Noordermeer 2007")
plt.plot(r_ma_n, vel_ma_n, 's-', label='HI from fig')
plt.ylim(0, 450)
plt.legend(loc='lower right');
Оно все конечно согласуется, т.к. это все одно и то же
Для случая бесконечного тонкого диска: $$\kappa=\frac{3}{R}\frac{d\Phi}{dR}+\frac{d^2\Phi}{dR^2}$$ где $\Phi$ - гравпотенциал, однако его знать не надо, т.к. есть проще формула: $$\kappa=\sqrt{2}\frac{\vartheta_c}{R}\sqrt{1+\frac{R}{\vartheta_c}\frac{d\vartheta_c}{dR}}$$
In [ ]:
#использовавшееся в дипломе приближение
ef = [104.55995900171898, 54.292488898429013, 34.856981586343409, 25.282657419650114, 19.107329429755854, 13.969314301465111,
10.21292602770659, 7.466641224935947, 5.458840202177822, 3.990942574473094, 2.9177667861366507, 2.133171013969345,
1.5595552723609254, 1.140186432133711, 0.8335870636080539]
fig = plt.figure(figsize=[12, 8])
plt.plot(test_points, [epicyclicFreq_real(gas_approx, x, scale) for x in test_points], '-')
plt.plot(test_points, [epicyclicFreq_real(spl_gas, x, scale) for x in test_points], '-')
# plt.plot(r_g_dens[1:], ef, '--')
plt.xlabel('$R, arcsec$')
plt.ylabel('$\kappa,\, km/s/kpc$', fontsize=15);
Плотность HI:
In [ ]:
r_g_dens, gas_dens = zip(*np.loadtxt("gas_density.dat", float))
plt.plot(r_g_dens, gas_dens, 's-');
В этой работе https://arxiv.org/pdf/astro-ph/0610453v1.pdf (B2 0258+35) есть $\Sigma T_a dv$ = 0.88.
In [ ]:
0.88*3.2
Полная масса $H_2$ получается для $R$
In [ ]:
2*math.pi * 24.20**2 * (scale * 1000.)**2 * 2.81/1e9
Не уверен, что это то, что нужно. Общая выведенная масса $H_2$ у них в работе - log M_H2 = 8.1, правда постоянаая Хаббла $H_0 = 100$.
In [ ]:
np.power(10., 8.1)/1e9
Если изменить постоянную Хаббла:
In [ ]:
(100./75.)**2 * np.power(10., 8.1)/1e9
Центральное значение:
In [ ]:
0.22*1e9/(2*math.pi*24.20**2 * (scale * 1000.)**2)
Мало конечно, но что делать.
Масса атомарного водорода у Ноордермеера $HI$ - $17.09\times10^9 M_{sun}$. Проверим, сколько у нас водорода:
In [ ]:
import scipy.integrate as integrate
import scipy.interpolate
tmp_ = scipy.interpolate.interp1d(r_g_dens, gas_dens)
result = integrate.quad(lambda l: 2*np.pi*l*tmp_(l), r_g_dens[0], r_g_dens[-1])
print (scale * 1000.)**2 * result[0]/1e9
В принципе похоже, оть и не должно так отличаться.
Из заметок:
Для справки: Использованная полная масса COCO для 1167 не согласуется с этой работой https://arxiv.org/pdf/1510.00729.pdf (у нас в 4 раза больше). Причем в этой https://arxiv.org/pdf/1408.7106.pdf работе для 1167 тоже в 4 раза меньше. Плюс там есть заметка:
NGC 1167 has a very large HI disk (Struve et al. 2010), and a large dust mass. The molecular gas derived from the detected CO emission is however modest. Given the usual exponential radial distribution of CO emission in galaxies (e.g., Young & Scoville 1991), it is not likely that we are missing some CO emission outside of the 7.8kpc beam. The galaxy is currently active, hosting the compact steep spectrum radio source B2 0258+35, and shows signs of previous AGN outbursts in the forms of large-scale (∼10′ / 200 kpc) low surface brightness lobes (Shulevski et al. 2012).
и правда, у нас меньше примерно в 1.5 раза, а эти две оценки согласуются. Не думаю, что это важно конечно.
In [ ]:
np.power(10., 8.5)/1e9, 3.3*1e8/1e9
Фотометрий две - в R и в B, которые вроде как согласуются (для максимального диска в R - 4)
In [ ]:
all_photometry = []
In [ ]:
Image('diplom_photom.png')
In [ ]:
Image('photometry_pdf.png')
In [ ]:
M_R = -22.94 #11.69 - это правильно? надо брать абсолютные? в дипломе были относительные, тут разница уже существенная
M_B = -21.54 #13.40
In [ ]:
print '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'))
In [ ]:
Image('noord_pics/n1167_photom.png')
In [ ]:
r_phot, mu_phot = zip(*np.loadtxt("n1167_noord_photoRB.dat", float, delimiter=','))
plt.plot(r_phot, mu_phot, 's')
plt.xlim(-10, 175)
plt.ylim(30, 15);
In [ ]:
# R-band
r_eff_R = 6.7
mu_eff_R = 19.89 # уточнить это ли число
n_R = 1.7
mu0d_R = 20.35 # и тут тоже
h_disc_R = 24.2
mu_eff_Rc = 19.40 # уточнить это ли число
mu0d_Rc = 20.12 # и тут тоже
In [ ]:
# B-band
r_eff_B = 6.7
mu_eff_B = 21.52 # уточнить это ли число
n_B = 1.7
mu0d_B = 22.24 # и тут тоже
h_disc_B = 27.5
mu_eff_Bc = 20.73 # уточнить это ли число
mu0d_Bc = 21.71 # и тут тоже
In [ ]:
p_ = np.arange(0., 150., 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.xlim(-10, 150)
plt.ylim(30, 15)
plt.legend();
Отлично, похоже на правду. Массовые модели:
In [ ]:
b_r_color = M_B-M_R
M_to_L_Rc = 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_Rc, band='R') for l in p_]
plt.plot(p_, surf_R, '-', label='R [M/L={:2.2f}] absmag'.format(M_to_L_Rc))
M_to_L_Bc = 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_Bc, band='B') for l in p_]
plt.plot(p_, surf_B, '-', label='B [M/L={:2.2f}] absmag'.format(M_to_L_Bc))
diploma_color = 13.40 - 11.69
M_to_L_R = 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_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(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_B, band='B') for l in p_]
plt.plot(p_, surf_B, '-', label='B [M/L={:2.2f}]'.format(M_to_L_B))
# максимальный диск
surf_R_max = [surf_density(mu=mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L=4., band='R') for l in p_]
plt.plot(p_, surf_R_max, '--', label='R [M/L={:2.2f}]'.format(4.))
plt.legend();
Так получились схожими и вполне не завышенными, M/L меньше ноордермееровского максимального.
В http://aas.aanda.org/articles/aas/pdf/2000/06/h1643.pdf она обозначена как 0258+35 и есть очень короткая $V$ фотометрия (3 секунды), соответственно диск оттуда не извлечь.
Тут http://adsabs.harvard.edu/cgi-bin/bib_query?2003AJ....125.2307A есть какая-то декомпозиция, но нет данных.
In [ ]:
all_photometry.append(('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(('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')))
#максимальный диск
all_photometry.append(('Noorder R_max', r_eff_R, mu_eff_Rc, n_R, mu0d_Rc, h_disc_R, 4.,
lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L=4., band='R')))
CALIFA $g$, $r$, $i$:
sdss http://skyserver.sdss.org/dr13/en/tools/chart/navi.aspx coords 45.42666667, 35.20583333
http://skyserver.sdss.org/dr13/en/tools/explore/Summary.aspx?id=1237670959481356507
Тут данные не исправлены за поглощение, как написал Сережа, надо взять FOREGROUND GALACTIC EXTINCTION из NED и исправить.
Раньше я брал цвета из SDSS таблицы, но это было неверно и к тому же там было предупреждение что данные не точны.
In [ ]:
A_lambda_g = 0.598
A_lambda_r = 0.414
A_lambda_i = 0.307
In [ ]:
MUE_g=21.46976
Re_g=7.57209
n_g=2.20557
MU0_g=21.89252-A_lambda_g
hi_g=24.4525
MUE_r=20.7038
Re_r=9.09017
n_r=2.6612
MU0_r=21.0106-A_lambda_r
hi_r=25.8222
MUE_i=20.6073
Re_i=11.4278
n_i=3.20625
MU0_i=20.624-A_lambda_i
hi_i=27.508
TODO: правильно ли исправлено за поглощение
TODO: правильно ли исправлены яркости
In [ ]:
M_g = disc_totmag(MU0_g, hi_g, scale)
M_r = disc_totmag(MU0_r, hi_r, scale)
M_i = disc_totmag(MU0_i, hi_i, scale)
fig = plt.figure(figsize=[16, 10])
sdss_surf = []
for color_desc, color in [('g-r', M_g-M_r), ('g-i', M_g-M_i) , ('r-i', M_r-M_i)]:
surf_g = [surf_density(mu=mu_disc(l, mu0=MU0_g, h=hi_g), M_to_L=bell_mass_to_light(color, 'g', color_desc), band='g') for l in p_]
plt.plot(p_, surf_g, '-', label='g [M/L={:2.2f}] {}'.format(bell_mass_to_light(color, 'g', color_desc), color_desc))
sdss_surf.append([surf_g[0], bell_mass_to_light(color, 'g', color_desc), 'g', color_desc])
surf_r = [surf_density(mu=mu_disc(l, mu0=MU0_r, h=hi_r), M_to_L=bell_mass_to_light(color, 'r', color_desc), band='r') for l in p_]
plt.plot(p_, surf_r, '-', label='r [M/L={:2.2f}] {}'.format(bell_mass_to_light(color, 'r', color_desc), color_desc))
sdss_surf.append([surf_r[0], bell_mass_to_light(color, 'r', color_desc), 'r', color_desc])
surf_i = [surf_density(mu=mu_disc(l, mu0=MU0_i, h=hi_i), M_to_L=bell_mass_to_light(color, 'i', color_desc), band='i') for l in p_]
plt.plot(p_, surf_i, '-', label='i [M/L={:2.2f}] {}'.format(bell_mass_to_light(color, 'i', color_desc), color_desc))
sdss_surf.append([surf_i[0], bell_mass_to_light(color, 'i', color_desc), 'i', color_desc])
plt.legend();
In [ ]:
sorted(sdss_surf)
Добавим экстремальные представители в таблицу (т.к. более-менее равномерно заполняет плоскость и чтобы были представлены все три модели):
In [ ]:
all_photometry.append(('califa g (g-i)', Re_g, MUE_g, n_g, MU0_g, hi_g, 5.59,
lambda l: surf_density(mu=mu_disc(l, mu0=MU0_g, h=hi_g), M_to_L=5.59, band='g')))
all_photometry.append(('califa r (g-r)', Re_r, MUE_r, n_r, MU0_r, hi_r, 3.88,
lambda l: surf_density(mu=mu_disc(l, mu0=MU0_r, h=hi_r), M_to_L=3.88, band='r'))) #модель из середины
all_photometry.append(('califa i (r-i)', Re_i, MUE_i, n_i, MU0_i, hi_i, 2.95,
lambda l: surf_density(mu=mu_disc(l, mu0=MU0_i, h=hi_i), M_to_L=2.95, band='i')))
In [ ]:
show_all_photometry_table(all_photometry, scale)
Можно провести тест-сравнение с кривой вращения тонкого диска при заданной фотометрии, если она слишком массивная - то не брать ее (это ограничение сверху).
In [ ]:
fig = plt.figure(figsize=[10,6])
plt.plot(test_points, spl_gas(test_points), '-', label='spl')
plt.plot(r_ma_n, vel_ma_n, 's-', label='HI from fig')
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, 450)
plt.xlim(0, 300)
plt.legend(loc='best');
Все более-менее одинаковое.
Spiral-like star-forming patterns in CALIFA early-type galaxies http://adsabs.harvard.edu/cgi-bin/bib_query?2016A%26A...585A..92G $H_{\alpha}$ (нижняя картинка) и SDSS + слабый SFR (верхняя), все это соответствует SFR 0.1-0.3 $M_{o}$ в год:
In [ ]:
Image('califa_2016_halpha_plus_sfr.png')
In [ ]:
# TODO: разобраться и сделать обоснованно
def plot_SF(ax):
# ax.plot([55., 65.], [0., 0.], '-', lw=7., color='red')
ax.plot([22., 38.], [0., 0.], '-', lw=7., color='red')
plot_SF(plt.gca())
plt.xlim(0, 350)
plt.ylim(0, 200)
In [ ]:
# from instabilities import *
In [ ]:
Image('diplom_results.png')
Устойчиво, когда > 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 [ ]:
sound_vel = 6. #скорость звука в газе, км/с
data_lim = max(max(r_sig_ma), max(r_sig_mi)) #где заканчиваются данные
In [ ]:
fig = plt.figure(figsize=[12, 6])
gd_data = zip(r_g_dens, gas_dens)
plt.plot(r_g_dens, [1./Qg(epicycl=l[0], sound_vel=sound_vel, gas_density=l[1]) for l in
zip([epicyclicFreq_real(gas_approx, x, scale) for x in r_g_dens], gas_dens)], 's-', label='$Q_{gas}$')
plt.plot(r_g_dens, [1./Qg(epicycl=l[0], sound_vel=4., gas_density=l[1]) for l in
zip([epicyclicFreq_real(gas_approx, x, scale) for x in r_g_dens], gas_dens)], 's-', label='$Q_{gas}$ & c=4')
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],
[sig_R_maj_max(_) for _ in r_g_dens],
[surf_density(l_, M_to_L_R, 'R') for l_ in [mu_disc(ll, mu0=mu0d_R, h=h_disc_R) for ll in r_g_dens]])], 's-', label='$Q_{star}^{max}$')
plt.plot(r_g_dens, [1./Qs(epicycl=l[0], sigma=l[1], star_density=l[2]) for l in
zip([epicyclicFreq_real(gas_approx, x, scale) for x in r_g_dens],
[sig_R_maj_min(_) for _ in r_g_dens],
[surf_density(l_, M_to_L_R, 'R') for l_ in [mu_disc(ll, mu0=mu0d_R, h=h_disc_R) for ll in r_g_dens]])], 's-', label='$Q_{star}^{min}$')
plt.plot(map(lambda l: l/scale, r_q_n), q_n, 'o-', label='noordermeer 1F')
plt.axhline(y=1, ls='--')
plt.legend()
plot_SF(plt.gca())
plt.ylabel('$Q^{-1}$', fontsize=15);
У Noordermeer значения больше, потому что похоже он использовал значение скорости звука меньше $c=4\, km/s$
НЕ ИСПРАВЛЕНО ЗА 1.6! Т.к. тут еще сравнение с звездным Q. И дальше тоже не будет исправлено, поскольку эту константу трудно определить.
Кинетическое приближение: $$\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 [ ]:
total_gas_data = zip(r_g_dens, map(lambda l: He_coeff*l, gas_dens))[1:]
disk_scales = [(l[5], l[0].split(' ')[1]) for l in all_photometry]
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), M_to_L_R, 'R'),
star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L_R, 'R'),
data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='R (not abs mag) maj max/min')
plt.ylim(0., 2.5)
plt.axhline(y=1., ls='-', color='grey')
plot_SF(ax)
plt.grid();
Сравним влияние наклона и как выглядит максимальный диск.
In [ ]:
fig = plt.figure(figsize=[10, 6])
ax = plt.gca()
incl = 49.0
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_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L_R, 'R'),
star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L_R, 'R'),
data_lim=data_lim, color='m', alpha=0.3, disk_scales=disk_scales, label='R (not abs mag) maj max/min i=49deg')
incl = 36.0
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_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L_R, 'R'),
star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L_R, 'R'),
data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='R (not abs mag) maj max/min i=36deg')
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), 4., 'R'),
star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), 4., 'R'),
data_lim=data_lim, color='b', alpha=0.3, disk_scales=disk_scales, label='R (max disk) maj max/min')
plt.ylim(0., 2.5)
plt.axhline(y=1., ls='-', color='grey')
plot_SF(ax)
plt.grid();
In [ ]:
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=all_photometry[i][-1],
star_density_min=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 [ ]:
# anim.save('..\\..\pics\\'+name+'.gif', writer='imagemagick', fps=1)
In [ ]:
# from IPython.display import HTML
# HTML(anim.to_html5_video())
Существует ограничение на максимальный диск в ~0.85 (изотермическое гало) и на субмаксимальный в 0.55-0.6 (NFW гало). Попробуем дотянуть фотметрию до максимальных дисков и посмотрим, как изменятся M/L (скорость зависит как корень из M/L):
In [ ]:
fig = plt.figure(figsize=[14,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 photom in all_photometry:
disc_max = 2.2*photom[5]
max_coeff = 0.85*spl_gas(disc_max)/disc_vel(disc_max, photom[7](0), photom[5], scale)
submax_coeff = 0.6*spl_gas(disc_max)/disc_vel(disc_max, photom[7](0), photom[5], scale)
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, 450)
plt.xlim(0, 300)
plt.legend(bbox_to_anchor=(1.15, 1.0));
Опять же, ничего сильно не выделяется.
In [ ]:
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=zip(r_g_dens, map(lambda l: He_coeff*l, gas_dens))[1:],
epicycl=epicyclicFreq_real, gas_approx=spl_gas, sound_vel=sound_vel, scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: max_coeffs[all_photometry[i][0]][0]*all_photometry[i][-1](l),
star_density_min=lambda l: max_coeffs[all_photometry[i][0]][0]*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=zip(r_g_dens, map(lambda l: He_coeff*l, gas_dens))[1:],
epicycl=epicyclicFreq_real, gas_approx=spl_gas, sound_vel=sound_vel, scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: max_coeffs[all_photometry[i][0]][1]*all_photometry[i][-1](l),
star_density_min=lambda l: max_coeffs[all_photometry[i][0]][1]*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)
ax.set_xlim(0., 160.)
plot_SF(ax)
return ax
anim = FuncAnimation(plt.gcf(), animate, frames=len(all_photometry), interval=1000)
In [ ]:
# anim.save('..\\..\pics\\'+name+'_MAXDISCS.gif', writer='imagemagick', fps=1)
In [ ]:
# from IPython.display import HTML
# HTML(anim.to_html5_video())
Есть очень неплохие модели.
In [ ]:
fig = plt.figure(figsize=[10, 4])
# def y_interp_(r, h):
# return 2.8*np.exp(-r/h)
def y_interp_(r, h):
return 0.22*1e9/(2*math.pi*h**2 * (scale * 1000.)**2) * np.exp(-r/h)
points = np.linspace(0.1, 100., 100.)
plt.plot(r_g_dens, gas_dens, 's-')
plt.plot(points, y_interp_(points, 24.), '--', alpha=0.5)
plt.plot(r_g_dens, map(lambda l: y_interp_(l[0], 24.) + l[1], zip(r_g_dens, gas_dens)), '--', alpha=0.5);
In [ ]:
fig = plt.figure(figsize=[10, 6])
total_gas_data_ = zip(r_g_dens, [He_coeff*(y_interp_(l[0], h_disc_R) + l[1]) for l in zip(r_g_dens, gas_dens)])
ax = plt.gca()
plot_2f_vs_1f(ax=ax, total_gas_data=total_gas_data_[1:], epicycl=epicyclicFreq_real,
gas_approx=spl_gas, sound_vel=sound_vel, scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), 4., 'R'),
star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), 4., 'R'),
data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='R (max disc) maj max/min')
plt.ylim(0., 2.5)
plt.axhline(y=1., ls='-', color='grey')
plot_SF(ax)
plt.grid();
Остальные фотометрии:
In [ ]:
from matplotlib.animation import FuncAnimation
fig = plt.gcf()
plt.figure(figsize=(10,6))
ax = plt.gca()
def animate(i):
ax.cla()
total_gas_data_ = zip(r_g_dens, [He_coeff*(y_interp_(l[0], all_photometry[i][-3]) + l[1]) for l in zip(r_g_dens, gas_dens)])
plot_2f_vs_1f(ax=ax, total_gas_data=total_gas_data_[1:], epicycl=epicyclicFreq_real, gas_approx=spl_gas, sound_vel=sound_vel, scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=all_photometry[i][-1],
star_density_min=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 [ ]:
# anim.save('..\\..\pics\\'+name+'_molec.gif', writer='imagemagick', fps=1)
In [ ]:
# from IPython.display import HTML
# HTML(anim.to_html5_video())
И оценки с максимальным диском:
In [ ]:
from matplotlib.animation import FuncAnimation
fig = plt.gcf()
plt.figure(figsize=(10,6))
ax = plt.gca()
def animate(i):
ax.cla()
total_gas_data_ = zip(r_g_dens, [He_coeff*(y_interp_(l[0], all_photometry[i][-3]) + l[1]) for l in zip(r_g_dens, gas_dens)])
plot_2f_vs_1f(ax=ax, total_gas_data=total_gas_data_[1:], epicycl=epicyclicFreq_real, gas_approx=spl_gas, sound_vel=sound_vel, scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: max_coeffs[all_photometry[i][0]][0]*all_photometry[i][-1](l),
star_density_min=lambda l: max_coeffs[all_photometry[i][0]][0]*all_photometry[i][-1](l),
data_lim=data_lim, color=np.random.rand(3), alpha=0.3, disk_scales=[(all_photometry[i][5], '')])
ax.axhline(y=1., ls='-', color='grey')
ax.set_title(all_photometry[i][0])
ax.set_ylim(0., 2.5)
return ax
anim = FuncAnimation(plt.gcf(), animate, frames=len(all_photometry), interval=1000)
In [ ]:
# anim.save('..\\..\pics\\'+name+'_molec_MAX.gif', writer='imagemagick', fps=1)
In [ ]:
# from IPython.display import HTML
# HTML(anim.to_html5_video())
In [ ]:
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('ngc1167_SDSS.jpeg'), aspect='auto')
axes[0].set_title(name)
try:
axes[1].errorbar(r_sig_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='.', mew=0, color='red', label='$\sigma_{los}^{maj}$')
axes[1].plot(points, map(sig_R_maj_min, points))
axes[1].plot(points, map(sig_R_maj_max, points))
axes[1].plot(points, map(sig_R_maj_maxmaxtrue, points))
except Exception:
pass
try:
axes[1].errorbar(r_sig_mi, sig_mi, yerr=e_sig_mi, fmt='.', marker='.', mew=0, color='b', label='$\sigma_{los}^{min}$')
axes[1].plot(points, map(sig_R_minor_min, points), '--')
axes[1].plot(points, map(sig_R_minor_max, points), '--')
except Exception:
pass
axes[1].set_ylim(0,250)
axes[1].set_xlim(0, 105)
axes[1].grid()
axes[1].legend()
axes[1].set_title('Dispersions')
for photom in all_photometry:
axes[2].plot(r_g_dens, map(photom[-1], r_g_dens), '-', label='{} [M/L={:2.2f}]'.format(photom[0], photom[-2]))
axes[2].set_xlim(0, 150)
axes[2].set_ylim(0, 300)
axes[2].set_title('Photometry')
axes[2].grid()
axes[2].legend()
axes[3].plot(r_g_dens, gas_dens, 'd-')
axes[3].plot(r_g_dens, [He_coeff*(y_interp_(l[0], h_disc_R) + l[1]) for l in zip(r_g_dens, gas_dens)], '*-')
axes[3].plot(r_g_dens, [y_interp_(l, h_disc_R) for l in r_g_dens], '--', label='H2 (R-photom)')
axes[3].set_title('Gas')
axes[3].grid()
axes[3].set_xlim(0, 200)
axes[3].legend()
#change this
plot_2f_vs_1f(ax=axes[4], total_gas_data=zip(r_g_dens, [He_coeff*(y_interp_(l[0], h_disc_R) + l[1]) for l in zip(r_g_dens, gas_dens)])[1:],
epicycl=epicyclicFreq_real,
gas_approx=spl_gas,
sound_vel=sound_vel,
scale=scale,
sigma_max=sig_R_maj_max,
sigma_min=sig_R_maj_min,
star_density_max=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), 5.53, 'R'),
star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), 5.53, 'R'),
data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='R maxdisc')
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)