In [1]:

    

from IPython.display import HTML
from IPython.display import Image
import os
from PIL import Image as ImagePIL

%pylab
%matplotlib inline
%run ../../../utils/load_notebook.py


    

    




Using matplotlib backend: Qt4Agg
Populating the interactive namespace from numpy and matplotlib

In [2]:

    

from photometry import *


    

    




importing Jupyter notebook from photometry.ipynb
Using matplotlib backend: Qt4Agg
Populating the interactive namespace from numpy and matplotlib

In [3]:

    

from instabilities import *


    

    




importing Jupyter notebook from instabilities.ipynb
Using matplotlib backend: Qt4Agg
Populating the interactive namespace from numpy and matplotlib

In [4]:

    

from utils import *


    

    




importing Jupyter notebook from utils.ipynb

In [5]:

    

name = 'N1167'
gtype = 'S0(A)' #CALIFA
incl = 36.0  #Zasov, in LEDA 49deg, CALIFA 42deg TODO: check, хотя разброс и не влияет на результат сильно
scale = 0.321 #kpc/arcsec from NED

data_path = '../../../data/ngc1167'
sin_i, cos_i = np.sin(incl*np.pi/180.), np.cos(incl*np.pi/180.)


    

In [6]:

    

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


    

In [7]:

    

os.chdir(data_path)


    

In [8]:

    

# Данные из NED
HTML('<iframe src=http://ned.ipac.caltech.edu/cgi-bin/objsearch?objname=ngc+1167&extend=no&hconst=73&omegam=0.27&omegav=0.73\
&corr_z=1&out_csys=Equatorial&out_equinox=J2000.0&obj_sort=RA+or+Longitude&of=pre_text&zv_breaker=30000.0&list_limit=5\
&img_stamp=YES width=1000 height=350></iframe>')


    

    
Out[8]:

In [9]:

    

# Данные из HYPERLEDA
HTML('<iframe src=http://leda.univ-lyon1.fr/ledacat.cgi?o=ngc1167 width=1000 height=350></iframe>')


    

    
Out[9]:

In [10]:

    

fig, subplot = plt.subplots(1, 2, figsize=[12, 12])
subplot[0].imshow(np.asarray(ImagePIL.open("ngc1167_JHK.jpg")))
subplot[0].set_title("2MASS image JHK")
subplot[0].set_xticks([]);subplot[0].set_yticks([])
subplot[1].imshow(np.asarray(ImagePIL.open("ngc1167_SDSS.jpeg")))
subplot[1].set_title("SDSS DR9 whole image")
subplot[1].set_xticks([]);subplot[1].set_yticks([]);


    

In [11]:

    

Image('noord_pics/instab.png')


    

    
Out[11]:

In [12]:

    

r_q_n, q_n = zip(*np.loadtxt("Q1F_from_webplotdigitizer.dat", float, delimiter=','))
fig = plt.figure(figsize=[8,2])
plt.plot(r_q_n, q_n, 'o-')
plt.xlim(0, 80)
plt.ylim(0, 1.2);


    

In [13]:

    

Image('noord_pics/rot_vel.png', width=400)


    

    
Out[13]:

In [14]:

    

r_ma_n, vel_ma_n = zip(*np.loadtxt("v_gas_from_webplotdigitizer.dat", float, delimiter=','))
plt.plot(r_ma_n, vel_ma_n, 's-')
plt.xlim(0, 300)
plt.ylim(0, 450);


    

In [15]:

    

Image('zasov_pics/rot_vel.png')


    

    
Out[15]:

In [16]:

    

r_ma_z, vel_ma_z = zip(*np.loadtxt("v_stars_from_webplotdigitizer_zasov.dat", float, delimiter=','))
plt.plot(r_ma_z, vel_ma_z, 's')
plt.xlim(0, 120)
plt.ylim(0, 420);


    

In [17]:

    

# Данные по звездной кинематике Засова 2012 вдоль большей полуоси, не исправленные за наклон 
r_ma, vel_ma, e_vel_ma, sig_ma, e_sig_ma = zip(*np.loadtxt("v_stars_maZ.dat", float))

# Данные по звездной кинематике Засова 2012 вдоль малой полуоси, не исправленные за наклон 
r_mi, vel_mi, e_vel_mi, sig_mi, e_sig_mi = zip(*np.loadtxt("v_stars_miZ.dat", float))

plt.plot(r_ma, vel_ma, '.-', label="Zasov 2008, maj")
plt.plot(r_mi, vel_mi, '.-', label="Zasov 2008, min")
plt.legend();


    

In [18]:

    

%run ../../utils/rotvelutils #import util functions for rotcurv correction


    

In [19]:

    

v0 = 4959.3 #TODO: вспомнить, как нормально называется эта величина
r_ma_b, vel_ma_b, e_vel_b = correct_rotation_curve(r_ma, vel_ma, e_vel_ma,  0.0, v0, incl)
r_mi_b, vel_mi_b, e_vel_mi_b = correct_rotation_curve(r_mi, vel_mi, e_vel_mi,  0.0, v0, incl)

fig = plt.figure(figsize=[10,8])
plt.errorbar(r_ma_b, vel_ma_b, yerr=e_vel_b, fmt='.', marker='.', mew=0, color='blue', label = 'Zasov star maj')
plt.errorbar(r_mi_b, vel_mi_b, yerr=e_vel_mi_b, fmt='.', marker='.', mew=0, color='green', label = 'Zasov star min')
plt.plot(r_ma_z, vel_ma_z, 'o', color='m', label='Zasov from fig')
plt.legend()
plt.ylim(0, 450);


    

In [20]:

    

fig = plt.figure(figsize=[10,8])
plt.errorbar(r_ma_b, vel_ma_b, yerr=e_vel_b, fmt='.', marker='.', mew=0, color='blue', label = 'Zasov star maj')
plt.plot(r_ma_z, vel_ma_z, 'o', color='m', label='Zasov from fig')

test_points = np.linspace(0.0, max(r_ma_b), 100)

poly_star = poly1d(polyfit(r_ma_z, vel_ma_z, deg=3))
plt.plot(test_points, poly_star(test_points), '-', label='poly deg=3*')

poly_star = poly1d(polyfit(r_ma_b, vel_ma_b, deg=3))
plt.plot(test_points, poly_star(test_points), '-', label='poly deg=3')

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

import scipy.interpolate as inter
spl = inter.UnivariateSpline(r_ma_b, vel_ma_b, k=3, s=100., w=w(e_vel_b))
plt.plot(test_points, spl(test_points), '-', label='spl s=100 w^2')

plt.legend(loc='lower right')
plt.ylim(0, 450);


    

In [21]:

    

star_approx = poly1d(polyfit(r_ma_b, vel_ma_b, deg=3))


    

In [22]:

    

#TODO: свидетельство того, что надо раздвигать

# Исправляем значения вдоль малой оси на синус угла:    
def correct_min(R):    
    return R / cos_i

plt.errorbar(r_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='.', mew=0, color='blue', label='$\sigma_{los}^{maj}$')
r_mi_extend = map(correct_min, r_mi)
plt.errorbar(r_mi_extend, sig_mi, yerr=e_sig_mi, fmt='.', marker='.', mew=0, color='black', label='$\sigma_{los}^{min}$')
plt.xlabel('$R$')
plt.ylabel('$\sigma$')
plt.legend();


    

In [23]:

    

bind_curve = lambda p: (abs(p[0]), abs(p[1]), p[2])
sig_maj_data = sorted(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) 

sig_min_data = sorted(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)

points = np.arange(0, max(radii_min), 0.1)
plt.errorbar(radii_maj, sig_maj_p, yerr=e_sig_maj_p, fmt='.', marker='.', mew=0, color='blue', label='$\sigma_{los}^{maj}$')
plt.errorbar(radii_min, sig_min_p, yerr=e_sig_min_p, fmt='.', marker='.', mew=0, color='red', label='$\sigma_{los}^{min}$')
plt.legend()
plt.ylim(0,250)
plt.xlim(0,70);


    

In [24]:

    

r_sig_ma, sig_ma, e_sig_ma = radii_maj, sig_maj_p, e_sig_maj_p

r_sig_mi, sig_mi, e_sig_mi = radii_min, sig_min_p, e_sig_min_p

spl_min = inter.UnivariateSpline(r_sig_mi, sig_mi, k=3, s=7000.)
sig_min_lim = max(r_sig_mi)

# hack for smooth at the end
spl_maj = inter.UnivariateSpline(r_sig_ma[:-2] + tuple(np.linspace(53, 60, 10)), sig_ma[:-2] + tuple(np.linspace(90, 80, 10)), k=3, s=10000.)
# spl_maj = inter.UnivariateSpline(r_sig_ma, sig_ma, k=3, s=10000.)
sig_maj_lim = max(r_sig_ma)

points = np.linspace(0.1, max(r_sig_ma)+15., 100)


    

In [25]:

    

@flat_end(sig_maj_lim)
def sig_R_maj_minmin(r, spl_maj=spl_maj):
    return spl_maj(r).item()

@flat_end(sig_maj_lim)
def sig_R_maj_min(r, spl_maj=spl_maj):
    return spl_maj(r).item()/sqrt(sin_i**2 + 0.49*cos_i**2)
    
@flat_end(sig_maj_lim)
def sig_R_maj_max(r, spl_maj=spl_maj):
    return spl_maj(r).item()/sqrt(0.5*sin_i**2 + 0.09*cos_i**2)

@flat_end(sig_maj_lim)
def sig_R_maj_maxmax(r, spl_maj=spl_maj):
    return spl_maj(r)*sqrt(2)/sin_i
    
@flat_end(sig_maj_lim)
def sig_R_maj_maxmaxtrue(r, spl_maj=spl_maj):
    return spl_maj(r)/sin_i/sqrt(sigPhi_to_sigR_real(r))

@flat_end(sig_min_lim)
def sig_R_minor_minmin(r, spl_min=spl_min):
    return spl_min(r).item()

@flat_end(sig_min_lim)
def sig_R_minor_min(r, spl_min=spl_min):
    return spl_min(r).item()/sqrt(sin_i**2 + 0.49*cos_i**2)
    
@flat_end(sig_min_lim)
def sig_R_minor_max(r, spl_min=spl_min):
    return spl_min(r).item()/sqrt(sin_i**2 + 0.09*cos_i**2)

@flat_end(sig_min_lim)
def sig_R_minor_maxmax(r, spl_min=spl_min):
    return spl_min(r)/sin_i


    

In [26]:

    

plt.errorbar(r_sig_mi, sig_mi, yerr=e_sig_mi, fmt='.', marker='.', mew=0, color='red', label='$\sigma_{los}^{min}$')

plt.plot(points, map(sig_R_minor_minmin, points), label = 'minmin')
plt.plot(points, map(sig_R_minor_min, points), label = 'min')
plt.plot(points, map(sig_R_minor_max, points), label = 'max')
plt.plot(points, map(sig_R_minor_maxmax, points), label = 'maxmax')

plt.legend()
plt.ylim(0,400)
plt.xlim(0,100);


    

In [27]:

    

def sigPhi_to_sigR_real(R):
    try:
        return 0.5 * (1 + R*star_approx.derivative()(R) / star_approx(R))
    except:
        return 0.5 * (1 + R*star_approx.deriv()(R) / star_approx(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);


    

In [28]:

    

plt.errorbar(r_sig_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='.', mew=0, color='red', label='$\sigma_{los}^{maj}$')
plt.plot(points, map(sig_R_maj_minmin, points), label = 'minmin')
plt.plot(points, map(sig_R_maj_min, points), label = 'min')
plt.plot(points, map(sig_R_maj_max, points), label = 'max')
plt.plot(points, map(sig_R_maj_maxmax, points), label = 'maxmax')
plt.plot(points, map(sig_R_maj_maxmaxtrue, points), label = 'maxmaxtrue')

plt.legend()
plt.ylim(0,400)
plt.xlim(0,100);


    

In [29]:

    

fig = plt.figure(figsize=[10, 7])

plt.plot(points, map(sig_R_maj_minmin, points), label = 'maj minmin')
plt.plot(points, map(sig_R_maj_min, points), label = 'maj min')
plt.plot(points, map(sig_R_maj_max, points), label = 'maj max')
plt.plot(points, map(sig_R_maj_maxmax, points), label = 'maj maxmax')
plt.plot(points, map(sig_R_maj_maxmaxtrue, points), label = 'maj maxmaxtrue')

plt.plot(points, map(sig_R_minor_minmin, points), '--', label = 'minor minmin')
plt.plot(points, map(sig_R_minor_min, points), '--', label = 'minor min')
plt.plot(points, map(sig_R_minor_max, points), '--', label = 'minor max')
plt.plot(points, map(sig_R_minor_maxmax, points), '--', label = 'minor maxmax')

plt.legend()
plt.ylim(50,400)
plt.xlim(0,100);


    

In [30]:

    

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.

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

fig = plt.figure(figsize=[10, 7])

points_ = np.arange(6.7, 30., 0.1)
# alpha = am1
# sig_R_0 = sm1
# h_kin = hkin1 
# plt.plot(points, [sigR_ger_exp(l, sm1, hkin1) 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.fill_between(points_, [sigR_ger_exp(l, sm1-dsm1, hkin1) for l in points_], [sigR_ger_exp(l, sm1+dsm1, hkin1) for l in points_], color='red', alpha=0.3)

points_ = np.arange(30., 60., 0.1)
# alpha = am2
# sig_R_0 = sm2
# h_kin = hkin2
# plt.plot(points, [sigR_ger_exp(l, sm2, hkin2) 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.fill_between(points_, [sigR_ger_exp(l, sm2-dsm2, hkin2) for l in points_], [sigR_ger_exp(l, sm2+dsm2, hkin2) for l in points_], color='red', alpha=0.3)

plt.plot(points, map(sig_R_maj_minmin, points), label = 'maj minmin')
plt.plot(points, map(sig_R_maj_min, points), label = 'maj min')
plt.plot(points, map(sig_R_maj_max, points), label = 'maj max')
plt.plot(points, map(sig_R_maj_maxmax, points), label = 'maj maxmax')
plt.plot(points, map(sig_R_maj_maxmaxtrue, points), label = 'maj maxmaxtrue')

plt.plot(points, map(sig_R_minor_minmin, points), '--', label = 'minor minmin')
plt.plot(points, map(sig_R_minor_min, points), '--', label = 'minor min')
plt.plot(points, map(sig_R_minor_max, points), '--', label = 'minor max')
plt.plot(points, map(sig_R_minor_maxmax, points), '--', label = 'minor maxmax')

plt.legend()
plt.ylim(50,400)
plt.xlim(0,100);


    

In [31]:

    

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

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

gas_approx = poly1d(polyfit(r_ma_n, vel_ma_n, deg=5))
test_points = np.linspace(min(r_ma_n), max(r_ma_n), 100)
plt.plot(test_points, gas_approx(test_points), '--', label='approx')

spl_gas = inter.UnivariateSpline(r_noord, vel_noord, k=3, s=10000.)
plt.plot(test_points, spl_gas(test_points), '-', label='spl')

plt.plot(r_wsrt, vel_wsrt, '.-', label="gas Struve")
plt.plot(r_noord, vel_noord, '.-', label="gas Noordermeer 2007")
plt.plot(r_ma_n, vel_ma_n, 's-', label='HI from fig')
plt.ylim(0, 450)
plt.legend(loc='lower right');


    

In [32]:

    

#использовавшееся в дипломе приближение
ef = [104.55995900171898, 54.292488898429013, 34.856981586343409, 25.282657419650114, 19.107329429755854, 13.969314301465111, 
      10.21292602770659, 7.466641224935947, 5.458840202177822, 3.990942574473094, 2.9177667861366507, 2.133171013969345, 
      1.5595552723609254, 1.140186432133711, 0.8335870636080539]

fig = plt.figure(figsize=[12, 8])
plt.plot(test_points, [epicyclicFreq_real(gas_approx, x, scale) for x in test_points], '-')
plt.plot(test_points, [epicyclicFreq_real(spl_gas, x, scale) for x in test_points], '-')
# plt.plot(r_g_dens[1:], ef, '--')
plt.xlabel('$R, arcsec$')
plt.ylabel('$\kappa,\, km/s/kpc$', fontsize=15);


    

In [33]:

    

r_g_dens, gas_dens = zip(*np.loadtxt("gas_density.dat", float))

plt.plot(r_g_dens, gas_dens, 's-');


    

In [34]:

    

0.88*3.2


    

    
Out[34]:





2.8160000000000003

In [35]:

    

all_photometry = []


    

In [36]:

    

Image('diplom_photom.png')


    

    
Out[36]:

In [37]:

    

Image('photometry_pdf.png')


    

    
Out[37]:

In [38]:

    

M_R = -22.94 #11.69 - это правильно? надо брать абсолютные? в дипломе были относительные, тут разница уже существенная
M_B = -21.54 #13.40


    

In [39]:

    

print 'B : {:2.2f}; R: {:2.2f}.'.format(bell_mass_to_light(M_B-M_R, 'B', 'B-R'), bell_mass_to_light(M_B-M_R, 'R', 'B-R'))


    

    




B : 3.80; R: 2.71.

In [40]:

    

Image('noord_pics/n1167_photom.png')


    

    
Out[40]:

In [41]:

    

r_phot, mu_phot = zip(*np.loadtxt("n1167_noord_photoRB.dat", float, delimiter=','))
plt.plot(r_phot, mu_phot, 's')
plt.xlim(-10, 175)
plt.ylim(30, 15);


    

In [42]:

    

# R-band
r_eff_R = 6.7
mu_eff_R = 19.89 # уточнить это ли число
n_R = 1.7
mu0d_R = 20.35 # и тут тоже
h_disc_R = 24.2

mu_eff_Rc = 19.40 # уточнить это ли число
mu0d_Rc = 20.12 # и тут тоже


    

In [43]:

    

# B-band
r_eff_B = 6.7
mu_eff_B = 21.52 # уточнить это ли число
n_B = 1.7
mu0d_B = 22.24 # и тут тоже
h_disc_B = 27.5

mu_eff_Bc = 20.73 # уточнить это ли число
mu0d_Bc = 21.71 # и тут тоже


    

In [44]:

    

p_ = np.arange(0., 150., 0.1)

fig = plt.figure(figsize=[8, 5])
plt.plot(r_phot, mu_phot, 's')
plt.plot(p_, total_mu_profile([mu_bulge(l, mu_eff=mu_eff_R, r_eff=r_eff_R, n=n_R) for l in p_], 
                              [mu_disc(l, mu0=mu0d_R, h=h_disc_R) for l in p_]), '-', label='sum R', color='#007700')
plt.plot(p_, total_mu_profile([mu_bulge(l, mu_eff=mu_eff_B, r_eff=r_eff_B, n=n_B) for l in p_], 
                              [mu_disc(l, mu0=mu0d_B, h=h_disc_B) for l in p_]), '-', label='sum B', color='#000077')
plt.xlim(-10, 150)
plt.ylim(30, 15)
plt.legend();


    

In [45]:

    

b_r_color = M_B-M_R

M_to_L_Rc = bell_mass_to_light(b_r_color, 'R', 'B-R')
surf_R = [surf_density(mu=mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L=M_to_L_Rc, band='R') for l in p_]
plt.plot(p_, surf_R, '-', label='R [M/L={:2.2f}] absmag'.format(M_to_L_Rc))

M_to_L_Bc = bell_mass_to_light(b_r_color, 'B', 'B-R')
surf_B = [surf_density(mu=mu_disc(l, mu0=mu0d_Bc, h=h_disc_B), M_to_L=M_to_L_Bc, band='B') for l in p_]
plt.plot(p_, surf_B, '-', label='B [M/L={:2.2f}] absmag'.format(M_to_L_Bc))

diploma_color = 13.40 - 11.69

M_to_L_R = bell_mass_to_light(diploma_color, 'R', 'B-R')
surf_R = [surf_density(mu=mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L=M_to_L_R, band='R') for l in p_]
plt.plot(p_, surf_R, '-', label='R [M/L={:2.2f}]'.format(M_to_L_R))

M_to_L_B = bell_mass_to_light(diploma_color, 'B', 'B-R')
surf_B = [surf_density(mu=mu_disc(l, mu0=mu0d_Bc, h=h_disc_B), M_to_L=M_to_L_B, band='B') for l in p_]
plt.plot(p_, surf_B, '-', label='B [M/L={:2.2f}]'.format(M_to_L_B))

# максимальный диск
surf_R_max = [surf_density(mu=mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L=4., band='R') for l in p_]
plt.plot(p_, surf_R_max, '--', label='R [M/L={:2.2f}]'.format(4.))

plt.legend();


    

In [46]:

    

all_photometry.append(('Noorder R', r_eff_R, mu_eff_Rc, n_R, mu0d_Rc, h_disc_R, M_to_L_Rc, 
                       lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L=M_to_L_Rc, band='R')))
all_photometry.append(('Noorder B', r_eff_B, mu_eff_Bc, n_B, mu0d_Bc, h_disc_B, M_to_L_Bc, 
                       lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_Bc, h=h_disc_B), M_to_L=M_to_L_Bc, band='B')))

#максимальный диск
all_photometry.append(('Noorder R_max', r_eff_R, mu_eff_Rc, n_R, mu0d_Rc, h_disc_R, 4., 
                       lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L=4., band='R')))


    

In [47]:

    

A_lambda_g = 0.598
A_lambda_r = 0.414
A_lambda_i = 0.307


    

In [48]:

    

MUE_g=21.46976
Re_g=7.57209
n_g=2.20557
MU0_g=21.89252-A_lambda_g
hi_g=24.4525

MUE_r=20.7038
Re_r=9.09017
n_r=2.6612
MU0_r=21.0106-A_lambda_r
hi_r=25.8222


MUE_i=20.6073
Re_i=11.4278
n_i=3.20625
MU0_i=20.624-A_lambda_i
hi_i=27.508


    

In [49]:

    

M_g = disc_totmag(MU0_g, hi_g, scale)
M_r = disc_totmag(MU0_r, hi_r, scale)
M_i = disc_totmag(MU0_i, hi_i, scale)

fig = plt.figure(figsize=[16, 10])

sdss_surf = []

for color_desc, color in [('g-r', M_g-M_r), ('g-i', M_g-M_i) , ('r-i', M_r-M_i)]:

    surf_g = [surf_density(mu=mu_disc(l, mu0=MU0_g, h=hi_g), M_to_L=bell_mass_to_light(color, 'g', color_desc), band='g') for l in p_]
    plt.plot(p_, surf_g, '-', label='g [M/L={:2.2f}] {}'.format(bell_mass_to_light(color, 'g', color_desc), color_desc))
    sdss_surf.append([surf_g[0], bell_mass_to_light(color, 'g', color_desc), 'g', color_desc])

    surf_r = [surf_density(mu=mu_disc(l, mu0=MU0_r, h=hi_r), M_to_L=bell_mass_to_light(color, 'r', color_desc), band='r') for l in p_]
    plt.plot(p_, surf_r, '-', label='r [M/L={:2.2f}] {}'.format(bell_mass_to_light(color, 'r', color_desc), color_desc))
    sdss_surf.append([surf_r[0], bell_mass_to_light(color, 'r', color_desc), 'r', color_desc])

    surf_i = [surf_density(mu=mu_disc(l, mu0=MU0_i, h=hi_i), M_to_L=bell_mass_to_light(color, 'i', color_desc), band='i') for l in p_]
    plt.plot(p_, surf_i, '-', label='i [M/L={:2.2f}] {}'.format(bell_mass_to_light(color, 'i', color_desc), color_desc))
    sdss_surf.append([surf_i[0], bell_mass_to_light(color, 'i', color_desc), 'i', color_desc])

plt.legend();


    

In [50]:

    

sorted(sdss_surf)


    

    
Out[50]:





[[581.51153235522474, 2.9543283874174473, 'i', 'r-i'],
 [598.90332964036588, 3.0426861887143533, 'i', 'g-r'],
 [603.80473161710029, 3.0675874162779109, 'i', 'g-i'],
 [680.58878059249344, 3.7551321144337044, 'r', 'r-i'],
 [704.21878564344934, 3.8855100950313668, 'r', 'g-r'],
 [713.51635423752532, 3.936809204580852, 'r', 'g-i'],
 [945.02782535508686, 5.2523424198628579, 'g', 'r-i'],
 [990.84446798666215, 5.506985393511914, 'g', 'g-r'],
 [1007.4771691477172, 5.5994277952288014, 'g', 'g-i']]

In [51]:

    

all_photometry.append(('califa g (g-i)', Re_g, MUE_g, n_g, MU0_g, hi_g, 5.59, 
                       lambda l: surf_density(mu=mu_disc(l, mu0=MU0_g, h=hi_g), M_to_L=5.59, band='g')))

all_photometry.append(('califa r (g-r)', Re_r, MUE_r, n_r, MU0_r, hi_r, 3.88, 
                       lambda l: surf_density(mu=mu_disc(l, mu0=MU0_r, h=hi_r), M_to_L=3.88, band='r'))) #модель из середины

all_photometry.append(('califa i (r-i)', Re_i, MUE_i, n_i, MU0_i, hi_i, 2.95, 
                       lambda l: surf_density(mu=mu_disc(l, mu0=MU0_i, h=hi_i), M_to_L=2.95, band='i')))


    

In [52]:

    

show_all_photometry_table(all_photometry, scale)


    

    




+------+----------------+---------+----------+------+---------+----------+-------+-------------+-----------+
|      | Name           |   r_eff |   mu_eff |    n |   mu0_d |   h_disc |   M/L | M_d/M_sun   |   Sigma_0 |
|------+----------------+---------+----------+------+---------+----------+-------+-------------+-----------|
| 0.00 | Noorder R      |    6.70 |    19.40 | 1.70 |   20.12 |    24.20 |  2.71 | 2.30E+11.   |       605 |
| 1.00 | Noorder B      |    6.70 |    20.73 | 1.70 |   21.71 |    27.50 |  3.80 | 2.55E+11.   |       520 |
| 2.00 | Noorder R_max  |    6.70 |    19.40 | 1.70 |   20.12 |    24.20 |  4.00 | 3.39E+11.   |       893 |
| 3.00 | califa g (g-i) |    7.57 |    21.47 | 2.21 |   21.29 |    24.45 |  5.59 | 3.89E+11.   |      1006 |
| 4.00 | califa r (g-r) |    9.09 |    20.70 | 2.66 |   20.60 |    25.82 |  3.88 | 3.04E+11.   |       703 |
| 5.00 | califa i (r-i) |   11.43 |    20.61 | 3.21 |   20.32 |    27.51 |  2.95 | 2.85E+11.   |       581 |
+------+----------------+---------+----------+------+---------+----------+-------+-------------+-----------+

In [53]:

    

fig = plt.figure(figsize=[10,6])

plt.plot(test_points, spl_gas(test_points), '-', label='spl')
plt.plot(r_ma_n, vel_ma_n, 's-', label='HI from fig')

for photom in all_photometry:
    plt.plot(test_points, map(lambda l: disc_vel(l, photom[7](0), photom[5], scale), test_points), '--', label=photom[0])


plt.ylim(0, 450)
plt.xlim(0, 300)
plt.legend(loc='best');


    

In [54]:

    

Image('califa_2016_halpha_plus_sfr.png')


    

    
Out[54]:

In [55]:

    

# TODO: разобраться и сделать обоснованно 
def plot_SF(ax):
#     ax.plot([55., 65.], [0., 0.], '-', lw=7., color='red')
    ax.plot([22., 38.], [0., 0.], '-', lw=7., color='red')
    
plot_SF(plt.gca())
plt.xlim(0, 350)
plt.ylim(0, 200)


    

    
Out[55]:





(0, 200)






    

In [56]:

    

# from instabilities import *


    

In [57]:

    

Image('diplom_results.png')


    

    
Out[57]:

In [58]:

    

sound_vel = 6.  #скорость звука в газе, км/с
data_lim = max(max(r_sig_ma), max(r_sig_mi)) #где заканчиваются данные


    

In [59]:

    

fig = plt.figure(figsize=[12, 6])
gd_data = zip(r_g_dens, gas_dens)

plt.plot(r_g_dens, [1./Qg(epicycl=l[0], sound_vel=sound_vel, gas_density=l[1]) for l in 
                    zip([epicyclicFreq_real(gas_approx, x, scale) for x in r_g_dens], gas_dens)], 's-', label='$Q_{gas}$')
plt.plot(r_g_dens, [1./Qg(epicycl=l[0], sound_vel=4., gas_density=l[1]) for l in 
                    zip([epicyclicFreq_real(gas_approx, x, scale) for x in r_g_dens], gas_dens)], 's-', label='$Q_{gas}$ & c=4')

plt.plot(r_g_dens, [1./Qs(epicycl=l[0], sigma=l[1], star_density=l[2]) for l in 
                    zip([epicyclicFreq_real(gas_approx, x, scale) for x in r_g_dens],
                        [sig_R_maj_max(_) for _ in r_g_dens], 
                        [surf_density(l_, M_to_L_R, 'R') for l_ in [mu_disc(ll, mu0=mu0d_R, h=h_disc_R) for ll in r_g_dens]])], 's-', label='$Q_{star}^{max}$')

plt.plot(r_g_dens, [1./Qs(epicycl=l[0], sigma=l[1], star_density=l[2]) for l in 
                    zip([epicyclicFreq_real(gas_approx, x, scale) for x in r_g_dens],
                        [sig_R_maj_min(_) for _ in r_g_dens], 
                        [surf_density(l_, M_to_L_R, 'R') for l_ in [mu_disc(ll, mu0=mu0d_R, h=h_disc_R) for ll in r_g_dens]])], 's-', label='$Q_{star}^{min}$')

plt.plot(map(lambda l: l/scale, r_q_n), q_n, 'o-', label='noordermeer 1F')

plt.axhline(y=1, ls='--')
plt.legend()
plot_SF(plt.gca())
plt.ylabel('$Q^{-1}$', fontsize=15);


    

In [60]:

    

total_gas_data = zip(r_g_dens, map(lambda l: He_coeff*l, gas_dens))[1:]
disk_scales = [(l[5], l[0].split(' ')[1]) for l in all_photometry]


fig = plt.figure(figsize=[10, 6])
ax = plt.gca()

plot_2f_vs_1f(ax=ax, total_gas_data=total_gas_data, epicycl=epicyclicFreq_real, 
              gas_approx=spl_gas, sound_vel=sound_vel, scale=scale, 
              sigma_max=sig_R_maj_max, 
              sigma_min=sig_R_maj_min, 
              star_density_max=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L_R, 'R'), 
              star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L_R, 'R'), 
              data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='R (not abs mag) maj max/min')

plt.ylim(0., 2.5)
plt.axhline(y=1., ls='-', color='grey')
plot_SF(ax)
plt.grid();


    

In [61]:

    

fig = plt.figure(figsize=[10, 6])
ax = plt.gca()

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

plot_2f_vs_1f(ax=ax, total_gas_data=total_gas_data, epicycl=epicyclicFreq_real, 
              gas_approx=spl_gas, sound_vel=sound_vel, scale=scale, 
              sigma_max=sig_R_maj_max, 
              sigma_min=sig_R_maj_min, 
              star_density_max=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L_R, 'R'), 
              star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L_R, 'R'), 
              data_lim=data_lim, color='m', alpha=0.3, disk_scales=disk_scales, label='R (not abs mag) maj max/min i=49deg')

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

plot_2f_vs_1f(ax=ax, total_gas_data=total_gas_data, epicycl=epicyclicFreq_real, 
              gas_approx=spl_gas, sound_vel=sound_vel, scale=scale, 
              sigma_max=sig_R_maj_max, 
              sigma_min=sig_R_maj_min, 
              star_density_max=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L_R, 'R'), 
              star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L_R, 'R'), 
              data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='R (not abs mag) maj max/min i=36deg')

plot_2f_vs_1f(ax=ax, total_gas_data=total_gas_data, epicycl=epicyclicFreq_real, 
              gas_approx=spl_gas, sound_vel=sound_vel, scale=scale, 
              sigma_max=sig_R_maj_max, 
              sigma_min=sig_R_maj_min, 
              star_density_max=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), 4., 'R'), 
              star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), 4., 'R'), 
              data_lim=data_lim, color='b', alpha=0.3, disk_scales=disk_scales, label='R (max disk) maj max/min')

plt.ylim(0., 2.5)
plt.axhline(y=1., ls='-', color='grey')
plot_SF(ax)
plt.grid();


    

In [62]:

    

from matplotlib.animation import FuncAnimation

fig = plt.gcf()
plt.figure(figsize=(10,6))

ax = plt.gca()

def animate(i):
    ax.cla()
    plot_2f_vs_1f(ax=ax, total_gas_data=total_gas_data, epicycl=epicyclicFreq_real, gas_approx=spl_gas, sound_vel=sound_vel, scale=scale, 
              sigma_max=sig_R_maj_max, 
              sigma_min=sig_R_maj_min, 
              star_density_max=all_photometry[i][-1], 
              star_density_min=all_photometry[i][-1],
              data_lim=data_lim, color=np.random.rand(3), alpha=0.3, disk_scales=[(all_photometry[i][5], '')])
    ax.axhline(y=1., ls='-', color='grey')
    ax.set_title(all_photometry[i][0])
    ax.set_ylim(0., 2.5)
    return ax
anim = FuncAnimation(plt.gcf(), animate, frames=len(all_photometry), interval=1000)


    

    




C:\Anaconda\lib\site-packages\matplotlib\axes\_axes.py:519: UserWarning: No labelled objects found. Use label='...' kwarg on individual plots.
  warnings.warn("No labelled objects found. "






    






<matplotlib.figure.Figure at 0xf4ef320>






    

In [63]:

    

anim.save('..\\..\pics\\'+name+'.gif', writer='imagemagick', fps=1)


    

In [64]:

    

from IPython.display import HTML
HTML(anim.to_html5_video())


    

    
Out[64]:

In [65]:

    

fig = plt.figure(figsize=[14,6])

plt.plot(test_points, spl_gas(test_points), '-', label='spline')
plt.plot(test_points, 0.85*spl_gas(test_points), '--', label='max disc')
# plt.plot(test_points, 0.6*spl_gas(test_points), '--', label='submax disc')
# plt.errorbar(r_wsrt, vel_wsrt, yerr=e_vel_wsrt, fmt='.', marker='.', mew=0, label = 'WSRT')
# plt.plot(r, vel, '.', label = 'Noord thesis')

max_coeffs = {}

for photom in all_photometry:
    disc_max = 2.2*photom[5]
    max_coeff = 0.85*spl_gas(disc_max)/disc_vel(disc_max, photom[7](0), photom[5], scale)
    submax_coeff = 0.6*spl_gas(disc_max)/disc_vel(disc_max, photom[7](0), photom[5], scale)
    
    plt.plot(test_points, map(lambda l: disc_vel(l, max_coeff**2 * photom[7](0), photom[5], scale), test_points), next(linecycler), label=photom[0] + '_MAX')
    
    print '{:15s}: M/L was {:2.2f} and for max it equal {:2.2f}, for submax equal {:2.2f}'.format(photom[0], photom[6], photom[6]*max_coeff**2, photom[6]*submax_coeff**2)
    max_coeffs[photom[0]] = [max_coeff**2, submax_coeff**2]


plt.ylim(0, 450)
plt.xlim(0, 300)
plt.legend(bbox_to_anchor=(1.15, 1.0));


    

    




Noorder R      : M/L was 2.71 and for max it equal 5.53, for submax equal 2.76
Noorder B      : M/L was 3.80 and for max it equal 7.93, for submax equal 3.95
Noorder R_max  : M/L was 4.00 and for max it equal 5.53, for submax equal 2.76
califa g (g-i) : M/L was 5.59 and for max it equal 6.80, for submax equal 3.39
califa r (g-r) : M/L was 3.88 and for max it equal 6.39, for submax equal 3.18
califa i (r-i) : M/L was 2.95 and for max it equal 5.52, for submax equal 2.75






    

In [66]:

    

from matplotlib.animation import FuncAnimation

fig = plt.gcf()
plt.figure(figsize=(10,6))

ax = plt.gca()

def animate(i):
    ax.cla()
    
    plot_2f_vs_1f(ax=ax, total_gas_data=zip(r_g_dens, map(lambda l: He_coeff*l, gas_dens))[1:],
              epicycl=epicyclicFreq_real, gas_approx=spl_gas, sound_vel=sound_vel, scale=scale, 
              sigma_max=sig_R_maj_max, 
              sigma_min=sig_R_maj_min, 
              star_density_max=lambda l: max_coeffs[all_photometry[i][0]][0]*all_photometry[i][-1](l), 
              star_density_min=lambda l: max_coeffs[all_photometry[i][0]][0]*all_photometry[i][-1](l),
              data_lim=data_lim, color=np.random.rand(3), alpha=0.2, disk_scales=[(all_photometry[i][5], '')])
    
    plot_2f_vs_1f(ax=ax, total_gas_data=zip(r_g_dens, map(lambda l: He_coeff*l, gas_dens))[1:], 
              epicycl=epicyclicFreq_real, gas_approx=spl_gas, sound_vel=sound_vel, scale=scale, 
              sigma_max=sig_R_maj_max, 
              sigma_min=sig_R_maj_min, 
              star_density_max=lambda l: max_coeffs[all_photometry[i][0]][1]*all_photometry[i][-1](l), 
              star_density_min=lambda l: max_coeffs[all_photometry[i][0]][1]*all_photometry[i][-1](l),
              data_lim=data_lim, color=np.random.rand(3), alpha=0.2, disk_scales=[(all_photometry[i][5], '')])
    
    ax.axhline(y=1., ls='-', color='grey')
    ax.set_title(all_photometry[i][0])
    ax.set_ylim(0., 2.5)
    ax.set_xlim(0., 160.)
    plot_SF(ax)
    return ax
anim = FuncAnimation(plt.gcf(), animate, frames=len(all_photometry), interval=1000)


    

    






<matplotlib.figure.Figure at 0xef7b978>






    

In [67]:

    

anim.save('..\\..\pics\\'+name+'_MAXDISCS.gif', writer='imagemagick', fps=1)


    

In [68]:

    

from IPython.display import HTML
HTML(anim.to_html5_video())


    

    
Out[68]:

In [69]:

    

fig = plt.figure(figsize=[10, 4])

def y_interp_(r, h):
    return 2.8*np.exp(-r/h)

points = np.linspace(0.1, 100., 100.)

plt.plot(r_g_dens, gas_dens, 's-')

plt.plot(points, y_interp_(points, 24.),  '--', alpha=0.5)
plt.plot(r_g_dens, map(lambda l: y_interp_(l[0], 24.) + l[1], zip(r_g_dens, gas_dens)),  '--', alpha=0.5);


    

    




C:\Anaconda\lib\site-packages\ipykernel\__main__.py:6: DeprecationWarning: object of type <type 'float'> cannot be safely interpreted as an integer.






    

In [70]:

    

fig = plt.figure(figsize=[10, 6])

total_gas_data_ = zip(r_g_dens, [He_coeff*(y_interp_(l[0], h_disc_R) + l[1]) for l in zip(r_g_dens, gas_dens)])

ax = plt.gca()

plot_2f_vs_1f(ax=ax, total_gas_data=total_gas_data_[1:], epicycl=epicyclicFreq_real, 
              gas_approx=spl_gas, sound_vel=sound_vel, scale=scale, 
              sigma_max=sig_R_maj_max, 
              sigma_min=sig_R_maj_min, 
              star_density_max=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), 4., 'R'), 
              star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), 4., 'R'), 
              data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='R (max disc) maj max/min')

plt.ylim(0., 2.5)
plt.axhline(y=1., ls='-', color='grey')
plot_SF(ax)
plt.grid();