In [1]:
%matplotlib inline
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 pandas as pd
import PIL as pil
import heapq
import matplotlib.pylab as pylab
pylab.rcParams['figure.figsize'] = 18, 14
# import seaborn as sns
# sns.set_palette("deep", desat=.6)

try:
    from PIL import Image
except:
    import Image
    
incl = 30.

sin_i2 = np.sin(incl*np.pi/180.)**2
cos_i2 = np.cos(incl*np.pi/180.)**2

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

plt.imshow(np.asarray(Image.open("shapiro_fg3a.png")))
plt.plot()


Out[2]:
[]

In [3]:
plt.imshow(np.asarray(Image.open("shapiro_fg4a.png")))
plt.plot()


Out[3]:
[]

In [4]:
pylab.rcParams['figure.figsize'] = 15, 5

#Velocity
lu = (215, 67)
lu_val = (30, 200)
rd = (1217, 463)
rd_val = (80, 25)

data = [(389,259), (411,257), (467,241), (513,249), (575,263), (589,245), (693,275), (699,261), (809,269), (871,293), (953,293)]

def extend_values(data, lu, lu_val, rd, rd_val):
    xscale = 1.0*(rd_val[0]-lu_val[0])/(rd[0]-lu[0])
    yscale = 1.0*(rd_val[1]-lu_val[1])/(lu[1]-rd[1])
    extended = []
    for d in data:
        extended.append((xscale*(d[0]-lu[0])+lu_val[0], yscale*(rd[1] - d[1])+rd_val[1]))
    return extended
        
e_data = extend_values(data, lu, lu_val, rd, rd_val)
plt.plot(zip(*e_data)[0],zip(*e_data)[1], 'or')

vel_fit = lambda l: 356.*np.power(l, -0.21)
plt.plot(np.linspace(30., 90., 100), map(vel_fit, np.linspace(30., 90., 100)), '--')
plt.ylim(25, 225)
plt.xlim(30, 95)
plt.show()



In [5]:
#Sig_maj
lu = (219, 509)
lu_val = (30, 120)
rd = (1215, 913)
rd_val = (80, 30)

data = [(387,689), (409,699), (463,721), (471,753), (515,731), (575,765), (587,783), (691,763), (699,813), (809,797), (873,819), (953,775)]

sig_maj_data = extend_values(data, lu, lu_val, rd, rd_val)
plt.plot(zip(*sig_maj_data)[0],zip(*sig_maj_data)[1], 'or')

sigR = lambda l: 213.*np.exp(-l/72.)
sigZ = lambda l: 124.*np.exp(-l/72.)
phi_to_R = lambda l: 0.5*(1 - 0.21)

sigR_min = lambda l: 193.*np.exp(-l/66.)
sigZ_min = lambda l: 115.*np.exp(-l/66.)
phi_to_R_min = lambda l: 0.5*(1 - 0.24)

sig_maj = lambda l: sqrt(sigR(l)**2 * (phi_to_R(l) * sin_i2 + sigZ(l)**2 * cos_i2/sigR(l)**2))
sig_maj_min = lambda l: sqrt(sigR_min(l)**2 * (phi_to_R_min(l) * sin_i2 + sigZ_min(l)**2 * cos_i2/sigR_min(l)**2))

plt.plot(np.linspace(30., 90., 100), map(sig_maj, np.linspace(30., 90., 100)), '--')
plt.plot(np.linspace(30., 90., 100), map(sig_maj_min, np.linspace(30., 90., 100)), '-')
plt.ylim(30, 130)
plt.xlim(30, 95)
plt.show()



In [6]:
#Sig_min
lu = (219, 957)
lu_val = (30, 120)
rd = (1217, 1359)
rd_val = (80, 30)

data = [(387,1071), (449,1083), (467,1085), (539,1103), (547,1133), (609,1149), (667,1187), (689,1191), (779,1193), (813,1215), 
        (875,1205), (963,1229), (971,1237), (1123,1225), (1131,1249), (1331,1269), (1365,1227)]

sig_min_data = extend_values(data, lu, lu_val, rd, rd_val)
plt.plot(zip(*sig_min_data)[0], zip(*sig_min_data)[1], 'or')

sig_min = lambda l: sqrt(sigR(l)**2 * sin_i2 + sigZ(l)**2 * cos_i2)
sig_min_min = lambda l: sqrt(sigR_min(l)**2 * sin_i2 + sigZ_min(l)**2 * cos_i2)

plt.plot(np.linspace(30., 90., 100), map(sig_min_min, np.linspace(30., 90., 100)), '-')
plt.plot(np.linspace(30., 90., 100), map(sig_min, np.linspace(30., 90., 100)), '--')
plt.ylim(30, 130)
plt.xlim(30, 95)
plt.show()


Картинка сравнения нашего приближения с Герсеновским:


In [24]:
import scipy.interpolate as inter

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=[10, 10], sharex=True)
fig.tight_layout()

radii_maj, sig_maj_p = zip(*sig_maj_data)
spl_maj = inter.UnivariateSpline(radii_maj, sig_maj_p, k=3, s=10000)
radii_min, sig_min_p = zip(*sig_min_data)
spl_min = inter.UnivariateSpline(radii_min, sig_min_p, k=3, s=10000)

sigR = lambda l: 213.*np.exp(-l/72.)
sigZ = lambda l: 124.*np.exp(-l/72.)
phi2_to_R2 = lambda l: 0.5*(1 - 0.21)

sig_maj = lambda l: sigR(l)*sqrt(phi2_to_R2(l) * sin_i2 + (124./213.)**2 * cos_i2)
sig_min = lambda l: sqrt(sigR(l)**2 * sin_i2 + sigZ(l)**2 * cos_i2)

ax1.plot(zip(*sig_maj_data)[0],zip(*sig_maj_data)[1], 'or')
ax1.set_xlim(30, 95)
ax1.set_ylim(30, 130)
ax1.set_ylabel('maj')
ax1.plot(np.linspace(30., 90., 100), map(sig_maj, np.linspace(30., 90., 100)), '--')
ax1.plot(np.linspace(35., 70., 100), map(spl_maj, np.linspace(35., 70., 100)), '-')

ax2.plot(zip(*sig_min_data)[0],zip(*sig_min_data)[1], 'or')
ax2.set_ylim(30, 130)
ax2.set_ylabel('min')
ax2.set_xlabel('R, arcsec')
ax2.plot(np.linspace(30., 90., 100), map(sig_min, np.linspace(30., 90., 100)), '--')
ax2.plot(np.linspace(30., 90., 100), map(spl_min, np.linspace(30., 90., 100)), '-')

plt.show()



In [8]:
h_kin = 72.
beta = 0.21

def sig_maj_exp(R):
    global alpha, sigR_0, h_kin, beta
    return np.exp(-R/h_kin)*sigR_0*sqrt(0.5*(1 - beta) * sin_i2 + alpha**2 * cos_i2)

def sig_min_exp(R):
    global alpha, sigR_0, h_kin
    return np.exp(-R/h_kin)*sigR_0*sqrt(sin_i2 + alpha**2 * cos_i2)

alphas = np.arange(0.25, 1., 0.01)
sigmas = np.arange(80.0, 250, 0.25)

def compute_chi2_maps(alph=(), sigm=()):
    '''Вычисляем все изображения, чтобы потом только настройки менять'''
    image_min = np.random.uniform(size=(len(sigm), len(alph)))
    image_maj = np.random.uniform(size=(len(sigm), len(alph)))
    image = np.random.uniform(size=(len(sigmas), len(alphas)))
    for i,si in enumerate(sigm):
        for j,al in enumerate(alph):
            global alpha, sigR_0
            alpha = al
            sigR_0 = si
            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)
            image_maj[i][j] = sqerr_maj
            image_min[i][j] = sqerr_min
    return image_maj, image_min

pics_path = "C:\\science\\2FInstability\\data\\ngc1068\\"
if not os.path.exists(pics_path):
    os.makedirs(pics_path)
if os.path.isfile(pics_path + 'chi2_map_maj.npy'):
    image_maj = np.load(pics_path + "chi2_map_maj.npy")
    image_min = np.load(pics_path + "chi2_map_min.npy")
else:
    image_maj, image_min = compute_chi2_maps(alph=alphas, sigm=sigmas)
    np.save(pics_path + 'chi2_map_maj', image_maj)
    np.save(pics_path + 'chi2_map_min', image_min)

In [9]:
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.):
    if image is not None:
        im = ax.imshow(image, cmap='jet', vmin=image.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)
        ax.set_title(title, size=20.)
        ax.set_ylabel('$\sigma_{R,0}$', size=20.)
        ax.set_xlabel(r'$\alpha$', size=20.)
        ax.grid(True)
        
min_sigmas = np.where(image_min < image_min.min() + 10.)
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)
poly_slice_min = poly1d(polyfit(slice_alph, slice_sig, deg=3))
maj_sigmas = np.where(image_maj < image_maj.min() + 10.)
slice_alph, slice_sig = maj_sigmas[1], maj_sigmas[0]
slice_alph = map(lambda l: alphas[0] + (alphas[-1] - alphas[0])*l/len(image_maj[0]) , slice_alph)
slice_sig = map(lambda l: sigmas[0] + (sigmas[-1] - sigmas[0])*l/len(image_maj), slice_sig)
poly_slice_maj = poly1d(polyfit(slice_alph, slice_sig, deg=3))
       
fig, axes = plt.subplots(nrows=1, ncols=2, sharex=False, sharey=True, figsize=[15,5])
plot_chi2_map(image_maj, axes[0], log_scale=False, title='$\chi^2_{maj}$', is_contour=False, vmax=300.)
plot_chi2_map(image_min, axes[1], log_scale=False, title='$\chi^2_{min}$', is_contour=False, vmax=300.)
axes[0].plot(0.58, 213., 'o', color='m')
axes[0].errorbar(0.58, 213., xerr=0.07, yerr=20., fmt='.', marker='.', mew=0, color='m')
axes[0].plot(alphas[20:], map(poly_slice_maj, alphas[20:]), '.-', label = 'maj -> maj', color= 'red')
axes[0].plot(alphas[20:], map(poly_slice_min, alphas[20:]), '.-', label = 'min -> maj', color= 'pink')
axes[1].plot(0.58, 213., 'o', color='m')
axes[1].errorbar(0.58, 213., xerr=0.07, yerr=20., fmt='.', marker='.', mew=0, color='m')
axes[1].plot(alphas[20:], map(poly_slice_min, alphas[20:]), '.-', label = 'min -> min', color='red')
axes[1].plot(alphas[20:], map(poly_slice_maj, alphas[20:]), '.-', label = 'maj -> min', color='pink') 
plt.show()

fig, axes = plt.subplots(nrows=1, ncols=2, sharex=False, sharey=False, figsize=[15,5])
plot_chi2_map((image_min + image_maj)/2, axes[0], log_scale=False, title='$\chi^2_{maj}+\chi^2_{min}$', is_contour=False, vmax=150.)
axes[0].plot(0.58, 213., 'o', color='m')
axes[0].errorbar(0.58, 213., xerr=0.07, yerr=20., fmt='.', marker='.', mew=0, color='m')

err_maj_1, err_maj_2 = [], []
err_min_1, err_min_2 = [], []
for al in alphas:
    global alpha, sigR_0
    alpha = al
    sigR_0 = poly_slice_min(alpha)
    err_maj_1.append(sum(power([sig_maj_exp(p[0]) - p[1] for p in sig_maj_data], 2))/len(sig_maj_data))
    err_min_1.append(sum(power([sig_min_exp(p[0]) - p[1] for p in sig_min_data], 2))/len(sig_min_data))
    sigR_0 = poly_slice_maj(alpha)
    err_maj_2.append(sum(power([sig_maj_exp(p[0]) - p[1] for p in sig_maj_data], 2))/len(sig_maj_data))
    err_min_2.append(sum(power([sig_min_exp(p[0]) - p[1] for p in sig_min_data], 2))/len(sig_min_data))
    
axes[1].plot(alphas, err_maj_1, '.-', label = 'min -> maj', color= 'red')
axes[1].plot(alphas, err_min_1, '.-', label = 'min -> min', color='orange')
axes[1].plot(alphas, err_maj_2, '.-', label = 'maj -> maj', color= 'green')
axes[1].plot(alphas, err_min_2, '.-', label = 'maj -> min', color='blue') 
axes[1].legend()

plt.show()