Сначала всякие настройки и импорты
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)
Выведем также для удобства данные из обоих файлов:
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)
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])
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)
Теперь восстановим исходные неизвестные - $\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)
Построим полученные профили дисперсий скоростей:
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()
In [34]:
# Перебор alpha
alphas = np.arange(0.1, 0.85, 0.15)
# Перебор sig_R_0
sigmas = np.arange(120., 180., 12.)
# Те картинки, на которые стоит обратить особое внимание
good_pics = []
def plot_ranges(sigmas_range, alphas_range, good_pics=[], calc_chi=False, best_err=3):
'''
Для всех предложенных вариантов sigR и alpha
рисует графики исходных и восстановленных дисперсий скоростей los.
Если calc_chi = True, то также считает ошибку по наблюдаемым точкам.
Если ошибка считается, то отмечаются best_err лучших (наименьших) результата.
Синий - для большой оси, красный - малой, зеленый - полусумма.
Изменяет глобальные значения sig_R_0 и alpha!'''
nrows = alphas.size
ncols = sigmas.size
fig, axes = plt.subplots(nrows=nrows, ncols=ncols, sharex=True, sharey=True, figsize=[16,12])
plt_index = 0
# Последнее - среднее геометрическое
sqerr_majs, sqerr_mins, sqerr_mean = [],[],[]
for al in alphas_range:
for si in sigmas_range:
global alpha, sig_R_0
alpha = al
sig_R_0 = si
ax = axes[plt_index/ncols, plt_index % ncols]
ax.set_title(r'$\alpha = %s, \sigma_{R,0}=%s$' % (al,si))
# ax.plot(points, poly_sig_maj(points), '-', color='blue')
ax.plot(points_obs, [sig_maj_exp(Rr, Robs) for (Rr, Robs) in zip(points_real, points_obs)],
'--', color='blue')
# ax.plot(points, poly_sig_min(points), '-', color='red')
ax.plot(points, [sig_min_exp(R) for R in points], '--', color='red')
ax.plot(radii_min, sig_min_p, 's', color='red', ms=1)
ax.plot(radii_maj, sig_maj_p, 's', color='blue', ms=1)
if calc_chi:
# sqerr_maj = sum(power([sig_maj_exp(Robs, p[0]) - p[1] for
# (Robs, p) in zip(r_ma_bind, sig_maj_data)], 2))/len(sig_maj_data)
# sqerr_min = sum(power([sig_min_exp(p[0]) - p[1] for p in sig_min_data], 2))/len(sig_min_data)
sqerr_maj = 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])
ax.text(1, 5, "$\chi^2_{maj}=%5.0f\, \chi^2_{min}=%5.0f$" % (sqerr_maj, sqerr_min), fontsize=12)
sqerr_majs.append(sqerr_maj);sqerr_mins.append(sqerr_min)
sqerr_mean.append(0.5*sqerr_maj+0.5*sqerr_min)
ax.set_ylim(0, 250)
ax.set_xlim(0, 80)
if (plt_index/ncols, plt_index % ncols) in good_pics:
ax.plot([40], [200], 'o', markersize=12., color=(0.2,1.0,0.))
plt_index = plt_index + 1
if calc_chi:
best_maj_err = heapq.nsmallest(best_err, sqerr_majs)
for b_maj in best_maj_err:
b_maj_ind = sqerr_majs.index(b_maj)
ax = axes[b_maj_ind/ncols, b_maj_ind % ncols]
#ax.plot([35], [200], 'o', markersize=12., color='b')
ax.text(70, 200, "%s" % (best_maj_err.index(b_maj)+1), fontsize=12, color='b',
bbox=dict(facecolor='none', edgecolor='b', boxstyle='round'))
best_min_err = heapq.nsmallest(best_err, sqerr_mins)
for b_min in best_min_err:
b_min_ind = sqerr_mins.index(b_min)
ax = axes[b_min_ind/ncols, b_min_ind % ncols]
#ax.plot([30], [200], 'o', markersize=12., color='r')
ax.text(60, 200, "%s" % (best_min_err.index(b_min)+1), fontsize=12, color='r',
bbox=dict(facecolor='none', edgecolor='r', boxstyle='round'))
best_mean_err = heapq.nsmallest(best_err, sqerr_mean)
for b_mean in best_mean_err:
b_mean_ind = sqerr_mean.index(b_mean)
ax = axes[b_mean_ind/ncols, b_mean_ind % ncols]
#ax.plot([25], [200], 's', markersize=12., color='g')
ax.text(50, 200, "%s" % (best_mean_err.index(b_mean)+1), fontsize=12, color='g',
bbox=dict(facecolor='none', edgecolor='g', boxstyle='round'))
plot_ranges(sigmas, alphas, good_pics=good_pics, calc_chi=True)
plt.show()
In [35]:
tex_vkr_dir = 'C:\\Users\\root\\Dropbox\\RotationCurves\\PhD\\VKR\\imgs\\'
tex_imgs_dir = "C:\\Users\\root\\Dropbox\\RotationCurves\\PhD\\paper1\\text\\imgs\\"
In [36]:
os.chdir("C:\\Users\\root\\Dropbox\\RotationCurves\\PhD\\paper1\\text\\imgs")
main_slice = lambda l: sig_min_0/sqrt(sin_i**2 + cos_i**2 * l**2)
alphas = np.arange(0.1, 1.2, 0.03)
sigmas = np.arange(10.0, 400, 3.)
import matplotlib.mlab as mlab
import matplotlib
fig, axes = plt.subplots(nrows=3, ncols=1, sharex=True, sharey=False, figsize=[8,16])
ax = axes[0]
# levels = np.linspace(start=image_min.min(), stop=20., num=5)
# levels = [2., 10., 50., 100.]
# im = ax.imshow(image_min, cmap='jet', vmin=image_min.min(), vmax=20., interpolation='spline16',
# origin="lower", aspect="auto")
# plt.show()
# levels = np.linspace(start=image_min.min()+1., stop=image_min.min()+20., num=5)
levels = np.linspace(start=image_min.min()*1.3, stop=image_min.min()*1.3+20, num=5)
cset=ax.contour(image_min, levels, colors = 'k', origin='lower', extent=[alphas[0],alphas[-1],sigmas[0],sigmas[-1]])
min_map_gutter = cset.collections[0].get_paths()
v1,v2 = min_map_gutter[1].vertices, min_map_gutter[0].vertices
x1,x2 = v1[:,0], v2[:,0]
y1,y2 = v1[:,1], v2[:,1]
plt.clabel(cset, inline=1, fontsize=10, fmt='%1.1f',)
# ax.text(0.87, 260, '$\chi^2_{min}$', size = 24.)
ax.set_ylabel('$\sigma_{R,0}$', size=28.)
xx = np.arange(0.25, 1.0, 0.01)
ax.plot(xx, map(main_slice, xx), '--', color='black')
ax.set_ylim(50, 200)
ax.set_xlim(0.25, 0.99)
ax.fill_between(x1, y1, 0, color='gray', alpha=0.3)
ax.fill_between(x2, y2, 0, color='white')
xlim, ylim = ax.get_xlim(), ax.get_ylim()
ax.text(0.85*(xlim[1]-xlim[0])+xlim[0], 0.85*(ylim[1]-ylim[0]) + ylim[0], '$\chi^2_{min}$', size = 34.)
min_sigmas = np.where(image_min < image_min.min() + 0.03)
slice_alph, slice_sig = min_sigmas[1], min_sigmas[0]
slice_alph = map(lambda l: alphas[0] + (alphas[-1] - alphas[0])*l/len(image_min[0]) , slice_alph)
slice_sig = map(lambda l: sigmas[0] + (sigmas[-1] - sigmas[0])*l/len(image_min), slice_sig)
# ax.plot(slice_alph, slice_sig, '.', color='pink')
poly_slice = poly1d(polyfit(slice_alph, slice_sig, deg=3))
# ax.plot(xx, poly_slice(xx), '.-', color='black')
ax = axes[1]
# levels = np.linspace(start=image_maj.min()-4.3, stop=10., num=10)
# levels = [7., 10., 50., 100.]
levels = np.linspace(start=image_maj.min()+1., stop=image_maj.min()+20., num=6)
cset=ax.contour(image_maj, levels, hold='on', colors = 'k', origin='lower', extent=[alphas[0],alphas[-1],sigmas[0],sigmas[-1]])
plt.clabel(cset, inline=1, fontsize=10, fmt='%1.1f',)
# ax.text(0.87, 260, '$\chi^2_{maj}$', size = 24.)
ax.set_ylabel('$\sigma_{R,0}$', size=28.)
xx = np.arange(0.25, 1.0, 0.01)
ax.plot(xx, map(main_slice, xx), '--', color='black')
ax.fill_between(x1, y1, 0, color='gray', alpha=0.3)
ax.fill_between(x2, y2, 0, color='white')
ax.set_ylim(50, 200)
xlim, ylim = ax.get_xlim(), ax.get_ylim()
ax.text(0.85*(xlim[1]-xlim[0])+xlim[0], 0.85*(ylim[1]-ylim[0]) + ylim[0], '$\chi^2_{maj}$', size = 34.)
ax = axes[2]
err_maj = []
for al in alphas:
global alpha, sig_R_0
alpha = al
sig_R_0 = main_slice(al)
# 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_maj = calc_chi2_normal(sig_maj_p, e_sig_maj_p, [sig_maj_exp(r, r) for r in radii_maj])
err_maj.append(sqerr_maj)
ax.plot(alphas, err_maj, '--', color='black')
err_maj1 = []
for pa in zip(x2,y2):
global alpha, sig_R_0
alpha = pa[0]
sig_R_0 = pa[1]
sqerr_maj = calc_chi2_normal(sig_maj_p, e_sig_maj_p, [sig_maj_exp(r, r) for r in radii_maj])
# 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)])
err_maj1.append(sqerr_maj)
# ax.plot(x2, err_maj1, '-', color='black')
err_maj2 = []
for pa in zip(x1,y1):
global alpha, sig_R_0
alpha = pa[0]
sig_R_0 = pa[1]
sqerr_maj = calc_chi2_normal(sig_maj_p, e_sig_maj_p, [sig_maj_exp(r, r) for r in radii_maj])
# 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)])
err_maj2.append(sqerr_maj)
# ax.plot(x1, err_maj2, '-', color='black')
ax.set_ylabel(r'$\chi^2$', size=28.)
ax.set_xlabel(r'$\alpha$', size=28.)
# ax.fill_between(x1, 0, err_maj2, color='grey', alpha=0.3)
# ax.fill_between(x2, 0, err_maj1, color='white')
# ax.set_ylim(20, 50)
import scipy.interpolate as sp
try:
f1 = sp.interp1d(x2, err_maj1, kind='linear')
ax.fill_between(x1, map(f1, x1), err_maj2, color='grey', alpha=0.3)
except Exception:
f2 = sp.interp1d(x1, err_maj2, kind='linear')
ax.fill_between(x2, map(f2, x2), err_maj1, color='grey', alpha=0.3)
ax.set_ylim(0, 15)
fig.subplots_adjust(hspace=0.)
axes[0].yaxis.get_major_ticks()[0].set_visible(False)
axes[1].yaxis.get_major_ticks()[0].set_visible(False)
ax.set_xlim(0.25, 0.99)
# plt.savefig('ngc3245_maps.eps', format='eps')
plt.savefig(tex_vkr_dir+'ngc3245_maps_large.png', format='png',bbox_inches='tight')
# plt.savefig('ngc3245_maps.pdf', format='pdf', dpi=150)
plt.show()
os.chdir("C:\\science\\2FInstability\\data\\ngc3245")
In [37]:
os.chdir("C:\\Users\\root\\Dropbox\\RotationCurves\\PhD\\paper1\\text\\imgs")
main_slice = lambda l: sig_min_0/sqrt(sin_i**2 + cos_i**2 * l**2)
alphas = np.arange(0.1, 1.2, 0.03)
sigmas = np.arange(10.0, 400, 3.)
import matplotlib.mlab as mlab
import matplotlib
fig, axes = plt.subplots(nrows=3, ncols=1, sharex=True, sharey=False, figsize=[8,16])
ax = axes[0]
# levels = np.linspace(start=image_min.min(), stop=20., num=5)
# levels = [2., 10., 50., 100.]
# im = ax.imshow(image_min, cmap='jet', vmin=image_min.min(), vmax=20., interpolation='spline16',
# origin="lower", aspect="auto")
# plt.show()
# levels = np.linspace(start=image_min.min()+1., stop=image_min.min()+20., num=5)
levels = np.linspace(start=image_min.min()*1.3, stop=image_min.min()*1.3+20, num=5)
cset=ax.contour(image_min, levels, colors = 'k', origin='lower', extent=[alphas[0],alphas[-1],sigmas[0],sigmas[-1]])
min_map_gutter = cset.collections[0].get_paths()
v1,v2 = min_map_gutter[1].vertices, min_map_gutter[0].vertices
x1,x2 = v1[:,0], v2[:,0]
y1,y2 = v1[:,1], v2[:,1]
plt.clabel(cset, inline=1, fontsize=10, fmt='%1.1f',)
# ax.text(0.87, 260, '$\chi^2_{min}$', size = 24.)
ax.set_ylabel('$\sigma_{R,0}$', size=20.)
xx = np.arange(0.25, 1.0, 0.01)
ax.plot(xx, map(main_slice, xx), '--', color='black')
ax.set_ylim(50, 200)
ax.set_xlim(0.25, 0.99)
ax.fill_between(x1, y1, 0, color='gray', alpha=0.3)
ax.fill_between(x2, y2, 0, color='white')
xlim, ylim = ax.get_xlim(), ax.get_ylim()
ax.text(0.85*(xlim[1]-xlim[0])+xlim[0], 0.85*(ylim[1]-ylim[0]) + ylim[0], '$\chi^2_{min}$', size = 24.)
min_sigmas = np.where(image_min < image_min.min() + 0.03)
slice_alph, slice_sig = min_sigmas[1], min_sigmas[0]
slice_alph = map(lambda l: alphas[0] + (alphas[-1] - alphas[0])*l/len(image_min[0]) , slice_alph)
slice_sig = map(lambda l: sigmas[0] + (sigmas[-1] - sigmas[0])*l/len(image_min), slice_sig)
# ax.plot(slice_alph, slice_sig, '.', color='pink')
poly_slice = poly1d(polyfit(slice_alph, slice_sig, deg=3))
# ax.plot(xx, poly_slice(xx), '.-', color='black')
ax = axes[1]
# levels = np.linspace(start=image_maj.min()-4.3, stop=10., num=10)
# levels = [7., 10., 50., 100.]
levels = np.linspace(start=image_maj.min()+1., stop=image_maj.min()+20., num=6)
cset=ax.contour(image_maj, levels, hold='on', colors = 'k', origin='lower', extent=[alphas[0],alphas[-1],sigmas[0],sigmas[-1]])
plt.clabel(cset, inline=1, fontsize=10, fmt='%1.1f',)
# ax.text(0.87, 260, '$\chi^2_{maj}$', size = 24.)
ax.set_ylabel('$\sigma_{R,0}$', size=20.)
xx = np.arange(0.25, 1.0, 0.01)
ax.plot(xx, map(main_slice, xx), '--', color='black')
ax.fill_between(x1, y1, 0, color='gray', alpha=0.3)
ax.fill_between(x2, y2, 0, color='white')
ax.set_ylim(50, 200)
xlim, ylim = ax.get_xlim(), ax.get_ylim()
ax.text(0.85*(xlim[1]-xlim[0])+xlim[0], 0.85*(ylim[1]-ylim[0]) + ylim[0], '$\chi^2_{maj}$', size = 24.)
ax = axes[2]
err_maj = []
for al in alphas:
global alpha, sig_R_0
alpha = al
sig_R_0 = main_slice(al)
# 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_maj = calc_chi2_normal(sig_maj_p, e_sig_maj_p, [sig_maj_exp(r, r) for r in radii_maj])
err_maj.append(sqerr_maj)
ax.plot(alphas, err_maj, '--', color='black')
err_maj1 = []
for pa in zip(x2,y2):
global alpha, sig_R_0
alpha = pa[0]
sig_R_0 = pa[1]
sqerr_maj = calc_chi2_normal(sig_maj_p, e_sig_maj_p, [sig_maj_exp(r, r) for r in radii_maj])
# 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)])
err_maj1.append(sqerr_maj)
# ax.plot(x2, err_maj1, '-', color='black')
err_maj2 = []
for pa in zip(x1,y1):
global alpha, sig_R_0
alpha = pa[0]
sig_R_0 = pa[1]
sqerr_maj = calc_chi2_normal(sig_maj_p, e_sig_maj_p, [sig_maj_exp(r, r) for r in radii_maj])
# 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)])
err_maj2.append(sqerr_maj)
# ax.plot(x1, err_maj2, '-', color='black')
ax.set_ylabel(r'$\chi^2$', size=20.)
ax.set_xlabel(r'$\alpha$', size=20.)
# ax.fill_between(x1, 0, err_maj2, color='grey', alpha=0.3)
# ax.fill_between(x2, 0, err_maj1, color='white')
# ax.set_ylim(20, 50)
import scipy.interpolate as sp
try:
f1 = sp.interp1d(x2, err_maj1, kind='linear')
ax.fill_between(x1, map(f1, x1), err_maj2, color='grey', alpha=0.3)
except Exception:
f2 = sp.interp1d(x1, err_maj2, kind='linear')
ax.fill_between(x2, map(f2, x2), err_maj1, color='grey', alpha=0.3)
ax.set_ylim(0, 15)
fig.subplots_adjust(hspace=0.)
axes[0].yaxis.get_major_ticks()[0].set_visible(False)
axes[1].yaxis.get_major_ticks()[0].set_visible(False)
ax.set_xlim(0.25, 0.99)
plt.savefig(tex_imgs_dir+'ngc3245_maps.eps', format='eps',bbox_inches='tight')
plt.savefig(tex_imgs_dir+'ngc3245_maps.png', format='png',bbox_inches='tight')
plt.savefig(tex_imgs_dir+'ngc3245_maps.pdf', format='pdf', dpi=150,bbox_inches='tight')
plt.show()
os.chdir("C:\\science\\2FInstability\\data\\ngc3245")
In [45]:
os.chdir(tex_imgs_dir)
alphas = np.arange(0.1, 1.2, 0.03)
sigmas = np.arange(10.0, 400, 3.)
fig = plt.figure(figsize=[10,10])
ax = fig.add_subplot(111)
corr_image = (image_min*len(sig_min_p) + image_maj*len(sig_maj_p)) / (len(sig_min_p) + len(sig_maj_p))
# plot_chi2_map(corr_image, ax, log_scale=False, title='$\chi^2$', is_contour=True, vmax=10.)
levels = np.linspace(start=corr_image.min()+0.4, stop=corr_image.min()*1.1+15, num=10)
levels = [3.3, 4.3, 5.1, 6.7, 8.4, 10., 13., 15., 19.]
cset=ax.contour(corr_image, levels, colors = 'k', origin='lower', extent=[alphas[0],alphas[-1],sigmas[0],sigmas[-1]])
min_map_gutter = cset.collections[0].get_paths()
plt.clabel(cset, inline=1, fontsize=10, fmt='%1.1f',)
ax.set_ylim(50, 200)
ax.set_xlim(0.25, 1.)
ax.set_ylabel('$\sigma_{R,0}$', size=22.)
ax.set_xlabel(r'$\alpha$', size=22.)
for tick in ax.yaxis.get_major_ticks():
tick.label.set_fontsize(12)
for tick in ax.xaxis.get_major_ticks():
tick.label.set_fontsize(12)
plt.savefig(tex_imgs_dir+'n3245_pure_chi.eps', format='eps', bbox_inches='tight')
plt.savefig(tex_imgs_dir+'n3245_pure_chi.png', format='png', bbox_inches='tight')
plt.savefig(tex_imgs_dir+'n3245_pure_chi.pdf', format='pdf', dpi=150, bbox_inches='tight')
plt.show()
os.chdir("C:\\science\\2FInstability\\data\\ngc3245")
In [ ]:
In [67]:
import scipy.optimize as opt
def chisqfunc((x_sig, x_alpha)):
global sig_R_0, alpha
sig_R_0 = x_sig
alpha = x_alpha
sqerr_ma = calc_chi2_normal(sig_maj_p, e_sig_maj_p, [sig_maj_exp(l, l) for l in radii_maj])
sqerr_mi = calc_chi2_normal(sig_min_p, e_sig_min_p, [sig_min_exp(l) for l in radii_min])
# print 'chisqf ', sqerr_mi*len(sig_min_p), sqerr_ma*len(sig_maj_p)
chisq = (sqerr_mi*len(sig_min_p) + sqerr_ma*len(sig_maj_p)) / (len(sig_min_p) + len(sig_maj_p) - 2)
return chisq
x0 = np.array([130., 0.7])
res = opt.minimize(chisqfunc, x0, bounds=[(sigmas[0], sigmas[-1]), (alphas[0], alphas[-1])], method='L-BFGS-B')
print res
In [68]:
def gen_next_normal(radii, sig, esig):
randomDelta = np.array([np.random.normal(0., derr/2, 1)[0] for derr in esig] )
randomdataY = sig + randomDelta
return zip(radii, randomdataY)
In [69]:
plt.plot(radii_maj, sig_maj_p, 's', label='$\sigma_{los}^{maj}$', color='blue')
plt.errorbar(radii_maj, sig_maj_p, yerr=e_sig_maj_p, fmt='o', marker='.', color='blue')
plt.plot(points, poly_sig_maj(points), label = '$\sigma_{los}^{maj}\, splinefit$', color='blue')
for i in range(3):
r, s = zip(*gen_next_normal(radii_maj, sig_maj_p, e_sig_maj_p))
plt.plot(r, s, 's', color='red')
plt.ylim(0., 250.)
plt.legend()
plt.show()
In [70]:
import time
# os.chdir("C:\\science\\2FInstability\\data\\ngc1068")
N = 10000
result = []
start_time = time.time()
fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, figsize=(16, 16))
radii_maj1, sig_maj_p1, e_sig_maj_p1 = radii_maj, sig_maj_p, e_sig_maj_p
radii_min1, sig_min_p1, e_sig_min_p1 = radii_min, sig_min_p, e_sig_min_p
for i in log_progress(range(N)):
global spl_maj, spl_min
global radii_min, radii_maj, sig_min_p, sig_maj_p, sig_min_0
r, s = zip(*gen_next_normal(radii_maj1, sig_maj_p1, e_sig_maj_p1))
spl_maj = inter.UnivariateSpline(r, s, k=4, s=10000.)
radii_maj, sig_maj_p = r, s
ax1.plot(points, spl_maj(points), label = '$\sigma_{los}^{maj}\, splinefit$', color='blue')
r, s = zip(*gen_next_normal(radii_min1, sig_min_p1, e_sig_min_p1))
spl_min = inter.UnivariateSpline(r, s, k=4, s=10000.)
sig_min_0 = spl_min(radii_min[0])
radii_min, sig_min_p = r, s
ax2.plot(points, spl_min(points), label = '$\sigma_{los}^{maj}\, splinefit$', color='red')
# res = opt.minimize(chisqfunc, x0, bounds=[(sigmas[0], sigmas[-1]), (alphas[0], alphas[-1])], method='L-BFGS-B')
res = opt.minimize(chisqfunc, x0, bounds=[(sigmas[0], sigmas[-1]), (0, alphas[-1])], method='L-BFGS-B')
# print res.fun
result.append(res.x)
# np.save(pics_path + 'monte_carlo', np.array(result))
print("--- %s seconds ---" % (time.time() - start_time))
ax1.errorbar(radii_maj, sig_maj_p, yerr=e_sig_maj_p, fmt='o', marker='.', color='red')
ax2.errorbar(radii_min, sig_min_p, yerr=e_sig_min_p, fmt='o', marker='.', color='blue')
ax1.set_ylim(0., 200.)
ax2.set_ylim(0., 200.)
plt.show()
radii_maj, sig_maj_p, e_sig_maj_p = radii_maj1, sig_maj_p1, e_sig_maj_p1
radii_min, sig_min_p, e_sig_min_p = radii_min1, sig_min_p1, e_sig_min_p1
In [71]:
s,a = zip(*result)
plt.plot(a, s, '.')
plt.plot(alphas, map(main_slice, alphas), '--')
plt.xlim(0.0, 0.99)
plt.ylim(100, 200)
plt.show()
In [72]:
from scipy.stats import norm
fig = plt.figure(figsize=(12,6))
ax = fig.add_subplot(111)
n, bins, patches = ax.hist(s, 20, normed=1, facecolor='green', alpha=0.75)
mu, std = norm.fit(s)
xmin, xmax = plt.xlim()
x = np.linspace(xmin, xmax, 100)
p = norm.pdf(x, mu, std)
plt.plot(x, p, 'k', linewidth=2)
ax.set_title('$\mu=%s,\ \sigma=%s$' % (mu, std), fontsize=18)
ax.grid(True)
plt.show()
In [73]:
from scipy.stats import norm
fig = plt.figure(figsize=(12,6))
ax = fig.add_subplot(111)
n, bins, patches = ax.hist(a, 20, normed=1, facecolor='green', alpha=0.75)
mu, std = norm.fit(a)
xmin, xmax = plt.xlim()
x = np.linspace(xmin, xmax, 100)
p = norm.pdf(x, mu, std)
plt.plot(x, p, 'k', linewidth=2)
ax.set_title('$\mu=%s,\ \sigma=%s$' % (mu, std), fontsize=18)
ax.grid(True)
plt.show()
In [ ]:
In [39]:
from scipy.optimize import curve_fit
sig_maj_data = zip(radii_maj, sig_maj_p, e_sig_maj_p)
sig_maj_data = filter(lambda l: l[0] > 20., 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] > 20., sig_min_data)
radii_min, sig_min_p, e_sig_min_p = zip(*sig_min_data)
def fuu(x, A, B):
return A*x + B
A,B = curve_fit(fuu, radii_min, map(np.log, sig_min_p))[0]
print A, B
In [40]:
def fuu1(x, A, B):
return np.exp(A*x+B)
In [41]:
h_kin = float("{0:.2f}".format(-1./A))
print h_kin
plt.plot(points, fuu1(points, A, B), 'x', color='red')
plt.plot(points, map(lambda l: np.exp(B)*np.exp(-l/h_kin), points), '+')
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, spl_maj(points), label = '$\sigma_{los}^{maj}\, splinefit$', color='blue')
plt.plot(radii_min, sig_min_p, 's', label='$\sigma_{los}^{min}$', color='red')
plt.errorbar(radii_min, sig_min_p, yerr=e_sig_min_p, fmt='.', marker='.', mew=0, color='red')
plt.plot(points, spl_min(points), label = '$\sigma_{los}^{min}\, splinefit$', color='red')
# plt.axvline(x=radii_min[14], color='black')
plt.legend()
plt.ylim(0, 250)
plt.show()
In [42]:
sig_min_0 = fuu1(radii_min[0], A, B)
def sigR_ger_exp(R):
return sig_R_0*exp(-R/h_kin)/exp(-radii_min[0]/h_kin)
def sigZ_ger_exp(R):
return sigR_ger_exp(R)*alpha
def sig_maj_exp(R, R1):
return sqrt(sigPhi_to_sigR_real(R) * sigR_ger_exp(R)**2 * sin_i**2 + sigZ_ger_exp(R)**2 * cos_i**2)
def sig_min_exp(R):
# return sig_R_0*exp(-R/h_kin)*sqrt(sin_i**2 + alpha**2 * cos_i**2)
return sqrt(sigR_ger_exp(R)**2 * sin_i**2 + sigZ_ger_exp(R)**2 * cos_i**2)
In [43]:
# sig_min_0 = fuu1(radii_min[0], A, B)/np.exp(radii_min[0]*A)
sig_min_0 = fuu1(radii_min[0], A, B)
print sig_min_0
sig_R_0 = 118.
alpha = 0.6
plt.plot(points, spl_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, spl_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()
In [44]:
alphas = np.arange(0.1, 1.2, 0.03)
sigmas = np.arange(10.0, 400, 3.)
radii_maj1, sig_maj_p1, e_sig_maj_p1 = radii_maj, sig_maj_p, e_sig_maj_p
radii_min1, sig_min_p1, e_sig_min_p1 = radii_min, sig_min_p, e_sig_min_p
image, image_maj, image_min = compute_chi2_maps(alphas=alphas, sigmas=sigmas)
radii_maj, sig_maj_p, e_sig_maj_p = radii_maj1, sig_maj_p1, e_sig_maj_p1
radii_min, sig_min_p, e_sig_min_p = radii_min1, sig_min_p1, e_sig_min_p1
In [45]:
alphas = np.arange(0.1, 1.2, 0.03)
sigmas = np.arange(10.0, 400, 3.)
fig, axes = plt.subplots(nrows=3, ncols=1, sharex=False, sharey=True, figsize=[12,24])
plot_chi2_map(image_maj, axes[0], log_scale=False, title='$\chi^2_{maj}$', is_contour=False, vmax=20.)
plot_chi2_map(image_min, axes[1], log_scale=False, title='$\chi^2_{min}$', is_contour=False, vmax=20.)
corr_image = (image_min*len(sig_min_p) + image_maj*len(sig_maj_p)) / (len(sig_min_p) + len(sig_maj_p))
print 'N1_maj={},\t N2_min={},\t chi^2_corr[0][0]={} (was {} and {})'.format(len(sig_maj_p), len(sig_min_p), corr_image[0][0],
image_min[0][0], image_maj[0][0])
plot_chi2_map(corr_image, axes[2], log_scale=False, title='$\chi^2$', is_contour=True, vmax=20.)
plt.show()
In [46]:
alphas = np.arange(0.3, 0.9, 0.11)
sigmas = np.arange(86., 166., 16.)
plot_ranges(sigmas, alphas, good_pics=good_pics, calc_chi=True)
plt.show()
In [47]:
alphas = np.arange(0.1, 1.2, 0.03)
sigmas = np.arange(10.0, 400, 3.)
import matplotlib.mlab as mlab
import matplotlib
fig, axes = plt.subplots(nrows=3, ncols=1, sharex=True, sharey=False, figsize=[8,16])
ax = axes[0]
levels = np.linspace(start=image_min.min()*1.1, stop=image_min.min()*1.1+4, num=5)
levels = np.linspace(start=image_min.min()*1.2, stop=image_min.min()*1.1+14, num=5)
cset=ax.contour(image_min, levels, colors = 'k', origin='lower', extent=[alphas[0],alphas[-1],sigmas[0],sigmas[-1]])
min_map_gutter = cset.collections[0].get_paths()
v1,v2 = min_map_gutter[1].vertices, min_map_gutter[0].vertices
x1,x2 = v1[:,0], v2[:,0]
y1,y2 = v1[:,1], v2[:,1]
plt.clabel(cset, inline=1, fontsize=10, fmt='%1.1f',)
# ax.text(0.87, 172, '$\chi^2_{min}$', size = 24.)
ax.text(0.87, 130, '$\chi^2_{min}$', size = 24.)
ax.set_ylabel('$\sigma_{R,0}$', size=20.)
xx = np.arange(0.25, 1.0, 0.01)
ax.plot(xx, map(main_slice, xx), '--', color='black')
ax.set_ylim(0, 150)
ax.set_ylim(40, 200)
min_sigmas = np.where(image_min < image_min.min() + 0.03)
slice_alph, slice_sig = min_sigmas[1], min_sigmas[0]
slice_alph = map(lambda l: alphas[0] + (alphas[-1] - alphas[0])*l/len(image_min[0]) , slice_alph)
slice_sig = map(lambda l: sigmas[0] + (sigmas[-1] - sigmas[0])*l/len(image_min), slice_sig)
# ax.plot(slice_alph, slice_sig, '.', color='pink')
poly_slice = poly1d(polyfit(slice_alph, slice_sig, deg=3))
# ax.plot(xx, poly_slice(xx), '.-', color='black')
ax.fill_between(x1, y1, 0, color='gray', alpha=0.3)
ax.fill_between(x2, y2, 0, color='white')
ax = axes[1]
# levels = np.append(np.linspace(start=image_maj.min()+0.1, stop=image_maj.min()+4.1, num=6), np.array([image_maj.min()+0.25]))
levels = np.linspace(start=image_maj.min()*1.1, stop=image_maj.min()*1.1+5, num=6)
# levels = [image_maj.min()+0.1, image_maj.min()+0.25, image_maj.min()+1.1, image_maj.min()+2.1, image_maj.min()+3.1,
# image_maj.min()+4.1]
levels = sorted(levels)
cset=ax.contour(image_maj, levels, hold='on', colors = 'k', origin='lower', extent=[alphas[0],alphas[-1],sigmas[0],sigmas[-1]])
plt.clabel(cset, inline=1, fontsize=10, fmt='%1.1f',)
# ax.text(0.87, 172, '$\chi^2_{maj}$', size = 24.)
ax.text(0.87, 110, '$\chi^2_{maj}$', size = 24.)
ax.set_ylabel('$\sigma_{R,0}$', size=20.)
xx = np.arange(0.25, 1.0, 0.01)
ax.plot(xx, map(main_slice, xx), '--', color='black')
ax.fill_between(x1, y1, 0, color='gray', alpha=0.3)
ax.fill_between(x2, y2, 0, color='white')
ax.set_ylim(0, 150)
ax.set_ylim(40, 200)
ax = axes[2]
err_maj = []
for al in alphas:
global alpha, sig_R_0
alpha = al
sig_R_0 = main_slice(al)
sqerr_maj = calc_chi2_normal(sig_maj_p, e_sig_maj_p, [sig_maj_exp(r, r) for r in radii_maj])
err_maj.append(sqerr_maj)
ax.plot(alphas, err_maj, '--', color='black')
err_maj1 = []
for pa in zip(x2,y2):
global alpha, sig_R_0
alpha = pa[0]
sig_R_0 = pa[1]
sqerr_maj = calc_chi2_normal(sig_maj_p, e_sig_maj_p, [sig_maj_exp(r, r) for r in radii_maj])
err_maj1.append(sqerr_maj)
# ax.plot(x2, err_maj1, '-', color='black')
err_maj2 = []
for pa in zip(x1,y1):
global alpha, sig_R_0
alpha = pa[0]
sig_R_0 = pa[1]
sqerr_maj = calc_chi2_normal(sig_maj_p, e_sig_maj_p, [sig_maj_exp(r, r) for r in radii_maj])
err_maj2.append(sqerr_maj)
# ax.plot(x1, err_maj2, '-', color='black')
ax.set_ylabel(r'$\chi^2$', size=20.)
ax.set_xlabel(r'$\alpha$', size=20.)
import scipy.interpolate as sp
try:
f1 = sp.interp1d(x2, err_maj1, kind='linear')
ax.fill_between(x1, map(f1, x1), err_maj2, color='grey', alpha=0.3)
except Exception:
f2 = sp.interp1d(x1, err_maj2, kind='linear')
ax.fill_between(x2, map(f2, x2), err_maj1, color='grey', alpha=0.3)
ax.set_ylabel(r'$\chi^2$', size=20.)
ax.set_xlabel(r'$\alpha$', size=20.)
# ax.set_ylim(0.5, 3.)
fig.subplots_adjust(hspace=0.)
axes[0].yaxis.get_major_ticks()[0].set_visible(False)
axes[1].yaxis.get_major_ticks()[0].set_visible(False)
ax.set_xlim(0.25, 0.99)
plt.show()
In [48]:
import time
N = 100
result = []
start_time = time.time()
fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, figsize=(16, 16))
radii_maj2, sig_maj_p2, e_sig_maj_p2 = radii_maj, sig_maj_p, e_sig_maj_p
radii_min2, sig_min_p2, e_sig_min_p2 = radii_min, sig_min_p, e_sig_min_p
for i in log_progress(range(N)):
global radii_min, radii_maj, sig_min_p, sig_maj_p
global A,B,h_kin
r, s = zip(*gen_next_normal(radii_maj2, sig_maj_p2, e_sig_maj_p2))
radii_maj, sig_maj_p = r, s
ax1.plot(points, spl_maj(points), label = '$\sigma_{los}^{maj}\, splinefit$', color='blue')
r, s = zip(*gen_next_normal(radii_min2, sig_min_p2, e_sig_min_p2))
A, B = curve_fit(fuu, r, map(np.log, s))[0]
# print A,B
h_kin = float("{0:.2f}".format(-1./A))
# print h_kin
radii_min, sig_min_p = r, s
ax2.plot(points, map(lambda l: np.exp(B)*np.exp(-l/h_kin), points), '-', color='red')
ax2.plot(points, fuu1(points, A, B), 'x', color='red')
res = opt.minimize(chisqfunc, x0, bounds=[(sigmas[0], sigmas[-1]), (alphas[0], alphas[-1])], method='L-BFGS-B')
result.append(res.x)
print("--- %s seconds ---" % (time.time() - start_time))
ax1.errorbar(radii_maj, sig_maj_p, yerr=e_sig_maj_p, fmt='o', marker='.', color='red')
ax2.errorbar(radii_min, sig_min_p, yerr=e_sig_min_p, fmt='o', marker='.', color='blue')
ax1.set_ylim(0., 210.)
ax2.set_ylim(0., 250.)
plt.show()
radii_maj, sig_maj_p, e_sig_maj_p = radii_maj1, sig_maj_p1, e_sig_maj_p1
radii_min, sig_min_p, e_sig_min_p = radii_min1, sig_min_p1, e_sig_min_p1
In [49]:
s,a = zip(*result)
plt.plot(a, s, '.')
plt.plot(alphas, map(main_slice, alphas), '--')
# plt.xlim(0.0, 0.99)
plt.ylim(0, 400)
plt.show()
In [ ]:
In [ ]:
In [50]:
from scipy.optimize import curve_fit
sig_maj_data = zip(radii_maj, sig_maj_p, e_sig_maj_p)
sig_maj_data = filter(lambda l: l[0] > 30. and l[0] < 45, 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] > 30. and l[0] < 45, sig_min_data)
radii_min, sig_min_p, e_sig_min_p = zip(*sig_min_data)
def fuu(x, A, B):
return A*x + B
A,B = curve_fit(fuu, radii_min, map(np.log, sig_min_p))[0]
print A, B
In [51]:
def fuu1(x, A, B):
return np.exp(A*x+B)
In [52]:
h_kin = float("{0:.2f}".format(-1./A))
print h_kin
plt.plot(points, fuu1(points, A, B), 'x', color='red')
plt.plot(points, map(lambda l: np.exp(B)*np.exp(-l/h_kin), points), '+')
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, spl_maj(points), label = '$\sigma_{los}^{maj}\, splinefit$', color='blue')
plt.plot(radii_min, sig_min_p, 's', label='$\sigma_{los}^{min}$', color='red')
plt.errorbar(radii_min, sig_min_p, yerr=e_sig_min_p, fmt='.', marker='.', mew=0, color='red')
plt.plot(points, spl_min(points), label = '$\sigma_{los}^{min}\, splinefit$', color='red')
# plt.axvline(x=radii_min[14], color='black')
plt.legend()
plt.ylim(0, 250)
plt.show()
In [53]:
sig_min_0 = fuu1(radii_min[0], A, B)
def sigR_ger_exp(R):
return sig_R_0*exp(-R/h_kin)/exp(-radii_min[0]/h_kin)
def sigZ_ger_exp(R):
return sigR_ger_exp(R)*alpha
def sig_maj_exp(R, R1):
return sqrt(sigPhi_to_sigR_real(R) * sigR_ger_exp(R)**2 * sin_i**2 + sigZ_ger_exp(R)**2 * cos_i**2)
def sig_min_exp(R):
# return sig_R_0*exp(-R/h_kin)*sqrt(sin_i**2 + alpha**2 * cos_i**2)
return sqrt(sigR_ger_exp(R)**2 * sin_i**2 + sigZ_ger_exp(R)**2 * cos_i**2)
In [54]:
sig_min_0 = fuu1(radii_min[0], A, B)
sig_R_0 = 118.
alpha = 0.6
plt.plot(points, spl_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, spl_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()
In [55]:
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)
alpha = 0.9
sig_R_0 = 90.
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()
In [56]:
alphas = np.arange(0.1, 1.2, 0.03)
sigmas = np.arange(10.0, 400, 3.)
radii_maj1, sig_maj_p1, e_sig_maj_p1 = radii_maj, sig_maj_p, e_sig_maj_p
radii_min1, sig_min_p1, e_sig_min_p1 = radii_min, sig_min_p, e_sig_min_p
image, image_maj, image_min = compute_chi2_maps(alphas=alphas, sigmas=sigmas)
radii_maj, sig_maj_p, e_sig_maj_p = radii_maj1, sig_maj_p1, e_sig_maj_p1
radii_min, sig_min_p, e_sig_min_p = radii_min1, sig_min_p1, e_sig_min_p1
In [57]:
alphas = np.arange(0.1, 1.2, 0.03)
sigmas = np.arange(10.0, 400, 3.)
fig, axes = plt.subplots(nrows=3, ncols=1, sharex=False, sharey=True, figsize=[12,24])
plot_chi2_map(image_maj, axes[0], log_scale=False, title='$\chi^2_{maj}$', is_contour=False, vmax=20.)
plot_chi2_map(image_min, axes[1], log_scale=False, title='$\chi^2_{min}$', is_contour=False, vmax=20.)
corr_image = (image_min*len(sig_min_p) + image_maj*len(sig_maj_p)) / (len(sig_min_p) + len(sig_maj_p))
print 'N1_maj={},\t N2_min={},\t chi^2_corr[0][0]={} (was {} and {})'.format(len(sig_maj_p), len(sig_min_p), corr_image[0][0],
image_min[0][0], image_maj[0][0])
plot_chi2_map(corr_image, axes[2], log_scale=False, title='$\chi^2$', is_contour=True, vmax=20.)
plt.show()
In [58]:
alphas = np.arange(0.3, 0.9, 0.11)
sigmas = np.arange(86., 166., 16.)
plot_ranges(sigmas, alphas, good_pics=good_pics, calc_chi=True)
plt.show()
In [59]:
alphas = np.arange(0.1, 1.2, 0.03)
sigmas = np.arange(10.0, 400, 3.)
import matplotlib.mlab as mlab
import matplotlib
fig, axes = plt.subplots(nrows=3, ncols=1, sharex=True, sharey=False, figsize=[8,16])
ax = axes[0]
levels = np.linspace(start=image_min.min()*1.1, stop=image_min.min()*1.1+4, num=5)
levels = np.linspace(start=image_min.min()*1.2, stop=image_min.min()*1.1+14, num=5)
cset=ax.contour(image_min, levels, colors = 'k', origin='lower', extent=[alphas[0],alphas[-1],sigmas[0],sigmas[-1]])
min_map_gutter = cset.collections[0].get_paths()
v1,v2 = min_map_gutter[1].vertices, min_map_gutter[0].vertices
x1,x2 = v1[:,0], v2[:,0]
y1,y2 = v1[:,1], v2[:,1]
plt.clabel(cset, inline=1, fontsize=10, fmt='%1.1f',)
# ax.text(0.87, 172, '$\chi^2_{min}$', size = 24.)
ax.text(0.87, 130, '$\chi^2_{min}$', size = 24.)
ax.set_ylabel('$\sigma_{R,0}$', size=20.)
xx = np.arange(0.25, 1.0, 0.01)
ax.plot(xx, map(main_slice, xx), '--', color='black')
ax.set_ylim(0, 150)
ax.set_ylim(40, 200)
min_sigmas = np.where(image_min < image_min.min() + 0.03)
slice_alph, slice_sig = min_sigmas[1], min_sigmas[0]
slice_alph = map(lambda l: alphas[0] + (alphas[-1] - alphas[0])*l/len(image_min[0]) , slice_alph)
slice_sig = map(lambda l: sigmas[0] + (sigmas[-1] - sigmas[0])*l/len(image_min), slice_sig)
# ax.plot(slice_alph, slice_sig, '.', color='pink')
poly_slice = poly1d(polyfit(slice_alph, slice_sig, deg=3))
# ax.plot(xx, poly_slice(xx), '.-', color='black')
ax.fill_between(x1, y1, 0, color='gray', alpha=0.3)
ax.fill_between(x2, y2, 0, color='white')
ax = axes[1]
# levels = np.append(np.linspace(start=image_maj.min()+0.1, stop=image_maj.min()+4.1, num=6), np.array([image_maj.min()+0.25]))
levels = np.linspace(start=image_maj.min()*1.1, stop=image_maj.min()*1.1+5, num=6)
# levels = [image_maj.min()+0.1, image_maj.min()+0.25, image_maj.min()+1.1, image_maj.min()+2.1, image_maj.min()+3.1,
# image_maj.min()+4.1]
levels = sorted(levels)
cset=ax.contour(image_maj, levels, hold='on', colors = 'k', origin='lower', extent=[alphas[0],alphas[-1],sigmas[0],sigmas[-1]])
plt.clabel(cset, inline=1, fontsize=10, fmt='%1.1f',)
# ax.text(0.87, 172, '$\chi^2_{maj}$', size = 24.)
ax.text(0.87, 110, '$\chi^2_{maj}$', size = 24.)
ax.set_ylabel('$\sigma_{R,0}$', size=20.)
xx = np.arange(0.25, 1.0, 0.01)
ax.plot(xx, map(main_slice, xx), '--', color='black')
ax.fill_between(x1, y1, 0, color='gray', alpha=0.3)
ax.fill_between(x2, y2, 0, color='white')
ax.set_ylim(0, 150)
ax.set_ylim(40, 200)
ax = axes[2]
err_maj = []
for al in alphas:
global alpha, sig_R_0
alpha = al
sig_R_0 = main_slice(al)
sqerr_maj = calc_chi2_normal(sig_maj_p, e_sig_maj_p, [sig_maj_exp(r, r) for r in radii_maj])
err_maj.append(sqerr_maj)
ax.plot(alphas, err_maj, '--', color='black')
err_maj1 = []
for pa in zip(x2,y2):
global alpha, sig_R_0
alpha = pa[0]
sig_R_0 = pa[1]
sqerr_maj = calc_chi2_normal(sig_maj_p, e_sig_maj_p, [sig_maj_exp(r, r) for r in radii_maj])
err_maj1.append(sqerr_maj)
# ax.plot(x2, err_maj1, '-', color='black')
err_maj2 = []
for pa in zip(x1,y1):
global alpha, sig_R_0
alpha = pa[0]
sig_R_0 = pa[1]
sqerr_maj = calc_chi2_normal(sig_maj_p, e_sig_maj_p, [sig_maj_exp(r, r) for r in radii_maj])
err_maj2.append(sqerr_maj)
# ax.plot(x1, err_maj2, '-', color='black')
ax.set_ylabel(r'$\chi^2$', size=20.)
ax.set_xlabel(r'$\alpha$', size=20.)
import scipy.interpolate as sp
try:
f1 = sp.interp1d(x2, err_maj1, kind='linear')
ax.fill_between(x1, map(f1, x1), err_maj2, color='grey', alpha=0.3)
except Exception:
f2 = sp.interp1d(x1, err_maj2, kind='linear')
ax.fill_between(x2, map(f2, x2), err_maj1, color='grey', alpha=0.3)
ax.set_ylabel(r'$\chi^2$', size=20.)
ax.set_xlabel(r'$\alpha$', size=20.)
# ax.set_ylim(0.5, 3.)
fig.subplots_adjust(hspace=0.)
axes[0].yaxis.get_major_ticks()[0].set_visible(False)
axes[1].yaxis.get_major_ticks()[0].set_visible(False)
ax.set_xlim(0.25, 0.99)
plt.show()
In [60]:
import time
N = 100
result = []
start_time = time.time()
fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, figsize=(16, 16))
radii_maj2, sig_maj_p2, e_sig_maj_p2 = radii_maj, sig_maj_p, e_sig_maj_p
radii_min2, sig_min_p2, e_sig_min_p2 = radii_min, sig_min_p, e_sig_min_p
for i in log_progress(range(N)):
global radii_min, radii_maj, sig_min_p, sig_maj_p
global A,B,h_kin
r, s = zip(*gen_next_normal(radii_maj2, sig_maj_p2, e_sig_maj_p2))
radii_maj, sig_maj_p = r, s
ax1.plot(points, spl_maj(points), label = '$\sigma_{los}^{maj}\, splinefit$', color='blue')
r, s = zip(*gen_next_normal(radii_min2, sig_min_p2, e_sig_min_p2))
A, B = curve_fit(fuu, r, map(np.log, s))[0]
# print A,B
h_kin = float("{0:.2f}".format(-1./A))
# print h_kin
radii_min, sig_min_p = r, s
ax2.plot(points, map(lambda l: np.exp(B)*np.exp(-l/h_kin), points), '-', color='red')
ax2.plot(points, fuu1(points, A, B), 'x', color='red')
res = opt.minimize(chisqfunc, x0, bounds=[(sigmas[0], sigmas[-1]), (alphas[0], alphas[-1])], method='L-BFGS-B')
result.append(res.x)
print("--- %s seconds ---" % (time.time() - start_time))
ax1.errorbar(radii_maj, sig_maj_p, yerr=e_sig_maj_p, fmt='o', marker='.', color='red')
ax2.errorbar(radii_min, sig_min_p, yerr=e_sig_min_p, fmt='o', marker='.', color='blue')
ax1.set_ylim(0., 210.)
ax2.set_ylim(0., 250.)
plt.show()
radii_maj, sig_maj_p, e_sig_maj_p = radii_maj1, sig_maj_p1, e_sig_maj_p1
radii_min, sig_min_p, e_sig_min_p = radii_min1, sig_min_p1, e_sig_min_p1
In [61]:
s,a = zip(*result)
plt.plot(a, s, '.')
plt.plot(alphas, map(main_slice, alphas), '--')
# plt.xlim(0.0, 0.99)
plt.ylim(0, 400)
plt.show()
In [ ]:
In [ ]:
In [ ]:
In [ ]:
Как мы помним из карт выше - тут есть вырождение. Попробуем применить герсеновскую методику, чтобы понять, есть ли вырождение у них. Для этого возьмем прямой участок и подгоним экспонентами: $$\sigma_R = \sigma_{R,0}e^{-R/h}$$ $$\sigma_Z = \sigma_{Z,0}e^{-R/h}$$ Подставим в уравнение для малой оси и возьмем логарифм: $$\log{\sigma_{los, min}} = \log{\sigma_{R, 0}} - R/h + 0.5\log{(\sin^2i + \sigma_{z,0}^2/\sigma_{R,0}^2 \cos^2i)}$$ И теперь если подогнать прямой $Y = AX + B$ то $A=-1/h$, $B = \log{\sigma_{R, 0}} + 0.5\log{(\sin^2i + \sigma_{z,0}^2/\sigma_{R,0}^2 \cos^2i)}$, то видно, где есть место для вырождения, можем нарисовать карту. Точки обрежем по 21".
In [62]:
plt.plot(points, [log(poly_sig_min(R)) for R in points], '.', color='m', label='$\sigma_{R, exp}\, polyfit$')
plt.plot(radii_min, map(np.log, sig_min_p), 's', label='$\sigma_{los}^{min}$', color='red')
poly_exp_line = poly1d(polyfit(radii_min[18:], map(np.log, sig_min_p[18:]), deg=1))
plt.plot(points, poly_exp_line(points), '-', label='%s' % poly_exp_line)
plt.axvline(x=21.)
plt.legend()
plt.show()
In [ ]:
plt.plot(radii_min, sig_min_p, 's', label='$\sigma_{los}^{min}$', color='red')
plt.plot(radii_maj, sig_maj_p, 's', label='$\sigma_{los}^{maj}$', color='blue')
plt.plot(points, [poly_sig_min(R) for R in points], '-', color='m', label='$\sigma_{R, exp} polyfit$')
plt.plot(points, map(np.exp, poly_exp_line(points)), '-')
plt.show()
In [ ]:
os.chdir("C:\\science\\2FInstability\\data\\ngc3245")
h_kin = -poly_exp_line[1]
alphas = np.arange(0.2, 0.9, 0.03)
sigmas = np.arange(10.0, 200, 3.)
sig_maj_gers_exp = lambda l: sig_R_0*exp(-l*h_kin)*sqrt(sigPhi_to_sigR_real(l) * sin_i**2 + alpha **2 * cos_i**2)
sig_min_gers_exp = lambda l: sig_R_0*exp(-l*h_kin)*sqrt(sin_i**2 + alpha **2 * cos_i**2)
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)))
global alpha, sig_R_0
cut_data = lambda l: l[0] > 21.
ma_data = filter(cut_data, zip(radii_maj, sig_maj_p))
mi_data = filter(cut_data, zip(radii_min, sig_min_p))
for i,si in enumerate(sigmas):
for j,al in enumerate(alphas):
alpha = al
sig_R_0 = si
sqerr_maj = sum(power([sig_maj_gers_exp(p[0]) - p[1] for p in ma_data], 2))/len(ma_data)
sqerr_min = sum(power([sig_min_gers_exp(p[0]) - p[1] for p in mi_data], 2))/len(mi_data)
image_maj[i][j] = sqerr_maj
image_min[i][j] = sqerr_min
image[i][j] = 0.5*(sqerr_maj+sqerr_min)
return image, image_maj, image_min
pics_path = '.\\pics\\'
if not os.path.exists(pics_path):
os.makedirs(pics_path)
if os.path.isfile(pics_path + 'gers_chi2_map_maj.npy'):
image_maj = np.load(pics_path + "gers_chi2_map_maj.npy")
image_min = np.load(pics_path + "gers_chi2_map_min.npy")
image = np.load(pics_path + "gers_chi2_map.npy")
else:
image, image_maj, image_min = compute_chi2_maps(alphas=alphas, sigmas=sigmas)
np.save(pics_path + 'gers_chi2_map_maj', image_maj)
np.save(pics_path + 'gers_chi2_map_min', image_min)
np.save(pics_path + 'gers_chi2_map', image)
Связь, по которой по идее должно идти вырождение, описывается как $\sigma_{R,0}=\frac{e^B}{\sqrt{\sin^2i + \alpha^2\cos^2i}}$
In [ ]:
def plot_chi2_map(image, ax, log_scale=False, title='$\chi^2$', is_contour=False, vmax=0.):
if image is not None:
image_log = image
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)
degenereation_slice = lambda l: exp(poly_exp_line[0])/ sqrt(sin_i**2 + l**2 * cos_i**2)
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_{mean}$', is_contour=False, vmax=300.)
plot_chi2_map(image_maj, axes[1], log_scale=False, title='$\chi^2_{maj}$', is_contour=False, vmax=300.)
plot_chi2_map(image_min, axes[2], log_scale=False, title='$\chi^2_{min}$', is_contour=False, vmax=300.)
axes[0].plot(xx, map(degenereation_slice, xx), '-', color='m')
axes[1].plot(xx, map(degenereation_slice, xx), '-', color='m')
axes[2].plot(xx, map(degenereation_slice, xx), '-', color='m')
plt.show()
In [ ]:
plt.figure(figsize=(8, 8))
for al in alphas:
global alpha, sig_R_0
cut_data = lambda l: l[0] > 21.
ma_data = filter(cut_data, zip(radii_maj, sig_maj_p))
mi_data = filter(cut_data, zip(radii_min, sig_min_p))
alpha = al
sig_R_0 = degenereation_slice(alpha)
sqerr_maj = sum(power([sig_maj_gers_exp(p[0]) - p[1] for p in ma_data], 2))/len(ma_data)
sqerr_min = sum(power([sig_min_gers_exp(p[0]) - p[1] for p in mi_data], 2))/len(mi_data)
plt.plot(alpha, sqerr_maj, 'x', color='blue')
plt.plot(alpha, sqerr_min, 'x', color='red')
plt.show()
In [ ]:
alphas = np.arange(0.25, 0.85, 0.1)
sigmas = np.arange(90., 150., 8.)
good_pics = []
def plot_ranges_gers(sigmas_range, alphas_range, good_pics=[], calc_chi=False, best_err=3):
nrows = alphas.size
ncols = sigmas.size
fig, axes = plt.subplots(nrows=nrows, ncols=ncols, sharex=True, sharey=True, figsize=[20,12])
plt_index = 0
# Последнее - среднее геометрическое
sqerr_majs, sqerr_mins, sqerr_mean = [],[],[]
for al in alphas_range:
for si in sigmas_range:
global alpha, sig_R_0
alpha = al
sig_R_0 = si
ax = axes[plt_index/ncols, plt_index % ncols]
ax.set_title(r'$\alpha = %s, \sigma_{R,0}=%s$' % (al,si))
ax.plot(points, poly_sig_maj(points), '-', color='blue')
ax.plot(points, [sig_maj_gers_exp(R) for R in points], '--', color='blue')
ax.plot(points, poly_sig_min(points), '-', color='red')
ax.plot(points, [sig_min_gers_exp(R) for R in points], '--', color='red')
ax.plot(radii_min, sig_min_p, 's', color='red', ms=1)
ax.plot(radii_maj, sig_maj_p, 's', color='blue', ms=1)
ax.set_ylim(0, 170)
ax.set_xlim(15, 70)
plt_index = plt_index + 1
plot_ranges_gers(sigmas, alphas, good_pics=good_pics, calc_chi=True)
plt.show()
Проверим в целях найти минимум еще две идеи - построим карту по честному отношению и по степенной кривой вращения:
In [ ]:
os.chdir("C:\\science\\2FInstability\\data\\ngc3245")
alphas = np.arange(0.2, 0.9, 0.03)
sigmas = np.arange(50.0, 200, 3.)
def sig_maj_gers_exp(l):
return sig_R_0*exp(-l*h_kin)*sqrt(sigPhi_to_sigR_real(l) * sin_i**2 + alpha **2 * cos_i**2)
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)))
global alpha, sig_R_0
cut_data = lambda l: l[0] > 21.
ma_data = filter(cut_data, zip(radii_maj, sig_maj_p))
mi_data = filter(cut_data, zip(radii_min, sig_min_p))
cut_data = lambda l: l[0] < 54.
ma_data = filter(cut_data, ma_data)
mi_data = filter(cut_data, mi_data)
for i,si in enumerate(sigmas):
for j,al in enumerate(alphas):
alpha = al
sig_R_0 = si
sqerr_maj = sum(power([sig_maj_gers_exp(p[0]) - p[1] for p in ma_data], 2))/len(ma_data)
sqerr_min = sum(power([sig_min_gers_exp(p[0]) - p[1] for p in mi_data], 2))/len(mi_data)
image_maj[i][j] = sqerr_maj
image_min[i][j] = sqerr_min
image[i][j] = 0.5*(sqerr_maj+sqerr_min)
return image, image_maj, image_min
pics_path = '.\\pics\\'
if not os.path.exists(pics_path):
os.makedirs(pics_path)
if os.path.isfile(pics_path + 'gers2_chi2_map_maj.npy'):
image_maj = np.load(pics_path + "gers2_chi2_map_maj.npy")
image_min = np.load(pics_path + "gers2_chi2_map_min.npy")
image = np.load(pics_path + "gers2_chi2_map.npy")
else:
image, image_maj, image_min = compute_chi2_maps(alphas=alphas, sigmas=sigmas)
np.save(pics_path + 'gers2_chi2_map_maj', image_maj)
np.save(pics_path + 'gers2_chi2_map_min', image_min)
np.save(pics_path + 'gers2_chi2_map', image)
In [ ]:
def plot_chi2_map(image, ax, log_scale=False, title='$\chi^2$', is_contour=False, vmax=0.):
if image is not None:
image_log = image
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)
degenereation_slice = lambda l: exp(poly_exp_line[0])/ sqrt(sin_i**2 + l**2 * cos_i**2)
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_{mean}$', is_contour=False, vmax=300.)
plot_chi2_map(image_maj, axes[1], log_scale=False, title='$\chi^2_{maj}$', is_contour=False, vmax=300.)
plot_chi2_map(image_min, axes[2], log_scale=False, title='$\chi^2_{min}$', is_contour=False, vmax=300.)
axes[0].plot(xx, map(degenereation_slice, xx), '-', color='m')
axes[1].plot(xx, map(degenereation_slice, xx), '-', color='m')
axes[2].plot(xx, map(degenereation_slice, xx), '-', color='m')
plt.show()
In [ ]:
plt.figure(figsize=(8, 8))
for al in alphas:
global alpha, sig_R_0
cut_data = lambda l: l[0] > 21.
ma_data = filter(cut_data, zip(radii_maj, sig_maj_p))
mi_data = filter(cut_data, zip(radii_min, sig_min_p))
cut_data = lambda l: l[0] < 54.
ma_data = filter(cut_data, ma_data)
mi_data = filter(cut_data, mi_data)
alpha = al
sig_R_0 = degenereation_slice(alpha)
sqerr_maj = sum(power([sig_maj_gers_exp(p[0]) - p[1] for p in ma_data], 2))/len(ma_data)
sqerr_min = sum(power([sig_min_gers_exp(p[0]) - p[1] for p in mi_data], 2))/len(mi_data)
plt.plot(alpha, sqerr_maj, 'x', color='blue')
plt.plot(alpha, sqerr_min, 'x', color='red')
plt.ylim(0, 500)
plt.show()
In [ ]:
from scipy.optimize import curve_fit
def powerfunc(x,a,b):
return a*np.power(x, b)
fitpars, covmat = curve_fit(powerfunc, r_ma_b[24:-3], vel_ma_b[24:-3])
variances = covmat.diagonal()
std_devs = np.sqrt(variances)
print fitpars,std_devs
print fitpars[0], fitpars[1]
plt.plot(r_ma_b, vel_ma_b, 'x-', color='blue', markersize=6)
plt.plot(r_ma_b[24:-3], vel_ma_b[24:-3], 'o', color='blue', markersize=6)
test_points = np.arange(0.0, max(r_ma_b), 0.1)
plt.plot(test_points, powerfunc(test_points, fitpars[0], fitpars[1]), 'r-')
plt.xlabel('$R$'); plt.ylim(0)
plt.ylabel('$V^{maj}_{\phi}(R)$')
plt.ylim(100, 250)
plt.show()
Теперь в ${\displaystyle \sigma_{\varphi}^{2}/\sigma_{R}^{2}=0.5\left(1+\frac{R}{\bar{v}_{\varphi}}\frac{d\bar{v}_{\varphi}}{dR}\right)}$ подставим производную для степенного закона $\bar{v}_{\varphi}(R)=V_0R^{\beta}$
In [ ]:
def sigPhi_to_sigR_power(R, a=fitpars[0], b=fitpars[1]):
return 0.5 * (1 + b)
plt.plot(test_points, [sigPhi_to_sigR_power(R) for R in test_points], 'd-', color='blue')
plt.axhline(y=0.5)
plt.xlabel('$R$')
plt.ylabel(r"$\sigma_{\varphi}^2/\sigma_{R}^2$")
plt.show()
In [ ]:
os.chdir("C:\\science\\2FInstability\\data\\ngc3245")
alphas = np.arange(0.2, 0.9, 0.03)
sigmas = np.arange(10.0, 200, 3.)
sig_maj_gers_exp = lambda l: sig_R_0*exp(-l*h_kin)*sqrt(sigPhi_to_sigR_power(l) * sin_i**2 + alpha **2 * cos_i**2)
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)))
global alpha, sig_R_0
cut_data = lambda l: l[0] > 21.
ma_data = filter(cut_data, zip(radii_maj, sig_maj_p))
mi_data = filter(cut_data, zip(radii_min, sig_min_p))
for i,si in enumerate(sigmas):
for j,al in enumerate(alphas):
alpha = al
sig_R_0 = si
sqerr_maj = sum(power([sig_maj_gers_exp(p[0]) - p[1] for p in ma_data], 2))/len(ma_data)
sqerr_min = sum(power([sig_min_gers_exp(p[0]) - p[1] for p in mi_data], 2))/len(mi_data)
image_maj[i][j] = sqerr_maj
image_min[i][j] = sqerr_min
image[i][j] = 0.5*(sqerr_maj+sqerr_min)
return image, image_maj, image_min
pics_path = '.\\pics\\'
if not os.path.exists(pics_path):
os.makedirs(pics_path)
if os.path.isfile(pics_path + 'gers3_chi2_map_maj.npy'):
image_maj = np.load(pics_path + "gers3_chi2_map_maj.npy")
image_min = np.load(pics_path + "gers3_chi2_map_min.npy")
image = np.load(pics_path + "gers3_chi2_map.npy")
else:
image, image_maj, image_min = compute_chi2_maps(alphas=alphas, sigmas=sigmas)
np.save(pics_path + 'gers3_chi2_map_maj', image_maj)
np.save(pics_path + 'gers3_chi2_map_min', image_min)
np.save(pics_path + 'gers3_chi2_map', image)
In [ ]:
def plot_chi2_map(image, ax, log_scale=False, title='$\chi^2$', is_contour=False, vmax=0.):
if image is not None:
image_log = image
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)
degenereation_slice = lambda l: exp(poly_exp_line[0])/ sqrt(sin_i**2 + l**2 * cos_i**2)
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_{mean}$', is_contour=False, vmax=300.)
plot_chi2_map(image_maj, axes[1], log_scale=False, title='$\chi^2_{maj}$', is_contour=False, vmax=300.)
plot_chi2_map(image_min, axes[2], log_scale=False, title='$\chi^2_{min}$', is_contour=False, vmax=300.)
axes[0].plot(xx, map(degenereation_slice, xx), '-', color='m')
axes[1].plot(xx, map(degenereation_slice, xx), '-', color='m')
axes[2].plot(xx, map(degenereation_slice, xx), '-', color='m')
plt.show()
In [ ]:
plt.figure(figsize=(8, 8))
for al in alphas:
global alpha, sig_R_0
cut_data = lambda l: l[0] > 21.
ma_data = filter(cut_data, zip(radii_maj, sig_maj_p))
mi_data = filter(cut_data, zip(radii_min, sig_min_p))
alpha = al
sig_R_0 = degenereation_slice(alpha)
sqerr_maj = sum(power([sig_maj_gers_exp(p[0]) - p[1] for p in ma_data], 2))/len(ma_data)
sqerr_min = sum(power([sig_min_gers_exp(p[0]) - p[1] for p in mi_data], 2))/len(mi_data)
plt.plot(alpha, sqerr_maj, 'x', color='blue')
plt.plot(alpha, sqerr_min, 'x', color='red')
plt.show()
In [ ]:
image_min = np.load(pics_path + "chi2_map_min.npy")
image_maj = np.load(pics_path + "chi2_map_maj.npy")
In [ ]:
os.chdir("C:\\Users\\root\\Dropbox\\RotationCurves\\PhD\\paper1\\text\\imgs")
alphas = np.arange(0.1, 1.2, 0.03)
sigmas = np.arange(100.0, 400, 3.)
import matplotlib.mlab as mlab
import matplotlib
fig, axes = plt.subplots(nrows=3, ncols=1, sharex=True, sharey=False, figsize=[8,16])
ax = axes[0]
# levels = np.linspace(start=image_min.min(), stop=20., num=5)
# levels = [2., 10., 50., 100.]
# im = ax.imshow(image_min, cmap='jet', vmin=image_min.min(), vmax=20., interpolation='spline16',
# origin="lower", aspect="auto")
# plt.show()
# levels = np.linspace(start=image_min.min()+1., stop=image_min.min()+20., num=5)
levels = np.linspace(start=image_min.min()*1.3, stop=image_min.min()*1.3+20, num=5)
cset=ax.contour(image_min, levels, colors = 'k', origin='lower', extent=[alphas[0],alphas[-1],sigmas[0],sigmas[-1]])
min_map_gutter = cset.collections[0].get_paths()
v1,v2 = min_map_gutter[1].vertices, min_map_gutter[0].vertices
x1,x2 = v1[:,0], v2[:,0]
y1,y2 = v1[:,1], v2[:,1]
plt.clabel(cset, inline=1, fontsize=10, fmt='%1.1f',)
ax.text(0.87, 260, '$\chi^2_{min}$', size = 24.)
ax.set_ylabel('$\sigma_{R,0}$', size=20.)
xx = np.arange(0.25, 1.0, 0.01)
ax.plot(xx, map(main_slice, xx), '--', color='black')
ax.set_ylim(170, 270)
ax.fill_between(x1, y1, 0, color='gray', alpha=0.3)
ax.fill_between(x2, y2, 0, color='white')
min_sigmas = np.where(image_min < image_min.min() + 0.03)
slice_alph, slice_sig = min_sigmas[1], min_sigmas[0]
slice_alph = map(lambda l: alphas[0] + (alphas[-1] - alphas[0])*l/len(image_min[0]) , slice_alph)
slice_sig = map(lambda l: sigmas[0] + (sigmas[-1] - sigmas[0])*l/len(image_min), slice_sig)
# ax.plot(slice_alph, slice_sig, '.', color='pink')
poly_slice = poly1d(polyfit(slice_alph, slice_sig, deg=3))
# ax.plot(xx, poly_slice(xx), '.-', color='black')
ax = axes[1]
# levels = np.linspace(start=image_maj.min()-4.3, stop=10., num=10)
# levels = [7., 10., 50., 100.]
levels = np.linspace(start=image_maj.min()+1., stop=image_maj.min()+60., num=6)
cset=ax.contour(image_maj, levels, hold='on', colors = 'k', origin='lower', extent=[alphas[0],alphas[-1],sigmas[0],sigmas[-1]])
plt.clabel(cset, inline=1, fontsize=10, fmt='%1.1f',)
ax.text(0.87, 260, '$\chi^2_{maj}$', size = 24.)
ax.set_ylabel('$\sigma_{R,0}$', size=20.)
xx = np.arange(0.25, 1.0, 0.01)
ax.plot(xx, map(main_slice, xx), '--', color='black')
ax.fill_between(x1, y1, 0, color='gray', alpha=0.3)
ax.fill_between(x2, y2, 0, color='white')
ax.set_ylim(170, 270)
ax = axes[2]
err_maj = []
for al in alphas:
global alpha, sig_R_0
alpha = al
sig_R_0 = main_slice(al)
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_maj = calc_chi2_normal(sig_maj_p, e_sig_maj_p, [sig_maj_exp(r) for r in radii_maj])
err_maj.append(sqerr_maj)
ax.plot(alphas, err_maj, '--', color='black')
err_maj1 = []
for pa in zip(x2,y2):
global alpha, sig_R_0
alpha = pa[0]
sig_R_0 = pa[1]
# sqerr_maj = calc_chi2_normal(sig_maj_p, e_sig_maj_p, [sig_maj_exp(r) for r in radii_maj])
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)])
err_maj1.append(sqerr_maj)
# ax.plot(x2, err_maj1, '-', color='black')
err_maj2 = []
for pa in zip(x1,y1):
global alpha, sig_R_0
alpha = pa[0]
sig_R_0 = pa[1]
# sqerr_maj = calc_chi2_normal(sig_maj_p, e_sig_maj_p, [sig_maj_exp(r) for r in radii_maj])
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)])
err_maj2.append(sqerr_maj)
# ax.plot(x1, err_maj2, '-', color='black')
ax.set_ylabel(r'$\chi^2$', size=20.)
ax.set_xlabel(r'$\alpha$', size=20.)
# ax.fill_between(x1, 0, err_maj2, color='grey', alpha=0.3)
# ax.fill_between(x2, 0, err_maj1, color='white')
# ax.set_ylim(20, 50)
import scipy.interpolate as sp
# f1 = sp.interp1d(x2, err_maj1, kind='linear')
# ax.fill_between(x1, map(f1, x1), err_maj2, color='grey', alpha=0.3)
f2 = sp.interp1d(x1, err_maj2, kind='linear')
ax.fill_between(x2, map(f2, x2), err_maj1, color='grey', alpha=0.3)
fig.subplots_adjust(hspace=0.)
axes[0].yaxis.get_major_ticks()[0].set_visible(False)
axes[1].yaxis.get_major_ticks()[0].set_visible(False)
ax.set_xlim(0.25, 0.99)
plt.savefig('ngc3245_maps.eps', format='eps')
plt.savefig('ngc3245_maps.png', format='png')
plt.savefig('ngc3245_maps.pdf', format='pdf', dpi=150)
plt.show()