In [1]:
import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize
from math import *
#Размер изображений
import matplotlib.pylab as pylab
pylab.rcParams['figure.figsize'] = 12, 8
#Наклон
incl = 15
Для начала посмотрим на галактику (моделируемый наклон $15^{\circ}$) плашмя и с ребра:
In [2]:
import Image
fig, subplot = subplots(1, 2)
subplot[0].imshow(np.asarray(Image.open("C:\\science\\2FInstability\\data\\model\\1\\rot15\\xy_plot0.png")))
subplot[1].imshow(np.asarray(Image.open("C:\\science\\2FInstability\\data\\model\\1\\rot15\\zy_plot0.png")))
plt.show()
Находим файл с данными:
In [3]:
import os
os.chdir("C:\\science\\2FInstability\\data\\model\\1")
Строим карту скоростей:
In [4]:
velocity_data_points = np.loadtxt("C:\\science\\2FInstability\\data\\model\\1\\rot15\\mvz_05.txt", float)
vel_values = np.array([p[2] for p in velocity_data_points])
# Пересчет координат [-10,10] -> [0,40] и обратно
def coord(x):
return int((x + 10)*2)
def reverse_coord(x):
return x/2.0 - 10.0
image = np.random.uniform(size=(40, 40))
for point in velocity_data_points:
image[coord(point[0])][coord(point[1])] = point[2]
plt.imshow(image, cmap="winter", vmin=vel_values.min(), vmax=vel_values.max(), interpolation='none', extent=[-10,10,-10,10])
plt.colorbar()
plt.show()
Вырезаем кривую вращения $V_{\phi}(R)$ вдоль большой ($maj$) и малой ($min$) оси:
In [5]:
def is_min_axis(map_point):
return coord(map_point[0]) == 20
radii = np.arange(-10.0, 10.0, 0.5)
vel_min = []
for min_axis in filter(is_min_axis, velocity_data_points):
vel_min.append(min_axis[2])
plt.plot(radii, vel_min)
plt.ylim(-0.15, 0.15)
plt.xlabel('$R$')
plt.ylabel('$V^{min}_{\phi}(R)$')
plt.show()
In [6]:
def is_maj_axis(map_point):
return coord(map_point[1]) == 20
vel_maj = []
for maj_axis in filter(is_maj_axis, velocity_data_points):
vel_maj.append(maj_axis[2])
plt.plot(radii, vel_maj)
plt.xlabel('$R$')
plt.ylabel('$V^{maj}_{\phi}(R)$')
plt.show()
Перегнутая кривая вращения:
In [7]:
data_together = zip(radii, vel_maj)
bind_curve = lambda p: (abs(p[0]), abs(p[1]))
binded_data = map(bind_curve, data_together)
binded_data.sort()
radii_b, vel_maj_b = zip(*binded_data)
plt.plot(radii_b, vel_maj_b, 'x-')
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}$.
Найдем характерный масштаб $R_{0}$. Для этого сначала приблизим $V_{\phi}^{maj}$ полиномом.
In [8]:
poly_star = poly1d(polyfit(radii_b, vel_maj_b, deg=9))
plt.plot(radii_b, vel_maj_b, 'x-', color='blue', markersize=6)
test_points = np.arange(0.0, max(radii_b), 0.1)
plt.plot(test_points, poly_star(test_points), '-', color='red')
plt.xlabel('$R$')
plt.ylabel('$V^{maj}_{\phi}(R)$')
plt.show()
Изобразим получившийся профиль $\sigma_{\varphi}^{2}/\sigma_{R}^{2}$, вычисляемый через производную полинома:
In [9]:
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='green')
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.show()
Найдем теперь характерный масштаб $f=0.5(1+e^{-R/R_{0}})$:
In [10]:
def f(R, Ro):
return 0.5*(1 + np.exp( -R/Ro ))
# Обрежем на правом краю по 8.0 и выкидываем первую точку
xdata = np.array(filter(lambda x: x < 8.0 and x > 0.001, 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.show()
Теперь знаем значение отношения $\sigma_{\varphi}^{2}/\sigma_{R}^{2}$ в любой точке, заведем соответствующую функцию:
In [11]:
def sigPhi_to_sigR(R):
return sqrt(f(R, Ro))
Теперь построим карту дисперсий скоростей:
In [12]:
dispersion_data_points = np.loadtxt("C:\\science\\2FInstability\\data\\model\\1\\rot15\\svz_05.txt", float)
disp_values = np.array([p[2] for p in dispersion_data_points])
disp_image = np.random.uniform(size=(40, 40))
for point in dispersion_data_points:
disp_image[coord(point[0])][coord(point[1])] = point[2]
plt.imshow(disp_image, cmap="winter", vmin=disp_values.min(), vmax=disp_values.max(), interpolation='none', extent=[-10,10,-10,10])
plt.colorbar()
plt.show()
Как и ранее, делаем разрезы вдоль большой и малой оси чтобы получить дисперсии скоростей на луче зрения ($\sigma_{los}^{maj}$ и $\sigma_{los}^{min}$):
In [13]:
sig_los_maj = []
for maj_axis in filter(is_maj_axis, dispersion_data_points):
sig_los_maj.append(maj_axis[2])
sig_los_min = []
for min_axis in filter(is_min_axis, dispersion_data_points):
sig_los_min.append(min_axis[2])
# Исправляем значения вдоль малой оси на синус угла:
def correct_min(R):
return R / cos(incl * pi / 180)
radii_m = map(correct_min, radii)
plt.plot(radii, sig_los_maj, 's-', label='$\sigma_{los}^{maj}$')
plt.plot(radii_m, sig_los_min, 's-', label='$\sigma_{los}^{min}$')
plt.xlabel('$R$')
plt.ylabel('$\sigma$')
plt.legend()
plt.show()
Перегнем и приблизим полиномами:
In [14]:
sig_maj_data = zip(radii, sig_los_maj)
sig_maj_data = map(bind_curve, sig_maj_data)
sig_maj_data.sort()
radii_maj, sig_maj_p = zip(*sig_maj_data)
# Добавляем лишние точки чтобы протянуть дальше
num_fake_points = 10; expscale = 6
fake_radii, fake_sig = zip(*[(10.2 + i, 0.05* exp(- (10.2+i) / expscale )) for i in range(1, num_fake_points+1)])
poly_sig_maj = poly1d(polyfit(radii_maj + fake_radii, sig_maj_p + fake_sig, deg=9))
sig_min_data = zip(radii_m, sig_los_min)
sig_min_data = map(bind_curve, sig_min_data)
sig_min_data.sort()
radii_min, 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.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.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.show()
Посмотрим теперь на настоящие профили дисперсий скоростей в радиальном ($\sigma_{R}(R)$), азимутальном ($\sigma_{\varphi}(R)$) и вертикальном ($\sigma_{z}(R)$) направлениях:
In [15]:
sigR_real_points = np.loadtxt("C:\\science\\2FInstability\\data\\model\\1\\real\\svR_R.txt", float)
sigR_real_radii, sigR_real = zip(*sigR_real_points)
sigPhi_real_points = np.loadtxt("C:\\science\\2FInstability\\data\\model\\1\\real\\svphi_R.txt", float)
sigPhi_real_radii, sigPhi_real = zip(*sigPhi_real_points)
sigZ_real_points = np.loadtxt("C:\\science\\2FInstability\\data\\model\\1\\real\\svZ_R.txt", float)
sigZ_real_radii, sigZ_real = zip(*sigZ_real_points)
plt.plot(sigR_real_radii, sigR_real, '.-', color='red', label='$\sigma_{R}$')
plt.plot(sigPhi_real_radii, sigPhi_real, '.-', color='blue', label=r'$\sigma_{\varphi}$')
plt.plot(sigZ_real_radii, sigZ_real, '.-', color='black', label='$\sigma_{Z}$')
plt.legend()
plt.show()
Построим график отношения $\sigma_Z/\sigma_R$:
In [23]:
plt.plot(sigR_real_radii, [sigZ_real[i]/(sigR_real[i]+0.001) for i in range(0, sigR_real.__len__())], 's-', color='red')
plt.ylim(0, 1.0)
plt.ylabel('$\sigma_Z/\sigma_R$')
plt.xlabel('$R$')
plt.show()
Проверим, что эти профили соответствуют вырезанным из поля дисперсий. Связь профилей описывается следующими уравнениями: $$\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 [17]:
def restore_sig_los_maj(sigPhi, sigZ):
return sqrt(sigPhi**2 * sin(incl*pi/180)**2 + sigZ**2 * cos(incl*pi/180)**2)
def restore_sig_los_min(sigR, sigZ):
return sqrt(sigR**2 * sin(incl*pi/180)**2 + sigZ**2 * cos(incl*pi/180)**2)
sig_los_maj_restored = map(restore_sig_los_maj, sigPhi_real, sigZ_real)
sig_los_min_restored = map(restore_sig_los_min, sigR_real, sigZ_real)
plt.plot(sigR_real_radii, sig_los_maj_restored, '.', color='blue', label='$\sigma_{los}^{maj} restored$')
plt.plot(points, poly_sig_maj(points), label = '$\sigma_{los}^{maj} polyfit$', color='blue')
plt.plot(sigR_real_radii, sig_los_min_restored, '.', color='red', label='$\sigma_{los}^{min} restored$')
plt.plot(points, poly_sig_min(points), label = '$\sigma_{los}^{min} polyfit$', color='red')
plt.legend()
plt.show()
In [18]:
poly_sig_maj = poly1d(polyfit(sigR_real_radii, sig_los_maj_restored, deg=9))
poly_sig_min = poly1d(polyfit(sigR_real_radii, sig_los_min_restored, deg=9))
# plt.plot(sigR_real_radii, sig_los_maj_restored, '.', color='blue', label='$\sigma_{los}^{maj} restored$')
plt.plot(points, poly_sig_maj(points), label = '$\sigma_{los}^{maj} polyfit$', color='blue')
# plt.plot(sigR_real_radii, sig_los_min_restored, '.', color='red', label='$\sigma_{los}^{min} restored$')
plt.plot(points, poly_sig_min(points), label = '$\sigma_{los}^{min} polyfit$', color='red')
plt.legend()
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 [19]:
#Обрезаем данные по x > r_eff
r_eff = 1.0
#Значение sig_los_min в 0
sig_min_0 = poly_sig_min(0)
#Количество точек N и сами точки
N = 100
radii_points = np.arange(r_eff, 20.0, 20.0/N)
#Как вычислять ошибку
def residuals(params, xdata, ydata):
return (ydata - numpy.dot(xdata, params))
#Начальное приближение (А,В)
x0 = [0.01, 0.01]
#Уравнения: одно для 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))
print 'Solution: A = %s, B = %s' % (A, B)
print 'Chi^2:', chi2
#Подставить в уравнение в какой-нибудь точке:
# rrr = 5.0
# print (poly_sig_maj(rrr)**2) * sig_min_0**2 - (poly_sig_min(rrr)**2) * (A*f(rrr, Ro)+B)
Теперь восстановим исходные неизвестные - $\alpha$ и $\sigma_{R, 0}$:
In [20]:
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 = 0.12
# alpha = 0.4
print "sig_R_0 = %s, alpha = %s" % (sig_R_0, alpha)
Сравним полученные профили дисперсий скоростей с настоящими (заданными в модели):
In [21]:
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(sigR_real_radii, sigR_real, '.', color='red', label='$\sigma_{R}$')
plt.plot(sigR_real_radii, [sigR_exp(R) for R in sigR_real_radii], '-', color='red', label='$\sigma_{R, exp}$')
plt.plot(sigPhi_real_radii, sigPhi_real, '.', color='blue', label=r'$\sigma_{\varphi}$')
plt.plot(sigPhi_real_radii, [sigPhi_exp(R) for R in sigPhi_real_radii], '-', color='blue', label=r'$\sigma_{\varphi, exp}$')
plt.plot(sigZ_real_radii, sigZ_real, '.', color='black', label='$\sigma_{Z}$')
plt.plot(sigZ_real_radii, [sigZ_exp(R) for R in sigZ_real_radii], '-', color='black', label='$\sigma_{Z, exp}$')
plt.ylim(0, 0.15)
plt.legend()
plt.show()
И восстановим профили $\sigma_{los}^{maj}$ и $\sigma_{los}^{min}$:
In [22]:
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)
plt.plot(points, poly_sig_maj(points), '-', label = '$\sigma_{los}^{maj} polyfit$', color='blue')
plt.plot(points, [sig_maj_exp(R) for R in points], '--', color='blue', label='$\sigma_{maj, exp}$')
plt.plot(points, poly_sig_min(points), '-', label = '$\sigma_{los}^{min} polyfit$', color='red')
plt.plot(points, [sig_min_exp(R) for R in points], '--', color='red', label='$\sigma_{min, exp}$')
plt.legend()
plt.show()