NGC7217-checkpoint


Эксперименты по восстановлению профилей дисперсий в трех направлениях для NGC7217 (UGC11914)

Сначала всякие настройки и импорты


In [3]:
import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize
from math import *
from IPython.display import HTML
from IPython.display import Image
import os
import pandas as pd

#Размер изображений
import matplotlib.pylab as pylab
pylab.rcParams['figure.figsize'] = 12, 12

#Наклон галактики по данным Сильченко
incl=30.0

# Масштаб пк/секунда из NED
scale=80

Всякие картинки и БД для большего удобства:


In [7]:
# Данные из SDSS DR9
HTML('<iframe src=http://skyserver.sdss3.org/dr9/en/tools/explore/obj.asp?ra=22:07:52.4&dec=+31:21:34 width=1000 height=350></iframe>')


Out[7]:

In [5]:
# Данные из HYPERLEDA
HTML('<iframe src=http://leda.univ-lyon1.fr/ledacat.cgi?o=ngc7217 width=1000 height=350></iframe>')


Out[5]:

In [6]:
# Данные из NED
HTML('<iframe src=http://ned.ipac.caltech.edu/cgi-bin/objsearch?objname=ngc+7217&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[6]:

In [8]:
os.chdir("C:\\science\\2FInstability\\data\\ngc7217")
try:
    from PIL import Image
except:
    import Image

fig, subplot = subplots(1, 2)
subplot[0].imshow(np.asarray(Image.open("ngc7217_JHK.jpg")))
subplot[0].set_title("2MASS image JHK")
subplot[1].imshow(np.asarray(Image.open("ngc7217_SDSS.jpeg")))
subplot[1].set_title("SDSS DR9 whole image")
plt.show()

plt.imshow(np.asarray(Image.open("ugc11914.png")))
plt.title("Image with HI surf. dens.")
plt.show()



In [10]:
#Выводим данные за header-а файла
for line in file("v_stars_ma.dat"):
    if line[0] == '#':
        print(line)


# UGC 11914 = NGC 7217

# Silchenko+2011

# PA=81

# SN15

# corrected for V_sys

# inclination = 30 deg

# data are not corrected for inclination

#  r"       v      dv      sig    dsig

Выведем также для удобства данные из обоих файлов:


In [11]:
ma = pd.read_csv('v_stars_ma.dat', skiprows=8, header=False)
HTML(ma.to_html())


Out[11]:
-95.9 -123.5 10.4 103.8 14.8
0 -72.0 -111.0 6.1 105.8 8.6
1 -62.6 -106.1 4.6 114.2 6.2
2 -55.1 -100.3 4.4 134.8 5.7
3 -49.7 -96.7 3.7 111.1 5.0
4 -46.0 -108.9 3.8 107.5 5.3
5 -42.8 -109.5 3.3 112.0 4.5
6 -39.6 -104.2 3.2 114.4 4.4
7 -36.9 -112.2 2.9 93.3 4.2
8 -34.8 -112.5 2.9 104.0 4.1
9 -32.6 -103.5 2.6 107.1 3.7
10 -30.5 -110.3 2.3 99.9 3.3
11 -28.4 -110.7 2.2 97.8 3.1
12 -26.8 -98.4 2.8 106.2 4.0
13 -25.7 -114.2 2.9 94.1 4.2
14 -24.6 -116.0 2.9 96.8 4.2
15 -23.6 -119.5 3.0 89.8 4.4
16 -22.5 -103.6 3.2 116.7 4.3
17 -21.4 -104.6 2.9 112.5 4.1
18 -20.3 -96.0 2.8 116.1 4.0
19 -19.3 -95.4 2.6 123.9 3.6
20 -18.2 -95.3 2.6 115.4 3.6
21 -17.1 -99.4 2.4 107.9 3.4
22 -16.1 -100.1 2.3 113.3 3.3
23 -15.0 -89.2 2.3 120.9 3.1
24 -13.9 -91.2 2.3 120.5 3.1
25 -12.9 -83.4 2.2 125.7 2.9
26 -11.8 -84.1 2.1 120.8 2.8
27 -10.7 -85.5 1.9 116.0 2.6
28 -9.6 -81.7 1.8 119.6 2.5
29 -8.6 -76.7 1.7 121.2 2.3
30 -7.5 -67.5 1.7 128.8 2.2
31 -6.4 -62.7 1.5 132.5 2.0
32 -5.4 -54.5 1.3 130.4 1.7
33 -4.3 -48.3 1.2 134.4 1.6
34 -3.2 -37.7 1.1 136.4 1.5
35 -2.1 -27.5 1.0 135.0 1.3
36 -1.1 -13.2 0.8 129.3 1.1
37 0.0 3.5 0.8 122.5 1.1
38 1.1 18.5 0.8 125.4 1.0
39 2.1 28.5 0.8 127.2 1.2
40 3.2 40.7 0.9 130.2 1.2
41 4.3 48.5 1.1 132.9 1.5
42 5.4 59.5 1.1 127.0 1.5
43 6.4 65.3 1.1 120.5 1.6
44 7.5 71.4 1.3 120.4 1.8
45 8.6 77.7 1.4 115.2 2.0
46 9.6 83.1 1.4 103.7 2.1
47 10.7 90.0 1.5 102.8 2.2
48 11.8 90.0 1.7 106.5 2.4
49 12.9 97.4 1.8 101.2 2.7
50 13.9 95.9 2.0 100.5 2.9
51 15.0 98.6 2.1 100.8 3.2
52 16.1 97.3 2.2 101.9 3.2
53 17.1 93.7 2.3 106.2 3.4
54 18.2 99.2 2.5 115.2 3.6
55 19.3 107.6 2.1 90.7 3.3
56 20.3 107.7 2.3 100.3 3.4
57 21.4 112.3 2.3 89.6 3.5
58 22.5 113.4 2.2 83.0 3.4
59 23.6 108.2 2.4 86.5 3.7
60 24.6 115.2 2.2 72.2 3.6
61 25.7 104.4 2.7 87.7 4.1
62 27.3 110.8 2.1 84.8 3.3
63 29.4 113.8 2.0 70.8 3.3
64 31.6 113.3 2.2 78.5 3.5
65 33.7 111.7 2.3 81.6 3.6
66 35.9 111.1 2.5 72.3 4.0
67 38.0 115.1 2.6 82.0 4.2
68 40.7 107.3 3.0 101.5 4.4
69 43.9 105.5 2.8 79.6 4.3
70 47.6 103.6 3.1 86.9 4.7
71 51.9 96.2 3.6 96.9 5.2
72 57.6 96.3 4.4 104.4 6.3
73 66.5 106.9 4.3 73.9 6.8
74 73.9 114.1 4.9 84.2 7.5
75 83.6 103.9 6.0 65.8 9.8
76 102.4 83.3 9.9 75.3 15.6

In [13]:
mi = pd.read_csv('v_stars_mi.dat', skiprows=8, header=False)
HTML(mi.to_html())


Out[13]:
-85.4 -40.2 16.4 167.1 20.2
0 -59.4 -19.0 6.2 113.8 8.6
1 -45.8 -11.8 3.9 91.9 5.8
2 -37.8 -10.0 3.2 93.4 4.6
3 -32.6 -4.8 3.3 111.2 4.5
4 -28.9 -8.6 2.9 101.6 4.1
5 -25.6 -8.6 2.5 99.1 3.7
6 -23.0 -5.7 2.7 101.2 3.9
7 -20.9 -8.8 2.7 106.1 3.8
8 -18.7 -7.0 2.4 113.3 3.4
9 -16.6 -6.8 2.1 97.9 3.1
10 -15.0 -15.3 2.5 94.7 3.8
11 -13.9 -11.6 2.6 93.6 3.8
12 -12.9 -2.5 2.3 100.2 3.3
13 -11.8 -2.7 2.1 97.8 3.1
14 -10.7 -1.9 2.1 105.0 3.1
15 -9.6 -4.3 1.9 98.8 2.8
16 -8.6 -3.2 1.9 103.7 2.8
17 -7.5 -2.7 1.8 113.7 2.5
18 -6.4 -2.6 1.7 120.0 2.4
19 -5.4 -1.6 1.6 128.9 2.2
20 -4.3 1.5 1.5 132.9 2.0
21 -3.2 1.1 1.3 132.3 1.8
22 -2.1 -0.7 1.2 134.2 1.6
23 -1.1 1.0 1.0 128.5 1.4
24 0.0 1.5 1.0 126.7 1.3
25 1.1 2.7 1.1 126.9 1.4
26 2.1 2.1 1.2 130.3 1.6
27 3.2 2.4 1.3 128.0 1.8
28 4.3 2.6 1.4 121.2 2.0
29 5.4 5.1 1.6 117.0 2.2
30 6.4 4.8 1.7 111.7 2.4
31 7.5 1.5 1.8 108.6 2.6
32 8.6 5.9 1.9 106.3 2.7
33 9.6 -0.1 2.0 103.5 3.0
34 10.7 3.0 2.3 103.6 3.4
35 11.8 2.0 2.3 96.9 3.5
36 12.9 2.5 2.3 91.9 3.5
37 14.4 5.7 2.0 93.2 3.0
38 16.6 3.2 2.1 98.8 3.1
39 18.7 4.0 2.3 87.7 3.5
40 20.9 0.1 2.8 104.4 4.1
41 23.5 -3.2 2.7 95.0 4.1
42 26.8 0.3 2.9 80.4 4.5
43 30.5 -0.7 2.6 88.2 4.0
44 35.2 -1.7 3.4 102.7 4.9
45 42.5 0.8 4.4 111.4 6.2
46 56.6 8.8 6.0 104.0 8.6
47 80.7 -12.9 17.1 164.0 21.7

Посмотрим теперь на все доступные данные по кривым вращения.


In [14]:
# Данные по звездной кинематике Silchenko+2011 вдоль большей полуоси, не исправленные за наклон 
silch_raw_data = np.loadtxt("v_stars_ma.dat", float)
r_ma, vel_ma, e_vel_ma, sig_ma, e_sig_ma = zip(*silch_raw_data)

# Данные по звездной кинематике Silchenko+2011 вдоль малой полуоси, не исправленные за наклон 
silch_raw_data = np.loadtxt("v_stars_mi.dat", float)
r_mi, vel_mi, e_vel_mi, sig_mi, e_sig_mi = zip(*silch_raw_data)

plt.plot(r_ma, vel_ma, '.-', label="Silchenko+2011, maj")
plt.plot(r_mi, vel_mi, '.-', label="Silchenko+2011, min")
plt.legend()
plt.plot()


Out[14]:
[]

In [18]:
def incline_velocity(v, angle):
    return v / sin(angle * pi / 180)

# Переносит центр в (r0,v0) и перегибает кривую вращения, 
# а также исправляет за наклон если необходимо
def correct_rotation_curve(rdata, vdata, dvdata, r0, v0, incl):
    rdata_tmp = [abs(r-r0) for r in rdata]
    vdata_tmp = [incline_velocity(abs(v-v0), incl) for v in vdata]
    data = zip(rdata_tmp, vdata_tmp, dvdata)
    data.sort()
    return zip(*data)

r_ma_b, vel_ma_b, e_vel_b = correct_rotation_curve(r_ma, vel_ma, e_vel_ma,  0.0, 2.5, incl)
r_mi_b, vel_mi_b, e_vel_mi_b = correct_rotation_curve(r_mi, vel_mi, e_vel_mi,  0.0, 2.5, incl)

plt.plot(r_ma_b, vel_ma_b, 'd', label = 'Silchenko star maj')
plt.errorbar(r_ma_b, vel_ma_b, yerr=e_vel_b, fmt='.', marker='.', mew=0, color='blue')
plt.plot(r_mi_b, vel_mi_b, '.', label = 'Silchenko star min', color='green')
plt.errorbar(r_mi_b, vel_mi_b, yerr=e_vel_mi_b, fmt='.', marker='.', mew=0, color='green')
plt.legend()
plt.plot()


Out[18]:
[]

В дальнейшем используем только данные по звездам по большой полуоси, приблизим их полиномом.


In [19]:
poly_star = poly1d(polyfit(r_ma_b, vel_ma_b, deg=7))

plt.plot(r_ma_b, vel_ma_b, 'x-', color='blue', markersize=6)
test_points = np.arange(0.0, max(r_ma_b), 0.1)
plt.plot(test_points, poly_star(test_points), '-', color='red')
plt.xlabel('$R$'); plt.ylim(0)
plt.ylabel('$V^{maj}_{\phi}(R)$')
plt.show()


Кривая вращения нам нужна для нахождения соотношения $\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) и приближается гладко функцией $f=0.5(1+e^{-R/R_{0}}),$ где $R_{0}$ --- характерный масштаб.

${\bf Примечание:}$ Такое приближение оправдано следующими соображениями. Для равновесного диска верно уравнение, описанное выше. Для твердотельного участка вращения в центральных областях выражение в скобках равно 2, а $\sigma_{\varphi}^{2}/\sigma_{R}^{2}=1$. На плоском участке кривой вращения на периферии диска $\sigma_{\varphi}^{2}/\sigma_{R}^{2}\thickapprox0.5$. Функция $f$ как раз аппроксимирует такое поведение отношения $\sigma_{\varphi}^{2}/\sigma_{R}^{2}$.

Изобразим получившийся профиль $\sigma_{\varphi}^{2}/\sigma_{R}^{2}$, вычисляемый через производную полинома:


In [20]:
def sigPhi_to_sigR_real(R):
        return 0.5 * (1 + R*poly_star.deriv()(R) / poly_star(R))

plt.plot(test_points, [sigPhi_to_sigR_real(R) for R in test_points], 'd-', color='blue')
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.2)
plt.show()


Найдем теперь характерный масштаб $f=0.5(1+e^{-R/R_{0}})$:


In [21]:
def f(R, Ro):
    return 0.5*(1 + np.exp( -R/Ro ))

xdata = test_points
ydata = sigPhi_to_sigR_real(xdata)

from scipy.optimize import curve_fit
popt, pcov = curve_fit(f, xdata, ydata, p0=[1.0])
Ro = popt[0]

plt.plot(xdata, ydata, 'x-')
plt.plot(xdata, [f(p, Ro) for p in xdata], 's')
plt.axhline(y=0.5)
plt.axhline(y=0.0)
plt.title('$R_{0} = %s $' % Ro)
plt.ylim(0, 2)
plt.show()


Теперь знаем значение отношения $\sigma_{\varphi}^{2}/\sigma_{R}^{2}$ в любой точке, заведем соответствующую функцию:


In [22]:
def sigPhi_to_sigR(R):
    return sqrt(f(R, Ro))

Построим графики дисперсий скоростей на луче зрения вдоль большой и малой оси ($\sigma_{los}^{maj}$ и $\sigma_{los}^{min}$):


In [23]:
# Исправляем значения вдоль малой оси на синус угла:    
def correct_min(R):    
    return R / cos(incl * pi / 180) 

r_mi_extend = map(correct_min, r_mi)
    
plt.plot(r_ma, sig_ma, 's-', label='$\sigma_{los}^{maj}$')
plt.errorbar(r_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='.', mew=0, color='blue')
plt.plot(r_mi_extend, sig_mi, 's-', label='$\sigma_{los}^{min}$')
plt.errorbar(r_mi_extend, sig_mi, yerr=e_sig_mi, fmt='.', marker='.', mew=0, color='black')
plt.xlabel('$R$')
plt.ylabel('$\sigma$')
plt.legend()
plt.show()


Перегнем и приблизим полиномами:


In [24]:
bind_curve = lambda p: (abs(p[0]), abs(p[1]), p[2])
sig_maj_data = 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) 

# Добавляем лишние точки чтобы протянуть дальше
num_fake_points = 10; expscale = 200.0
# fake_radii, fake_sig = zip(*[(31.0 + i, 115*exp(- i / expscale )) for i in range(1, num_fake_points+1)])
fake_radii, fake_sig = (),()

poly_sig_maj = poly1d(polyfit(radii_maj + fake_radii, sig_maj_p + fake_sig, deg=9))

sig_min_data = 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) 

poly_sig_min = poly1d(polyfit(radii_min, sig_min_p, deg=9))

points = np.arange(0, max(radii_min), 0.1)
plt.plot(radii_maj, sig_maj_p, 's', label='$\sigma_{los}^{maj}$', color='blue')
plt.errorbar(radii_maj, sig_maj_p, yerr=e_sig_maj_p, fmt='.', marker='.', mew=0, color='blue')
plt.plot(points, poly_sig_maj(points), label = '$\sigma_{los}^{maj} polyfit$', color='blue')
plt.plot(radii_min, sig_min_p, 's', label='$\sigma_{los}^{min}$', color='red')
plt.errorbar(radii_min, sig_min_p, yerr=e_sig_min_p, fmt='.', marker='.', mew=0, color='red')
plt.plot(points, poly_sig_min(points), label = '$\sigma_{los}^{min} polyfit$', color='red')
plt.plot(fake_radii, fake_sig, 'bs', color='green', label='$fake points$')
plt.legend()
plt.ylim(0,150)
# plt.xlim(0,55)
plt.show()


Методика восстановления профилей $\sigma_{R}(R)$, $\sigma_{\varphi}(R)$ и $\sigma_{z}(R)$ следующая. Представим, что $\sigma_{Z}/\sigma_{R} \equiv \alpha \equiv const$. Тогда, зная значения $\sigma_{\varphi}^{2}/\sigma_{R}^{2}=f(R)$ в каждой точке, получаем из уравнений, описанных выше: $$\sigma_{los,maj}^2=\sigma_R^2[f\sin^2i+\alpha^2\cos^2i]$$ $$\sigma_{los,min}^2=\sigma_R^2[\sin^2i+\alpha^2\cos^2i]$$ Представим теперь $\sigma_R(R)=\sigma_{R,0}\times F(R)$, где $F(0)=1$. Значение в квадратных скобках для $\sigma_{los,min}$ равно константе и, следуя предположению, получаем представление для дисперсии вдоль луча зрения для малой оси как $\sigma_{los,min}(R)=\sigma_{min,0}\times F(R)$. Очевидно $\sigma_{min,0} = \sigma_{los,min}(0)$, а значит мы знаем в каждой точке значение $F(R)=\sigma_{los,min(R)}/\sigma_{min,0}$. Описанная выше система уравнений вырождается в следующую: $$\sigma_{los,maj}^2(R)=\frac{\sigma_{R,0}^2\sigma_{los,min}^2(R)[f\sin^2i+\alpha^2\cos^2i]}{\sigma_{min,0}^2}$$ $$\sigma_{min,0}^2=\sigma_{R,0}^2\sin^2i+\sigma_{R,0}^2\alpha^2\cos^2i$$ Сделаем замену: $\sigma_{R,0}^2\sin^2 i \equiv A,\ \sigma_{R,0}^2\cos^2 i\times \alpha^2 \equiv B$. Окончательно, имеем $N+1$ линейное уравнение для $N$ точек, которые можем решить МНК: $$\left\{ \begin{array}{lr} \sigma_{los,maj}^2(R_j)\times \sigma_{min,0}^2 =\sigma_{los,min}^2(R_j)[Af(R_j)+B]\\ \sigma_{min,0}^2=A+B \end{array} \right. $$


In [38]:
#Обрезаем данные по x > r_eff
r_eff = 5.0

#Значение sig_los_min в 0
sig_min_0 = poly_sig_min(0)

#Правая граница
r_max = 60.0

#Количество точек N и сами точки
N = 100
radii_points = np.arange(r_eff, r_max, (r_max-r_eff)/N)

#Как вычислять ошибку 
def residuals(params, xdata, ydata):
            return (ydata - numpy.dot(xdata, params))

#Начальное приближение (А,В)
x0 = [1000, 100]

#Уравнения:  одно для min и N для maj
#Левая часть:
eq_left = np.concatenate( ([sig_min_0**2], 
                              [poly_sig_maj(p)**2 * sig_min_0**2  for p in radii_points]) )
#Правая часть:
eq_right = np.transpose(
            np.array([
                np.concatenate(([1.0],
                                   [poly_sig_min(R)**2 * sigPhi_to_sigR(R)**2 for R in radii_points])),
                np.concatenate(([1.0],
                                   [poly_sig_min(R)**2 for R in radii_points]))]))

# Что будет, если выкинуть уравнение для min:
# eq_left = np.array([poly_sig_maj(p)**2 * sig_min_0**2  for p in radii_points])
# eq_right = np.transpose(
#             np.array([[poly_sig_min(R)**2 * sigPhi_to_sigR(R)**2 for R in radii_points],
#                                    [poly_sig_min(R)**2 for R in radii_points]]))

# МНК для получившихся уравнений:
solution = scipy.optimize.leastsq(residuals, x0, args=(eq_right, eq_left))[0]
A, B = solution[0], solution[1]

chi2 = sum(power(residuals(solution, eq_right, eq_left), 2))/N
        
print 'Solution: A = %s, B = %s' % (A, B)
print 'Chi^2:', chi2

#Подставить в уравнение в какой-нибудь точке:
# rrr = 25.0
# print (poly_sig_maj(rrr)**2) * sig_min_0**2 - (poly_sig_min(rrr)**2) * (A*f(rrr, Ro)+B)


Solution: A = 7329.88300374, B = 11606.2840983
Chi^2: 5.60375340903e+14

Теперь восстановим исходные неизвестные - $\alpha$ и $\sigma_{R, 0}$:


In [28]:
sig_R_0 = round( sqrt(A) / sin(incl*pi/180), 3)
alpha = round( sqrt(B)/ (cos(incl*pi/180) * sig_R_0), 3)

# sig_R_0 = 160
# alpha = 0.55

print "sig_R_0 = %s, alpha = %s" % (sig_R_0, alpha)


sig_R_0 = 171.229, alpha = 0.727

Построим полученные профили дисперсий скоростей:


In [30]:
def sigR_exp(R):
    return sig_R_0*poly_sig_min(R)/sig_min_0

def sigZ_exp(R):
    return alpha * sigR_exp(R)

def sigPhi_exp(R):
    return sigPhi_to_sigR(R) * sigR_exp(R)


plt.plot(points, [sigR_exp(R) for R in points], '-', color='red', label='$\sigma_{R, exp}$')
plt.plot(points, [sigPhi_exp(