In [1]:

    

from IPython.display import HTML
from IPython.display import Image
import os

%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 = 'N4258'
gtype = 'SABb' #LEDA, 'SBbc' from Heraudeau98
incl = 68.3  #LEDA, 56 здесь http://mnras.oxfordjournals.org/content/312/1/2.full.pdf
scale = 0.043 #kpc/arcsec according to NED

data_path = '../../../data/n4258_u7353'
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)

# Данные из NED
HTML('<iframe src=http://ned.ipac.caltech.edu/cgi-bin/objsearch?objname=ngc+4258&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[7]:

In [8]:

    

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


    

    
Out[8]:

In [9]:

    

#SDSS
Image('n4258_SDSS.jpg', width=300)


    

    
Out[9]:

In [10]:

    

Image('n4258_SDSS_labeled.jpg', width=500)


    

    
Out[10]:

In [11]:

    

#JHK
Image('n4258_2MASS.jpg', width=300)


    

    
Out[11]:

In [12]:

    

Image('http://www.spitzer.caltech.edu/uploaded_files/graphics/high_definition_graphics/0008/2016/ssc2007-06a_Ti.jpg', width=300)


    

    
Out[12]:

In [13]:

    

Image('http://www.spitzer.caltech.edu/uploaded_files/graphics/high_definition_graphics/0010/1167/sig14-020a_Inline.jpg', width=300)


    

    
Out[13]:

In [14]:

    

Image('http://www.spitzer.caltech.edu/uploaded_files/graphics/high_definition_graphics/0010/1185/sig14-020_Inline.jpg', width=300)


    

    
Out[14]:

In [15]:

    

Image('GCC_and_R25.png')


    

    
Out[15]:

In [16]:

    

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

fig = plt.figure(figsize=[8,5])
plt.errorbar(r_ma, vel_ma, e_vel_ma, fmt='.', marker='.', mew=0, label="Heraudeau 1998 stars maj")
plt.legend()
plt.ylim(-350., 350.);


    

In [17]:

    

Image('her98_rot.png') #оригинал


    

    
Out[17]:

In [18]:

    

r_ma_b, vel_ma_b, e_vel_b = zip(*sorted(zip(np.abs(r_ma), np.abs(vel_ma), e_vel_ma)))


    

In [19]:

    

# %%time
fig = plt.figure(figsize=[8,4])
plt.errorbar(r_ma_b, vel_ma_b, yerr=e_vel_b, fmt='.', marker='.', mew=0, color='blue', label = 'Her98 star maj')

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

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

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=10000., w=w(e_vel_b))
plt.plot(test_points, spl(test_points), '-', label='spl s=10000 w^2')

spl = inter.UnivariateSpline(r_ma_b, vel_ma_b, k=3, s=1800.)
plt.plot(test_points, spl(test_points), '-', label='spl s=1800')

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


    

In [20]:

    

star_approx = spl


    

In [21]:

    

r_sig_ma = r_ma #Heraudeau+1998

fig = plt.figure(figsize=[6., 4.])
plt.errorbar(r_sig_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='.', mew=0, color='blue', label=r'$maj\, Heraudeau $')

plt.xlabel('$R$')
plt.ylabel('$\sigma$')
plt.ylim(0, 350)
plt.legend();


    

In [22]:

    

Image('her98_sig.png') #из статьи


    

    
Out[22]:

In [23]:

    

r_sig_ma = map(abs, r_sig_ma)

fig = plt.figure(figsize=[6., 4.])
plt.errorbar(r_sig_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='.', mew=0, color='blue', label=r'$maj\, Heraudeau $')
plt.xlabel('$R$')
plt.ylabel('$\sigma$')
plt.ylim(0, 200)
plt.legend();


    

In [24]:

    

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_ma)+15., 100)


    

In [25]:

    

# TODO: find way to move to external file
@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))


    

In [26]:

    

def sigPhi_to_sigR_real(R):
        return 0.5 * (1 + R*star_approx.derivative()(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 [27]:

    

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,250)
plt.xlim(0,100);


    

In [28]:

    

Image('CO_rot.png') #from Sawada-Satoh 2007, правда [P.A.] = 86


    

    
Out[28]:

In [29]:

    

# Данные по кинематике газа Sawada-Satoh 2007 в CO
r_co, vel_co = zip(*np.loadtxt("CO_rot.dat", float, delimiter=','))
plt.plot(r_co, vel_co, 'o', label='CO Sawada-Satoh 2007')
plt.ylim(0, 250)
plt.xlim(0, 2250.)


    

    
Out[29]:





(0, 2250.0)






    

In [30]:

    

Image('HI_rot.png') #HI WHISP from van Eymeren 2011, PA=331(-180=151, все ок), R25=19.46kpc, разные цвета - это две стороны, черная - усредненная


    

    
Out[30]:

In [31]:

    

R25 = 19.46

# Данные по кинематике газа van Eymeren 2011 в HI
r_hi, vel_hi = zip(*np.loadtxt("HI_rot.dat", float, delimiter=','))
plt.plot(r_hi, vel_hi, 'd', label='HI Eymeren+2011')
plt.ylim(0, 250)
plt.xlim(-0.5, 2.)


    

    
Out[31]:





(-0.5, 2.0)






    

In [32]:

    

Image('CO_HI_rot.png') # CO+HI Sofue 97


    

    
Out[32]:

In [33]:

    

plt.plot([l/1000./scale for l in r_co], vel_co, '.', label='CO Sawada-Satoh 2007')
plt.plot([l*R25/scale for l in r_hi], vel_hi, 'd', label='HI Eymeren+2011')
plt.legend()


    

    
Out[33]:





<matplotlib.legend.Legend at 0xd3cc6d8>






    

In [34]:

    

r_hi = [l*R25/scale for l in r_hi]
r_co = [l/1000./scale for l in r_co]


    

In [35]:

    

fig = plt.figure(figsize=[10,6])
_1,_2, = [0.0,],[0.0,]
_1.extend(r_hi)
_2.extend(vel_hi)
_1,_2 = zip(*sorted(zip(_1,_2)))

gas_approx = poly1d(polyfit(_1, _2, deg=9))
test_points = np.linspace(0, max(r_hi), 100)
plt.plot(test_points, gas_approx(test_points), '--', label='poly approx')

spl_gas = inter.UnivariateSpline(_1, _2, k=3, s=2000.)
plt.plot(test_points, spl_gas(test_points), '-', label='spline')

plt.plot(r_hi, vel_hi, 'd', label='HI Eymeren+2011')

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


    

In [36]:

    

_1,_2 = zip(*filter(lambda l: l[0] > 100., zip(r_hi, vel_hi)))
_1 = list(r_co) + list(_1)
_2 = list(vel_co) + list(_2)

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

gas_approx = poly1d(polyfit(_1, _2, deg=9))
test_points = np.linspace(0, max(r_hi), 100)
plt.plot(test_points, gas_approx(test_points), '--', label='poly approx')

spl_gas = inter.UnivariateSpline(_1, _2, k=3, s=2000.)
plt.plot(test_points, spl_gas(test_points), '-', label='spline')

# plt.plot(r_co, vel_co, 'd', label='CO')
plt.plot(_1, _2, '.', label='CO+HI')

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


    

In [37]:

    

fig = plt.figure(figsize=[8, 4])
plt.plot(test_points, [epicyclicFreq_real(gas_approx, x, scale) for x in test_points], '-', label='poly')
plt.plot(test_points, [epicyclicFreq_real(spl_gas, x, scale) for x in test_points], '-', label='spline')
plt.xlabel('$R, arcsec$')
plt.ylabel('$\kappa,\, km/s/kpc$', fontsize=15)
plt.ylim(0, 1000)
plt.legend();


    

In [38]:

    

Image('u7353_gas_dens.png')


    

    
Out[38]:

In [39]:

    

r_g_dens, gas_dens = zip(*np.loadtxt("gas_data.dat", float, delimiter=','))

plt.semilogy([l/scale for l in r_g_dens], gas_dens, 's');


    

In [40]:

    

plt.semilogy([l/scale for l in r_g_dens[48:56]], gas_dens[48:56], 'o', color='r')
plt.semilogy([l/scale for l in r_g_dens[56:]], gas_dens[56:], 'o');


    

In [41]:

    

plt.plot([l/scale for l in r_g_dens[48:56]], gas_dens[48:56], '-o', color='r')
plt.plot([l/scale for l in r_g_dens[56:]], gas_dens[56:], '-o') #только сумма
plt.xlim(0, 200);


    

In [42]:

    

r_HI_dens, HI_dens = [l/scale for l in r_g_dens[:48]], gas_dens[:48]
r_mol_dens, mol_dens = [l/scale for l in r_g_dens[48:56]], gas_dens[48:56]
r_g_dens, gas_dens = [l/scale for l in r_g_dens[56:]], gas_dens[56:] #используем только полный газ


    

In [43]:

    

Image('VRI_photom_profiles.png', width=400)


    

    
Out[43]:

In [44]:

    

r_phot, mu_phot = zip(*np.loadtxt("VRI_photom_profiles.dat", float, delimiter=','))

fig = plt.figure(figsize=[8, 5])
plt.plot(map(lambda l: l/scale, r_phot), mu_phot, 's') 
plt.ylim(22, 14)
plt.xlim(0, 100.);


    

In [45]:

    

r_phot = map(lambda l: l/scale, r_phot) # но не факт, что там такое расстояние - т.к. 
# All distances and linear sizes have been calculated from the radial
# velocities listed in Table 2, assuming H0 =ˆ 55 km s21 Mpc21 and q0=0 , 
# скорость у нее 488 км/с (я прикинул - вроде похоже, у нас расст. 8.97 МПк и из Хаббла похоже)


    

In [46]:

    

r_photI, mu_photI = r_phot[:43], mu_phot[:43]
r_photR, mu_photR = r_phot[43:86], mu_phot[43:86]
r_photV, mu_photV = r_phot[86:], mu_phot[86:]

fig = plt.figure(figsize=[8, 5])
plt.plot(r_photI, mu_photI, 's')
plt.plot(r_photR, mu_photR, 's')
plt.plot(r_photV, mu_photV, 's')
plt.ylim(22, 14)
plt.xlim(0, 100.);


    

In [47]:

    

from scipy.optimize import curve_fit

p_ = np.arange(0.1, 100., 0.1)

def func(x, a, b):
    return a * x + b

popt, pcov = curve_fit(func, r_photV[25:], mu_photV[25:])

fig = plt.figure(figsize=[16, 5])
plt.subplot(121)
plt.plot(r_photV, mu_photV, 's', label='V')
plt.plot(r_photV, map(lambda l: func(l, *popt), r_photV), '-', label="Fitted Curve")
plt.ylim(22, 14)
plt.xlim(0, 100.)
plt.legend()

plt.subplot(122)
M_to_L_V = 1.74
plt.plot(p_, [surf_density(mu=m, M_to_L=M_to_L_V, band='V') for m in [mu_disc(l, mu0=popt[1], h=1./popt[0]) for l in p_]], '-', label='V [M/L={:2.2f}]'.format(M_to_L_V))
plt.title('mu0={:2.1f}mag, h_d={:2.0f}"'.format(popt[1], 1./popt[0]))
plt.legend();


    

In [48]:

    

popt, pcov = curve_fit(func, r_photI[25:], mu_photI[25:])

fig = plt.figure(figsize=[16, 5])
plt.subplot(121)
plt.plot(r_photI, mu_photI, 's', label='I')
plt.plot(r_photI, map(lambda l: func(l, *popt), r_photI), '-', label="Fitted Curve")
plt.ylim(22, 14)
plt.xlim(0, 100.)
plt.legend()

plt.subplot(122)
M_to_L_I = 0.97
plt.plot(p_, [surf_density(mu=m, M_to_L=M_to_L_I, band='I') for m in [mu_disc(l, mu0=popt[1], h=1./popt[0]) for l in p_]], '-', label='I [M/L={:2.2f}]'.format(M_to_L_I))
plt.title('mu0={:2.1f}mag, h_d={:2.0f}"'.format(popt[1], 1./popt[0]))
plt.legend();


    

In [49]:

    

all_photometry = []


    

In [50]:

    

# for I band
r_eff_I = 14.89
mu_eff_I = 18.48
n_I = 1.5
mu0d_I = 18.26
h_disc_I = 74.24


    

In [51]:

    

p_ = np.arange(0.1, 300., 0.1)

r_phot, mu_phot = zip(*np.loadtxt("photom_I.dat", float, delimiter=','))

fig = plt.figure(figsize=[8, 5])
plt.plot(map(abs,r_phot), mu_phot, 's', label='this research')
plt.plot(r_photI, mu_photI, 's', label='prev. research H=55')
plt.plot(map(lambda l:l*9./(488./75.), r_photI), mu_photI, 's', label='prev. research H=75')
plt.ylim(23, 14)
plt.xlim(0, 100.)
plt.legend();


    

In [52]:

    

bulge_I = [mu_bulge(l, mu_eff=mu_eff_I, r_eff=r_eff_I, n=n_I) for l in p_]
disc_I = [mu_disc(l, mu0=mu0d_I, h=h_disc_I) for l in p_]
total_profile = total_mu_profile(bulge_I, disc_I)

fig = plt.figure(figsize=[8, 5])
plt.plot(r_phot, mu_phot, 's')
plt.plot(p_, bulge_I, '-', label='bulge', color='red')
plt.plot(p_, disc_I, '-', label='disk', color='green')
plt.plot(p_, total_profile, '-', label='sum', color='m')
plt.ylim(23, 14)
plt.legend();


    

In [53]:

    

M_to_L_I = 0.97
surf_I = [surf_density(mu=m, M_to_L=M_to_L_I, band='I') for m in disc_I]

plt.plot(p_, surf_I, '-', label='YOSHINO I [M/L={:2.2f}]'.format(M_to_L_I))
plt.legend()


    

    
Out[53]:





<matplotlib.legend.Legend at 0xf3a3978>






    

In [54]:

    

# for V band
r_eff_V = 10.34
mu_eff_V = 18.57
n_V = 1.5
mu0d_V = 19.54
h_disc_V = 78.55
M_to_L_V = 1.74
disc_V = [mu_disc(l, mu0=mu0d_V, h=h_disc_V) for l in p_]


    

In [55]:

    

# for J band
r_eff_J = 7.31
mu_eff_J = 16.43
n_J = 1.5
mu0d_J = 16.37
h_disc_J = 38.95
M_to_L_J = 0.95
disc_J = [mu_disc(l, mu0=mu0d_J, h=h_disc_J) for l in p_]


    

In [56]:

    

surf_V = [surf_density(mu=m, M_to_L=M_to_L_V, band='V') for m in disc_V]
surf_J = [surf_density(mu=m, M_to_L=M_to_L_J, band='J') for m in disc_J]

plt.plot(p_, surf_I, '-', label='YOSHINO I [M/L={:2.2f}]'.format(M_to_L_I))
plt.plot(p_, surf_V, '-', label='YOSHINO V [M/L={:2.2f}]'.format(M_to_L_V))
plt.plot(p_, surf_J, '-', label='YOSHINO J [M/L={:2.2f}]'.format(M_to_L_J))
plt.ylabel('$M_{sun}/{pc}^2$', fontsize=15.)
plt.legend()


    

    
Out[56]:





<matplotlib.legend.Legend at 0xd9273c8>






    

In [57]:

    

all_photometry.append(('YOSHINO I', r_eff_I, mu_eff_I, n_I, mu0d_I, h_disc_I, M_to_L_I, 
                       lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_I, h=h_disc_I), M_to_L=M_to_L_I, band='I')))

all_photometry.append(('YOSHINO V', r_eff_V, mu_eff_V, n_V, mu0d_V, h_disc_V, M_to_L_V, 
                       lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_V, h=h_disc_V), M_to_L=M_to_L_V, band='V')))

all_photometry.append(('YOSHINO J', r_eff_J, mu_eff_J, n_J, mu0d_J, h_disc_J, M_to_L_J, 
                       lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_J, h=h_disc_J), M_to_L=M_to_L_J, band='J')))


    

In [58]:

    

# for 3.6 band
dist_36 =  8.17 #Mpc
r_eff_36 = np.power(10., 2.77)/scale/1000./(dist_36/8.97) #pc
mu_eff_36 = 17.41
n_36 = 2.8
mu0d_36 = 18.82
h_disc_36 = np.power(10., 3.5)/scale/1000./(dist_36/8.97) #because original is Log(), pc


    

In [59]:

    

M_to_L_s4g = s4g_mass_to_light(-21.328, -20.881)
M_to_L_s4g


    

    
Out[59]:





0.65393261738835051

In [60]:

    

surf_36 = [s4g_surf_density(m, M_to_L=M_to_L_s4g) for m in [mu_disc(l, mu0=mu0d_36, h=h_disc_36) for l in p_]]

plt.plot(p_, surf_36, '-', label='3.6 [M/L={:2.2f}]'.format(0.6))
plt.legend()


    

    
Out[60]:





<matplotlib.legend.Legend at 0x10419278>






    

In [61]:

    

bv_color = 0.77
mcgaugh_mass_to_light(bv_color, 'mu36')


    

    
Out[61]:





array([ 0.47,  0.58,  0.7 ,  0.62])

In [62]:

    

all_photometry.append(('SPITZER 3.6', r_eff_36, mu_eff_36, n_36, mu0d_36, h_disc_36, M_to_L_s4g, 
                       lambda l: s4g_surf_density(mu_disc(l, mu0=mu0d_36, h=h_disc_36), M_to_L=M_to_L_s4g)))

print mu_face_on(mu0d_36, cos_i)

all_photometry.append(('SPITZER 3.6 faceon', r_eff_36, mu_eff_36, n_36, mu_face_on(mu0d_36, cos_i), h_disc_36, M_to_L_s4g, 
                       lambda l: s4g_surf_density(mu_disc(l, mu0=mu_face_on(mu0d_36, cos_i), h=h_disc_36), M_to_L=M_to_L_s4g)))


    

    




19.9002390653

In [63]:

    

Image('https://s3.amazonaws.com/aasie/images/0004-637X/780/1/69/apj488463f25_hr.jpg', width=400)


    

    
Out[63]:

In [64]:

    

# используем вторую модель, т.к. она лучше; 
# также надо пересчитывать из светимости яркость - возьму поверхн яркость с картинки - примерно 17.8
m_b_K = 9.24 
r_eff_K = 6.27
n_K = 3.26 
m_d_K = 6.00
h_disc_K = 146.


    

In [65]:

    

M_to_L_K = bell_mass_to_light(bv_color, 'K', 'B-V') #оценку берем из пред. работы
M_to_L_K


    

    
Out[65]:





0.79058760299863684

In [66]:

    

surf_K = [surf_density(mu=m, M_to_L=M_to_L_K, band='K') for m in [mu_disc(l, mu0=17.8, h=h_disc_K) for l in p_]]

plt.plot(p_, surf_K, '-', label='K [M/L={:2.2f}]'.format(M_to_L_K))
plt.legend()


    

    
Out[66]:





<matplotlib.legend.Legend at 0xd527b38>






    

In [67]:

    

all_photometry.append(('GALFIT K', r_eff_K, None, n_K, 17.8, h_disc_K, M_to_L_K, 
                       lambda l: surf_density(mu=mu_disc(l, mu0=17.8, h=h_disc_K), M_to_L=M_to_L_K, band='K')))


    

In [68]:

    

h_disc_s4g = 238.3528*0.75
mu0d_s4g =  3.24 - (log10((9.9858 + 5*log10(h_disc_s4g*scale) + 38.57)/4.255) - 8)/0.4
h_disc_s4g, mu0d_s4g


    

    
Out[68]:





(178.7646, 20.50187578788826)

In [69]:

    

p_ = np.arange(0.1, 200., 0.1)

surf_s4g = [s4g_surf_density(mu_disc(l, mu0=mu0d_s4g, h=h_disc_s4g), M_to_L_s4g) for l in p_]
plt.plot(p_, surf_s4g, '-', label='S4G [M/L={:2.2f}]'.format(M_to_L_s4g))
plt.legend();


    

In [70]:

    

all_photometry.append(('S4G 3.6', None, None, None, mu0d_s4g, h_disc_s4g, M_to_L_s4g, 
                       lambda l: s4g_surf_density(mu_disc(l, mu0=mu0d_s4g, h=h_disc_s4g), M_to_L_s4g)))


    

In [71]:

    

for photom in all_photometry:
    plt.plot(p_, map(photom[-1], p_), '-', label=photom[0].format(M_to_L_I))
plt.ylabel('$M_{sun}/{pc}^2$', fontsize=15.)
plt.legend(loc='best')
plt.legend()


    

    
Out[71]:





<matplotlib.legend.Legend at 0xd9027f0>






    

In [72]:

    

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 | YOSHINO I          |   14.89 |    18.48 |   1.50 |   18.26 |    74.24 |  0.97 | 5.62E+10.   |       878 |
| 1.00 | YOSHINO V          |   10.34 |    18.57 |   1.50 |   19.54 |    78.55 |  1.74 | 6.93E+10.   |       967 |
| 2.00 | YOSHINO J          |    7.31 |    16.43 |   1.50 |   16.37 |    38.95 |  0.95 | 5.76E+10.   |      3271 |
| 3.00 | SPITZER 3.6        |   15.03 |    17.41 |   2.80 |   18.82 |    80.74 |  0.65 | 8.45E+10.   |      1116 |
| 4.00 | SPITZER 3.6 faceon |   15.03 |    17.41 |   2.80 |   19.90 |    80.74 |  0.65 | 3.12E+10.   |       413 |
| 5.00 | GALFIT K           |    6.27 |   nan    |   3.26 |   17.80 |   146.00 |  0.79 | 1.30E+11.   |       523 |
| 6.00 | S4G 3.6            |  nan    |   nan    | nan    |   20.50 |   178.76 |  0.65 | 8.80E+10.   |       237 |
+------+--------------------+---------+----------+--------+---------+----------+-------+-------------+-----------+

In [73]:

    

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

plt.plot(test_points, spl_gas(test_points), '-', label='spline')

plt.plot(_1, _2, '.', label='CO+HI')

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, 300)
plt.xlim(0, 300)
plt.legend(loc='best');


    

In [74]:

    

Image('halpha.png', width=1700) #две картинки - Halpha и UV


    

    
Out[74]:

In [75]:

    

Image('halpha_dist.png', width=600)


    

    
Out[75]:

In [76]:

    

Image('large-fg16_19.jpeg')


    

    
Out[76]:

In [77]:

    

Image('large-fg16_19_dist.png', width=500)


    

    
Out[77]:

In [78]:

    

def plot_SF(ax):
    ax.plot([0., 86./2/(175./600.)], [0., 0.], '-', lw=7., color='midnightblue')
#     ax.plot([0., 123./2/(175./600.)], [0., 0.], '-', lw=7., color='midnightblue') # глазами видно, что там больше, но там именно спирали маленькие
    ax.plot([0., 150*58./180.], [0., 0.], '-', lw=7., color='r')
    ax.plot([250./2/(175./600.), 260./2/(175./600.)], [0., 0.], '-', lw=7., color='b') #внешняя спираль
    ax.plot([220./2/(175./600.), 310./2/(175./600.)], [0., 0.], '-', lw=7., color='b') #внешняя спираль, на глаз
    
plot_SF(plt.gca())
plt.xlim(0, 650)
plt.ylim(0, 200)


    

    
Out[78]:





(0, 200)






    

In [79]:

    

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


    

In [80]:

    

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],
                        map(sig_R_maj_max, r_g_dens), 
                        [surf_density(l_, M_to_L_I, 'I') for l_ in [mu_disc(ll, mu0=mu0d_I, h=h_disc_I) 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],
                        map(sig_R_maj_min, r_g_dens), 
                        [surf_density(l_, M_to_L_I, 'I') for l_ in [mu_disc(ll, mu0=mu0d_I, h=h_disc_I) for ll in r_g_dens]])], 's-', label='$Q_{star}^{min}$')

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


    

In [81]:

    

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

plt.plot(r_g_dens[:-14], [1./Qg(epicycl=l[0], sound_vel=sound_vel, gas_density=l[1]) for l in 
                    zip([epicyclicFreq_real(spl_gas, x, scale) for x in r_g_dens[:-14]], gas_dens[:-14])], 's-', label='$Q_{g}^{-1}$')

plt.axhline(y=1, ls='--')
plt.legend()
plot_SF(plt.gca())
plt.ylabel('$Q^{-1}$', fontsize=15)
plt.title(name)
plt.xlim(0, 650)
plot_Q_levels(plt.gca(), [1., 1.5, 2., 3.])
plt.savefig('..\\..\pics\\instab_spirals\\'+name+'_spiral'+'.png', format='png', bbox_inches='tight');


    

In [82]:

    

total_gas_data = zip(r_g_dens, map(lambda l: l, gas_dens))[:20]
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_J, h=h_disc_J), M_to_L_J, 'J'), 
              star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_I, h=h_disc_I), M_to_L_I, 'I'), 
              data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='J/I maj max/min')

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_I, h=h_disc_I), M_to_L_I, 'I'), 
              star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_I, h=h_disc_I), M_to_L_I, 'I'), 
              data_lim=data_lim, color='m', alpha=0.3, disk_scales=disk_scales, label='I maj max/min')

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_V, h=h_disc_V), M_to_L_V, 'V'), 
                  star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_V, h=h_disc_V), M_to_L_V, 'V'), 
                  data_lim=data_lim, color='r', alpha=0.2, disk_scales=disk_scales, label='V maj max/min')

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


    

In [83]:

    

gas_data = zip(r_g_dens, gas_dens)[:7]
invQg, invQs, invQeff = zip(*get_invQeff_from_data(gas_data=gas_data, 
                                epicycl=epicyclicFreq_real, 
                                gas_approx=spl_gas,
                                sound_vel=sound_vel, 
                                scale=scale,
                                sigma=sig_R_maj_min,
                                star_density=lambda l: surf_density(mu_disc(l, mu0=mu0d_I, h=h_disc_I), M_to_L_I, 'I'), verbose=True))


    

    




r = 9.88; gas_d = 170.26; epicycl = 828.67; sig = 119.96; star_d = 768.96
	Qs = 8.91; Qg = 2.15; Qeff = 2.10
r = 28.64; gas_d = 90.39; epicycl = 272.92; sig = 109.81; star_d = 597.24
	Qs = 3.46; Qg = 1.33; Qeff = 1.28
r = 45.43; gas_d = 35.24; epicycl = 212.11; sig = 93.27; star_d = 476.37
	Qs = 2.86; Qg = 2.66; Qeff = 2.36
r = 63.20; gas_d = 16.97; epicycl = 91.21; sig = 81.75; star_d = 374.95
	Qs = 1.37; Qg = 2.38; Qeff = 1.19
r = 80.98; gas_d = 14.67; epicycl = 82.19; sig = 81.75; star_d = 295.11
	Qs = 1.57; Qg = 2.48; Qeff = 1.35
r = 98.75; gas_d = 10.43; epicycl = 75.31; sig = 81.75; star_d = 232.28
	Qs = 1.83; Qg = 3.19; Qeff = 1.58
r = 116.53; gas_d = 10.96; epicycl = 63.68; sig = 81.75; star_d = 182.82
	Qs = 1.96; Qg = 2.57; Qeff = 1.66

In [84]:

    

from matplotlib.animation import FuncAnimation

total_gas_data = zip(r_g_dens, map(lambda l: l, gas_dens))[:20]

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)
    ax.set_xlim(0., 265.)
    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 0xf03da20>






    

In [85]:

    

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


    

In [86]:

    

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


    

    
Out[86]:

In [87]:

    

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')
plt.plot(_1, _2, '.', label='CO+HI')

max_coeffs = {}

for ind, photom in enumerate(all_photometry):
    if type(photom[5]) == tuple:
        disc_v = lambda l: disc_vel(l, photom[7][0](0), photom[5][0], scale, Sigma0_2=photom[7][1](0), h_2=photom[5][1])
    else:
        disc_v = lambda l: disc_vel(l, photom[7](0), photom[5], scale)
    
    max_coeff = 100.
    submax_coeff = 100.
    for d1, d2 in zip(_1, _2)[20:]:
        max_coeff = min(0.85*d2/disc_v(d1), max_coeff)
        submax_coeff = min(0.6*d2/disc_v(d1), submax_coeff)
        
              
#     values = map(disc_v, test_points[1:])
#     disc_max = test_points[values.index(max(values))]
        
#     max_coeff = 0.85*spl_gas(disc_max)/disc_v(disc_max)
#     submax_coeff = 0.6*spl_gas(disc_max)/disc_v(disc_max)
    
    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 '{:25s}: 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, 300)
plt.xlim(0, 500)
plt.legend(bbox_to_anchor=(1.15, 1.0));


    

    




YOSHINO I                : M/L was 0.97 and for max it equal 0.82, for submax equal 0.41
YOSHINO V                : M/L was 1.74 and for max it equal 1.22, for submax equal 0.61
YOSHINO J                : M/L was 0.95 and for max it equal 0.53, for submax equal 0.27
SPITZER 3.6              : M/L was 0.65 and for max it equal 0.38, for submax equal 0.19
SPITZER 3.6 faceon       : M/L was 0.65 and for max it equal 1.03, for submax equal 0.51
GALFIT K                 : M/L was 0.79 and for max it equal 0.51, for submax equal 0.26
S4G 3.6                  : M/L was 0.65 and for max it equal 0.80, for submax equal 0.40






    

In [88]:

    

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=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=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: 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 0xf2ea828>






    

In [89]:

    

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