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


In [1]:
import matplotlib.pyplot as plt
import numpy as np
from numpy import poly1d, polyfit, power
import scipy.optimize
from math import *
from IPython.display import HTML
from IPython.display import Image
import os
import PIL as pil
import heapq
from matplotlib.ticker import MultipleLocator, FormatStrFormatter
import matplotlib.cm as cm
import scipy.interpolate as inter
%matplotlib inline

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

#Наклон галактики по данным Засова
incl=36.0

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

#Эффективный радиус балджа
r_eb = 6.7

In [2]:
os.chdir("C:\\science\\2FInstability\\data\\ngc1167")

In [3]:
# Данные по звездной кинематике Засова 2012 вдоль большей полуоси, не исправленные за наклон 
zasov_raw_data = np.loadtxt("v_stars_maZ.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_miZ.dat", float)
r_mi, vel_mi, e_vel_mi, sig_mi, e_sig_mi = zip(*zasov_raw_data)

# Данные по кинематике газа Struve, WSRT (не исправлено за наклон)
wsrt_raw_data = np.loadtxt("v_gas_WSRT.dat", float)
r_wsrt, vel_wsrt, e_vel_wsrt = zip(*wsrt_raw_data)

# Данные по кинематике газа Noordermee 2007, WSRT (не исправлено за наклон?)
noord_raw_data = np.loadtxt("v_gas_noord.dat", float)
r_noord, vel_noord, e_vel_noord = zip(*noord_raw_data)

plt.plot(r_ma, vel_ma, '.-', label="Zasov 2008, maj")
plt.plot(r_mi, vel_mi, '.-', label="Zasov 2008, min")
plt.plot(r_wsrt, vel_wsrt, '.-', label="gas Struve")
plt.plot(r_noord, vel_noord, '.-', label="gas Noordermeer 2007")
plt.legend()
plt.plot()


Out[3]:
[]

In [4]:
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, 4959.3, incl)
r_mi_b, vel_mi_b, e_vel_mi_b = correct_rotation_curve(r_mi, vel_mi, e_vel_mi,  0.0, 4959.3, incl)

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.plot()


Out[4]:
[]

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


In [5]:
poly_star = poly1d(polyfit(r_ma_b, vel_ma_b, deg=3))

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 [6]:
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 [7]:
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, 2)
plt.show()


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


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

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


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


Перегнем и приблизим полиномами:


In [10]:
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=9))

points = np.arange(0, max(radii_min), 0.1)
plt.plot(radii_maj, sig_maj_p, 's', label='$\sigma_{los}^{maj}$', color='blue')
plt.errorbar(radii_maj, sig_maj_p, yerr=e_sig_maj_p, fmt='.', marker='.', mew=0, color='blue')
plt.plot(points, poly_sig_maj(points), label = '$\sigma_{los}^{maj} polyfit$', color='blue')
plt.plot(radii_min, sig_min_p, 's', label='$\sigma_{los}^{min}$', color='red')
plt.errorbar(radii_min, sig_min_p, yerr=e_sig_min_p, fmt='.', marker='.', mew=0, color='red')
plt.plot(points, poly_sig_min(points), label = '$\sigma_{los}^{min} polyfit$', color='red')
plt.plot(fake_radii, fake_sig, 'bs', color='green', label='$fake points$')
plt.legend()
plt.ylim(0,250)
plt.xlim(0,55)
plt.show()


Table of Contents


In [11]:
%%javascript 
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')



In [ ]:

0. 0, sig_R_0=poly(0)


In [12]:
spl_maj = poly_sig_maj
spl_min = poly_sig_min

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.legend()
plt.ylim(0,250)
plt.xlim(0,55)
plt.show()

poly_sig_maj = spl_maj
poly_sig_min = spl_min



In [13]:
#Значение sig_los_min в 0
sig_min_0 = poly_sig_min(0)
print sig_min_0


185.454960214

И восстановим профили $\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 [14]:
# 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)


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


def sig_maj_exp(R):
    tmp = sigPhi_to_sigR_real(R) * sin_i**2 + alpha**2 * cos_i**2
    if tmp > 0:
        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)
    else:
        return -1000000
#     return sig_R_0*spl_min(R)/sig_min_0 * sqrt(sigPhi_to_sigR(R)**2 * sin_i**2 + alpha**2 * cos_i**2)
#     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 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)

Теперь то, с чего надо было начинать - построим картинки для разных значений $\alpha$ и $\sigma_{R,0}$. Для того, чтобы найти где минимум, попробуем построить просто двумерные карты $\chi^2$ для разных $\sigma_{R,0}$ $\alpha$: (это очень долго, так что пересчитывать в крайнем случае)


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

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) 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])
            if alpha >  0.39 and alpha < 0.41 and sig_R_0 > 135. and sig_R_0 < 143.:
                print alpha, sig_R_0, sqerr_maj, sqerr_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)


0.4 136.0 54.7919796696 18.9764586125
0.4 139.0 52.6711782658 18.1831330033
0.4 142.0 50.5935380703 17.4069431301

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



In [17]:
# Перебор alpha
alphas = np.arange(0.1, 0.7, 0.11)

# Перебор sig_R_0
sigmas = np.arange(270., 350., 16.)

# Те картинки, на которые стоит обратить особое внимание
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, [sig_maj_exp(Rr) for Rr in points], '--', 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(p[0]) - p[1] for p in 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)
                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, 50)
            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(35, 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(30, 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.text(25, 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 [18]:
main_slice = lambda l: sig_min_0/sqrt(sin_i**2 + cos_i**2 * l**2)

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)

# 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 = [100., 125., 150., 175., 200.]
# levels = [image_min.min()+0.02, image_min.min()+0.4, image_min.min()+1.1, image_min.min()+2., 
#           image_min.min()+3.1, image_min.min()+4.1]
# levels = np.linspace(start=image_min.min()+0.1, stop=image_min.min()+4.1, num=5)
levels = np.linspace(start=image_min.min()*1.1, stop=image_min.min()*1.1+4, num=5)
# im = ax.imshow(image_min, cmap='jet', vmin=image_min.min(), vmax=20., interpolation='spline16', 
#                    origin="lower", aspect="auto")
# plt.show()
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, 280, '$\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(180, 300)
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 = [image_maj.min()+0.2, image_maj.min()+0.7, image_maj.min()+1.1, image_maj.min()+2.1, image_maj.min()+3.1, 
#           image_maj.min()+4.1]
levels = np.linspace(start=image_maj.min()+0.3, stop=image_maj.min()+4.1, num=5)
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, 280, '$\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(150, 320)

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) 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])
    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])
    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(3.1, 3.5)


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('ngc1167_maps.eps', format='eps')
# plt.savefig('ngc1167_maps.png', format='png')
# plt.savefig('ngc1167_maps.pdf', format='pdf', dpi=150)

plt.show()


C:\Anaconda\lib\site-packages\matplotlib\text.py:52: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal
  if rotation in ('horizontal', None):
C:\Anaconda\lib\site-packages\matplotlib\text.py:54: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal
  elif rotation == 'vertical':

In [19]:
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=10.)
plot_chi2_map(image_min, axes[1], log_scale=False, title='$\chi^2_{min}$', is_contour=False, vmax=10.)
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=5.)
plt.show()


N1_maj=68,	 N2_min=55,	 chi^2_corr[0][0]=68.7042177524 (was 33.8471376298 and 96.8974443221)

In [ ]:


In [20]:
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(r) for r in radii_maj])
    sqerr_mi = calc_chi2_normal(sig_min_p, e_sig_min_p, [sig_min_exp(r) for r in radii_min])
    chisq = (sqerr_mi*len(sig_min_p) + sqerr_ma*len(sig_maj_p)) / (len(sig_min_p) + len(sig_maj_p))
    return chisq

x0 = np.array([100., 0.5])

res = opt.minimize(chisqfunc, x0, bounds=[(sigmas[0], sigmas[-1]), (alphas[0], alphas[-1])], method='L-BFGS-B')
print res


  status: 0
 success: True
    nfev: 78
     fun: 0.94009584882108832
       x: array([ 206.97301776,    0.81899892])
 message: 'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
     jac: array([ -3.33066907e-08,  -4.10782519e-06])
     nit: 18

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

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



In [22]:
import time

N = 300

result = []

start_time = time.time()

for i in log_progress(range(N)):
    global poly_sig_maj, poly_sig_min
    r, s = zip(*gen_next_normal(radii_maj, sig_maj_p, e_sig_maj_p))
#     poly_sig_maj = inter.UnivariateSpline(r, s, k=3, s=10000., w=w(e_sig_maj_p))
    poly_sig_maj = poly1d(polyfit(r, s,  deg=9))
    
    r, s = zip(*gen_next_normal(radii_min, sig_min_p, e_sig_min_p))
#     poly_sig_min = inter.UnivariateSpline(r, s, k=3, s=10000., w=w(e_sig_min_p))
    poly_sig_min = poly1d(polyfit(r, s,  deg=9))
    
#     result.append((opt.minimize(chisqfunc, x0, bounds=[(sigmas[0], sigmas[-1]), (alphas[0], alphas[-1])], method='L-BFGS-B').x,
#                   poly_sig_maj.get_coeffs(), poly_sig_min.get_coeffs()))
    result.append((opt.minimize(chisqfunc, x0, bounds=[(sigmas[0], sigmas[-1]), (alphas[0], alphas[-1])], method='L-BFGS-B').x,
                  poly_sig_maj.coeffs, poly_sig_min.coeffs))
print("--- %s seconds ---" % (time.time() - start_time))


---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-22-1f0220689313> in <module>()
      7 start_time = time.time()
      8 
----> 9 for i in log_progress(range(N)):
     10     global poly_sig_maj, poly_sig_min
     11     r, s = zip(*gen_next_normal(radii_maj, sig_maj_p, e_sig_maj_p))

C:\Users\root\.ipython\profile_default\startup\log_progress.py in log_progress(sequence, every, size)
      1 def log_progress(sequence, every=None, size=None):
----> 2     from ipywidgets import IntProgress, HTML, VBox
      3     from IPython.display import display
      4 
      5     is_iterator = False

C:\Anaconda\lib\site-packages\ipywidgets\__init__.py in <module>()
      3 from IPython import get_ipython
      4 from ._version import version_info, __version__
----> 5 from .widgets import *
      6 
      7 

C:\Anaconda\lib\site-packages\ipywidgets\widgets\__init__.py in <module>()
----> 1 from .widget import Widget, DOMWidget, CallbackDispatcher, register, widget_serialization
      2 
      3 from .trait_types import Color, EventfulDict, EventfulList
      4 
      5 from .widget_bool import Checkbox, ToggleButton, Valid

C:\Anaconda\lib\site-packages\ipywidgets\widgets\widget.py in <module>()
     11 
     12 from IPython.core.getipython import get_ipython
---> 13 from ipykernel.comm import Comm
     14 from traitlets.config import LoggingConfigurable
     15 from ipython_genutils.importstring import import_item

C:\Anaconda\lib\site-packages\ipykernel\__init__.py in <module>()
      1 from ._version import version_info, __version__, kernel_protocol_version_info, kernel_protocol_version
----> 2 from .connect import *

C:\Anaconda\lib\site-packages\ipykernel\connect.py in <module>()
     11 
     12 from IPython.core.profiledir import ProfileDir
---> 13 from IPython.paths import get_ipython_dir
     14 from ipython_genutils.path import filefind
     15 from ipython_genutils.py3compat import str_to_bytes

ImportError: No module named paths

In [ ]:
len(result)

In [ ]:
s, _, _ = zip(*result)
s,a = zip(*s)
plt.plot(a, s, '.')
plt.plot(alphas, map(main_slice, alphas), '--')
# plt.xlim(0.0, 0.99)
plt.ylim(0, 420)
plt.show()

In [ ]:


In [23]:
sig_maj_data = zip(r_ma[:-1], sig_ma[:-1], e_sig_ma[:-1])
sig_maj_data = map(bind_curve, sig_maj_data)
sig_maj_data.sort()
radii_maj1, sig_maj_p1, e_sig_maj_p1 = zip(*sig_maj_data) 

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_min1, sig_min_p1, e_sig_min_p1 = zip(*sig_min_data) 

points = np.arange(0, max(radii_min), 0.1)

1. r_ef, sig_R_0=spline(0)


In [24]:
# Граница. по которой обрезаем
cutted = r_eb


sig_maj_data = zip(radii_maj1, sig_maj_p1, e_sig_maj_p1)
sig_maj_data = filter(lambda l: l[0] > cutted, sig_maj_data)
radii_maj, sig_maj_p, e_sig_maj_p = zip(*sig_maj_data) 

sig_min_data = zip(radii_min1, sig_min_p1, e_sig_min_p1)
sig_min_data = filter(lambda l: l[0] > cutted, sig_min_data)
radii_min, sig_min_p, e_sig_min_p = zip(*sig_min_data) 


def w(arr):
    return map(lambda l: 1/(1. + l**2), arr)

spl_maj = inter.UnivariateSpline(radii_maj[::-1], sig_maj_p[::-1], k=3, s=10000., w=w(e_sig_maj_p))
spl_min = inter.UnivariateSpline(radii_min[::-1], sig_min_p[::-1], k=3, s=10000., w=w(e_sig_min_p))

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=cutted, color='black')
plt.legend()
plt.ylim(0,250)
plt.xlim(0,55)
plt.show()

poly_sig_maj = spl_maj
poly_sig_min = spl_min



In [25]:
#Значение sig_los_min в 0
sig_min_0 = poly_sig_min(0)
print sig_min_0


255.840463841

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

Monte-Carlo


In [26]:
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(r) for r in radii_maj])
    sqerr_mi = calc_chi2_normal(sig_min_p, e_sig_min_p, [sig_min_exp(r) for r in radii_min])
    chisq = (sqerr_mi*len(sig_min_p) + sqerr_ma*len(sig_maj_p)) / (len(sig_min_p) + len(sig_maj_p))
#     print sig_R_0, alpha, chisq
    return chisq

x0 = np.array([100., 0.5])

res = opt.minimize(chisqfunc, x0, bounds=[(sigmas[0], sigmas[-1]), (alphas[0], alphas[-1])], method='L-BFGS-B')
print res


  status: 0
 success: True
    nfev: 66
     fun: 0.79275721664809851
       x: array([  3.97000000e+02,   3.18968178e-01])
 message: 'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
     jac: array([-0.00048679,  0.        ])
     nit: 15

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

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



In [28]:
os.chdir("C:\\science\\2FInstability\\data\\ngc1167")
pics_path = '.cutted\\pics\\'
import time

sig_min_0 = 200.
N = 2000

result = []

start_time = time.time()
# fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, figsize=(16, 16))

if not os.path.exists(pics_path):
    os.makedirs(pics_path)
if os.path.isfile(pics_path + 'monte_carlo.npy'):
    result = np.load(pics_path + "monte_carlo.npy")
else:
    for i in log_progress(range(N)):
        global poly_sig_maj, poly_sig_min
        r, s = zip(*gen_next_normal(radii_maj, sig_maj_p, e_sig_maj_p))
        poly_sig_maj = inter.UnivariateSpline(r, s, k=3, s=10000., w=w(e_sig_maj_p))
#         poly_sig_maj = inter.UnivariateSpline(r, s, k=3, s=10000.)
#         ax1.plot(points, poly_sig_maj(points), label = '$\sigma_{los}^{maj}\, splinefit$', color='blue')
        
        r, s = zip(*gen_next_normal(radii_min, sig_min_p, e_sig_min_p))
        poly_sig_min = inter.UnivariateSpline(r, s, k=3, s=10000., w=w(e_sig_min_p))
#         poly_sig_min = inter.UnivariateSpline(r, s, k=3, s=10000.)
        
#         ax2.plot(points, poly_sig_min(points), label = '$\sigma_{los}^{maj}\, splinefit$', color='red')
        
        result.append((opt.minimize(chisqfunc, x0, bounds=[(sigmas[0], sigmas[-1]), (alphas[0], alphas[-1])], method='L-BFGS-B').x,
                      poly_sig_maj.get_coeffs(), poly_sig_min.get_coeffs()))
    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., 250.)
# ax2.set_ylim(0., 250.)
# plt.show()

In [29]:
len(result)


Out[29]:
2000

In [30]:
s, _, _ = zip(*result)
s,a = zip(*s)
plt.plot(a, s, '.')
plt.plot(alphas, map(main_slice, alphas), '--')
# plt.xlim(0.0, 0.99)
plt.ylim(0, 420)
plt.show()



In [31]:
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 [32]:
min_n = len(radii_min)

def func(x, alph, sig):
    global alpha, sig_R_0
    alpha = alph
    sig_R_0 = sig
    return [sig_min_exp(x[l]) if l > min_n else sig_maj_exp(x[l]) for l in range(len(x))]

In [33]:
rr = radii_maj + radii_min
print len(rr), min_n
sgs = sig_maj_p + sig_min_p
esgs = e_sig_maj_p + e_sig_min_p
esgs = [l/2 for l in esgs]


88 40

In [34]:
popt, pcov = opt.curve_fit(func, rr, sgs, sigma=esgs, absolute_sigma=True)
print popt, pcov

s_sq = np.array([((np.array(func(rr, popt[0], popt[1]))-np.array(sgs))**2)[l]/esgs[l] for l in range(len(rr))]).sum()/(len(rr)-2)
pcov = pcov * s_sq
for i in range(len(pcov)):
    print sqrt(pcov[i][i])


[   0.37590811  298.67179456] [[  7.12036159e-03  -1.37287302e+00]
 [ -1.37287302e+00   2.66950705e+02]]
0.492040282182
95.2719755719

In [35]:
def err_pcov(popt, pcov):
    s_sq = np.array([((np.array(func(rr, popt[0], popt[1]))-np.array(sgs))**2)[l]
                     /esgs[l]**2 for l in range(len(rr))]).sum()/(len(rr)-2)
    pcov = pcov * s_sq
    return (sqrt(pcov[0][0]), sqrt(pcov[1][1]))

In [36]:
import time

N = 100

result1 = []
start_time = time.time()

for i in log_progress(range(N)):
    global spl_maj, spl_min
    r, s = zip(*gen_next_normal(radii_maj, sig_maj_p, e_sig_maj_p))
    r1, s1 = zip(*gen_next_normal(radii_min, sig_min_p, e_sig_min_p))
    rr = r + r1
    sgs = s + s1
    popt, pcov = opt.curve_fit(func, rr, sgs, sigma=esgs, absolute_sigma=True)
    err = err_pcov(popt, pcov)
    result1.append((popt[0], popt[1], err[0], err[1]))
print("--- %s seconds ---" % (time.time() - start_time))


---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-36-654c1e57385f> in <module>()
      6 start_time = time.time()
      7 
----> 8 for i in log_progress(range(N)):
      9     global spl_maj, spl_min
     10     r, s = zip(*gen_next_normal(radii_maj, sig_maj_p, e_sig_maj_p))

C:\Users\root\.ipython\profile_default\startup\log_progress.py in log_progress(sequence, every, size)
      1 def log_progress(sequence, every=None, size=None):
----> 2     from ipywidgets import IntProgress, HTML, VBox
      3     from IPython.display import display
      4 
      5     is_iterator = False

C:\Anaconda\lib\site-packages\ipywidgets\__init__.py in <module>()
      3 from IPython import get_ipython
      4 from ._version import version_info, __version__
----> 5 from .widgets import *
      6 
      7 

C:\Anaconda\lib\site-packages\ipywidgets\widgets\__init__.py in <module>()
----> 1 from .widget import Widget, DOMWidget, CallbackDispatcher, register, widget_serialization
      2 
      3 from .trait_types import Color, EventfulDict, EventfulList
      4 
      5 from .widget_bool import Checkbox, ToggleButton, Valid

C:\Anaconda\lib\site-packages\ipywidgets\widgets\widget.py in <module>()
     11 
     12 from IPython.core.getipython import get_ipython
---> 13 from ipykernel.comm import Comm
     14 from traitlets.config import LoggingConfigurable
     15 from ipython_genutils.importstring import import_item

C:\Anaconda\lib\site-packages\ipykernel\__init__.py in <module>()
      1 from ._version import version_info, __version__, kernel_protocol_version_info, kernel_protocol_version
----> 2 from .connect import *

C:\Anaconda\lib\site-packages\ipykernel\connect.py in <module>()
     11 
     12 from IPython.core.profiledir import ProfileDir
---> 13 from IPython.paths import get_ipython_dir
     14 from ipython_genutils.path import filefind
     15 from ipython_genutils.py3compat import str_to_bytes

ImportError: No module named paths

In [ ]:
a,s,erra,errs = zip(*result1)
# plt.plot(a, s, '.')
plt.errorbar(a, s, yerr=errs, xerr=erra, fmt='o', marker='.', color='red')
plt.plot(alphas, map(main_slice, alphas), '--')
plt.xlim(0.0, 0.99)
plt.ylim(0, 400)
plt.show()

In [ ]:
plt.errorbar(a, s, yerr=errs, fmt='o', marker='.', color='red')
plt.plot(alphas, map(main_slice, alphas), '--')
plt.xlim(0.0, 0.99)
plt.ylim(0, 400)
plt.show()

Exponent on the edge


In [37]:
cutted = 15.0
# cutted = r_eb

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

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()
sig_min_data = filter(lambda l: l[0] > cutted, 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

# ind = [True if radii_min[l] > 15. and radii_min[l] < 30. else False for l in range(len(radii_min))]
#from 15'' to 30''
A,B = curve_fit(fuu, radii_min[:14], map(np.log, sig_min_p[:14]))[0]
# A,B = curve_fit(fuu, radii_min[:28], map(np.log, sig_min_p[:28]))[0]
print A, B


-0.0190221634649 5.4363013152

In [38]:
def fuu1(x, A, B):
    return np.exp(A*x+B)

In [39]:
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), '+')

A1, B1 = A, B

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.axvline(x=radii_min[28], color='black')
plt.legend()
plt.ylim(0, 250)
plt.show()


52.57

In [40]:
def sigR_ger_exp(R):
    return sig_R_0*exp(-R/h_kin)

def sigZ_ger_exp(R):
    return sigR_ger_exp(R)*alpha

def sigPhi_ger_exp(R):
    return sigPhi_to_sigR_real(R) * sigR_ger_exp(R)

def sig_maj_exp(R):
    return sqrt(sigPhi_ger_exp(R)**2 * sin_i**2 + sigZ_ger_exp(R)**2 * cos_i**2)

def sig_min_exp(R):
    return sqrt(sigR_ger_exp(R)**2 * sin_i**2 + sigZ_ger_exp(R)**2 * cos_i**2)

In [41]:
alphas = np.arange(0.1, 1.2, 0.03)
sigmas = np.arange(100.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

sig_maj_data = zip(radii_maj, sig_maj_p, e_sig_maj_p)
sig_maj_data = filter(lambda l: l[0] < radii_min[14], sig_maj_data)
radii_maj, sig_maj_p, e_sig_maj_p = zip(*sig_maj_data) 
radii_min, sig_min_p, e_sig_min_p = radii_min[:14], sig_min_p[:14], e_sig_min_p[:14]


# sig_maj_data = zip(radii_maj, sig_maj_p, e_sig_maj_p)
# sig_maj_data = filter(lambda l: l[0] < radii_min[28], sig_maj_data)
# radii_maj, sig_maj_p, e_sig_maj_p = zip(*sig_maj_data) 
# radii_min, sig_min_p, e_sig_min_p = radii_min[:28], sig_min_p[:28], e_sig_min_p[:28]

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


0.4 136.0 45.5117270776 26.9278000765
0.4 139.0 44.4114203479 26.1585567278
0.4 142.0 43.3248906001 25.4005672883

In [42]:
alphas = np.arange(0.1, 1.2, 0.03)
sigmas = np.arange(100.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=True, vmax=10.)
plot_chi2_map(image_min, axes[1], log_scale=False, title='$\chi^2_{min}$', is_contour=True, vmax=10.)
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=10.)
plt.show()


N1_maj=33,	 N2_min=26,	 chi^2_corr[0][0]=55.8507334797 (was 40.6304165742 and 67.8424983143)

In [43]:
alphas = np.arange(0.3, 0.9, 0.11)
sigmas = np.arange(286., 366., 16.)

plot_ranges(sigmas, alphas, good_pics=good_pics, calc_chi=True)
plt.show()



In [44]:
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()*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(100, 400)
# 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(100, 400)
# ax.set_ylim(50, 120)

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) 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])
    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])
    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 [45]:
import time

N = 10000

result = []
start_time = time.time()

fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, figsize=(16, 16))

sig_maj_data = zip(radii_maj, sig_maj_p, e_sig_maj_p)
sig_maj_data = filter(lambda l: l[0] < radii_min[14], sig_maj_data)
radii_maj, sig_maj_p, e_sig_maj_p = zip(*sig_maj_data) 
radii_min, sig_min_p, e_sig_min_p = radii_min[:14], sig_min_p[:14], e_sig_min_p[:14]

# sig_maj_data = zip(radii_maj, sig_maj_p, e_sig_maj_p)
# sig_maj_data = filter(lambda l: l[0] < radii_min[28], sig_maj_data)
# radii_maj, sig_maj_p, e_sig_maj_p = zip(*sig_maj_data) 
# radii_min, sig_min_p, e_sig_min_p = radii_min[:28], sig_min_p[:28], e_sig_min_p[:28]


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


if not os.path.exists(pics_path):
    os.makedirs(pics_path)
if os.path.isfile(pics_path + 'monte_carlo_exp1.npy'):
    result = np.load(pics_path + "monte_carlo_exp1.npy")
else:
    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)
    np.save(pics_path + 'monte_carlo_exp1', 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., 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 [46]:
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()
s1, a1 = s,a



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


Exponent on the edge 2

Возьмем прямой участок и подгоним экспонентами: $$\sigma_R = \sigma_{R,0}e^{-R/h}$$ $$\sigma_Z = \sigma_{Z,0}e^{-R/h}$$


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

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()
sig_min_data = filter(lambda l: l[0] > 30. and l[0] < 62. and l[1] > 100., 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

#from 30'' to 60''
A,B = curve_fit(fuu, radii_min, map(np.log, sig_min_p))[0]
print A, B


-0.000744089214043 4.93524296143

In [50]:
def fuu1(x, A, B):
    return np.exp(A*x+B)

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


1343.92

In [52]:
def sigR_ger_exp(R):
    return sig_R_0*exp(-R/h_kin)

def sigZ_ger_exp(R):
    return sigR_ger_exp(R)*alpha

def sigPhi_ger_exp(R):
    return sigPhi_to_sigR_real(R) * sigR_ger_exp(R)

def sig_maj_exp(R):
    return sqrt(sigPhi_ger_exp(R)**2 * sin_i**2 + sigZ_ger_exp(R)**2 * cos_i**2)

def sig_min_exp(R):
    return sqrt(sigR_ger_exp(R)**2 * sin_i**2 + sigZ_ger_exp(R)**2 * cos_i**2)

In [53]:
alpha = 0.25
sig_R_0 = main_slice(alpha)
print sig_min_exp(0.0), fuu1(0.0, A, B)


200.0 139.106936047

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


0.4 136.0 2.46425302393 4.2073380908
0.4 139.0 2.30123783993 3.88565695725
0.4 142.0 2.14634066236 3.57780231262

In [55]:
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=True, vmax=10.)
plot_chi2_map(image_min, axes[1], log_scale=False, title='$\chi^2_{min}$', is_contour=True, vmax=10.)
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=10.)
plt.show()


N1_maj=15,	 N2_min=11,	 chi^2_corr[0][0]=22.7588833694 (was 30.5550739095 and 17.0416769734)

In [56]:
alphas = np.arange(0.3, 0.9, 0.11)
sigmas = np.arange(186., 266., 16.)

plot_ranges(sigmas, alphas, good_pics=good_pics, calc_chi=True)
plt.show()



In [57]:
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(10, 400)
# ax.set_ylim(40, 200)

ax.axvline(x=0.4)
ax.axhline(y=139)
ax.axhline(y=161)


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.axvline(x=0.4)
ax.axhline(y=139)
ax.axhline(y=161)


ax.fill_between(x1, y1, 0, color='gray', alpha=0.3)
ax.fill_between(x2, y2, 0, color='white')
ax.set_ylim(10, 400)
# ax.set_ylim(50, 120)

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) 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])
    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])
    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:
    f2 = sp.interp1d(x1, err_maj2, kind='linear')
    ax.fill_between(x2, map(f2, x2), err_maj1, color='grey', alpha=0.3)
except Exception:
    f1 = sp.interp1d(x2, err_maj1, kind='linear')
    ax.fill_between(x1, map(f1, x1), err_maj2, 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 [58]:
import time


alphas = np.arange(0.001, 1.2, 0.03)

N = 10000

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


if not os.path.exists(pics_path):
    os.makedirs(pics_path)
if os.path.isfile(pics_path + 'monte_carlo_exp2.npy'):
    result = np.load(pics_path + "monte_carlo_exp2.npy")
else:
    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)
    np.save(pics_path + 'monte_carlo_exp2', 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., 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 [59]:
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 [60]:
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 [61]:
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()


Profile reconstruction

15-30'': h_kin = 52.6, sig = 281+/-28, alpha=0.72+/-0.09

30-60'': h_kin = 1343, sig = 222+/-35, alpha=0.30+/-0.08


In [62]:
def plot_exp_profiles(points, sig, alph, hkin, stl):
    global alpha, sig_R_0, h_kin
    alpha = alph
    sig_R_0 = sig
    h_kin = hkin
    plt.plot(points, [sig_maj_exp(l) for l in points], ls=stl, color='blue')
    plt.plot(points, [sig_min_exp(l) for l in points], ls=stl, color='red')

In [63]:
cutted = r_eb

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

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()
sig_min_data = filter(lambda l: l[0] > cutted, sig_min_data)
radii_min, sig_min_p, e_sig_min_p = zip(*sig_min_data)

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(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')


plot_exp_profiles(np.arange(15., 30., 0.1), 281., 0.72, 52.6, '-')
plot_exp_profiles(np.arange(15., 30., 0.1), 281.+28., 0.72+0.09, 52.6, '--')
plot_exp_profiles(np.arange(15., 30., 0.1), 281.-28., 0.72-0.09, 52.6, '--')


plot_exp_profiles(np.arange(30., 60., 0.1), 222., 0.3, 1343., '-')
plot_exp_profiles(np.arange(30., 60., 0.1), 222.+35., 0.3+0.08, 1343., '--')
plot_exp_profiles(np.arange(30., 60., 0.1), 222.-35., 0.3-0.08, 1343., '--')

# plt.legend()
plt.ylim(0, 250)
plt.show()


Картинки для статьи


In [84]:
tex_vkr_dir = 'C:\\Users\\root\\Dropbox\\RotationCurves\\PhD\\VKR\\imgs\\'
tex_imgs_dir = "C:\\Users\\root\\Dropbox\\RotationCurves\\PhD\\paper1\\text\\imgs\\"

In [98]:
def sigR_ger_exp(R):
    return sig_R_0*exp(-R/h_kin)

def sigZ_ger_exp(R):
    return sigR_ger_exp(R)*alpha

def sigPhi_ger_exp(R):
    return sigPhi_to_sigR_real(R) * sigR_ger_exp(R)

def sig_maj_exp(R):
    return sqrt(sigPhi_ger_exp(R)**2 * sin_i**2 + sigZ_ger_exp(R)**2 * cos_i**2)

def sig_min_exp(R):
    return sqrt(sigR_ger_exp(R)**2 * sin_i**2 + sigZ_ger_exp(R)**2 * cos_i**2)

In [99]:
# cutted = 15.0
cutted = r_eb

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

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()
sig_min_data = filter(lambda l: l[0] > cutted, sig_min_data)
radii_min, sig_min_p, e_sig_min_p = zip(*sig_min_data)

In [100]:
from matplotlib import rc

rc('text', usetex=True)
rc('text.latex',unicode=True)
rc('text.latex',preamble='\usepackage[russian]{babel}')

In [101]:
def plot_exp_profiles(ax, points, sig, alph, hkin, stl, mi=False):
    global alpha, sig_R_0, h_kin
    alpha = alph
    sig_R_0 = sig
    h_kin = hkin
    _dt = None 
    if mi:
        _dt = [sig_min_exp(l) for l in points]
        ax.plot(points, [sig_min_exp(l) for l in points], ls=stl, color='black')
    else:
        _dt = [sig_maj_exp(l) for l in points]
        ax.plot(points, [sig_maj_exp(l) for l in points], ls=stl, color='black')
    return _dt 

fig, axes = plt.subplots(figsize=[16,10], ncols=1, nrows=2)
ax2, ax3 = axes

radii_maj1, sig_maj_p1, e_sig_maj_p1 = zip(*filter(lambda l: l[0] < 30., zip(radii_maj, sig_maj_p, e_sig_maj_p)))
ax2.errorbar(radii_maj1, sig_maj_p1, yerr=e_sig_maj_p1, fmt='.', marker='.', mew=1, color='red')
radii_maj2, sig_maj_p2, e_sig_maj_p2 = zip(*filter(lambda l: l[0] < 72. and l[0]>30., zip(radii_maj, sig_maj_p, e_sig_maj_p)))
ax2.errorbar(radii_maj2, sig_maj_p2, yerr=e_sig_maj_p2, fmt='.', marker='.', mew=1, color='blue')

am1 = 0.72
dam1 = 0.09
sm1 = 281.
dsm1 = 28.
hkin1 = 52.6

am2 = 0.3
dam2 = 0.08
sm2 = 222.
dsm2 = 35.
hkin2 = 1343.

plot_exp_profiles(ax2, np.arange(r_eb, 30., 0.1), sm1, am1, hkin1, '-')
_dt1 = plot_exp_profiles(ax2, np.arange(r_eb, 30., 0.1), sm1+dsm1, am1+dam1, hkin1, '--')
_dt2 = plot_exp_profiles(ax2, np.arange(r_eb, 30., 0.1), sm1-dsm1, am1-dam1, hkin1, '--')
ax2.fill_between(np.arange(r_eb, 30., 0.1), _dt1, _dt2, color='gray', alpha=0.3)
plot_exp_profiles(ax2, np.arange(30., 60., 0.1), sm2, am2, hkin2, '-')
_dt1 = plot_exp_profiles(ax2, np.arange(30., 60., 0.1), sm2+dsm2, am2+dam2, hkin2, '--')
_dt2 = plot_exp_profiles(ax2, np.arange(30., 60., 0.1), sm2-dsm2, am2-dam2, hkin2, '--')
ax2.fill_between(np.arange(30., 60., 0.1), _dt1, _dt2, color='gray', alpha=0.3)
ax2.set_ylim(0, 250)
ax2.set_xlim(0, 70)
# ax2.text(60., 200., '$\sigma_{los}^{maj}$', size = 24.)
ax2.set_xticklabels([])
ax2.set_ylabel(u'$\sigma_{\\rm{los}}^{\\rm{maj}},\,$км/с', fontsize=30)

radii_min1, sig_min_p1, e_sig_min_p1 = zip(*filter(lambda l: l[0] < 30., zip(radii_min, sig_min_p, e_sig_min_p)))
ax3.errorbar(radii_min1, sig_min_p1, yerr=e_sig_min_p1, fmt='.', marker='.', mew=1, color='red')
radii_min2, sig_min_p2, e_sig_min_p2 = zip(*filter(lambda l: l[0] < 72. and l[0]>30., zip(radii_min, sig_min_p, e_sig_min_p)))
ax3.errorbar(radii_min2, sig_min_p2, yerr=e_sig_min_p2, fmt='.', marker='.', mew=1, color='blue')

plot_exp_profiles(ax3, np.arange(r_eb, 30., 0.1), sm1, am1, hkin1, '-', mi=True)
_dt1 = plot_exp_profiles(ax3, np.arange(r_eb, 30., 0.1), sm1+dsm1, am1+dam1, hkin1, '--', mi=True)
_dt2 = plot_exp_profiles(ax3, np.arange(r_eb, 30., 0.1), sm1-dsm1, am1-dam1, hkin1, '--', mi=True)
ax3.fill_between(np.arange(r_eb, 30., 0.1), _dt1, _dt2, color='gray', alpha=0.3)
plot_exp_profiles(ax3, np.arange(30., 62., 0.1), sm2, am2, hkin2, '-', mi=True)
_dt1 = plot_exp_profiles(ax3, np.arange(30., 62., 0.1), sm2+dsm2, am2+dam2, hkin2, '--', mi=True)
_dt2 = plot_exp_profiles(ax3, np.arange(30., 62., 0.1), sm2-dsm2, am2-dam2, hkin2, '--', mi=True)
ax3.fill_between(np.arange(30., 62., 0.1), _dt1, _dt2, color='gray', alpha=0.3)
ax3.set_ylim(0, 260)
ax3.set_xlim(0, 70)
# ax3.text(60., 200., '$\sigma_{los}^{min}$', size = 24.)
ax3.set_ylabel(u'$\sigma_{\\rm{los}}^{\\rm{min}},\,$км/с', fontsize=30)
ax3.set_xlabel(r'$R,\, ^{\prime\prime}$', fontsize=30)


# ax3.plot(points, fuu1(points, A, B), 'x', color='m')
# ax3.plot(points, fuu1(points, A1, B1), 'x', color='m')
for tick in ax3.yaxis.get_major_ticks():
    tick.label.set_fontsize(12)
for tick in ax3.xaxis.get_major_ticks():
    tick.label.set_fontsize(12)
for tick in ax2.yaxis.get_major_ticks():
    tick.label.set_fontsize(12)
for tick in ax2.xaxis.get_major_ticks():
    tick.label.set_fontsize(12)

fig.subplots_adjust(wspace=0.0, hspace=0.0)
# plt.savefig('1167_exponential_sve.eps', format='eps')
plt.savefig(tex_vkr_dir+'1167_exponential_sve_large.png', format='png', bbox_inches='tight')
# plt.savefig('1167_exponential_sve.pdf', format='pdf', dpi=150)
plt.show()



In [102]:
import matplotlib as mpl
mpl.rcParams.update(mpl.rcParamsDefault)

In [103]:
def plot_exp_profiles(ax, points, sig, alph, hkin, stl, mi=False):
    global alpha, sig_R_0, h_kin
    alpha = alph
    sig_R_0 = sig
    h_kin = hkin
    _dt = None 
    if mi:
        _dt = [sig_min_exp(l) for l in points]
        ax.plot(points, [sig_min_exp(l) for l in points], ls=stl, color='black')
    else:
        _dt = [sig_maj_exp(l) for l in points]
        ax.plot(points, [sig_maj_exp(l) for l in points], ls=stl, color='black')
    return _dt 

fig, axes = plt.subplots(figsize=[16,10], ncols=1, nrows=2)
ax2, ax3 = axes

radii_maj1, sig_maj_p1, e_sig_maj_p1 = zip(*filter(lambda l: l[0] < 30., zip(radii_maj, sig_maj_p, e_sig_maj_p)))
ax2.errorbar(radii_maj1, sig_maj_p1, yerr=e_sig_maj_p1, fmt='s', marker='s', mew=0, color='red', ms=8)
radii_maj2, sig_maj_p2, e_sig_maj_p2 = zip(*filter(lambda l: l[0] < 72. and l[0]>30., zip(radii_maj, sig_maj_p, e_sig_maj_p)))
ax2.errorbar(radii_maj2, sig_maj_p2, yerr=e_sig_maj_p2, fmt='o', marker='o', mew=0, color='blue', ms=8)

am1 = 0.72
dam1 = 0.09
sm1 = 281.
dsm1 = 28.
hkin1 = 52.6

am2 = 0.3
dam2 = 0.08
sm2 = 222.
dsm2 = 35.
hkin2 = 1343.

plot_exp_profiles(ax2, np.arange(r_eb, 30., 0.1), sm1, am1, hkin1, '-')
_dt1 = plot_exp_profiles(ax2, np.arange(r_eb, 30., 0.1), sm1+dsm1, am1+dam1, hkin1, '--')
_dt2 = plot_exp_profiles(ax2, np.arange(r_eb, 30., 0.1), sm1-dsm1, am1-dam1, hkin1, '--')
ax2.fill_between(np.arange(r_eb, 30., 0.1), _dt1, _dt2, color='none', alpha=0.7, hatch="\\", edgecolor='grey')
plot_exp_profiles(ax2, np.arange(30., 60., 0.1), sm2, am2, hkin2, '-')
_dt1 = plot_exp_profiles(ax2, np.arange(30., 60., 0.1), sm2+dsm2, am2+dam2, hkin2, '--')
_dt2 = plot_exp_profiles(ax2, np.arange(30., 60., 0.1), sm2-dsm2, am2-dam2, hkin2, '--')
ax2.fill_between(np.arange(30., 60., 0.1), _dt1, _dt2, color='none', alpha=0.7, hatch="/", edgecolor='grey')
ax2.set_ylim(0, 250)
ax2.set_xlim(0, 70)
# ax2.text(60., 200., '$\sigma_{los}^{maj}$', size = 24.)
ax2.set_xticklabels([])
ax2.set_ylabel(u'$\sigma_{\\rm{los}}^{\\rm{maj}},\,$km/s', fontsize=24)

radii_min1, sig_min_p1, e_sig_min_p1 = zip(*filter(lambda l: l[0] < 30., zip(radii_min, sig_min_p, e_sig_min_p)))
ax3.errorbar(radii_min1, sig_min_p1, yerr=e_sig_min_p1, fmt='s', marker='s', mew=0, color='red', ms=8)
radii_min2, sig_min_p2, e_sig_min_p2 = zip(*filter(lambda l: l[0] < 72. and l[0]>30., zip(radii_min, sig_min_p, e_sig_min_p)))
ax3.errorbar(radii_min2, sig_min_p2, yerr=e_sig_min_p2, fmt='o', marker='o', mew=0, color='blue', ms=8)

plot_exp_profiles(ax3, np.arange(r_eb, 30., 0.1), sm1, am1, hkin1, '-', mi=True)
_dt1 = plot_exp_profiles(ax3, np.arange(r_eb, 30., 0.1), sm1+dsm1, am1+dam1, hkin1, '--', mi=True)
_dt2 = plot_exp_profiles(ax3, np.arange(r_eb, 30., 0.1), sm1-dsm1, am1-dam1, hkin1, '--', mi=True)
ax3.fill_between(np.arange(r_eb, 30., 0.1), _dt1, _dt2, color='none', alpha=0.7, hatch="\\", edgecolor='grey')
plot_exp_profiles(ax3, np.arange(30., 62., 0.1), sm2, am2, hkin2, '-', mi=True)
_dt1 = plot_exp_profiles(ax3, np.arange(30., 62., 0.1), sm2+dsm2, am2+dam2, hkin2, '--', mi=True)
_dt2 = plot_exp_profiles(ax3, np.arange(30., 62., 0.1), sm2-dsm2, am2-dam2, hkin2, '--', mi=True)
ax3.fill_between(np.arange(30., 62., 0.1), _dt1, _dt2, color='none', alpha=0.7, hatch="/", edgecolor='grey')
ax3.set_ylim(0, 260)
ax3.set_xlim(0, 70)
# ax3.text(60., 200., '$\sigma_{los}^{min}$', size = 24.)
ax3.set_ylabel(u'$\sigma_{\\rm{los}}^{\\rm{min}},\,$km/s', fontsize=24)
ax3.set_xlabel(r'$R,\, ^{\prime\prime}$', fontsize=24)


# ax3.plot(points, fuu1(points, A, B), 'x', color='m')
# ax3.plot(points, fuu1(points, A1, B1), 'x', color='m')
for tick in ax3.yaxis.get_major_ticks():
    tick.label.set_fontsize(12)
for tick in ax3.xaxis.get_major_ticks():
    tick.label.set_fontsize(12)
for tick in ax2.yaxis.get_major_ticks():
    tick.label.set_fontsize(12)
for tick in ax2.xaxis.get_major_ticks():
    tick.label.set_fontsize(12)

fig.subplots_adjust(wspace=0.0, hspace=0.0)
plt.savefig(tex_imgs_dir+'1167_exponential_sve.eps', format='eps', bbox_inches='tight')
plt.savefig(tex_imgs_dir+'1167_exponential_sve.png', format='png', bbox_inches='tight')
plt.savefig(tex_imgs_dir+'1167_exponential_sve.pdf', format='pdf', dpi=150, bbox_inches='tight')
plt.show()



In [113]:
points = np.arange(r_eb, 30., 0.1)
alpha = am1
sig_R_0 = sm1
h_kin = hkin1 
plt.plot(points, [sigR_ger_exp(l) for l in points], ls='-', color='red')

alpha = am1-dam1
sig_R_0 = sm1-dsm1
plt.plot(points, [sigR_ger_exp(l) for l in points], ls='--', color='red')

alpha = am1+dam1
sig_R_0 = sm1+dsm1
plt.plot(points, [sigR_ger_exp(l) for l in points], ls='--', color='red')

plt.plot(points, [sigPhi_ger_exp(l) for l in points], ls='-', color='green')
plt.plot(points, [sigZ_ger_exp(l) for l in points], ls='-', color='blue')
plt.plot(points, [sig_min_exp(l) for l in points], ls='--', color='black')

points = np.arange(30., 60., 0.1)
alpha = am2
sig_R_0 = sm2
h_kin = hkin2
plt.plot(points, [sigR_ger_exp(l) for l in points], ls='-', color='red')

alpha = am2-dam2
sig_R_0 = sm2-dsm2
plt.plot(points, [sigR_ger_exp(l) for l in points], ls='--', color='red')

alpha = am2+dam2
sig_R_0 = sm2+dsm2
plt.plot(points, [sigR_ger_exp(l) for l in points], ls='--', color='red')

plt.plot(points, [sigPhi_ger_exp(l) for l in points], ls='-', color='green')
plt.plot(points, [sigZ_ger_exp(l) for l in points], ls='-', color='blue')
plt.plot(points, [sig_min_exp(l) for l in points], ls='--', color='black')

plt.ylim(0, 300)


Out[113]:
(0, 300)

In [92]:
spl_maj = inter.UnivariateSpline(radii_maj[::-1], sig_maj_p[::-1], k=3, s=10000., w=w(e_sig_maj_p))
spl_min = inter.UnivariateSpline(radii_min[::-1], sig_min_p[::-1], k=3, s=10000., w=w(e_sig_min_p))
poly_sig_maj = spl_maj
poly_sig_min = spl_min

In [93]:
os.chdir(tex_imgs_dir)


def sigR_ger_exp(R):
    return sig_R_0*exp(-R/h_kin)

def sigZ_ger_exp(R):
    return sigR_ger_exp(R)*alpha

def sigPhi_ger_exp(R):
    return sigPhi_to_sigR_real(R) * sigR_ger_exp(R)

def sig_maj_exp(R):
    return sqrt(sigPhi_ger_exp(R)**2 * sin_i**2 + sigZ_ger_exp(R)**2 * cos_i**2)

def sig_min_exp(R):
    return sqrt(sigR_ger_exp(R)**2 * sin_i**2 + sigZ_ger_exp(R)**2 * cos_i**2)



def plot_exp_profiles(ax, points, sig, alph, hkin, stl, mi=False):
    global alpha, sig_R_0, h_kin
    alpha = alph
    sig_R_0 = sig
    h_kin = hkin
    _dt = None 
    if mi:
        _dt = [sig_min_exp(l) for l in points]
        ax.plot(points, [sig_min_exp(l) for l in points], ls=stl, color='black')
    else:
        _dt = [sig_maj_exp(l) for l in points]
        ax.plot(points, [sig_maj_exp(l) for l in points], ls=stl, color='black')
    return _dt 

fig, axes = plt.subplots(figsize=[16,10], ncols=1, nrows=2)
ax2, ax3 = axes

radii_maj1, sig_maj_p1, e_sig_maj_p1 = zip(*filter(lambda l: l[0] < 30., zip(radii_maj, sig_maj_p, e_sig_maj_p)))
ax2.errorbar(radii_maj1, sig_maj_p1, yerr=e_sig_maj_p1, fmt='.', marker='.', mew=1, color='red')
radii_maj2, sig_maj_p2, e_sig_maj_p2 = zip(*filter(lambda l: l[0] < 72. and l[0]>30., zip(radii_maj, sig_maj_p, e_sig_maj_p)))
ax2.errorbar(radii_maj2, sig_maj_p2, yerr=e_sig_maj_p2, fmt='.', marker='.', mew=1, color='blue')

am1 = 0.72
dam1 = 0.09
sm1 = 281.
dsm1 = 28.
hkin1 = 52.6

am2 = 0.3
dam2 = 0.08
sm2 = 222.
dsm2 = 35.
hkin2 = 1343.

plot_exp_profiles(ax2, np.arange(r_eb, 30., 0.1), sm1, am1, hkin1, '-')
_dt1 = plot_exp_profiles(ax2, np.arange(r_eb, 30., 0.1), sm1+dsm1, am1+dam1, hkin1, '--')
_dt2 = plot_exp_profiles(ax2, np.arange(r_eb, 30., 0.1), sm1-dsm1, am1-dam1, hkin1, '--')
ax2.fill_between(np.arange(r_eb, 30., 0.1), _dt1, _dt2, color='gray', alpha=0.3)
plot_exp_profiles(ax2, np.arange(30., 60., 0.1), sm2, am2, hkin2, '-')
_dt1 = plot_exp_profiles(ax2, np.arange(30., 60., 0.1), sm2+dsm2, am2+dam2, hkin2, '--')
_dt2 = plot_exp_profiles(ax2, np.arange(30., 60., 0.1), sm2-dsm2, am2-dam2, hkin2, '--')
ax2.fill_between(np.arange(30., 60., 0.1), _dt1, _dt2, color='gray', alpha=0.3)
ax2.set_ylim(0, 250)
ax2.set_xlim(0, 70)
# ax2.text(60., 200., '$\sigma_{los}^{maj}$', size = 24.)
ax2.set_xticklabels([])
ax2.set_ylabel(r'$\sigma_{\rm{los}}^{\rm{maj}}$', fontsize=20)

radii_min1, sig_min_p1, e_sig_min_p1 = zip(*filter(lambda l: l[0] < 30., zip(radii_min, sig_min_p, e_sig_min_p)))
ax3.errorbar(radii_min1, sig_min_p1, yerr=e_sig_min_p1, fmt='.', marker='.', mew=1, color='red')
radii_min2, sig_min_p2, e_sig_min_p2 = zip(*filter(lambda l: l[0] < 72. and l[0]>30., zip(radii_min, sig_min_p, e_sig_min_p)))
ax3.errorbar(radii_min2, sig_min_p2, yerr=e_sig_min_p2, fmt='.', marker='.', mew=1, color='blue')

plot_exp_profiles(ax3, np.arange(r_eb, 30., 0.1), sm1, am1, hkin1, '-', mi=True)
_dt1 = plot_exp_profiles(ax3, np.arange(r_eb, 30., 0.1), sm1+dsm1, am1+dam1, hkin1, '--', mi=True)
_dt2 = plot_exp_profiles(ax3, np.arange(r_eb, 30., 0.1), sm1-dsm1, am1-dam1, hkin1, '--', mi=True)
ax3.fill_between(np.arange(r_eb, 30., 0.1), _dt1, _dt2, color='gray', alpha=0.3)
plot_exp_profiles(ax3, np.arange(30., 62., 0.1), sm2, am2, hkin2, '-', mi=True)
_dt1 = plot_exp_profiles(ax3, np.arange(30., 62., 0.1), sm2+dsm2, am2+dam2, hkin2, '--', mi=True)
_dt2 = plot_exp_profiles(ax3, np.arange(30., 62., 0.1), sm2-dsm2, am2-dam2, hkin2, '--', mi=True)
ax3.fill_between(np.arange(30., 62., 0.1), _dt1, _dt2, color='gray', alpha=0.3)
ax3.set_ylim(0, 260)
ax3.set_xlim(0, 70)
# ax3.text(60., 200., '$\sigma_{los}^{min}$', size = 24.)
ax3.set_ylabel(r'$\sigma_{\rm{los}}^{\rm{min}}$', fontsize=20)
ax3.set_xlabel(r'$R,\, \rm{arcsec}$', fontsize=20)


# ax3.plot(points, fuu1(points, A, B), 'x', color='m')
# ax3.plot(points, fuu1(points, A1, B1), 'x', color='m')


def sig_maj_exp(R):
    tmp = sigPhi_to_sigR_real(R) * sin_i**2 + alpha**2 * cos_i**2
    if tmp > 0:
        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)
    else:
        return -1000000
#     return sig_R_0*spl_min(R)/sig_min_0 * sqrt(sigPhi_to_sigR(R)**2 * sin_i**2 + alpha**2 * cos_i**2)
#     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):
    if R >= cutted:
        return sig_R_0*poly_sig_min(R)/sig_min_0 * sqrt(sin_i**2 + alpha**2 * cos_i**2)
    else:
        return -1000000

    
sig_min_0 = 210.

alpha = am1
sig_R_0 = sm1
ax2.plot(np.arange(r_eb, 62., 0.1), [sig_maj_exp(R) for R in np.arange(r_eb, 62., 0.1)], '--', lw=2, color='m')
alpha = am2
sig_R_0 = sm2
ax2.plot(np.arange(r_eb, 62., 0.1), [sig_maj_exp(R) for R in np.arange(r_eb, 62., 0.1)], '--', lw=2, color='g')


alpha = am1
sig_R_0 = sm1
ax3.plot(np.arange(r_eb, 62., 0.1), [sig_min_exp(R) for R in np.arange(r_eb, 62., 0.1)], '--', lw=2, color='m')
alpha = am2
sig_R_0 = sm2
ax3.plot(np.arange(r_eb, 62., 0.1), [sig_min_exp(R) for R in np.arange(r_eb, 62., 0.1)], '--', lw=2, color='g')


for tick in ax3.yaxis.get_major_ticks():
    tick.label.set_fontsize(12)
for tick in ax3.xaxis.get_major_ticks():
    tick.label.set_fontsize(12)
for tick in ax2.yaxis.get_major_ticks():
    tick.label.set_fontsize(12)
for tick in ax2.xaxis.get_major_ticks():
    tick.label.set_fontsize(12)


fig.subplots_adjust(wspace=0.0, hspace=0.0)
# plt.savefig('1167_exponential_sve.eps', format='eps')
# plt.savefig('1167_exponential_sve.png', format='png')
# plt.savefig('1167_exponential_sve.pdf', format='pdf', dpi=150)
plt.show()

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



In [95]:
import matplotlib as mpl
mpl.rcParams.update(mpl.rcParamsDefault)

In [96]:
os.chdir(tex_imgs_dir)

def plot_exp_profiles(ax, points, sig, alph, hkin, stl, mi=False):
    global alpha, sig_R_0, h_kin
    alpha = alph
    sig_R_0 = sig
    h_kin = hkin
    _dt = None 
    if mi:
        _dt = [sig_min_exp(l) for l in points]
        ax.plot(points, [sig_min_exp(l) for l in points], ls=stl, color='blue')
    else:
        _dt = [sig_maj_exp(l) for l in points]
        ax.plot(points, [sig_maj_exp(l) for l in points], ls=stl, color='red')
    return _dt 

fig = plt.figure(figsize=[16,10])
ax2 = plt.subplot2grid((5,5), (0,0), colspan=4, rowspan=1)
ax1 = plt.subplot2grid(shape=(5,5), loc=(1,0), colspan=4, rowspan=4)
ax3 = plt.subplot2grid((5,5), (1,4), colspan=1, rowspan=4)


ax1.plot(a, s, '.', ms=1., color='blue')
ax1.plot(a1, s1, '.', ms=1., color='red')
ax1.set_ylim(100, 410)
ax1.set_xlim(0., 1.15)
ax1.set_ylabel(r'$\sigma_{R,0}$', size=30.)
ax1.set_xlabel(r'$\alpha$', size=30.)
ax1.set_xticks(np.arange(0.0, 1.15, 0.1))


n, bins, patches = ax2.hist(a, 20, normed=True, facecolor='blue', alpha=1.0)
mua, stda = norm.fit(a)
xmin, xmax = ax2.get_xlim()
x = np.linspace(xmin, xmax, 100)
p = norm.pdf(x, mua, stda)
ax2.plot(x, p, 'k', linewidth=2)

n, bins, patches = ax2.hist(a1, 20, normed=True, facecolor='red', alpha=0.75)
mua1, stda1 = norm.fit(a1)
xmin, xmax = ax2.get_xlim()
x = np.linspace(xmin, xmax, 100)
p = norm.pdf(x, mua1, stda1)
ax2.plot(x, p, 'k', linewidth=2)

# ax2.grid(True)
ax2.set_xlim(0., 1.15)
ax2.set_yticks([])
ax2.set_xticks([])

n, bins, patches = ax3.hist(s, 20, normed=True, facecolor='blue', alpha=0.75, orientation=u'horizontal')
mus, stds = norm.fit(s)
xmin, xmax = ax3.get_ylim()
x = np.linspace(xmin, xmax, 100)
p = norm.pdf(x, mus, stds)
ax3.plot(p, x, 'k', linewidth=2)

n, bins, patches = ax3.hist(s1, 20, normed=True, facecolor='red', alpha=1.0, orientation=u'horizontal')
mus1, stds1 = norm.fit(s1)
xmin, xmax = ax3.get_ylim()
x = np.linspace(xmin, xmax, 100)
p = norm.pdf(x, mus1, stds1)
ax3.plot(p, x, 'k', linewidth=2)


# ax3.grid(True)
ax3.set_ylim(100, 410)
ax3.set_yticks([])
ax3.set_xticks([])

def plot_cov_ellipse(cov, pos, volume=.5, ax=None, fc='none', ec=[0,0,0], a=1, lw=2):
    import numpy as np
    from scipy.stats import chi2
    import matplotlib.pyplot as plt
    from matplotlib.patches import Ellipse

    def eigsorted(cov):
        vals, vecs = np.linalg.eigh(cov)
        order = vals.argsort()[::-1]
        return vals[order], vecs[:,order]

    if ax is None:
        ax = plt.gca()

    vals, vecs = eigsorted(cov)
    theta = np.degrees(np.arctan2(*vecs[:,0][::-1]))

    kwrg = {'facecolor':fc, 'edgecolor':ec, 'alpha':a, 'linewidth':lw}

    # Width and height are "full" widths, not radius
    width, height = 2 * np.sqrt(chi2.ppf(volume,2)) * np.sqrt(vals)
    ellip = Ellipse(xy=pos, width=width, height=height, angle=theta, **kwrg)

    ax.add_artist(ellip)

plot_cov_ellipse(np.cov(a1, s1), [mua1, mus1], ax=ax1, volume=0.68)
plot_cov_ellipse(np.cov(a, s), [mua, mus], ax=ax1, volume=0.68)

for tick in ax1.yaxis.get_major_ticks():
    tick.label.set_fontsize(12)
for tick in ax1.xaxis.get_major_ticks():
    tick.label.set_fontsize(12)


fig.subplots_adjust(wspace=0.0, hspace=0.0)
# plt.savefig('1167_mk.eps', format='eps')
plt.savefig(tex_vkr_dir+'1167_mk_large.png', format='png', bbox_inches='tight')
# plt.savefig('1167_mk.pdf', format='pdf', dpi=150)
plt.show()

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



In [97]:
os.chdir(tex_imgs_dir)

def plot_exp_profiles(ax, points, sig, alph, hkin, stl, mi=False):
    global alpha, sig_R_0, h_kin
    alpha = alph
    sig_R_0 = sig
    h_kin = hkin
    _dt = None 
    if mi:
        _dt = [sig_min_exp(l) for l in points]
        ax.plot(points, [sig_min_exp(l) for l in points], ls=stl, color='blue')
    else:
        _dt = [sig_maj_exp(l) for l in points]
        ax.plot(points, [sig_maj_exp(l) for l in points], ls=stl, color='red')
    return _dt 

fig = plt.figure(figsize=[16,10])
ax2 = plt.subplot2grid((5,5), (0,0), colspan=4, rowspan=1)
ax1 = plt.subplot2grid(shape=(5,5), loc=(1,0), colspan=4, rowspan=4)
ax3 = plt.subplot2grid((5,5), (1,4), colspan=1, rowspan=4)


ax1.plot(a, s, '.', ms=1., color='blue')
ax1.plot(a1, s1, '.', ms=1., color='red')
ax1.set_ylim(100, 410)
ax1.set_xlim(0., 1.15)
ax1.set_ylabel(r'$\sigma_{R,0}$', size=24.)
ax1.set_xlabel(r'$\alpha$', size=24.)
ax1.set_xticks(np.arange(0.0, 1.15, 0.1))


n, bins, patches = ax2.hist(a, 20, normed=True, facecolor='blue', alpha=1.0)
mua, stda = norm.fit(a)
xmin, xmax = ax2.get_xlim()
x = np.linspace(xmin, xmax, 100)
p = norm.pdf(x, mua, stda)
ax2.plot(x, p, 'k', linewidth=2)

n, bins, patches = ax2.hist(a1, 20, normed=True, facecolor='red', alpha=0.75)
mua1, stda1 = norm.fit(a1)
xmin, xmax = ax2.get_xlim()
x = np.linspace(xmin, xmax, 100)
p = norm.pdf(x, mua1, stda1)
ax2.plot(x, p, 'k', linewidth=2)

# ax2.grid(True)
ax2.set_xlim(0., 1.15)
ax2.set_yticks([])
ax2.set_xticks([])

n, bins, patches = ax3.hist(s, 20, normed=True, facecolor='blue', alpha=0.75, orientation=u'horizontal')
mus, stds = norm.fit(s)
xmin, xmax = ax3.get_ylim()
x = np.linspace(xmin, xmax, 100)
p = norm.pdf(x, mus, stds)
ax3.plot(p, x, 'k', linewidth=2)

n, bins, patches = ax3.hist(s1, 20, normed=True, facecolor='red', alpha=1.0, orientation=u'horizontal')
mus1, stds1 = norm.fit(s1)
xmin, xmax = ax3.get_ylim()
x = np.linspace(xmin, xmax, 100)
p = norm.pdf(x, mus1, stds1)
ax3.plot(p, x, 'k', linewidth=2)


# ax3.grid(True)
ax3.set_ylim(100, 410)
ax3.set_yticks([])
ax3.set_xticks([])

def plot_cov_ellipse(cov, pos, volume=.5, ax=None, fc='none', ec=[0,0,0], a=1, lw=2):
    import numpy as np
    from scipy.stats import chi2
    import matplotlib.pyplot as plt
    from matplotlib.patches import Ellipse

    def eigsorted(cov):
        vals, vecs = np.linalg.eigh(cov)
        order = vals.argsort()[::-1]
        return vals[order], vecs[:,order]

    if ax is None:
        ax = plt.gca()

    vals, vecs = eigsorted(cov)
    theta = np.degrees(np.arctan2(*vecs[:,0][::-1]))

    kwrg = {'facecolor':fc, 'edgecolor':ec, 'alpha':a, 'linewidth':lw}

    # Width and height are "full" widths, not radius
    width, height = 2 * np.sqrt(chi2.ppf(volume,2)) * np.sqrt(vals)
    ellip = Ellipse(xy=pos, width=width, height=height, angle=theta, **kwrg)

    ax.add_artist(ellip)

plot_cov_ellipse(np.cov(a1, s1), [mua1, mus1], ax=ax1, volume=0.68)
plot_cov_ellipse(np.cov(a, s), [mua, mus], ax=ax1, volume=0.68)

for tick in ax1.yaxis.get_major_ticks():
    tick.label.set_fontsize(12)
for tick in ax1.xaxis.get_major_ticks():
    tick.label.set_fontsize(12)


fig.subplots_adjust(wspace=0.0, hspace=0.0)
plt.savefig(tex_imgs_dir+'1167_mk.eps', format='eps', bbox_inches='tight')
plt.savefig(tex_imgs_dir+'1167_mk.png', format='png', bbox_inches='tight')
plt.savefig(tex_imgs_dir+'1167_mk.pdf', format='pdf', dpi=150, bbox_inches='tight')
plt.show()

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



In [ ]: