Сначала всякие настройки и импорты
In [26]:
import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize
from math import *
from IPython.display import HTML
from IPython.display import Image
import os
import pandas as pd
#Размер изображений
import matplotlib.pylab as pylab
pylab.rcParams['figure.figsize'] = 12, 12
#Наклон галактики по данным NED
incl=65.7
# Масштаб пк/секунда из NED
scale=284
Всякие картинки и БД для большего удобства:
In [27]:
# Данные из SDSS DR9
HTML('<iframe src=http://skyserver.sdss3.org/dr9/en/tools/explore/obj.asp?ra=14:03:01.0&dec=+34:45:28 width=1000 height=350></iframe>')
Out[27]:
In [28]:
# Данные из HYPERLEDA
HTML('<iframe src=http://leda.univ-lyon1.fr/ledacat.cgi?o=ngc5440 width=1000 height=350></iframe>')
Out[28]:
In [29]:
# Данные из NED
HTML('<iframe src=http://ned.ipac.caltech.edu/cgi-bin/objsearch?objname=ngc+5440&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[29]:
In [30]:
os.chdir("C:\\science\\2FInstability\\data\\ngc5440")
#Выводим данные за header-а файла
for line in file("v_stars_ma.dat"):
if line[0] == '#':
print(line)
Выведем также для удобства данные из обоих файлов:
In [31]:
ma = pd.read_csv('v_stars_ma.dat', skiprows=9, sep=' ', header=False)
# display(ma)
In [32]:
mi = pd.read_csv('v_stars_mi.dat', skiprows=9, sep=' ', header=False)
# display(mi)
Посмотрим теперь на все доступные даные по кривым вращения. Т.к. некоторые из них скорректированы, а некоторые нет, разобьем этот этап на несколько шагов.
Данные Коткова по газу вдоль главной оси:
In [33]:
n5440pa230_gas = pd.read_csv('.\\gas\\n5440pa230_gas_corr.txt', sep=' ')
display(n5440pa230_gas)
In [34]:
plt.plot(n5440pa230_gas['r,"'], n5440pa230_gas['v_Hbeta'], 's')
plt.plot(n5440pa230_gas['r,"'], n5440pa230_gas['v_OIII'], 'x')
plt.show()
In [35]:
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 = n5440pa230_gas['r,"']
vel_g_H = n5440pa230_gas['v_Hbeta']
e_vel_g_H = n5440pa230_gas['e_V']
# Кислород
r_g_O = n5440pa230_gas['r,"']
vel_g_O = n5440pa230_gas['v_OIII']
e_vel_g_O = n5440pa230_gas['e_v']
fig, subplot = subplots(2, 1)
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].plot(r_ma, vel_ma, '.-', label="Zasov 2012, maj")
subplot[0].plot(r_mi, vel_mi, '.-', label="Zasov 2012, min")
subplot[0].legend(loc = 'lower left')
subplot[0].set_ylim(3300, 4100)
subplot[1].imshow(np.asarray(Image.open("zasov2012_fig4_part.png")))
subplot[1].set_title("From Z2012 fig.4 part")
plt.plot()
Out[35]:
Попробуем как-нибудь автоматически вычистить значения по газу:
In [36]:
def bad_gas_data(data, diff=100):
pos = []
neg = []
level = None
for d in filter(lambda l: l[0] >= 0, data):
if level is None:
level = d[1]
elif abs(level - d[1]) < diff:
level = d[1]
else:
pos.append(d)
level = None
for d in filter(lambda l: l[0] < 0, data):
if level is None:
level = d[1]
elif abs(level - d[1]) < diff:
level = d[1]
else:
neg.append(d)
return neg + pos
H_data = zip(r_g_H, vel_g_H)
bad_H_r, bad_H_v = zip(*bad_gas_data(H_data))
O_data = filter(lambda l: l[0] < -32. or l[0] > 0., zip(r_g_O, vel_g_O))
bad_O_r, bad_O_v = zip(*bad_gas_data(O_data, diff=220.))
fig, subplot = subplots(2, 1)
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].plot(r_ma, vel_ma, '.', label="Zasov 2012, maj")
subplot[0].plot(r_mi, vel_mi, '.-', label="Zasov 2012, min")
subplot[0].plot(bad_H_r, bad_H_v, '^m')
subplot[0].plot(bad_O_r, bad_O_v, 'sm')
subplot[0].legend(loc = 'lower left')
# move legend outside of plot
box = subplot[0].get_position()
subplot[0].set_position([box.x0, box.y0, box.width * 0.8, box.height])
subplot[0].legend(loc='center left', bbox_to_anchor=(1, 0.5))
subplot[1].imshow(np.asarray(Image.open("zasov2012_fig4_part.png")))
subplot[1].set_title("From Z2012 fig.4 part")
plt.plot()
Out[36]:
In [37]:
# Вычищаем и дочищаем:
good_H_data = zip(r_g_H, vel_g_H)
good_H_data = [q for q in good_H_data if q not in zip(bad_H_r, bad_H_v)]
good_H_data = filter(lambda l: l[0] < 51., good_H_data)
good_H_r, good_H_v = zip(*good_H_data)
good_O_data = zip(r_g_O, vel_g_O)
good_O_data = [q for q in good_O_data if q not in zip(bad_O_r, bad_O_v)]
good_O_data = filter(lambda l: l[0] < 51., good_O_data)
good_O_data = filter(lambda l: (l[0] < -61. and abs(l[1]-4000) < 20.) or l[0] > -61., good_O_data)
good_O_r, good_O_v = zip(*good_O_data)
fig, subplot = subplots(2, 1)
subplot[0].plot(good_H_r, good_H_v, '^', label="gas Hbeta Zasov 2012", mfc='none')
subplot[0].plot(good_O_r, good_O_v, 's', label="gas OIII Zasov 2012", mfc='none')
subplot[0].plot(r_ma, vel_ma, '.', label="Zasov 2012, maj")
subplot[0].plot(r_mi, vel_mi, '.-', label="Zasov 2012, min")
subplot[0].legend(loc = 'lower left')
subplot[0].set_ylim(3300, 4100)
subplot[0].set_xlim(-100, 60)
subplot[1].imshow(np.asarray(Image.open("zasov2012_fig4_part.png")))
subplot[1].set_title("From Z2012 fig.4 part")
plt.plot()
Out[37]:
Теперь построим график дисперсий скоростей на луче зрения вдоль большой и малой осей по данным Засова:
In [38]:
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}_{min}$ и ${PA}_{maj}$ были определены правильно, просто щель для малой оси оказалась немного смещена. Данные и по скоростям, и по дисперсиям необходимо исправить, методика исправления описана ниже.
In [39]:
from mpl_toolkits.axes_grid.axislines import SubplotZero
import matplotlib.pyplot as plt
import numpy as np
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.11, -0.03, "$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.22, -0.3, "$P_{obs}$", fontsize=20)
ax.arrow(0.2, 0., 0., -0.28, head_width=0.02, head_length=0.02, fc='k', ec='k', width=0.005)
ax.plot((0., 0.2), (0., -0.3), '--', lw=0.3)
ax.text(0.06, -0.06, r"$\theta$", fontsize=20)
ax.text(0.2, -0.16, "$R_{obs}$", fontsize=20)
ax.plot(0.2, -0.3, '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=0^{\circ}$. Угол $\theta$ и его тригонометрические функции в свою очередь считаются как $$\theta = \arctan \frac{R_{obs}}{x_0}, \sin \theta = \frac{R_{obs}}{\sqrt{R_{obs}^2 + x_0^2}}, \cos \theta = \frac{x_0}{\sqrt{R_{obs}^2 + x_0^2}}$$ После восстановления скоростей исправлять за наклон галактики профиль уже не нужно (случай $\theta=0^{\circ}$).
Остается лишь вопрос, как найти $y_0$. В рассматриваемом случае это сделать просто - достаточно взять максимальное значение скорости вдоль малой оси и найти схожее значениее на большой.
In [40]:
r_ma_b, vel_ma_b, e_vel_b = r_ma, map(lambda l: abs(l-3715.0), vel_ma), e_vel_ma
r_mi_b, vel_mi_b, e_vel_mi_b = r_mi, map(lambda l: abs(l-3715.0), vel_mi), e_vel_mi
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.axvline(x = 9.)
plt.axvline(x = -7.)
plt.axhline(y =182.)
plt.xlim(-50., 50.)
plt.text(-40, 150, r"$x_0 = \frac{9+7}{2}=8$", fontsize=25)
plt.legend()
plt.plot()
Out[40]:
In [41]:
x_0 = 8.*cos(incl * pi / 180)
def sin_theta(R_obs):
return R_obs/sqrt(R_obs**2 + x_0**2)
def cos_theta(R_obs):
return x_0/sqrt(R_obs**2 + x_0**2)
cos_i, sin_i = cos(incl * pi / 180), sin(incl * pi / 180)
def correct_radii_min(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))
In [42]:
r_mi1 = map(correct_radii_min, r_mi)
In [43]:
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., 3715., incl)
r_mi_b, vel_mi_b, e_vel_mi_b = correct_rotation_curve(r_mi1, vel_mi, e_vel_mi, 0., 3715., 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.legend()
plt.ylim(0)
plt.plot()
Out[43]:
В дальнейшем используем только засовские данные по звездам по большой полуоси, приблизим их полиномом.
In [44]:
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()
Кривая вращения нам нужна для нахождения соотношения $\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 [45]:
def sigPhi_to_sigR_real(R):
return 0.5 * (1 + R*poly_star.deriv()(R) / poly_star(R))
plt.plot(test_points, [sigPhi_to_sigR_real(R) for R in test_points], 'd-', color='blue')
plt.axhline(y=0.5)
plt.axhline(y=0.0)
plt.xlabel('$R$')
plt.ylabel(r"$\sigma_{\varphi}^2/\sigma_{R}^2$")
plt.ylim(0)
plt.show()
Найдем теперь характерный масштаб $f=0.5(1+e^{-R/R_{0}})$:
In [46]:
def f(R, Ro):
return 0.5*(1 + np.exp( -R/Ro ))
xdata = test_points
ydata = sigPhi_to_sigR_real(xdata)
from scipy.optimize import curve_fit
popt, pcov = curve_fit(f, xdata, ydata, p0=[1.0])
Ro = popt[0]
plt.plot(xdata, ydata, 'x-')
plt.plot(xdata, [f(p, Ro) for p in xdata], 's')
plt.axhline(y=0.5)
plt.axhline(y=0.0)
plt.title('$R_{0} = %s $' % Ro)
plt.ylim(0)
plt.show()
Теперь знаем значение отношения $\sigma_{\varphi}^{2}/\sigma_{R}^{2}$ в любой точке, заведем соответствующую функцию:
In [47]:
def sigPhi_to_sigR(R):
return sqrt(f(R, Ro))
Построим графики дисперсий скоростей на луче зрения вдоль большой и малой оси ($\sigma_{los}^{maj}$ и $\sigma_{los}^{min}$), пересчитав расстояния для малой оси как было описано ранее:
In [48]:
# Исправляем значения вдоль малой оси на синус угла:
def correct_min(R):
return R / cos(incl * pi / 180)
# # r_mi_extend = map(correct_min, r_mi)
r_mi_extend = map(correct_radii_min, 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$')
plt.legend()
plt.xlim(-50,50)
plt.show()
In [49]:
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_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 = 200.0
# fake_radii, fake_sig = zip(*[(31.0 + i, 115*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=3))
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(70,250)
plt.xlim(0,50)
plt.show()
Посчитаем величину невязок для полученного приближения:
In [51]:
sqerr_maj = sum(power([poly_sig_maj(p[0]) - p[1] for p in sig_maj_data], 2))
sqerr_min = sum(power([poly_sig_min(p[0]) - p[1] for p in sig_min_data], 2))
print "Poly chi^2 for maj = %s, mean = %s" % (sqerr_maj, sqerr_maj / sig_maj_p.__len__())
print "Poly chi^2 for min = %s, mean = %s" % (sqerr_min, sqerr_min / sig_min_p.__len__())
Методика восстановления профилей $\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,min}^2(R_{obs})=\sigma_R^2(R)[\cos^2\theta+f(R)\sin^2\theta]\sin^2i + \alpha^2\sigma_R^2(R)\cos^2i,$$ Соответственно неправомерно делать одно из ключевых предположений - $\sigma_{los,min}(R)=\sigma_{R}(R)\times const$ и систему нельзя записать как линейную. Но можно поступить следующим образом: сосчитаем разность $$\sigma_{los,min}^2(R_{obs}) - \sigma_{los,maj}^2\sin^2\theta = \sigma_R^2\cos^2\theta[\sin^2i + \alpha^2\cos^2i]$$ Исходя из этого выражения можем легко построить график невязок для большой оси для разных $\alpha$. (НЕ УВЕРЕH ЧТО ТАМ ДЕЛАТЬ С Robs!!)
In [74]:
alpha = 0.4
# Возьмем набор каких-то точек
points_obs = np.arange(8., 30., 0.5)
# Скорректируем
points_real = map(correct_radii_min, points_obs)
# Углы
cos2_theta = map(np.square, map(cos_theta, points_obs))
sin2_theta = map(np.square, map(sin_theta, points_obs))
# Вот не уверен я тут в obs или real - вроде бы мы уже исправили малую
sigR2_difr = [(poly_sig_min(p[0])**2 - poly_sig_maj(p[0])**2 * p[2]) / (p[1] * (sin_i**2 + alpha**2 * cos_i**2))
for p in zip(points_real, cos2_theta, sin2_theta)]
# plt.plot(points_real, map(sqrt, sigR2_difr), 'o')
# plt.plot(points_real, sigR2_difr, 'o')
# plt.plot(points_real, [(poly_sig_min(p[0])**2 - poly_sig_maj(p[0])**2 * p[1]) for p in zip(points_real, sin2_theta)], 'o')
# plt.plot(points_real, [(poly_sig_min(p[0])**2 - poly_sig_maj(p[1])**2 * p[2]) for p in zip(points_obs, points_real, sin2_theta)], 'o')
# 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(points_real, map(poly_sig_maj, points_real), 'x')
plt.plot(points_real, [sqrt(p[1]*(sigPhi_to_sigR(p[0])**2 * sin_i**2 + alpha**2 + cos_i**2)) for p in zip(points_real, sigR2_difr)], 'o')
plt.show()
Выклядит все это удручающе, но посчитать график мы все же можем.
In [93]:
alphas = np.arange(0.1, 0.9, 0.01)
sig_maj_difr_errs = []
for al in alphas:
sigR2_difr = [(poly_sig_min(p[0])**2 - poly_sig_maj(p[0])**2 * p[2]) / (p[1] * (sin_i**2 + al**2 * cos_i**2))
for p in zip(points_real, cos2_theta, sin2_theta)]
sig2_maj_difr = [sqrt(p[1]*(sigPhi_to_sigR(p[0])**2 * sin_i**2 + al**2 + cos_i**2)) for p in zip(points_real, sigR2_difr)]
sqerr_maj = sum(power([poly_sig_maj(p[0]) - p[1] for p in zip(points_real, sig2_maj_difr)], 2))/len(points_real)
sig_maj_difr_errs.append(sqerr_maj)
plt.plot(alphas, sig_maj_difr_errs, '^')
plt.xlabel(r'$\alpha$')
plt.ylabel('$err$')
plt.show()
In [26]:
#Новая функция f'
def f_prime(R_real, R_obs):
return f(R_real, Ro)*sin_theta(R_obs)**2 + cos_theta(R_obs)**2
#ВОТ НЕ ПОНЯТНО - ВРОДЕ poly_sig_min(R_real)
def sigR_exp(R_real, R_obs, alp):
return sqrt(poly_sig_min(R_real)/(f_prime(R_real, R_obs)*sin_i**2 + alp**2*cos_i**2))
def sigZ_exp(R_real, R_obs, alp):
return alp * sigR_exp(R_real, R_obs, alp)
def sigPhi_exp(R_real, R_obs, alp):
return sigPhi_to_sigR(R_real) * sigR_exp(R_real, R_obs, alp)
def sig_min_exp(R_real, R_obs, alp):
return sqrt(sigR_exp(R_real, R_obs, alp)**2*f_prime(R_real, R_obs)*sin_i**2 + sigZ_exp(R_real, R_obs, alp)**2 * cos_i**2)
def sig_maj_exp(R_real, R_obs, alp):
return sqrt(sigPhi_exp(R_real, R_obs, alp)**2 * sin_i**2 + sigZ_exp(R_real, R_obs, alp)**2 * cos_i**2)
Проверим, что все правильно. Построим профили для $\alpha=0.5$:
In [27]:
alpha = 0.5
points_obs = np.arange(0, max(radii_min), 0.1)
points_real = map(correct_radii_min, points_obs)
plt.plot(points_real, [sigR_exp(Rr, Ro, alpha) for (Rr, Ro) in zip(points_real, points_obs)], '-', color='red', label='$\sigma_{R, exp}$')
plt.plot(points_real, [sigPhi_exp(Rr, Ro, alpha) for (Rr, Ro) in zip(points_real, points_obs)],
'-', color='blue', label=r'$\sigma_{\varphi, exp}$')
plt.plot(points_real, [sigZ_exp(Rr, Ro, alpha) for (Rr, Ro) in zip(points_real, points_obs)],
'-', color='black', label='$\sigma_{Z, exp}$')
plt.legend()
plt.show()
In [28]:
plt.plot(points, poly_sig_maj(points), '-', label = '$\sigma_{los}^{maj} polyfit$', color='blue')
plt.plot(points, [sig_maj_exp(Rr, Robs, alpha) 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_real, [sig_min_exp(Rr, Robs, alpha) for (Rr, Robs) in zip(points_real, points_obs)],
'--', color='red', label='$\sigma_{min, exp}$')
plt.legend()
# plt.ylim(0, 300)
plt.show()
In [28]:
In [28]:
In [28]:
Теперь восстановим исходные неизвестные - $\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 = 160
# alpha = 0.55
print "sig_R_0 = %s, alpha = %s" % (sig_R_0, alpha)
Построим полученные профили дисперсий скоростей:
In [ ]:
def sigR_exp(R):
return sig_R_0*poly_sig_min(R)/sig_min_0
def sigZ_exp(R):
return alpha * sigR_exp(R)
def sigPhi_exp(R):
return sigPhi_to_sigR(R) * sigR_exp(R)
plt.plot(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.ylim(0, 350)
plt.xlim(0, 50)
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$$
In [ ]:
def sig_maj_exp(R):
return sqrt(sigPhi_exp(R)**2 * sin(incl*pi/180)**2 + sigZ_exp(R)**2 * cos(incl*pi/180)**2)
def sig_min_exp(R):
return sqrt(sigR_exp(R)**2 * sin(incl*pi/180)**2 + sigZ_exp(R)**2 * cos(incl*pi/180)**2)
plt.plot(points, poly_sig_maj(points), '-', label = '$\sigma_{los}^{maj} polyfit$', color='blue')
plt.plot(points, [sig_maj_exp(R) for R in points], '--', color='blue', label='$\sigma_{maj, exp}$')
plt.plot(points, poly_sig_min(points), '-', label = '$\sigma_{los}^{min} polyfit$', color='red')
plt.plot(points, [sig_min_exp(R) for R in points], '--', color='red', label='$\sigma_{min, exp}$')
plt.legend()
plt.ylim(0, 350)
plt.xlim(0, 50)
plt.show()