Malin1-checkpoint


Malin1


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 seaborn as sns

import warnings
warnings.filterwarnings('ignore')

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

%matplotlib inline

# %install_ext http://raw.github.com/jrjohansson/version_information/master/version_information.py
%load_ext version_information
%reload_ext version_information

%version_information numpy, scipy, matplotlib


C:\Anaconda\lib\site-packages\matplotlib\__init__.py:872: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.
  warnings.warn(self.msg_depr % (key, alt_key))
Out[1]:
SoftwareVersion
Python2.7.12 64bit [MSC v.1500 64 bit (AMD64)]
IPython5.0.0
OSWindows 7 6.1.7601 SP1
numpy1.10.4
scipy0.17.0
matplotlib1.5.1
Wed Jul 20 12:42:02 2016 RTZ 2 (зима)

In [2]:
incl=38.0
cos_i, sin_i = cos(incl * pi / 180), sin(incl * pi / 180)

1. Наблюдательные данные


In [3]:
os.chdir("..\\data\\malin1\\data-moisav")

In [4]:
gas_raw_data = np.loadtxt("malin1_pa200_gas.txt", float)

In [5]:
r_gas, v_Hbeta, e_VHbeta, sig_Hb, e_sigHb, v_OIII, e_vOIII, sig_OIII, e_sigOIII, _, _ = zip(*gas_raw_data)

In [6]:
fig = plt.figure(figsize=[10, 8])
plt.errorbar(r_gas, v_Hbeta, yerr=e_VHbeta, label=r'$\rm{H}_{\beta}$', fmt='s-')
plt.errorbar(r_gas, v_OIII, yerr=e_vOIII, label=r'$\rm{OIII}$', fmt='o-')
plt.legend(fontsize=20)
plt.show()



In [7]:
fig = plt.figure(figsize=[10, 8])
plt.errorbar(r_gas, sig_Hb, yerr=e_sigHb, label=r'$\rm{H}_{\beta}$')
plt.errorbar(r_gas, sig_OIII, yerr=e_sigOIII, label=r'$\rm{OIII}$')
plt.legend(fontsize=20)
plt.show()



In [8]:
star_raw_data = np.loadtxt("malin1_pa235_sn5.txt", float)
r_pa235_sn5, v_pa235_sn5, ev_pa235_sn5, sig_pa235_sn5, esig_pa235_sn5, _, _, _, _ = zip(*star_raw_data)

In [9]:
star_raw_data = np.loadtxt("malin1_pa235_sn8.txt", float)
r_pa235_sn8, v_pa235_sn8, ev_pa235_sn8, sig_pa235_sn8, esig_pa235_sn8, _, _, _, _ = zip(*star_raw_data)

In [10]:
star_raw_data = np.loadtxt("malin1_pa333_sn10.txt", float)
r_pa333_sn10, v_pa333_sn10, ev_pa333_sn10, sig_pa333_sn10, esig_pa333_sn10, _, _, _, _ = zip(*star_raw_data)

In [11]:
star_raw_data = np.loadtxt("malin1_pa333_sn5.txt", float)
r_pa333_sn5, v_pa333_sn5, ev_pa333_sn5, sig_pa333_sn5, esig_pa333_sn5, _, _, _, _ = zip(*star_raw_data)

In [12]:
fig = plt.figure(figsize=[10,10])
plt.errorbar(r_pa235_sn5, v_pa235_sn5, ev_pa235_sn5, label=r'$PA=235^{\rm{o}} sn=5$', fmt='s-')
plt.errorbar(r_pa235_sn8, v_pa235_sn8, ev_pa235_sn8, label=r'$PA=235^{\rm{o}} sn=8$', fmt='o-')
plt.legend(fontsize=20)
plt.show()



In [13]:
fig = plt.figure(figsize=[10,10])
plt.errorbar(r_pa333_sn5, v_pa333_sn5, ev_pa333_sn5, label=r'$PA=333^{\rm{o}} sn=5$', fmt='s-')
plt.errorbar(r_pa333_sn10, v_pa333_sn10, ev_pa333_sn10, label=r'$PA=333^{\rm{o}} sn=10$', fmt='o-')
plt.legend(fontsize=20)
plt.show()



In [14]:
fig = plt.figure(figsize=[10,10])
plt.errorbar(r_pa235_sn5, sig_pa235_sn5, esig_pa235_sn5, label=r'$PA=235^{\rm{o}} sn=5$', fmt='s-')
plt.errorbar(r_pa235_sn8, sig_pa235_sn8, esig_pa235_sn8, label=r'$PA=235^{\rm{o}} sn=8$', fmt='o-')
plt.errorbar(r_pa333_sn5, sig_pa333_sn5, esig_pa333_sn5, label=r'$PA=333^{\rm{o}} sn=5$', fmt='s-')
plt.errorbar(r_pa333_sn10, sig_pa333_sn10, esig_pa333_sn10, label=r'$PA=333^{\rm{o}} sn=10$', fmt='o-')
plt.legend(fontsize=20)
plt.show()



In [36]:
#see slide 19
fig = plt.figure(figsize=[10,10])
# plt.errorbar(r_pa235_sn5, sig_pa235_sn5, esig_pa235_sn5, label=r'$PA=235^{\rm{o}} sn=5$', fmt='s')
plt.errorbar(r_pa235_sn8, sig_pa235_sn8, esig_pa235_sn8, label=r'$PA=235^{\rm{o}} sn=8$', fmt='o')
plt.errorbar(r_pa333_sn5, sig_pa333_sn5, esig_pa333_sn5, label=r'$PA=333^{\rm{o}} sn=5$', fmt='s')
# plt.errorbar(r_pa333_sn10, sig_pa333_sn10, esig_pa333_sn10, label=r'$PA=333^{\rm{o}} sn=10$', fmt='o')
plt.legend(fontsize=20)
plt.show()



In [16]:
fig = plt.figure(figsize=[10,10])
plt.errorbar(r_pa235_sn8, sig_pa235_sn8, esig_pa235_sn8, label=r'$PA=235^{\rm{o}} sn=8$', fmt='o-')
plt.errorbar(r_pa333_sn10, sig_pa333_sn10, esig_pa333_sn10, label=r'$PA=333^{\rm{o}} sn=10$', fmt='o-')
plt.legend(fontsize=20)
plt.xlim(-10,10)
plt.show()



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

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

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

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

235 - малая, 333 - большая


In [18]:
fig = plt.figure(figsize=[10,10])
plt.errorbar(r_pa333_sn5, v_pa333_sn5, ev_pa333_sn5, label=r'$PA=333^{\rm{o}} sn=5$', fmt='s-')
plt.errorbar(r_pa333_sn10, v_pa333_sn10, ev_pa333_sn10, label=r'$PA=333^{\rm{o}} sn=10$', fmt='o-')
plt.errorbar(r_pa235_sn5, v_pa235_sn5, ev_pa235_sn5, label=r'$PA=235^{\rm{o}} sn=5$', fmt='s-')
plt.errorbar(r_pa235_sn8, v_pa235_sn8, ev_pa235_sn8, label=r'$PA=235^{\rm{o}} sn=8$', fmt='o-')
plt.legend(fontsize=20)
plt.show()



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

b_vel = 24713

r_ma_b, vel_ma_b, e_vel_b = correct_rotation_curve(r_pa333_sn5, v_pa333_sn5, ev_pa333_sn5,  0.0, b_vel, incl)

plt.errorbar(r_ma_b, vel_ma_b, yerr=e_vel_b, fmt='.', marker='.', mew=0, color='blue')
# plt.xlim(0, 80.0)
plt.ylim(0)
plt.legend()
plt.plot()


Out[37]:
[]

In [38]:
import scipy.interpolate as inter
poly_star = poly1d(polyfit(r_ma_b, vel_ma_b, deg=3))
# poly_star = inter.UnivariateSpline (r_ma_b, vel_ma_b, k=3, s=10000.)

plt.plot(r_ma_b, vel_ma_b, 'o', color='blue', markersize=6)
test_points = np.arange(0.0, max(r_ma_b), 0.1)
plt.plot(test_points, poly_star(test_points), '-', color='red')
plt.xlabel('$R$'); plt.ylim(0)
plt.ylabel('$V^{maj}_{\phi}(R)$')
plt.show()



In [39]:
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='green')
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()



In [55]:
# def correct_min(R):    
#     return R / cos(incl * pi / 180) 

r_mi, sig_mi, e_sig_mi = r_pa235_sn8[:-3], sig_pa235_sn8[:-3], esig_pa235_sn8[:-3]
r_ma, sig_ma, e_sig_ma = r_pa333_sn5, sig_pa333_sn5, esig_pa333_sn5

# r_mi_extend = map(correct_min, r_mi)

plt.errorbar(r_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='x', mew=1, color='blue', ms=15)
plt.errorbar(r_mi, sig_mi, yerr=e_sig_mi, fmt='.', marker='x', mew=1, color='green', ms=15)
plt.xlabel('$R$')
plt.ylabel('$\sigma$')
plt.legend()
plt.xlim(-6, 6)
plt.show()



In [57]:
os.chdir("..\\olddata")

In [58]:
star_raw_data = np.loadtxt("MalinI_3.txt", float)
oldr_pa235, oldv_pa235, oldsig_pa235, oldev_pa235, oldesig_pa235 = zip(*star_raw_data)

In [59]:
star_raw_data = np.loadtxt("MalinI_2.txt", float)
oldr_pa333, oldv_pa333, oldsig_pa333, oldev_pa333, oldesig_pa333 = zip(*star_raw_data)

In [61]:
fig = plt.figure(figsize=[10,10])
plt.errorbar(oldr_pa235, oldsig_pa235, oldesig_pa235, label=r'$PA=235^{\rm{o}}$', fmt='s-')
plt.errorbar(oldr_pa333, oldsig_pa333, oldesig_pa333, label=r'$PA=333^{\rm{o}}$', fmt='o-')
plt.legend(fontsize=20)
plt.show()



In [62]:
fig = plt.figure(figsize=[10,10])
plt.errorbar(oldr_pa333, oldv_pa333, oldev_pa333, label=r'$PA=333^{\rm{o}}$', fmt='s-')
plt.errorbar(oldr_pa235, oldv_pa235, oldev_pa235, label=r'$PA=235^{\rm{o}}$', fmt='o-')
plt.legend(fontsize=20)
plt.show()



In [64]:
fig = plt.figure(figsize=[10,10])
plt.errorbar(oldr_pa235, oldsig_pa235, oldesig_pa235, label=r'$PA=235^{\rm{o}}$', fmt='s-', color='green')
plt.errorbar(oldr_pa333, oldsig_pa333, oldesig_pa333, label=r'$PA=333^{\rm{o}}$', fmt='o-', color='blue')
plt.errorbar(r_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='x', mew=1, color='blue', ms=15)
plt.errorbar(r_mi, sig_mi, yerr=e_sig_mi, fmt='.', marker='x', mew=1, color='green', ms=15)

plt.legend(fontsize=20)
plt.show()



In [72]:
plt.errorbar(r_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='x', mew=1, color='blue', ms=15)
plt.errorbar(r_mi, sig_mi, yerr=e_sig_mi, fmt='.', marker='x', mew=1, color='green', ms=15)
plt.xlabel('$R$')
plt.ylabel('$\sigma$')
plt.legend()
plt.xlim(-6, 6)
plt.show()



In [79]:
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=5))

sig_min_data = zip(r_mi, 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) 

poly_sig_min = poly1d(polyfit(radii_min, sig_min_p, deg=5))

points = np.arange(0, max(radii_maj), 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.legend()
plt.ylim(0, 250)
plt.show()



In [80]:
cutted = r_eb

sig_maj_data = zip(radii_maj, sig_maj_p, e_sig_maj_p)
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_min, sig_min_p, e_sig_min_p)
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(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.axvline(x=cutted, color='black')
plt.legend()
plt.ylim(0, 250)
plt.show()



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

points = filter(lambda l: l > cutted, points)

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, sig_min_p, k=2, 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.show()



In [101]:
#Значение sig_los_min в cutted
# sig_min_0 = spl_min(cutted)
sig_min_0 = spl_min(0)

#Значение sig_R в cutted
sig_R_0 = 80.
# sig_R_0 = sig_min_0

alpha = 0.5

def sigR_exp(R):
    return sig_R_0*spl_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)

In [103]:
def sig_maj_exp(R):
    return sig_R_0*spl_min(R)/sig_min_0 * sqrt(sigPhi_to_sigR_real(R) * sin_i**2 + alpha**2 * cos_i**2)

def sig_min_exp(R):
    return sig_R_0*spl_min(R)/sig_min_0 * sqrt(sin_i**2 + alpha**2 * cos_i**2)

In [108]:
alphas = np.arange(0.25, 1., 0.01)
sigmas = np.arange(100.0, 500, 0.25)

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])
            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)
# pics_path = '.cutted\\pics\\'
# if not os.path.exists(pics_path):
#     os.makedirs(pics_path)
# if os.path.isfile(pics_path + 'chi2_map.npy'):
#     image = np.load(pics_path + "chi2_map.npy")
#     image_maj = np.load(pics_path + "chi2_map_maj.npy")
#     image_min = np.load(pics_path + "chi2_map_min.npy")
# else:
#     image, image_maj, image_min = compute_chi2_maps(alphas=alphas, sigmas=sigmas)
#     np.save(pics_path + 'chi2_map', image)
#     np.save(pics_path + 'chi2_map_maj', image_maj)
#     np.save(pics_path + 'chi2_map_min', image_min)

In [109]:
from mpl_toolkits.axes_grid1 import make_axes_locatable
from matplotlib import cm

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