NGC3245-Copy2-checkpoint


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

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


In [1]:
import matplotlib.pyplot as plt
import numpy as np
from numpy import poly1d, polyfit, power
import math
import scipy.optimize
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
from IPython.display import display
%matplotlib inline

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

#Наклон галактики по данным NED
incl=61.9

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

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


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


Out[2]:

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


Out[3]:

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

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


# UGC 5663 = NGC 3245

# Zasov+2012

# PA = 355 deg; PA_maj = 177.2 deg (from Mendez-Abreu+2008)

# not corrected for V_sys

# inclination = 61.9 deg (NED) (q_b = 0.52 - Mendez-Abreu+2008)

# data ARE NOT corrected for inclination

# r        vel       d_vel   sig     d_sig

# "        km/s      km/s    km/s    km/s

#(1)       (2)       (3)     (4)     (5)

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


In [6]:
pd.set_option('display.notebook_repr_html', True)
ma = pd.read_csv('v_stars_ma.dat', skiprows=9, sep='   ', engine='python')
# display(ma)

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

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

Данные Коткова по газу вдоль главной оси:


In [8]:
n3245pa355_gas = pd.read_csv('.\\gas\\n3245pa355_gas_corr.txt', sep=' ')
display(n3245pa355_gas)


r," v_Hbeta e_V max_Hb e_max sig_Hb e_sig v_OIII e_v max_OIII_1 e_max.1 sig_OIII e_sig.1 [OIII]/Hbeta e_oiii_hb S/N(HB) S/N(OIII)
0 -32.85 1518.62 18.8 777.49 116.0 116.87 21.80 1472.71 14.90 413.76 112.0 64.52 20.90 -0.10 0.2110 2.37 2.68
1 -21.65 1510.07 21.8 583.89 82.9 137.73 24.40 1471.74 10.90 559.51 215.0 43.05 18.70 -0.09 0.2700 2.74 3.87
2 -15.71 1426.74 20.0 622.31 182.0 74.43 26.70 1428.93 8.12 654.17 97.0 64.61 11.40 0.39 0.2240 2.57 6.36
3 -12.19 1415.81 20.0 552.01 168.0 72.51 26.90 1422.73 7.64 622.88 94.0 61.45 11.00 0.41 0.2310 2.75 7.45
4 -9.69 1424.72 13.8 1161.15 -0.0 34.67 5.02 1423.40 7.10 702.56 86.8 66.97 9.82 0.50 0.1040 4.38 7.46
5 -7.51 1414.65 21.0 895.11 -0.0 26.88 6.08 1420.45 8.43 553.72 50.2 94.05 10.30 0.77 0.1160 2.14 6.75
6 -6.04 1460.12 16.5 1209.44 -0.0 21.45 3.86 1414.45 9.20 444.38 40.3 100.74 11.00 0.67 0.0997 3.31 6.48
7 -5.17 1459.18 16.7 779.44 -0.0 32.58 5.75 1437.37 11.70 321.17 25.8 135.06 13.00 0.67 0.0941 3.37 6.68
8 -4.45 1431.82 14.8 678.43 174.0 66.60 20.50 1406.03 11.00 378.10 24.7 153.58 12.00 0.54 0.1800 3.93 6.23
9 -3.75 1433.98 13.8 786.38 187.0 66.84 19.00 1388.74 10.50 429.39 25.3 162.20 11.40 0.56 0.1660 5.08 6.46
10 -3.21 1386.35 25.9 306.33 57.9 127.45 29.40 1386.46 14.30 239.67 17.5 176.53 15.30 0.47 0.1390 3.95 5.84
11 -2.86 1354.83 27.5 307.29 52.2 146.54 30.30 1376.24 13.80 254.53 18.6 169.64 14.80 0.41 0.1260 3.35 6.27
12 -2.50 1312.62 32.7 282.93 44.8 181.42 34.90 1369.13 12.90 279.69 19.6 166.83 13.90 0.39 0.1180 2.68 7.37
13 -2.14 1300.70 32.0 306.16 43.8 194.76 33.70 1352.33 13.20 287.77 18.9 181.00 14.10 0.37 0.1070 2.59 7.30
14 -1.78 1287.60 24.6 406.14 46.0 190.16 26.00 1327.75 12.80 311.89 18.5 192.13 13.50 0.32 0.0869 3.44 7.65
15 -1.43 1277.23 21.1 477.60 48.1 184.04 22.40 1320.92 14.20 296.78 17.5 213.98 14.80 0.29 0.0791 4.05 6.90
16 -1.07 1228.67 17.9 563.88 50.2 178.81 19.20 1319.14 13.20 322.41 17.9 211.64 13.80 0.26 0.0710 4.64 7.54
17 -0.71 1235.10 18.8 553.77 48.8 188.12 20.00 1303.31 13.20 327.83 17.9 215.06 13.80 0.26 0.0702 4.46 8.03
18 -0.36 1245.86 17.8 589.83 49.6 186.25 18.90 1294.89 12.20 351.97 18.6 206.07 12.80 0.25 0.0673 4.97 8.18
19 0.00 1199.11 14.5 687.68 57.0 157.90 15.80 1291.62 13.10 337.72 17.7 220.62 13.60 0.27 0.0664 5.91 8.94
20 0.36 1183.38 13.8 722.69 56.7 158.35 14.90 1279.18 12.60 342.49 18.5 207.72 13.20 0.23 0.0645 5.66 8.90
21 0.71 1179.53 13.3 732.64 58.1 152.65 14.50 1250.51 11.90 345.76 19.8 186.44 12.70 0.19 0.0662 5.91 8.25
22 1.07 1153.74 12.0 782.75 63.7 136.21 13.30 1207.70 12.40 331.94 19.1 192.48 13.10 0.21 0.0675 7.31 8.49
23 1.43 1143.83 12.7 726.16 61.1 139.23 14.10 1208.31 11.40 353.22 19.0 189.27 12.00 0.25 0.0677 6.39 8.78
24 1.78 1142.40 14.5 616.37 62.0 133.54 16.20 1209.10 11.20 345.03 19.3 180.38 11.90 0.31 0.0782 5.51 8.22
25 2.14 1136.37 12.3 695.72 68.9 117.51 14.10 1209.77 12.90 291.30 18.6 181.74 13.80 0.24 0.0801 7.27 6.82
26 2.50 1144.32 12.7 650.97 70.4 111.99 14.70 1220.33 14.10 259.48 18.1 181.00 15.00 0.24 0.0876 6.48 7.20
27 2.86 1122.76 13.8 575.27 76.9 101.35 16.50 1215.14 14.20 244.82 18.3 171.14 15.30 0.29 0.1040 6.47 7.10
28 3.21 1115.41 12.6 608.64 81.8 94.11 15.40 1182.77 14.00 237.23 18.5 163.22 15.20 0.26 0.1060 6.83 6.86
29 3.74 1107.50 12.2 838.54 103.0 98.62 14.70 1189.17 10.70 410.56 26.2 152.66 11.70 0.31 0.0946 5.69 6.70
30 4.46 1108.55 16.0 595.43 91.6 101.86 19.10 1183.51 14.50 287.58 22.9 164.12 15.60 0.32 0.1180 5.62 5.73
31 5.15 1086.38 18.2 493.95 123.0 78.35 23.60 1196.13 15.00 255.74 22.3 155.97 16.30 0.45 0.1800 3.76 4.87
32 6.09 1119.45 16.4 1290.67 -0.0 19.61 3.53 1189.58 14.20 306.22 24.2 163.00 15.40 0.73 0.0948 2.74 5.31
33 7.87 1146.32 24.0 583.28 199.0 76.12 31.60 1208.87 10.90 529.97 43.1 125.56 12.30 0.61 0.2400 2.25 6.08
34 11.45 1138.77 26.4 953.58 -0.0 23.37 6.64 1226.82 15.30 379.64 38.2 140.48 17.00 0.81 0.1410 1.77 3.70
35 15.82 1144.73 18.3 877.53 699.0 38.22 34.10 1220.65 13.40 375.41 49.6 100.65 16.10 0.48 0.5270 2.57 3.90
36 26.11 1128.28 20.7 764.22 170.0 93.15 25.70 1176.62 16.70 403.16 117.0 67.22 23.30 0.01 0.2500 2.63 2.71

In [9]:
try:
    from PIL import Image
except:
    import Image

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

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

# Данные по кинематике газа Засова 2012 вдоль большой полуоси, не исправленные за наклон (они же Катков)
# Водород
r_g_H = n3245pa355_gas['r,"']
vel_g_H = n3245pa355_gas['v_Hbeta']
e_vel_g_H = n3245pa355_gas['e_V']

# Кислород
r_g_O = n3245pa355_gas['r,"']
vel_g_O = n3245pa355_gas['v_OIII']
e_vel_g_O = n3245pa355_gas['e_v']


fig, subplot = plt.subplots(2, 1)
subplot[0].plot(r_ma, vel_ma, '.-', label="Zasov 2012, maj")
subplot[0].plot(r_mi, vel_mi, '.-', label="Zasov 2012, min")
subplot[0].plot(r_g_H, vel_g_H, '^', label="gas Hbeta Zasov 2012", mfc='none')
subplot[0].plot(r_g_O, vel_g_O, 's', label="gas OIII Zasov 2012", mfc='none')
subplot[0].legend()

subplot[1].imshow(np.asarray(Image.open("zasov2012_fig3_part.png")))
subplot[1].set_title("From Z2012 fig.3 part")
plt.plot()


Out[9]:
[]

Теперь построим график дисперсий скоростей на луче зрения вдоль большой и малой осей по данным Засова:


In [10]:
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, sig_mi, 's-', label='$\sigma_{los}^{min}$')
plt.errorbar(r_mi, sig_mi, yerr=e_sig_mi, fmt='.', marker='.', mew=0, color='black')
plt.xlabel('$R$')
plt.ylabel('$\sigma$')
plt.legend()
plt.show()


Что мы видим? Даже без исправления значения вдоль малой оси на косинус угла видно, что центры двух кривых совсем не совпали, а это значит, что щели спектрометра не были точно соорентированы на центр и/или по осям галактики. Однако рассматриваемый случай оказывается достаточно простым - кривая вращения вдоль малой оси почти везде ноль, вдоль большой - симметрична, значит ${PA}_{maj}$ было определено правильно, просто щель для большой оси оказалась немного смещена. Данные и по скоростям, и по дисперсиям необходимо исправить, методика исправления описана ниже.


In [11]:
from mpl_toolkits.axes_grid.axislines import SubplotZero

fig = plt.figure(1)
ax = SubplotZero(fig, 111)
fig.add_subplot(ax)

for direction in ["xzero", "yzero"]:
    ax.axis[direction].set_axisline_style("-|>")
    ax.axis[direction].set_visible(True)

for direction in ["left", "right", "bottom", "top"]:
    ax.axis[direction].set_visible(False)

x = np.linspace(-1., 1., 100)
y = np.zeros(100)
ax.plot(x, y)

from pylab import figure, show, rand
from matplotlib.patches import Ellipse
e = Ellipse(xy=(0,0), width=1.9, height=1.0, angle=0.0)

ax.add_artist(e)
e.set_clip_box(ax.bbox)
e.set_alpha(0.1)
e.set_facecolor('blue')

ax.set_ylim(-0.75, 0.75)
ax.axes.get_xaxis().set_ticks([])
ax.axes.get_yaxis().set_ticks([])

ax.text(0.98, 0.02, "$x$", fontsize=20)
ax.text(0.02, 0.73, "$y$", fontsize=20)
ax.plot((0.0, 0.2), (0.0, 0.2), '--', lw=0.3)
ax.text(0.05, 0.015, r"$\varphi$", fontsize=20)
ax.plot((0.2, 0.2), (-0.49, 0.49), '-', lw=1.5, color='r')
ax.plot((-0.87, 0.87), (0.2, 0.2), '-', lw=1.5, color='r')
ax.plot(0.2, 0.2, 's')
ax.text(0.01, 0.11, "$y_0$", fontsize=20)
ax.text(0.07, -0.05, "$x_0$", fontsize=20)
ax.plot((0.0, 0.2), (0.0, 0.0), '-', lw=2.0, color='k')
ax.plot((0.0, 0.0), (0.0, 0.2), '-', lw=2.0, color='k')
ax.text(-0.3, 0.22, "$P_{obs}$", fontsize=20)
ax.arrow(0, 0.2, -0.27, 0.0, head_width=0.02, head_length=0.02, fc='k', ec='k', width=0.005)
ax.plot((0., -0.3), (0., 0.2), '--', lw=0.3)
ax.text(-0.06, 0.06, r"$\theta$", fontsize=20)
ax.text(-0.13, 0.16, "$R_{obs}$", fontsize=20)
ax.plot(-0.3, 0.2, 's')

plt.show()


Пусть центр координат находится в истинном центре галактики, углы осей ${PA}_{maj}$ и ${PA}_{min}$ определены правильно, но щели спектрометра пересекаются в точке $(x_0,y_0)$. Угол $\varphi$ очевидно равен $$\varphi=\arctan \frac{y_0}{x_0}$$ Тогда расстояния и скорости для точки $P_{obs}$ для соответствующей оси (на рисунке выше - для большой) пересчитываются по следующему правилу: $$R_{real} = R_{obs}\times\frac{\sqrt{\sin^2\theta\cos^2i+\cos^2\theta}}{\cos i\sin \theta},$$ $$V_{real} = V_{obs}\times\frac{\sqrt{\sin^2\theta\cos^2i+\cos^2\theta}}{\sin i\cos i\sin \theta},$$ где $\theta$ есть угол на соответствующую точку.

В нашем случае из-за того, что $\sigma_{los, maj}(0) < \sigma_{los, min}(0)$, а также по причинам, которые были указаны выше, смещена щель только вдоль большой оси, угол $\varphi=90^{\circ}$. Угол $\theta$ и его тригонометрические функции в свою очередь считаются как $$\theta = \arctan \frac{R_{obs}}{y_0}, \sin \theta = \frac{R_{obs}}{\sqrt{R_{obs}^2 + y_0^2}}, \cos \theta = \frac{y_0}{\sqrt{R_{obs}^2 + y_0^2}}$$ После восстановления скоростей исправлять за наклон галактики профиль уже не нужно (случай $\theta=0^{\circ}$).

Остается лишь вопрос, как найти $y_0$. В рассматриваемом случае это сделать просто (но не совсем честно) - достаточно взять исправленный профиль $\sigma_{los}^{min}(R)$ и посмотреть, на каком расстоянии значения будут равны значениям $\sigma_{los}^{maj}$ в центре. Иными словами, $y_0$ определяется условием $$\sigma_{los}^{min}(y_0) = \sigma_{los}^{maj}(0)$$ Найдем это значение, после чего исправим соответствующим образом кривую вращения.


In [12]:
# Исправляем значения вдоль малой оси на синус угла:    
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$')

#Подбираем значения y_0
plt.axhline(y=sig_ma[sig_ma.__len__()/2-1])
plt.axvline(x=-6.0)
plt.axvline(x=4.2)
plt.text(-60, 150, r"$y_0 = \frac{6+4.2}{2}=5.1$", fontsize=25)
plt.legend()
plt.show()


Также не забудем, что полученное значение $y_0$ надо исправить за $\cos i$, ибо мы считали его по малой оси:


In [13]:
y_0 = 5.1*cos(incl * pi / 180)

def sin_theta(R_obs):
    return R_obs/sqrt(R_obs**2 + y_0**2)

def cos_theta(R_obs):
    return y_0/sqrt(R_obs**2 + y_0**2)

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

def correct_radii_maj(R_obs):
    return R_obs * sqrt(sin_theta(R_obs)**2 * cos_i**2 + cos_theta(R_obs)**2) / cos_i /  abs(sin_theta(R_obs))

def correct_velocity_maj(V_obs, R_obs):
    return V_obs * sqrt(sin_theta(R_obs)**2 * cos_i**2 + cos_theta(R_obs)**2) / (sin_i * cos_i * abs(sin_theta(R_obs)))

In [14]:
r_ma1 = map(correct_radii_maj, r_ma)
vel_ma1 = map(correct_velocity_maj, map(lambda l: abs(l-1300.5), vel_ma), r_ma)
e_vel_ma1 = map(correct_velocity_maj, e_vel_ma, r_ma)

r_g_H1 = map(correct_radii_maj, r_g_H)
vel_g_H1 = map(correct_velocity_maj, map(lambda l: abs(l-1300.5), vel_g_H), r_g_H)
e_vel_g_H1 = map(correct_velocity_maj, e_vel_g_H, r_g_H)

r_g_O1 = map(correct_radii_maj, r_g_O)
vel_g_O1 = map(correct_velocity_maj, map(lambda l: abs(l-1300.5), vel_g_O), r_g_O)
e_vel_g_O1 = map(correct_velocity_maj, e_vel_g_O, r_g_O)

In [15]:
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, 1300.5, 90.)
r_ma_b, vel_ma_b, e_vel_b = correct_rotation_curve(r_ma1, vel_ma1, e_vel_ma1,  0., 0., 90.)
r_mi_b, vel_mi_b, e_vel_mi_b = correct_rotation_curve(r_mi, vel_mi, e_vel_mi,  0., 1276.5, 90.)

r_g_H_b, vel_g_H_b, e_vel_g_H_b = correct_rotation_curve(r_g_H1, vel_g_H1, e_vel_g_H1,  0., 0., 90.)
r_g_O_b, vel_g_O_b, e_vel_g_O_b = correct_rotation_curve(r_g_O1, vel_g_O1, e_vel_g_O1,  0., 0., 90.)

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.plot(r_g_H_b, vel_g_H_b, '^', label = 'Zasov gas Hbeta')
plt.plot(r_g_O_b, vel_g_O_b, 's', label = 'Zasov gas OIII')
plt.ylim(-10, 500)
plt.legend()
plt.plot()


Out[15]:
[]

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


In [16]:
left_bord = 7.
def cut_left(x):
    return filter(lambda l: l > left_bord, x)

start_ind = r_ma_b.index(cut_left(r_ma_b)[0])
r_ma_b = (0.0,) + r_ma_b[start_ind:]
vel_ma_b = (0.0,) + vel_ma_b[start_ind:]
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()



In [17]:
# os.chdir("C:\\Users\\root\\Dropbox\\RotationCurves\\PhD\\paper1\\text\\imgs")

# np.save('n3245_maj_rot', zip(r_ma_b, vel_ma_b, e_vel_b))
# np.save('n3245_rot_poly', zip(test_points, poly_star(test_points)))

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

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


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


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


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


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

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


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

# r_mi_extend = map(correct_min, r_mi)
    
plt.plot(r_ma1, sig_ma, 's-', label='$\sigma_{los}^{maj}$')
plt.errorbar(r_ma1, 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 [22]:
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) 

# poly_sig_maj = poly1d(polyfit(radii_maj, sig_maj_p, deg=9))

sig_maj_data = zip(r_ma1, 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) 

poly_sig_maj = poly1d(polyfit(radii_maj, sig_maj_p, 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) 

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

poly_sig_min = poly1d(polyfit(radii_min + fake_radii, sig_min_p + fake_sig, 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, 250)
plt.show()



In [23]:
plt.plot(radii_min, e_sig_min_p, 'o')
plt.plot(radii_maj, e_sig_maj_p, 'v')


Out[23]:
[<matplotlib.lines.Line2D at 0xc740c18>]

In [23]:
os.chdir("C:\\Users\\root\\Dropbox\\RotationCurves\\PhD\\paper1\\text\\imgs")

np.save('n3245_maj', zip(radii_maj, sig_maj_p, e_sig_maj_p))
np.save('n3245_min', zip(radii_min, sig_min_p, e_sig_min_p))
np.save('n3245_min_poly', zip(points, poly_sig_min(points)))
np.save('n3245_maj_poly', zip(points, poly_sig_maj(points)))

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

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


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

print "Poly chi^2 for maj = %s, mean = %s" % (sqerr_maj, sqerr_maj / sig_maj_p.__len__())
print "Poly chi^2 for min = %s, mean = %s" % (sqerr_min, sqerr_min / sig_min_p.__len__())


Poly chi^2 for maj = 945.654732111, mean = 12.2812302872
Poly chi^2 for min = 1045.379644, mean = 24.8899915237

Методика восстановления профилей $\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. $$

Однако это все было верно для случая несмещенной щели. Для нашего случая исходное уравнение для дисперсий вдоль большой оси меняется на $$\sigma_{los,maj}^2(R_{obs})=\sigma_R^2(R)[f(R)\sin^2\theta+\cos^2\theta]\sin^2i + \alpha^2\sigma_R^2(R)\cos^2i,$$ Если ввести новую функцию $f^{\prime}(R, R_{obs})=f(R)\sin^2\theta+\cos^2\theta$, то в тех же переменных, что и выше система для МНК будет записана как $$\left\{ \begin{array}{lr} \sigma_{los,maj}^2(R_{obs}(R_j))\times \sigma_{min,0}^2 =\sigma_{los,min}^2(R_j)[Af^{\prime}(R_j)+B]\\ \sigma_{min,0}^2=A+B \end{array} \right. $$

Попробуем поискать точки, в которых МНК решается лучше всего:


In [25]:
#Новая функция f'
def f_prime(R_real, R_obs):
    return f(R_real, Ro)*sin_theta(R_obs)**2 + cos_theta(R_obs)**2
    
#Как вычислять ошибку 
def residuals(params, xdata, ydata):
    return (ydata - np.dot(xdata, params))

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

#Восстановление первоначальных неизвестных
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 = [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_observ = np.arange(left, right, (right-left)/n)
            r_points_real = map(correct_radii_maj, r_points_observ)
            eq_left = np.concatenate( ([sig_min_0**2], 
                              [poly_sig_maj(p)**2 * sig_min_0**2  for p in r_points_observ]) )
            eq_right = np.transpose(np.array([
                            np.concatenate(([1.0],
                                               [poly_sig_min(R)**2 * f_prime(R, Rg) for (R, Rg) 
                                                in zip(r_points_real, r_points_observ)])),
                            np.concatenate(([1.0],
                                                [poly_sig_min(R)**2 for R in r_points_real]))]))
            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])


For N=50 top-10 best results:
	1 place: chi2 = 4.15897799029e+13 in range [12.0:32.0]; sig_R_0 = 200.73 alpha = 0.547
	2 place: chi2 = 4.34299841055e+13 in range [14.0:34.0]; sig_R_0 = 192.783 alpha = 0.69
	3 place: chi2 = 5.24955320833e+13 in range [12.0:34.0]; sig_R_0 = 208.448 alpha = 0.339
	4 place: chi2 = 6.89990143099e+13 in range [14.0:36.0]; sig_R_0 = 207.262 alpha = 0.348
For N=75 top-10 best results:
	1 place: chi2 = 4.03934849956e+13 in range [12.0:32.0]; sig_R_0 = 199.957 alpha = 0.563
	2 place: chi2 = 4.43850620931e+13 in range [14.0:34.0]; sig_R_0 = 193.204 alpha = 0.682
	3 place: chi2 = 5.19888202772e+13 in range [12.0:34.0]; sig_R_0 = 207.881 alpha = 0.355
	4 place: chi2 = 7.02488154073e+13 in range [14.0:36.0]; sig_R_0 = 207.885 alpha = 0.329
For N=100 top-10 best results:
	1 place: chi2 = 3.98148534633e+13 in range [12.0:32.0]; sig_R_0 = 199.571 alpha = 0.571
	2 place: chi2 = 4.4865003186e+13 in range [14.0:34.0]; sig_R_0 = 193.424 alpha = 0.678
	3 place: chi2 = 5.1758751359e+13 in range [12.0:34.0]; sig_R_0 = 207.6 alpha = 0.363
	4 place: chi2 = 7.08704846831e+13 in range [14.0:36.0]; sig_R_0 = 208.207 alpha = 0.318
Best of the best: N=100 chi2 = 3.98148534633e+13 on range[12.0:32.0]
C:\Anaconda\lib\site-packages\scipy\optimize\minpack.py:427: RuntimeWarning: Number of calls to function has reached maxfev = 600.
  warnings.warn(errors[info][0], RuntimeWarning)

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

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

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

#Пересчитываем точки в настоящие, чтобы уметь считать углы
radii_points_real = map(correct_radii_maj, radii_points_observ)

#Уравнения:  одно для 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_observ]) )
#Правая часть:
eq_right = np.transpose(
            np.array([
                np.concatenate(([1.0],
                                   [poly_sig_min(R)**2 * f_prime(R, Rg) for (R, Rg) 
                                    in zip(radii_points_real, radii_points_observ)])),
                np.concatenate(([1.0],
                                    [poly_sig_min(R)**2 for R in radii_points_real]))]))

# МНК для получившихся уравнений:
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 = 30992.5264638, B = 2876.59870463
Chi^2: 3.98148534633e+13

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


In [27]:
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 = 255.
# alpha = 0.01

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


sig_R_0 = 199.571, alpha = 0.571

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


In [28]:
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):
    if sigPhi_to_sigR_real(R) > 0:
        return sqrt(sigPhi_to_sigR_real(R)) * sigR_exp(R)
    else:
        return 0


plt.plot(points, [sigR_exp(R) for R in points], '-', color='red', label='$\sigma_{R, exp}$')
plt.plot(points, [sigPhi_exp(R) for R in points], '-', color='blue', label=r'$\sigma_{\varphi, exp}$')
plt.plot(points, [sigZ_exp(R) for R in points], '-', color='black', label='$\sigma_{Z, exp}$')
plt.legend()
plt.show()


И восстановим профили $\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$$

Однако надо не забыть исправить $\sigma^{maj}_{los}$: $$\sigma_{los,maj}^2(R_{obs})=[\sigma_R^2(R)\cos^2\theta + \sigma^2_{\varphi}(R)\sin^2\theta]\sin^2i + \sigma^2_{Z}(R)\cos^2i$$


In [29]:
# def sig_min_exp(R):
#     return sqrt(sigR_exp(R)**2 * sin_i**2 + sigZ_exp(R)**2 * cos_i**2)

def sig_maj_exp(R_real, R_obs):
    return sqrt((sigR_exp(R_real)**2 * cos_theta(R_obs)**2 + sigPhi_exp(R_real)**2 * sin_theta(R_obs)**2) * sin_i**2 + 
                sigZ_exp(R_real)**2 * cos_i**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)

points_obs = np.arange(0, max(radii_min), 0.1)
points_real = map(correct_radii_maj, points_obs)

# points_original = map(lambda l: g(l)*l, points)

plt.plot(points, poly_sig_maj(points), '-', label = '$\sigma_{los}^{maj} polyfit$', color='blue')
plt.plot(points_obs, [sig_maj_exp(Rr, Robs) for (Rr, Robs) in zip(points_real, points_obs)], 
         '--', 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.ylim(0, 300)
plt.show()


И посчитаем невязки для восстановленного профиля:


In [30]:
ppp =  map(lambda l: abs(l), r_ma)
ppp.sort()

sqerr_maj_final = sum(power([sig_maj_exp(p[0], p[1]) - p[2] for p in zip(radii_maj, ppp, sig_maj_p)], 2))
sqerr_min_final = sum(power([sig_min_exp(p[0]) - p[1] for p in sig_min_data], 2))

print "Poly chi^2 for maj = %s, mean = %s" % (sqerr_maj_final, sqerr_maj_final / sig_maj_p.__len__())
print "Poly chi^2 for min = %s, mean = %s" % (sqerr_min_final, sqerr_min_final / sig_min_p.__len__())


Poly chi^2 for maj = 14819.6594602, mean = 192.463109872
Poly chi^2 for min = 16480.8888173, mean = 392.402114698

Заметим, что можно было не решать систему МНК, а честно разрешить систему из двух уравнений $$\left\{ \begin{array}{lr} \sigma_{los,maj}^2(R_{obs}(R_j))\times \sigma_{min,0}^2 =\sigma_{los,min}^2(R_j)[Af^{\prime}(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^{\prime}(R)-1}\times(P(R, R_{obs})-1)$$ $$\alpha^2 = \tan^2 i\frac{f^{\prime}(R) - P(R, R_{obs})}{P(R, R_{obs})-1},$$ $$P(R, R_{obs})=\frac{\sigma_{los,maj}^2(R_{obs})}{\sigma_{los, min}^2(R)}$$ Имеет смысл также искать не $\alpha$, а $\alpha\cdot\sigma_{R,0}=\sigma_{Z,0}$.


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

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

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

def direct_find_sig_R_0(R, R_obs):
    return direct_solve_A(R, R_obs)

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

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


In [32]:
p_r = radii_maj
p_obs = map(lambda l: abs(l), r_ma)
p_obs.sort()

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

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, 320)
plt.show()


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


In [33]:
#Обрежем по 15
q=15.
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 = 284.608074765

Попробуем найти $\alpha$ как и раньше - решением избыточной линейной системы. Для этого вспомним уравнения: $$\sigma_{los,maj}^2(R_{obs})=\frac{\sigma_{R,0}^2\sigma_{los,min}^2(R)[f^{\prime}\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$$ Как мы видим, получается очевидная линейная система относительно $\alpha^2$.


In [34]:
sig_R_0 = poly_q[0]

def residuals(param, xdata, ydata):
            return (ydata - xdata*param)

x0 = [1.]
r_eff = 10.0
r_max = 40.0
N = 30
radii_points_observ = np.arange(r_eff, r_max, (r_max-r_eff)/N)
radii_points_real = map(correct_radii_maj, radii_points_observ)

#Уравнения - N штук для maj
#Левая часть:
eq_left = np.array([poly_sig_maj(Rg)**2 * sig_min_0**2 - 
                               poly_sig_min(R)**2 * f_prime(R, Rg) * sig_R_0**2 * sin_i**2
                               for (R, Rg) in zip(radii_points_real, radii_points_observ)])
#Правая часть:
eq_right = np.array([poly_sig_min(R)**2 * sig_R_0**2 * cos_i**2 for R in radii_points_real])


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

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


Solution: alpha^2 = -0.934488238496
Chi^2: 2.81292202956e+14

Не решается.


In [35]:
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, [sig_min_exp(R) for R in points], '--', color='red', label='$\sigma_{min, exp}$')

#Строим полученный на основе среднего профиль 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, 320)
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 [36]:
def sig_R_upper_lim(R):
    """Оценка сверху на sigR(R)"""
    return poly_sig_min(R)/sqrt(sin_i**2 + 0.09*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) for R in points], label = '$\sigma_R^{up}$', color='g')

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


Как видим, значения действительно оказались переоценены минимум на 40 км/c. Давайте ради интереса попробуем восстановить эллипсоид скоростей и исходные профили, исходя из условия об маржинальной устойчивости диска относительно изгибных возмущений:


In [37]:
sig_R_0 = 245.
alpha = 0.3

plt.plot(points, [sigR_exp(R) for R in points], '-', color='red', label='$\sigma_{R, exp}$')
plt.plot(points, [sigPhi_exp(R) for R in points], '-', color='blue', label=r'$\sigma_{\varphi, exp}$')
plt.plot(points, [sigZ_exp(R) for R in points], '-', color='black', label='$\sigma_{Z, exp}$')
plt.legend()
plt.show()