Эксперименты по восстановлению профилей дисперсий в трех направлениях для 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 0xcdab6a0>]

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