Галактика из диплома.
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 = 38.0 #mean value, in Zasov 36 deg, in LEDA 49deg, CALIFA 42deg TODO: check, хотя разброс и не влияет на результат сильно
distance = 67.4 #Mpc, Noordermeer
scale = 0.327 #kpc/arcsec according to Noordermeer distance
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 [18]:
%run ../../utils/rotvelutils #import util functions for rotcurv correction
In [19]:
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 [20]:
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 [21]:
star_approx = poly1d(polyfit(r_ma_b, vel_ma_b, deg=3))
In [22]:
#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 [23]:
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 [24]:
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 [25]:
def w(arr):
return map(lambda l: 1/(1. + l**2), arr)
for _ in range(1, 10, 3):
spl_maj = inter.UnivariateSpline(r_sig_ma[_::1]+ tuple(np.linspace(53, 60, 5)), sig_ma[_::1]+ tuple(np.linspace(55, 50, 5)), k=3, s=10000.,
w=w(e_sig_maj_p[_::1]+(5.,5,5,5,5,)))
plt.plot(points, spl_maj(points), '-')
plt.errorbar(r_sig_ma, sig_ma, yerr=e_sig_maj_p, fmt='.', marker='.', mew=0, color='blue', label='$\sigma_{los}^{maj}$')
plt.legend()
plt.ylim(0,250)
plt.xlim(0,70);
In [26]:
@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 [27]:
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 [28]:
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 [29]:
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 [30]:
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);