NGC3245-Copy1-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
import matplotlib.cm as cm
%matplotlib inline

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

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

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

r_eb = 15.1

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


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]:
from scipy.interpolate import UnivariateSpline

spl = UnivariateSpline(r_ma_b, vel_ma_b, k=4, s=10000)

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.plot(test_points, spl(test_points), '-', color='m')
plt.xlabel('$R$'); plt.ylim(0)
plt.ylabel('$V^{maj}_{\phi}(R)$')
plt.show()



In [18]:
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, spl(test_points)))

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

In [19]:
plt.plot(test_points, spl.derivative()(test_points), '.-')
plt.show()


Кривая вращения нам нужна для нахождения соотношения $\sigma_{\varphi}^{2}/\sigma_{R}^{2}$, которое описывается уравнением ${\displaystyle \sigma_{\varphi}^{2}/\sigma_{R}^{2}=0.5\left(1+\frac{R}{\bar{v}_{\varphi}}\frac{d\bar{v}_{\varphi}}{dR}\right)}$ (Binney & Tremaine, 1987) и приближается гладко функцией $f=0.5(1+e^{-R/R_{0}}),$ где $R_{0}$ --- характерный масштаб.

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

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


In [20]:
def sigPhi_to_sigR_real(R):
        return 0.5 * (1 + R*spl.derivative()(R) / spl(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.plot(r_ma_b, np.array(vel_ma_b)/250, 'x-', color='blue', markersize=6)
plt.plot(test_points, spl(test_points)/250, '-', color='m')

plt.show()


Имеет смысл обрезать данные по 60''

Построим графики дисперсий скоростей на луче зрения вдоль большой и малой оси ($\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]:
import scipy.interpolate as inter
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))
poly_sig_maj = inter.UnivariateSpline (radii_maj, sig_maj_p, k=4, s=10000)

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))
poly_sig_min = inter.UnivariateSpline (radii_min + fake_radii, sig_min_p + fake_sig, k=4, s=10000)

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]:
sig_maj_data = zip(radii_maj, sig_maj_p, e_sig_maj_p)
sig_maj_data = filter(lambda l: l[0] < 60. and l[0] > r_eb, sig_maj_data)
radii_maj, sig_maj_p, e_sig_maj_p = zip(*sig_maj_data) 

sig_min_data = zip(radii_min, sig_min_p, e_sig_min_p)
sig_min_data = filter(lambda l: l[0] < 60. and l[0] > r_eb, sig_min_data)
radii_min, sig_min_p, e_sig_min_p = zip(*sig_min_data) 

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 [24]:
plt.hist(e_sig_maj_p, color='blue')
plt.hist(e_sig_min_p, alpha=0.5, color='red')
plt.show()



In [25]:
plt.hist(map(lambda l: l[0]/l[1], zip(e_sig_maj_p, sig_maj_p)), color='blue')
plt.hist(map(lambda l: l[0]/l[1], zip(e_sig_min_p, sig_min_p)), alpha=0.5, color='red')
plt.show()



In [26]:
# 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")

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

def f(R, Ro):
    return sigPhi_to_sigR_real(R)

#Новая функция 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(radii_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 = 1.40000092904e+13 in range [26.0:46.0]; sig_R_0 = 56.164 alpha = 3.01
	2 place: chi2 = 1.50971179446e+13 in range [24.0:46.0]; sig_R_0 = 16.274 alpha = 11.213
	3 place: chi2 = 1.81023713002e+13 in range [16.0:44.0]; sig_R_0 = 30.66 alpha = 5.865
	4 place: chi2 = 1.81280726426e+13 in range [16.0:42.0]; sig_R_0 = 19.413 alpha = 9.441
	5 place: chi2 = 1.8443232415e+13 in range [28.0:48.0]; sig_R_0 = 110.447 alpha = 1.046
	6 place: chi2 = 1.9346697932e+13 in range [16.0:40.0]; sig_R_0 = 15.203 alpha = 12.112
	7 place: chi2 = 2.05217158105e+13 in range [16.0:38.0]; sig_R_0 = 23.2 alpha = 7.853
	8 place: chi2 = 2.05741150747e+13 in range [16.0:46.0]; sig_R_0 = 43.274 alpha = 4.02
	9 place: chi2 = 2.07617051378e+13 in range [16.0:36.0]; sig_R_0 = 36.349 alpha = 4.875
	10 place: chi2 = 2.25670586433e+13 in range [26.0:48.0]; sig_R_0 = 85.285 alpha = 1.702
For N=75 top-10 best results:
	1 place: chi2 = 1.42010256967e+13 in range [26.0:46.0]; sig_R_0 = 58.395 alpha = 2.873
	2 place: chi2 = 1.54093278912e+13 in range [24.0:46.0]; sig_R_0 = 21.708 alpha = 8.359
	3 place: chi2 = 1.73849315905e+13 in range [16.0:42.0]; sig_R_0 = 15.483 alpha = 11.882
	4 place: chi2 = 1.7506698834e+13 in range [16.0:44.0]; sig_R_0 = 28.833 alpha = 6.255
	5 place: chi2 = 1.85468629959e+13 in range [28.0:48.0]; sig_R_0 = 111.993 alpha = 1.009
	6 place: chi2 = 1.85724658179e+13 in range [16.0:40.0]; sig_R_0 = 8.383 alpha = 22.066
	7 place: chi2 = 1.98070177456e+13 in range [16.0:38.0]; sig_R_0 = 18.954 alpha = 9.664
	8 place: chi2 = 2.01624591184e+13 in range [16.0:36.0]; sig_R_0 = 33.62 alpha = 5.303
	9 place: chi2 = 2.0247482781e+13 in range [16.0:46.0]; sig_R_0 = 42.379 alpha = 4.113
	10 place: chi2 = 2.28047201626e+13 in range [26.0:48.0]; sig_R_0 = 87.078 alpha = 1.646
For N=100 top-10 best results:
	1 place: chi2 = 1.42991956329e+13 in range [26.0:46.0]; sig_R_0 = 59.49 alpha = 2.809
	2 place: chi2 = 1.55636435571e+13 in range [24.0:46.0]; sig_R_0 = 23.995 alpha = 7.541
	3 place: chi2 = 1.70177845244e+13 in range [16.0:42.0]; sig_R_0 = 13.092 alpha = 14.078
	4 place: chi2 = 1.72164237737e+13 in range [16.0:44.0]; sig_R_0 = 27.884 alpha = 6.477
	5 place: chi2 = 1.85946247453e+13 in range [28.0:48.0]; sig_R_0 = 112.763 alpha = 0.991
	6 place: chi2 = 1.94491755438e+13 in range [16.0:38.0]; sig_R_0 = 16.423 alpha = 11.183
	7 place: chi2 = 1.98610166019e+13 in range [16.0:36.0]; sig_R_0 = 32.167 alpha = 5.559
	8 place: chi2 = 2.00946267746e+13 in range [16.0:46.0]; sig_R_0 = 41.933 alpha = 4.161
	9 place: chi2 = 2.29189721451e+13 in range [26.0:48.0]; sig_R_0 = 87.969 alpha = 1.619
	10 place: chi2 = 2.35019378976e+13 in range [18.0:48.0]; sig_R_0 = 37.326 alpha = 4.678
Best of the best: N=50 chi2 = 1.40000092904e+13 on range[26.0:46.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 [28]:
#Обрезаем данные по 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 = 13473.368213, B = -704.000504481
Chi^2: 9.87693568894e+13

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


In [29]:
# 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 = 224.
alpha = 0.4

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


sig_R_0 = 224.0, alpha = 0.4

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


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

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

def sigPhi_exp(R):
    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 [31]:
# 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_maj_exp(R, R_obs):
    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)

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.errorbar(radii_maj, sig_maj_p, yerr=e_sig_maj_p, fmt='.', marker='.', mew=0, color='blue')
plt.errorbar(radii_min, sig_min_p, yerr=e_sig_min_p, fmt='.', marker='.', mew=0, color='red')
plt.legend()
plt.ylim(0, 300)
plt.show()


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


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

r_ma_bind = [abs(R) for R in r_ma]
r_ma_bind.sort()

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 = calc_chi2_normal(sig_maj_p, e_sig_maj_p, [sig_maj_exp(r, 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_maj = calc_chi2_normal(sig_maj_p, e_sig_maj_p, [sig_maj_exp(Robs, p[0]) for 
#                                    (Robs, p) in zip(r_ma_bind, sig_maj_data)])
#             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


image, image_maj, image_min = compute_chi2_maps(alphas=alphas, sigmas=sigmas)

In [33]:
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 = cm.colors.Normalize(vmax=image.max(), vmin=-image.max())
            cmap = cm.PRGn
            levels = np.linspace(start=image_log.min(), stop=vmax, num=10)
            cset=ax.contour(image_log, levels, hold='on', colors = 'k', origin='lower', extent=[alphas[0],alphas[-1],sigmas[0],sigmas[-1]])
        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=40.)
plot_chi2_map(image_maj, axes[1], log_scale=False, title='$\chi^2_{maj}$', is_contour=False, vmax=40.)
plot_chi2_map(image_min, axes[2], log_scale=False, title='$\chi^2_{min}$', is_contour=False, vmax=10.)
plt.show()