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

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


In [1]:
import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize
from numpy import poly1d, polyfit, power
import math
from math import *
from IPython.display import HTML
from IPython.display import Image
import os
import pandas as pd
import PIL as pil
import heapq
%matplotlib inline

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

#Наклон галактики по данным Zasov08, Noordermeer07
incl = 56

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

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


In [2]:
# Данные из SDSS DR9
HTML('<iframe src=http://skyserver.sdss3.org/dr9/en/tools/explore/obj.asp?ra=12:10:33.6&dec=+30:24:06 width=1000 height=350></iframe>')


Out[2]:

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


Out[3]:

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

In [5]:
os.chdir("C:\\science\\2FInstability\\data\\ngc4150")

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


# UGC 7165 = NGC 4150

# Zasov+08

# PA = 146 deg

# not corrected for V_sys

# inclination = 56 deg (adopted by Noordermeer+2007)

# data ARE NOT corrected for inclination

# SSBS_FUL: Cross correlation with star hd125560.ln

# Original spectra:n4150ma.ln

#  r"      v     dv  sig(TD)  dsig

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


In [6]:
# pd.set_option('display.notebook_repr_html', True)
# ma = pd.read_csv('v_stars_maZ.dat', skiprows=9, sep='  ', header=False)
# display(ma)

In [7]:
# mi = pd.read_csv('v_stars_mi.dat', skiprows=9, sep='   ', header=False)
# display(mi)

Посмотрим теперь на все доступные даные по кривым вращения. Т.к. некоторые из них скорректированы, а некоторые нет, разобьем этот этап на несколько шагов.


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

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

plt.plot(r_ma, vel_ma, '.-', label="Zasov 08, maj")
plt.plot(r_mi, vel_mi, '.-', label="Zasov 08, min")
plt.legend()
plt.plot()


Out[8]:
[]

In [9]:
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, 242., incl)
r_mi_b, vel_mi_b, e_vel_mi_b = correct_rotation_curve(r_mi, vel_mi, e_vel_mi,  0.0, 242., incl)

plt.plot(r_ma_b, vel_ma_b, 'd', label = 'Zasov 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 = 'Zasov 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[9]:
[]

Как видим, сигнал вдоль малой оси существенно не нулевой. Это означает, что смещена малая полуось и что по-хорошему наш подход тут не применим. Однако, как будет видно ниже, дисперсии скоростей относительно нормальные.


In [10]:
poly_star = poly1d(polyfit(r_ma_b, vel_ma_b, deg=5))

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



In [11]:
tex_imgs_dir = "C:\\Users\\root\\Dropbox\\RotationCurves\\PhD\\paper1\\text\\imgs"
try: 
    os.chdir(tex_imgs_dir)
except:
    tex_imgs_dir = "C:\\Users\\Alex March\\Dropbox\\RotationCurves\\PhD\\paper1\\text\\imgs"
    
os.chdir(tex_imgs_dir)

np.save('n4150_maj_rot', zip(r_ma_b, vel_ma_b, e_vel_b))
np.save('n4150_rot_poly', zip(test_points, poly_star(test_points)))

os.chdir("C:\\science\\2FInstability\\data\\ngc4150")

Кривая вращения нам нужна для нахождения соотношения $\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 [12]:
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, 5)
plt.show()


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


In [13]:
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, 5)
plt.show()


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


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

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


In [15]:
# Исправляем значения вдоль малой оси на синус угла:    
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.ylim(-50, 100)
plt.show()


Данные конечно же совершенно не вычещенные:


In [16]:
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) 

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) 

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='o', marker='.', 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='o', marker='.', color='red')
plt.legend()
plt.ylim(-5.,100)
plt.xlim(0,85)
plt.axhline(y=0., c='black')
plt.show()


Уберем лишние точки (нулевые, отступающие), отрежем центральную часть (до 8 секунд) и


In [17]:
# Граница. по которой обрезаем
cutted = 8.

os.chdir(tex_imgs_dir)

s_maj_data = zip(radii_maj, sig_maj_p, e_sig_maj_p)
#Обрезаем данные по 35" и убираем нулевые значения
s_maj_data = filter(lambda x: x[1] > 30 and x[0] < 35, s_maj_data)
#Убираем два отскока:
s_maj_data = filter(lambda x: x[0] < 20 or (x[0] > 20 and x[1] < 52), s_maj_data)
np.save('n4150_maj', s_maj_data)
s_maj_data = filter(lambda l: l[0] > cutted, s_maj_data)    
    
s_min_data = zip(radii_min, sig_min_p, e_sig_min_p)
#Обрезаем данные по 50" и убираем нулевые значения
s_min_data = filter(lambda x: x[0] < 50, s_min_data)
np.save('n4150_min', s_min_data)
s_min_data = filter(lambda l: l[0] > cutted, s_min_data)

os.chdir("C:\\science\\2FInstability\\data\\ngc4150")

plt.plot(radii_maj, sig_maj_p, 's', label='$\sigma_{los}^{maj} excl$', color='g')
plt.errorbar(radii_maj, sig_maj_p, yerr=e_sig_maj_p, fmt='o', marker='.', color='g')
radii_maj, sig_maj_p, e_sig_maj_p = zip(*s_maj_data)
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='o', marker='.', color='blue')

plt.plot(radii_min, sig_min_p, 's', label='$\sigma_{los}^{min} excl$', color='pink')
plt.errorbar(radii_min, sig_min_p, yerr=e_sig_min_p, fmt='o', marker='.', color='pink')
radii_min, sig_min_p, e_sig_min_p = zip(*s_min_data)
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='o', marker='.', color='red')

plt.axvline(x=cutted, color='black')
plt.ylim(0., 100.)
plt.legend()
plt.show()


Приблизим сплайнами ($k=3$), а весовую функцию возьмем $w(err)=\frac{1}{1+{err}^2}$:


In [18]:
def w(arr):
    return map(lambda l: 1/(1. + l**2), arr)

import scipy.interpolate as inter
points = np.arange(cutted, max(radii_min), 0.1)

poly_sig_maj = inter.UnivariateSpline (radii_maj[1:], sig_maj_p[1:], k=3, s=10000., w=w(e_sig_maj_p[1:]))
poly_sig_min = inter.UnivariateSpline (radii_min[1:-1:1], sig_min_p[1:-1:1], k=3, s=10000., w=w(e_sig_min_p[1:-1: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='o', marker='.', color='blue')
plt.plot(points, poly_sig_maj(points), label = '$\sigma_{los}^{maj}\, splinefit$', 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='o', marker='.', color='red')
plt.plot(points, poly_sig_min(points), label = '$\sigma_{los}^{min}\, splinefit$', color='red')
plt.ylim(0., 80.)
plt.legend()
plt.show()



In [19]:
os.chdir(tex_imgs_dir)

# np.save('n4150_maj', zip(radii_maj, sig_maj_p, e_sig_maj_p))
# np.save('n4150_min', zip(radii_min, sig_min_p, e_sig_min_p))
np.save('n4150_min_poly', zip(points, poly_sig_min(points)))
np.save('n4150_maj_poly', zip(points, poly_sig_maj(points)))

os.chdir("C:\\science\\2FInstability\\data\\ngc4150")

Посчитаем величину невязок для полученного приближения:


In [57]:
sqerr_maj = sum(power([poly_sig_maj(p[0]) - p[1] for p in sig_maj_data], 2))
sqerr_min = sum(power([poly_sig_min(p[0]) - p[1] for p in sig_min_data], 2))

chi2_maj = sqerr_maj / len(sig_maj_p)
chi2_min = sqerr_min / len(sig_min_p)

print "Poly chi^2 for maj full = %s, mean = %s" % (sqerr_maj, chi2_maj)
print "Poly chi^2 for min full = %s, mean = %s" % (sqerr_min, chi2_min)


Poly chi^2 for maj full = 17840.5701601, mean = 557.517817504
Poly chi^2 for min full = 14191.6395828, mean = 308.513903974

Методика восстановления профилей $\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 [58]:
#Как вычислять ошибку 
def residuals(params, xdata, ydata):
            return (ydata - np.dot(xdata, params))

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

In [59]:
#Восстановление первоначальных неизвестных
def physical_unknowns(A, B):
    sig_R_0 = round( sqrt(A) / sin(incl*pi/180), 3)
    alpha = round( sqrt(B)/ (cos(incl*pi/180) * sig_R_0), 3)
    return (sig_R_0, alpha)

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

Ns = [20, 50, 75, 100]

right_max = 60.

#Шаг сетки по расстоянию
dx = 2.

#Минимальный размер отрезка, на котором ищем ответ
min_size = 20.

lefts = np.arange(0, right_max, dx)

best = []

for n in Ns:
    result_for_N = []
    for left in lefts:
        rights = np.arange(left+min_size, right_max, dx)
        for right in rights:
            r_points = np.arange(left, right, (right-left)/n)
            eq_left = np.concatenate( ([sig_min_0**2], 
                              [poly_sig_maj(p)**2 * sig_min_0**2  for p in r_points]) )
            eq_right = np.transpose(
                        np.array([
                            np.concatenate(([1.0],
                                               [poly_sig_min(R)**2 * sigPhi_to_sigR(R)**2 for R in r_points])),
                            np.concatenate(([1.0],
                                               [poly_sig_min(R)**2 for R in r_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
            if A > 0 and B > 0:
                result_for_N.append([chi2, left, right-left, A, B])
    if result_for_N.__len__() > 0:
        result_for_N.sort()
        Z,X,Y,A,B = zip(*result_for_N)
        print "For N=%s top-10 best results:" % n
        for ind in range(0, min(10, Z.__len__())):
            prin = (ind+1, Z[ind], X[ind], X[ind]+Y[ind]) + physical_unknowns(A[ind], B[ind])
            print "\t%s place: chi2 = %s in range [%s:%s]; sig_R_0 = %s alpha = %s" % prin
        if best.__len__() == 0 or best[0] > Z[0]:
            best = [Z[0], X[0], Y[0], n]

if best.__len__() > 0:
    print "Best of the best: N=%s chi2 = %s on range[%s:%s]" % (best[3], best[0], best[1], best[2]+best[1])
else:
    print "No even one solution."


For N=20 top-10 best results:
	1 place: chi2 = 31742039451.7 in range [20.0:40.0]; sig_R_0 = 42.376 alpha = 2.709
	2 place: chi2 = 51924968726.0 in range [18.0:38.0]; sig_R_0 = 53.891 alpha = 2.016
	3 place: chi2 = 54849990714.0 in range [18.0:40.0]; sig_R_0 = 58.157 alpha = 1.814
	4 place: chi2 = 86878127070.9 in range [20.0:42.0]; sig_R_0 = 57.937 alpha = 1.814
	5 place: chi2 = 91925835176.4 in range [18.0:42.0]; sig_R_0 = 66.472 alpha = 1.481
	6 place: chi2 = 97809365055.0 in range [22.0:42.0]; sig_R_0 = 56.582 alpha = 1.874
	7 place: chi2 = 123106863707.0 in range [16.0:40.0]; sig_R_0 = 74.666 alpha = 1.219
	8 place: chi2 = 130642742902.0 in range [16.0:38.0]; sig_R_0 = 73.685 alpha = 1.25
	9 place: chi2 = 134920913272.0 in range [16.0:36.0]; sig_R_0 = 75.702 alpha = 1.185
	10 place: chi2 = 138276221375.0 in range [16.0:42.0]; sig_R_0 = 78.724 alpha = 1.093
For N=50 top-10 best results:
	1 place: chi2 = 36578096283.5 in range [20.0:40.0]; sig_R_0 = 43.083 alpha = 2.654
	2 place: chi2 = 45330543131.4 in range [18.0:38.0]; sig_R_0 = 51.422 alpha = 2.142
	3 place: chi2 = 53398890745.9 in range [18.0:40.0]; sig_R_0 = 57.069 alpha = 1.861
	4 place: chi2 = 105933337862.0 in range [18.0:42.0]; sig_R_0 = 66.897 alpha = 1.464
	5 place: chi2 = 106664774480.0 in range [20.0:42.0]; sig_R_0 = 60.292 alpha = 1.713
	6 place: chi2 = 109589534171.0 in range [16.0:40.0]; sig_R_0 = 72.649 alpha = 1.281
	7 place: chi2 = 114333833727.0 in range [16.0:38.0]; sig_R_0 = 70.891 alpha = 1.34
	8 place: chi2 = 118114541671.0 in range [22.0:42.0]; sig_R_0 = 60.902 alpha = 1.688
	9 place: chi2 = 121297952540.0 in range [16.0:36.0]; sig_R_0 = 72.529 alpha = 1.285
	10 place: chi2 = 139000563708.0 in range [16.0:42.0]; sig_R_0 = 77.776 alpha = 1.119
For N=75 top-10 best results:
	1 place: chi2 = 37827279914.4 in range [20.0:40.0]; sig_R_0 = 43.287 alpha = 2.639
	2 place: chi2 = 43928177300.5 in range [18.0:38.0]; sig_R_0 = 50.871 alpha = 2.172
	3 place: chi2 = 53307642174.9 in range [18.0:40.0]; sig_R_0 = 56.852 alpha = 1.871
	4 place: chi2 = 106809943234.0 in range [16.0:40.0]; sig_R_0 = 72.205 alpha = 1.295
	5 place: chi2 = 109535631955.0 in range [18.0:42.0]; sig_R_0 = 67.031 alpha = 1.459
	6 place: chi2 = 110696473126.0 in range [16.0:38.0]; sig_R_0 = 70.256 alpha = 1.361
	7 place: chi2 = 111389295659.0 in range [20.0:42.0]; sig_R_0 = 60.858 alpha = 1.689
	8 place: chi2 = 118140646139.0 in range [16.0:36.0]; sig_R_0 = 71.799 alpha = 1.308
	9 place: chi2 = 122799340123.0 in range [22.0:42.0]; sig_R_0 = 61.882 alpha = 1.649
	10 place: chi2 = 139767804190.0 in range [16.0:42.0]; sig_R_0 = 77.587 alpha = 1.124
For N=100 top-10 best results:
	1 place: chi2 = 38475425310.5 in range [20.0:40.0]; sig_R_0 = 43.395 alpha = 2.631
	2 place: chi2 = 43237896224.5 in range [18.0:38.0]; sig_R_0 = 50.595 alpha = 2.187
	3 place: chi2 = 53295392374.5 in range [18.0:40.0]; sig_R_0 = 56.747 alpha = 1.875
	4 place: chi2 = 105456547052.0 in range [16.0:40.0]; sig_R_0 = 71.983 alpha = 1.302
	5 place: chi2 = 108881344529.0 in range [16.0:38.0]; sig_R_0 = 69.937 alpha = 1.371
	6 place: chi2 = 111403323010.0 in range [18.0:42.0]; sig_R_0 = 67.104 alpha = 1.456
	7 place: chi2 = 113793667232.0 in range [20.0:42.0]; sig_R_0 = 61.147 alpha = 1.678
	8 place: chi2 = 116548361340.0 in range [16.0:36.0]; sig_R_0 = 71.431 alpha = 1.32
	9 place: chi2 = 125161995732.0 in range [22.0:42.0]; sig_R_0 = 62.375 alpha = 1.629
	10 place: chi2 = 140238683802.0 in range [16.0:42.0]; sig_R_0 = 77.495 alpha = 1.127
Best of the best: N=20 chi2 = 31742039451.7 on range[20.0:40.0]

Заведе функции для построения профилей дисперсий скоростей:


In [60]:
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)

И для профилей $\sigma_{los}^{maj}$ и $\sigma_{los}^{min}$. Связь профилей описывается следующими уравнениями: $$\sigma_{los,maj}^2=\sigma_{\varphi}^2\sin^2i+\sigma_Z^2\cos^2i$$ $$\sigma_{los,min}^2=\sigma_R^2\sin^2i+\sigma_Z^2\cos^2i$$


In [61]:
def sig_maj_exp(R):
    return sqrt(sigPhi_exp(R)**2 * sin(incl*pi/180)**2 + sigZ_exp(R)**2 * cos(incl*pi/180)**2)

def sig_min_exp(R):
    return sqrt(sigR_exp(R)**2 * sin(incl*pi/180)**2 + sigZ_exp(R)**2 * cos(incl*pi/180)**2)


cos_i, sin_i = cos(incl * pi / 180), sin(incl * pi / 180)

def sig_maj_exp(R):
    tmp = sigPhi_to_sigR_real(R) * sin_i**2 + alpha**2 * cos_i**2
    if tmp > 0:
        return sig_R_0*poly_sig_min(R)/sig_min_0 * sqrt(sigPhi_to_sigR_real(R) * sin_i**2 + alpha**2 * cos_i**2)
    else:
        return -100
#     return sig_R_0*poly_sig_min(R)/sig_min_0 * sqrt(sigPhi_to_sigR_real(R) * sin_i**2 + alpha**2 * cos_i**2)
#     return sig_R_0*spl_min(R)/sig_min_0 * sqrt(sigPhi_to_sigR(R)**2 * sin_i**2 + alpha**2 * cos_i**2)
#     return sqrt(sigPhi_exp(R)**2 * sin(incl*pi/180)**2 + sigZ_exp(R)**2 * cos(incl*pi/180)**2)

def sig_min_exp(R):
    return sig_R_0*poly_sig_min(R)/sig_min_0 * sqrt(sin_i**2 + alpha**2 * cos_i**2)
#     return sqrt(sigR_exp(R)**2 * sin(incl*pi/180)**2 + sigZ_exp(R)**2 * cos(incl*pi/180)**2)

Заметим, что можно было не решать систему МНК, а честно разрешить систему из двух уравнений $$\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. $$ относительно $A$ и $B$ для почти любого $R_j$ (а лучше даже относительно начальных неизвестных - $\sigma_{R,0}$ и $\alpha$). Решение: $$\sigma_{R,0}^2 = \frac{\sigma_{min,0}^2}{\sin^2 i}\times\frac{1}{f(R)-1}\times(P(R)-1)$$ $$\alpha^2 = \tan^2 i\frac{f(R) - P(R)}{P(R)-1},$$ $$P(R)=\frac{\sigma_{los,maj}^2(R)}{\sigma_{los, min}^2(R)}$$ Имеет смысл также искать не $\alpha$, а $\alpha\cdot\sigma_{R,0}=\sigma_{Z,0}$.


In [62]:
cos_i, sin_i = cos(incl * pi / 180), sin(incl * pi / 180)

def P(R):
    """Отношение maj к min, как описано выше"""
    return (poly_sig_maj(R)/poly_sig_min(R))**2

def direct_solve_A(R):
    """Аналитически находим значение sig_R_0 для уравнения в точке R"""
    res = sig_min_0**2 * (P(R) - 1) / (sin_i**2 * (sigPhi_to_sigR(R)**2 - 1))
    return sqrt(res) if res > 0 else 0

def direct_solve_B(R):
    """Аналитически находим значение alpha для уравнения в точке R"""
    res = (sigPhi_to_sigR(R)**2 - P(R))/(P(R) - 1) * (sin_i/cos_i)**2
    return sqrt(res) if res > 0 else 0

def direct_find_sig_R_0(R):
    return direct_solve_A(R)

def direct_find_sig_Z_0(R):
    return direct_solve_A(R) * direct_solve_B(R)

Найдем значения $\sigma_{R,0}$ и $\sigma_{Z,0}$ для всех точек на большой оси:


In [63]:
p_r = radii_maj

direct_sigR0 = map(direct_find_sig_R_0, p_r)
direct_sigZ0 = map(direct_find_sig_Z_0, p_r)

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(p_r, direct_sigR0, 'o', color='r', label='$\sigma_{R,0}$')
plt.plot(p_r, direct_sigZ0, 'o', color='k', label='$\sigma_{Z,0}$')

plt.legend()
plt.ylim(-10, 200)
plt.show()


Как видно - хорошо восстанавливаются значения $\sigma_{R,0}$ - лежат относительно одной прямой, найдем их среднее:


In [64]:
#Обрежем по 12
q=12.
ind_q = p_r.index(filter(lambda l: l > q, p_r)[0])

poly_q = poly1d(polyfit(p_r[ind_q:], direct_sigR0[ind_q:], deg=0))
print "sig_R poly mean = %s" % poly_q[0]


sig_R poly mean = 131.394205838

In [65]:
sig_R_0 = poly_q[0]

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(p_r, direct_sigR0, 'o', color='r', label='$\sigma_{R,0}$')
plt.axhline(y=sig_R_0)

#Строим полученный на основе среднего профиль sig_R
plt.plot(points, poly_sig_min(points)*sig_R_0/sig_min_0, label = '$\sigma_R$', color='m')

plt.legend()
plt.ylim(-10, 200)
plt.show()


Как мы видим, значения профиля дисперсии $\sigma_R(R)$ восстанавливаются довольно надежно, однако оказываются по-видимому больше реальных и поэтому не получается восстановить профиль в вертикальном направлении. Попробуем оценить, насколько этот вклад оказывается переоценен в данном случае. Известно, что отношение $\sigma_Z/\sigma_R$ не может быть меньше некоего порогового значения, в противном случае галактика будет неустойчива к осесимметричным изгибным возмущениям плотности. Многие авторы оценивали эту величину, в том числе Засов, однако последния статья Сотниковой и Радионова "Bending instability in galactic discs. Advocacy of the linear theory" (2013, http://arxiv.org/abs/1306.5975) продемонстрировала, что искомое попроговое значение близко к таковому, полученному из линейной теории Тумре в 1966 год и равно примерно 0.3. Чем это ценно для нас? Если мы примем, что $\frac{\sigma_Z}{\sigma_R} \gtrsim 0.3,$ то, исходя из уравнения $\sigma_{los,min}^2=\sigma_R^2\sin^2i+\sigma_Z^2\cos^2i$ можем получить оценку сверху на значения радиальной дисперсии: $$\frac{\sigma_{los,min}}{\sqrt{\sin^2i+0.09\cos^2i}} \gtrsim \sigma_R$$ Также аналогичную оценку, только чуть более сложную, можно сделать и для данных вдоль большой оси.


In [66]:
def sig_R_upper_lim(R, alpha):
    """Оценка сверху на sigR(R)"""
    return poly_sig_min(R)/sqrt(sin_i**2 + alpha**2 * cos_i**2)

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(p_r, direct_sigR0, 'o', color='r', label='$\sigma_{R,0}$')
plt.axhline(y=sig_R_0)
plt.plot(points, poly_sig_min(points)*sig_R_0/sig_min_0, label = '$\sigma_R$', color='m')
plt.plot(points, [sig_R_upper_lim(R, 0.3) for R in points], label = '$\sigma_R^{up}$', color='g')

plt.legend()
plt.ylim(-10, 200)
plt.show()


Как видим, значения действительно оказались сильно переоценены.

Теперь то, с чего надо было начинать - построим картинки для разных значений $\alpha$ и $\sigma_{R,0}$. Для того, чтобы найти где минимум, попробуем построить просто двумерные карты $\chi^2$ для разных $\sigma_{R,0}$ $\alpha$: (это очень долго, так что пересчитывать в крайнем случае)


In [67]:
alphas = np.arange(0.1, 1.2, 0.03)
sigmas = np.arange(10.0, 300, 3.)

def calc_chi2_normal(obs, obserr, predicted):
    return sum([(o-p)**2/err**2 for (o,p,err) in zip(obs, predicted, obserr)])/len(obs)

def compute_chi2_maps(alphas=(), sigmas=()):
    '''Вычисляем все изображения, чтобы потом только настройки менять'''
    image_min = np.random.uniform(size=(len(sigmas), len(alphas)))
    image_maj = np.random.uniform(size=(len(sigmas), len(alphas)))
    image = np.random.uniform(size=(len(sigmas), len(alphas)))
    for i,si in enumerate(sigmas):
        for j,al in enumerate(alphas):
            global alpha, sig_R_0
            alpha = al
            sig_R_0 = si
#             sqerr_maj = sum(power([sig_maj_exp(p[0]) - p[1] for p in sig_maj_data], 2))/len(sig_maj_data)
#             sqerr_min = sum(power([sig_min_exp(p[0]) - p[1] for p in sig_min_data], 2))/len(sig_min_data)
#             sqerr_maj = sum([(sig_maj_exp(p[0]) - p[1])**2/p[2]**2 for p in sig_maj_data])
#             sqerr_min = sum([(sig_min_exp(p[0]) - p[1])**2/p[2]**2 for p in sig_min_data])
            sqerr_maj = calc_chi2_normal(sig_maj_p, e_sig_maj_p, [sig_maj_exp(r) for r in radii_maj])
            sqerr_min = calc_chi2_normal(sig_min_p, e_sig_min_p, [sig_min_exp(r) for r in radii_min])
            sqerr_sum = 0.5*sqerr_maj+0.5*sqerr_min
            image[i][j] = sqerr_sum
            image_maj[i][j] = sqerr_maj
            image_min[i][j] = sqerr_min
    return image, image_maj, image_min

pics_path = '.cutted\\pics\\'
if not os.path.exists(pics_path):
    os.makedirs(pics_path)
if os.path.isfile(pics_path + 'chi2_map.npy'):
    image = np.load(pics_path + "chi2_map.npy")
    image_maj = np.load(pics_path + "chi2_map_maj.npy")
    image_min = np.load(pics_path + "chi2_map_min.npy")
else:
    image, image_maj, image_min = compute_chi2_maps(alphas=alphas, sigmas=sigmas)
    np.save(pics_path + 'chi2_map', image)
    np.save(pics_path + 'chi2_map_maj', image_maj)
    np.save(pics_path + 'chi2_map_min', image_min)

In [68]:
from mpl_toolkits.axes_grid1 import make_axes_locatable

def plot_chi2_map(image, ax, log_scale=False, title='$\chi^2$', is_contour=False, vmax=0.):
    '''Рисуем получившиеся карты.
    Colormaps: http://wiki.scipy.org/Cookbook/Matplotlib/Show_colormaps'''
    if image is not None:
        if log_scale:
            image_log = np.apply_along_axis(np.log, 1, image)
            vmax = image_log.max()
        else:
            image_log = image
        if is_contour:
            norm = plt.cm.colors.Normalize(vmax=image.max(), vmin=-image.max())
            cmap = plt.cm.PRGn
            levels = np.concatenate([np.array([image_log.min()*1.1,]), np.linspace(start=image_log.min(), stop=vmax, num=10)])
            levels = sorted(levels)
            cset=ax.contour(image_log, levels, hold='on', colors = 'k', origin='lower', 
                            extent=[alphas[0],alphas[-1],sigmas[0],sigmas[-1]])
            ax.clabel(cset, inline=1, fontsize=10, fmt='%1.1f',)
        im = ax.imshow(image_log, cmap='jet', vmin=image_log.min(), vmax=vmax, interpolation='spline16', 
                   origin="lower", extent=[alphas[0], alphas[-1],sigmas[0],sigmas[-1]], aspect="auto")
        divider = make_axes_locatable(ax)
        cax = divider.append_axes("right", size="5%", pad=0.05)
        plt.colorbar(im, cax=cax)
        min_sigma = sigmas[int(np.where(image == image.min())[0])]        
        ax.set_title(title + '$,\ \sigma(min)=%s$' % min_sigma, size=20.)
        ax.set_ylabel('$\sigma_{R,0}$', size=20.)
        ax.set_xlabel(r'$\alpha$', size=20.)
        ax.grid(True)
 
fig, axes = plt.subplots(nrows=3, ncols=1, sharex=False, sharey=True, figsize=[16,16])
plot_chi2_map(image, axes[0], log_scale=False, title='$\chi^2 = (\chi^2_{maj} + \chi^2_{min})/2$', is_contour=False, vmax=20.)
plot_chi2_map(image_maj, axes[1], log_scale=False, title='$\chi^2_{maj}$', is_contour=False, vmax=20.)
plot_chi2_map(image_min, axes[2], log_scale=False, title='$\chi^2_{min}$', is_contour=False, vmax=20.)
plt.show()