Галактика из диплома.
In [1]:
from IPython.display import HTML
from IPython.display import Image
import os
%pylab
%matplotlib inline
%run ../../../utils/load_notebook.py
In [2]:
from photometry import *
In [3]:
from instabilities import *
In [4]:
from utils import *
In [5]:
name = 'N3898'
gtype = 'SA(s)ab'
incl = 70. #(adopted by Epinat+2008)
scale = 0.092 #kpc/arcsec according to ApJ 142 145(31pp) 2011
data_path = '../../../data/ngc3898'
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+3898&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=ngc3898 width=1000 height=350></iframe>')
Out[8]:
In [9]:
#SDSS
Image('ngc3898_SDSS.jpeg', width=300)
Out[9]:
Картинка из выборки http://cosmo.nyu.edu/hogg/rc3/ и там есть с маштабом изображение:
In [10]:
Image('ngc3898_SDSS_labeled.jpeg', width=500)
Out[10]:
In [11]:
#JHK
Image('ngc3898_JHK.jpg', width=300)
Out[11]:
Есть пара картинок от HST, например вот, но не очень хороших.
Noordermeer, thesis, p. 113
UGC 6787 (NGC 3898) suffers from the same problem as UGC 6787: the luminous bulge distorts the isophotes of the disk out to large radii and the standard method to derive the inclination (section 3.5.1) cannot be used. Only after the model image of the bulge was subtracted could the disk ellipticity be determined. The final disk photometric profile and the intrinsic axis ratio of the bulge were derived using an inclination of $61^\circ$ . Note that this is still somewhat lower than the kinematic inclination derived from the HI velocity field (chapter 4)
Noordermeer, thesis, pp. 188-189
The rotation curve of UGC 6787 (NGC 3898) is well resolved and shows some characteristic ‘wiggles’ with an amplitude of 30 – 50 km/s. The kinematics in the central parts are only barely resolved in the optical spectrum, and due to the high inclination angle and resulting line-of-sight integration effects, the central line-profiles are strongly broadened. After the initial rise of the rotation curve to the peak velocity of 270 km/s, the rotation velocities drop to approximately 220 km/s at a radius of 3000 (∼ 2.75 kpc), after which they gradually rise again to 250 km/s at R ≈ 10000 (∼9 kpc). The rotation velocities then drop again to 220 km/s, after which they rise again to reach a more or less flat plateau at 250 km/s. Although there are clear indications that the gas disk of UGC 6787 is warped, the locations of the ‘wiggles’ in the rotation curve do not coincide with the radii where the position angle and the inclination change and the variations in the rotation velocity seem to be real. This is further confirmed by the fact that the variations occur symmetrically at all position angles over the velocity field. The discrepancy between the kinematic inclination angle derived here and the optical inclination from chapter 3 can be explained by the dominance of the bulge in the optical image. As was noted in chapter 3, the optical image of this galaxy is dominated by the spheroidal bulge out to large radii, which makes it impossible to obtain a reliable estimate for the inclination. The kinematical inclination derived here is free of such effects and thus reflects the true orientation of this galaxy more accurately.
In [12]:
Image('noord_pics/HI_data.png')
Out[12]:
In [13]:
Image('noord_pics/photometry.png')
Out[13]:
In [14]:
Image('noord_pics/rot_vel_all.png')
Out[14]:
In [15]:
#точки, снятые с предыдущей картинки с помощью WebPlotDigitizer
r_ma, vel_ma = zip(*np.loadtxt("v_gas_from_webplotdigitizer.dat", float, delimiter=','))
plt.plot(r_ma, vel_ma, '--')
plt.xlim(0, 100)
plt.ylim(0, 500);
In [16]:
Image('noord_pics/rot_vel_HI.png')
Out[16]:
In [17]:
Image('noord_pics/rot_vel_center.png', width=300)
Out[17]:
In [18]:
Image('noord_pics/rot_vel_large.png', width=300)
Out[18]:
Снимем отсюда точки - уже по наблюдательным данным (отдельно в центре и отдельно на периферии):
In [19]:
r_ma2, vel_ma2 = zip(*np.loadtxt("v_gas_from_webplotdigitizer2.dat", float, delimiter=','))
plt.plot(r_ma2, vel_ma2, '--')
plt.plot(map(lambda l: l/scale, r_ma), vel_ma, '--')
plt.xlim(0, 400)
plt.ylim(0, 320);
Видно, что снял я достаточно точно (с двух картинок данные совпали).
In [20]:
Image('noord_pics/isothermal_halo.png')
Out[20]:
In [21]:
Image('noord_pics/NFW.png')
Out[21]:
Ценное - одножидкостный критерий на сравнение:
In [22]:
Image('noord_pics/instab.png')
Out[22]:
In [23]:
r_q_n, q_n = zip(*np.loadtxt("noord_pics/Q1F.dat", float, delimiter=','))
fig = plt.figure(figsize=[8,2])
plt.plot(r_q_n, q_n, 'o-')
plt.xlim(0, 30)
plt.ylim(0, 1.2);
Наконец последнее - все табличные данные из диссертации:
In [24]:
Image('noord_pics/tables.png')
Out[24]:
TODO: подчистить данные из Noordermeera, убрать ненужное
TODO: добавить pignatelli
In [25]:
# Данные по звездной кинематике Noordermeer вдоль большей полуоси (не исправленные за наклон?)
r_ma3, vel_ma3, e_vel_ma3, sig_ma3, e_sig_ma3 = zip(*np.loadtxt("v_stars_ma.dat", float))
# Данные по звездной кинематике Heraudeau+1999 вдоль большой полуоси (не исправленные за наклон?)
r_ma4, vel_ma4, e_vel_ma4, sig_ma4, e_sig_ma4 = zip(*np.loadtxt("v_stars_HSM.dat", float))
fig = plt.figure(figsize=[10,8])
plt.errorbar(r_ma3, vel_ma3, e_vel_ma3, fmt='-.', marker='.', mew=0, label="Noord stars maj")
plt.errorbar(r_ma4, vel_ma4, e_vel_ma4, fmt='-.', marker='.', mew=0, label="Heraudeau stars maj")
plt.legend()
plt.show()
TODO: проверить, как хорошо были сняты звездные данные (обе кривые)
Определенно надо было минусы проставить в Heraudeau. Перегнем и сравним со снятыми данными:
In [26]:
r_ma3_, vel_ma3_, e_vel_ma3_ = zip(*sorted(zip(np.abs(r_ma3), np.abs(vel_ma3), e_vel_ma3)))
r_ma4_, vel_ma4_, e_vel_ma4_ = zip(*sorted(zip(np.abs(r_ma4), np.abs(vel_ma4), e_vel_ma4)))
In [27]:
fig = plt.figure(figsize=[10,8])
plt.errorbar(r_ma3_, vel_ma3_, e_vel_ma3_, fmt='-.', marker='.', mew=0, label="Noord stars maj")
plt.errorbar(r_ma4_, vel_ma4_, e_vel_ma4_, fmt='-.', marker='.', mew=0, label="Heraudeau stars maj")
plt.plot(r_ma2, vel_ma2, '--')
plt.plot(map(lambda l: l/scale, r_ma), vel_ma, '--')
plt.xlim(0, 400)
plt.ylim(0, 320)
plt.legend()
plt.xlabel('R, arcsec')
plt.show()
Видно, что в статье 2008 года кривая протянута приемлимо далеко, однако газ имеет сильно большие значения, чем звездная кривая вращения.
In [28]:
incl_corr = lambda l: l/sin_i
r_ma_b, vel_ma_b, e_vel_b = r_ma3_, map(incl_corr, vel_ma3_), map(incl_corr, e_vel_ma3_)
In [29]:
%%time
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 = 'Noord star maj')
test_points = np.linspace(0.0, max(r_ma_b), 100)
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=500.)
plt.plot(test_points, spl(test_points), '-', label='spl s=500')
spl1 = inter.UnivariateSpline(r_ma_b, vel_ma_b, k=3, s=10000., w=w(e_vel_b))
plt.plot(test_points, spl1(test_points), '-', label='spl s=500 w^2')
plt.legend(loc='lower right')
plt.ylim(0, 250)
plt.show()
C весами плохо получается, берем обычный достаточно гладкий сплайн.
In [30]:
star_approx = spl
Дисперсии из Pignatelli http://mnras.oxfordjournals.org/content/323/1/188.full.pdf:
In [31]:
Image('pignatelli_2001_kinem.png', width=400)
Out[31]:
In [32]:
r_p, v_p, dv_p, sig_p, dsig_p, _, _, _, _ = zip(*np.loadtxt("pignatelli_2001_kinem.dat", float))
fig = plt.figure(figsize=[8, 3])
plt.errorbar(r_p, sig_p, yerr=dsig_p, fmt='.', marker='.', mew=0, color='blue', label=r'$\sigma_{los}^{maj} Pignatelli$')
plt.xlabel('$R$')
plt.ylabel('$\sigma$')
plt.xlim(-40., 40.)
plt.ylim(0, 320)
plt.legend();
In [33]:
# Исправляем значения вдоль малой оси на синус угла:
def correct_min(R):
return R / cos_i
r_sig_ma, sig_ma, e_sig_ma = zip(*np.loadtxt("s_stars_maN.dat", float))
r_sig_mi, sig_mi, e_sig_mi = zip(*np.loadtxt("s_stars_miN.dat", float))
fig = plt.figure(figsize=[8, 8])
plt.errorbar(r_sig_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='.', mew=0, color='blue', label=r'$\sigma_{los}^{maj} Noord$')
plt.errorbar(map(abs, r_ma4), sig_ma4, yerr=e_sig_ma4, fmt='.', marker='.', mew=0, color='red', label=r'$maj Heraudeau $')
plt.errorbar(map(abs, r_p), sig_p, yerr=dsig_p, fmt='.', marker='.', mew=0, color='m', label=r'$\sigma_{los}^{maj} Pignatelli$')
r_sig_mi = map(correct_min, r_sig_mi)
plt.errorbar(r_sig_mi, sig_mi, yerr=e_sig_mi, fmt='.', marker='.', mew=0, color='black', label='$\sigma_{los}^{min} Noord$')
plt.xlabel('$R$')
plt.ylabel('$\sigma$')
plt.legend();
TODO: с малой осью явно что-то не то
In [34]:
spl_min = inter.UnivariateSpline(r_sig_mi, sig_mi, k=3, s=10000.)
sig_min_lim = max(r_sig_mi)
spl_maj = inter.UnivariateSpline(r_sig_ma, sig_ma, k=3, s=100.)
sig_maj_lim = max(r_sig_ma)
points = np.linspace(0.1, max(r_sig_ma)+15., 100)
In [35]:
@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 [36]:
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);
Используем соотношение $\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)
In [37]:
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 [38]:
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);
Сравним major vs minor оценки:
In [39]:
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,300)
plt.xlim(0,100);
"Настоящая" оценка сверху на этот раз имеет какую-то особенность.
Газовая кривая, нужна для эпициклического приближения:
In [40]:
fig = plt.figure(figsize=[10,8])
# Данные по кинематике WSRT
r_wsrt, vel_wsrt, e_vel_wsrt = zip(*np.loadtxt("v_gas_WSRT.dat", float))
# Данные по кинематике газа Epinat+2008 в Halpha
r_ha, _,_,_, vel_ha, e_vel_ha, _,_ = zip(*np.loadtxt("v_gasHa.dat", str))
r_ha, vel_ha, e_vel_ha = np.array(r_ha, dtype='float'), np.array(vel_ha, dtype='float'), np.array(e_vel_ha, dtype='float')
plt.plot(r_ma2, vel_ma2, '--', label='Noord thesis #1')
plt.plot(map(lambda l: l/scale, r_ma), vel_ma, '--', label='Noord thesis #2')
plt.errorbar(r_wsrt, vel_wsrt, yerr=e_vel_wsrt, fmt='.', marker='.', mew=0, label = 'WSRT')
plt.errorbar(r_ha, vel_ha, yerr=e_vel_ha, fmt='.', marker='d', mew=0, label = 'Halpha')
plt.ylim(0, 300)
plt.legend(loc='lower right');
Согласие достаточно хорошее, приблизим:
In [41]:
fig = plt.figure(figsize=[10,6])
gas_approx = poly1d(polyfit(r_wsrt, vel_wsrt, deg=15))
test_points = np.linspace(min(r_wsrt), max(r_wsrt), 100)
plt.plot(test_points, gas_approx(test_points), '--', label='poly approx')
spl_gas = inter.UnivariateSpline(r_wsrt, vel_wsrt, k=3, s=2400.)
plt.plot(test_points, spl_gas(test_points), '-', label='spline')
plt.errorbar(r_wsrt, vel_wsrt, yerr=e_vel_wsrt, fmt='.', marker='.', mew=0, label = 'WSRT')
plt.ylim(0, 300)
plt.legend(loc='lower right');
Для случая бесконечного тонкого диска: $$\kappa=\frac{3}{R}\frac{d\Phi}{dR}+\frac{d^2\Phi}{dR^2}$$ где $\Phi$ - гравпотенциал, однако его знать не надо, т.к. есть проще формула: $$\kappa=\sqrt{2}\frac{\vartheta_c}{R}\sqrt{1+\frac{R}{\vartheta_c}\frac{d\vartheta_c}{dR}}$$
In [42]:
fig = plt.figure(figsize=[12, 8])
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()
plt.show()
Плотность HI:
In [43]:
r_g_dens, gas_dens = zip(*np.loadtxt("gas_density.dat", float))
plt.plot(r_g_dens, gas_dens, 's-');
В этой давней работе http://adsabs.harvard.edu/abs/1987ApJ...322...64S есть указание что плотность постоянна, равна $CO$ 1.74 на протяжении примерно 30 секунд от центра.
В этой работе http://adsabs.harvard.edu/abs/2014A%26A...564A..65B (из трех частей) вроде как есть в последней части точка, в первой работе написано что она измерялась, в таблице интенсивности для нее не даны и надо вычислять по формуле $I(CO) = 5\sigma(W_{HI}\delta V_{CO})^{1/2}K km\, s−1$
Насколько я нашел все параметры, интенсивность $I(CO)$ получается
In [44]:
5*2.00*sqrt(463*15.7)/1000.
Out[44]:
А центральные плотности
In [45]:
5*2.00*sqrt(463*15.7)/1000. * 3.2
Out[45]:
По порядкам адекватно.
TODO: добавить
In [46]:
all_photometry = []
Фотометрий четыре - одна в J, две в R и одна в B.
In [47]:
Image('diplom_cite_p38.png')
Out[47]:
In [48]:
from wand.image import Image as WImage
img = WImage(filename='ngc3898.pdf', resolution=200)
img[:, 150:1800]
Out[48]:
In [49]:
img = WImage(filename='ngc3898.pdf[1]', resolution=200)
img[:, 150:1300]
Out[49]:
На этой странице ошибка - фотометрия не в J, она в R!
In [50]:
Image('noord_pics/photom_table.png')
Out[50]:
In [51]:
Image('photom_table2.png')
Out[51]:
In [52]:
Image('noord_pics/n3898_photom.png')
Out[52]:
In [53]:
r_phot, mu_phot = zip(*np.loadtxt("noord_pics/n3898_noord_photoRB.dat", float, delimiter=','))
plt.plot(r_phot, mu_phot, 's')
plt.xlim(-10, 250)
plt.ylim(30, 15);
Посмотрим, как ведет себя ноордермееровская фотометрия в $R$:
In [54]:
# for R band
r_eff_R = 8.8
mu_eff_R = 18.37
n_R = 2.3
mu0d_R = 20.49
h_disc_R = 36.2
In [55]:
p_ = np.arange(0.1, 220., 0.1)
bulge = [mu_bulge(l, mu_eff=mu_eff_R, r_eff=r_eff_R, n=n_R) for l in p_]
disc = [mu_disc(l, mu0=mu0d_R, h=h_disc_R) for l in p_]
total_profile = total_mu_profile(bulge, disc)
fig = plt.figure(figsize=[8, 5])
plt.plot(r_phot, mu_phot, 's')
plt.plot(p_, bulge, '-', label='bulge', color='red')
plt.plot(p_, disc, '-', label='disk', color='green')
plt.plot(p_, total_profile, '-', label='sum', color='m')
plt.xlim(-10, 175)
plt.ylim(30, 15)
plt.legend();
В $B$:
In [56]:
# for B band
r_eff_B = 8.8
mu_eff_B = 19.80
n_B = 2.3
mu0d_B = 22.0
h_disc_B = 42.9
In [57]:
bulge_B = [mu_bulge(l, mu_eff=mu_eff_B, r_eff=r_eff_B, n=n_B) for l in p_]
disc_B = [mu_disc(l, mu0=mu0d_B, h=h_disc_B) for l in p_]
total_profile_B = total_mu_profile(bulge_B, disc_B)
fig = plt.figure(figsize=[8, 5])
plt.plot(r_phot, mu_phot, 's')
plt.plot(p_, bulge_B, '-', label='bulge', color='red')
plt.plot(p_, disc_B, '-', label='disk', color='green')
plt.plot(p_, total_profile_B, '-', label='sum', color='m')
plt.xlim(-10, 175)
plt.ylim(30, 15)
plt.legend();
In [58]:
# TODO: для калибровки нужны светимости диска видимо?
M_R = -20.74
M_B = -19.55 # а у Гутиереза -20.5, что достаточно критично
M_to_L_R = bell_mass_to_light(M_B - M_R, 'R', 'B-R')
M_to_L_R
Out[58]:
In [59]:
# TODO: почему балдж не учитывается?
surf_R = [surf_density(mu=m, M_to_L=M_to_L_R, band='R') for m in disc]
In [60]:
M_to_L_B = bell_mass_to_light(M_B - M_R, 'B', 'B-R')
surf_B = [surf_density(mu=m, M_to_L=M_to_L_B, band='B') for m in disc_B]
plt.plot(p_, surf_R, '-', label='Noord R [M/L={:2.2f}]'.format(M_to_L_R))
plt.plot(p_, surf_B, '-', label='Noord B [M/L={:2.2f}]'.format(M_to_L_B))
plt.legend()
Out[60]:
In [61]:
all_photometry.append(('Noorder R', r_eff_R, mu_eff_R, n_R, mu0d_R, h_disc_R, M_to_L_R,
lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_R, h=h_disc_R), M_to_L=M_to_L_R, band='R')))
all_photometry.append(('Noorder B', r_eff_B, mu_eff_B, n_B, mu0d_B, h_disc_B, M_to_L_B,
lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_B, h=h_disc_B), M_to_L=M_to_L_B, band='B')))
$J$ из Mendez-Abreu:
In [62]:
# for J band
r_eff_J = 11.9
mu_eff_J = 18.13
n_J = 3.75
mu0d_J = 19.07
h_disc_J = 29.2
In [63]:
disc_J = [mu_disc(l, mu0=mu0d_J, h=h_disc_J) for l in p_]
M_to_L_J = bell_mass_to_light(M_B - M_R, 'J', 'B-R')
surf_J = [surf_density(mu=m, M_to_L=M_to_L_J, band='J') for m in disc_J]
plt.plot(p_, surf_R, '-', label='Noord R [M/L={:2.2f}]'.format(M_to_L_R))
plt.plot(p_, surf_B, '-', label='Noord B [M/L={:2.2f}]'.format(M_to_L_B))
plt.plot(p_, surf_J, '-', label='Mendez-Abreu J [M/L={:2.2f}]'.format(M_to_L_J))
plt.legend()
Out[63]:
Тоже весьма похоже.
In [64]:
all_photometry.append(('Mendez-Abreu 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')))
Последнее - Gutierrez в R:
In [65]:
Image('gutierrez_cite_p24.png')
Out[65]:
Тут два диска, видимо надо брать оба:
In [66]:
# данные из базы в Vizier
_, _, r_phot2, mu_phot2 = zip(*np.loadtxt("gutierrez_photomR.dat", str, delimiter='\t'))
r_phot2, mu_phot2 = map(float, r_phot2), map(float, mu_phot2)
plt.plot(r_phot2, mu_phot2, 's')
plt.ylim(28, 15);
In [67]:
h_inner = 30.
h_out = 59.9
mu0_inner = 19.54
mu0_out = 21.53
In [68]:
Image('gutierrez_photomR.png' , width=500)
Out[68]:
In [69]:
fig = plt.figure(figsize=[5, 5])
plt.plot(p_, [mu_disc(l, mu0=mu0_inner, h=h_inner) for l in p_], '--', label='in', color='#007700')
plt.plot(p_, [mu_disc(l, mu0=mu0_out, h=h_out) for l in p_], '--', label='out', color='#000077')
plt.plot(p_, total_mu_profile([mu_disc(l, mu0=mu0_inner, h=h_inner) for l in p_],
[mu_disc(l, mu0=mu0_out, h=h_out) for l in p_]), '-', label='sum', color='#770000')
plt.plot(r_phot2, mu_phot2, '.', ms=2)
plt.xlim(0., 400)
plt.ylim(28, 15)
plt.legend()
plt.axhline(y=23.3, ls='--', alpha=0.2)
plt.axvline(x = 111., linestyle='-.', alpha=0.2);
In [70]:
disc_1 = [mu_disc(l, mu0=mu0_inner, h=h_inner) for l in p_]
disc_2 = [mu_disc(l, mu0=mu0_out, h=h_out) for l in p_]
total_profile_2 = total_mu_profile(disc_1, disc_2)
fig = plt.figure(figsize=[8, 5])
plt.plot(r_phot2, mu_phot2, '.')
plt.plot(p_, disc_1, '-', label='disk inner', color='green')
plt.plot(p_, disc_2, '-', label='disk out', color='green')
plt.plot(p_, total_profile_2, '-', label='sum', color='m')
plt.ylim(28, 15)
plt.legend();
Сумма выглядит достаточно странно, ибо она больше очевидно чем сам профиль (и в самой работе такого много на самом деле).
In [71]:
surf_R2 = [surf_density(mu=m, M_to_L=M_to_L_R, band='R') for m in total_profile_2]
plt.plot(p_, surf_R2, '-', label='R2 [M/L={:2.2f}]'.format(M_to_L_R))
surf_R2_ = [surf_density(mu=mu_disc(l, mu0=mu0_inner, h=h_inner), M_to_L=M_to_L_R, band='R') for l in p_]
plt.plot(p_, surf_R2_, '-', label='inner only'.format(M_to_L_R)) #только внутренний диск
gutierrez_surf = lambda l: surf_density(total_mu_profile([mu_disc(l, mu0_inner, h_inner)], [mu_disc(l, mu0_out, h_out)]), M_to_L_R, 'R')
plt.legend();
Ожидаемо в принципе, слишком массивно получилось.
Добавлять не будем, явно переоцененная модель.
In [72]:
# all_photometry.append(('Gutierrez R (in)', None, None, None, mu0_inner, h_inner, M_to_L_R,
# lambda l: surf_density(mu=mu_disc(l, mu0=mu0_inner, h=h_inner), M_to_L=M_to_L_R, band='R')))
# all_photometry.append(('Gutierrez R (out)', None, None, None, mu0_out, h_out, M_to_L_R,
# lambda l: surf_density(mu=mu_disc(l, mu0=mu0_out, h=h_out), M_to_L=M_to_L_R, band='R')))
Попробую решить проблему с выступающим профилем:
In [73]:
residual_profile = [-2.5*np.log10(np.power(10, -np.array(l[1])/2.5) - np.power(10, -np.array(mu_disc(l[0], mu0=mu0_out, h=h_out))/2.5)) for l in zip(r_phot2, mu_phot2)]
fig = plt.figure(figsize=[8, 5])
plt.plot(r_phot2, residual_profile, '.')
plt.plot(r_phot2[:161], residual_profile[:161], '.')
plt.ylim(32, 15)
plt.legend();
In [74]:
fig = plt.figure(figsize=[8, 5])
plt.plot(r_phot2, residual_profile, '.')
plt.plot(r_phot2[100:161], residual_profile[100:161], '.')
inner_approx = poly1d(polyfit(r_phot2[100:161], residual_profile[100:161], deg=1))
plt.plot(p_, inner_approx(p_), '--')
plt.ylim(32, 15)
plt.legend();
In [75]:
inner_approx
Out[75]:
Т.к. для диска mu = mu0 + 1.0857*r/h, то
In [76]:
h_approx = 1.0857/inner_approx[1]
h_approx
Out[76]:
Т.е. у меня получилась короче шкала и массивнее внутренний диск.
Проверяем:
In [77]:
fig = plt.figure(figsize=[8, 5])
plt.plot(p_, [mu_disc(l, mu0=inner_approx[0], h=h_approx) for l in p_], '--', label='in', color='#007700')
plt.plot(p_, [mu_disc(l, mu0=mu0_out, h=h_out) for l in p_], '--', label='out', color='#000077')
plt.plot(p_, total_mu_profile([mu_disc(l, mu0=inner_approx[0], h=h_approx) for l in p_],
[mu_disc(l, mu0=mu0_out, h=h_out) for l in p_]), '-', label='sum', color='#770000')
plt.plot(r_phot2, mu_phot2, '.', ms=2)
plt.xlim(0., 400)
plt.ylim(28, 15)
plt.legend();
In [78]:
disc_1_ = [mu_disc(l, mu0=inner_approx[0], h=h_approx) for l in p_]
disc_2 = [mu_disc(l, mu0=mu0_out, h=h_out) for l in p_]
total_profile_2_ = total_mu_profile(disc_1_, disc_2)
surf_R2_ = [surf_density(mu=m, M_to_L=M_to_L_R, band='R') for m in total_profile_2_]
plt.plot(p_, surf_R2_, '-', label='R2 [M/L={:2.2f}]'.format(M_to_L_R))
gutierrez_surf_ = lambda l: surf_density(total_mu_profile([mu_disc(l, inner_approx[0], h_approx)], [mu_disc(l, mu0_out, h_out)]), M_to_L_R, 'R')
plt.legend();
Массивнее стало за счет того, что было $19.5^m$, а стало $19^m$, т.е. примерно в полтора раза увеличилось.
In [79]:
all_photometry.append(('Gutierrez R (new)', None, None, None, (inner_approx[0],mu0_out), (h_approx,h_out), M_to_L_R,
(lambda l: surf_density(mu=mu_disc(l, mu0=inner_approx[0], h=h_approx), M_to_L=M_to_L_R, band='R'),
lambda l: surf_density(mu=mu_disc(l, mu0=mu0_out, h=h_out), M_to_L=M_to_L_R, band='R'))))
Я нашел еще две работы в $V$ - Driel 1994 и Pignatelli 2013, можно использовать другой цвет в калибровке. В первой нет нормальной декомпозиции, возбмем более новую (цвет $B-V$ совпадает примерно 0.85 в первой и 0.90 во второй). Еще ценно, что там есть ссылки на кучу другой фотометрии:
TODO: почему тут написано Pignatelli 2013, хотя картинки из работы 2001 года?
In [80]:
Image('pignatelli_cite_p3.png', width=500)
Out[80]:
In [81]:
# Pignatelli 2013
# у них там две модели, непонятно, какую лучше брать
mu_eff_V = 20.9
r_eff_V = 25.6
n_V = 4 #они использовали Вокулера
mu0d_V = 21.6
h_disc_V = 45.9
In [82]:
Image('pignatelli_photomV.png')
Out[82]:
In [83]:
# данные из базы в Vizier
r_photV, mu_photV = zip(*np.loadtxt("pignatelli_photomV.dat", float, delimiter=','))
plt.plot(r_photV, mu_photV, 's')
plt.ylim(26, 16.5);
In [84]:
bulge_V = [mu_bulge(l, mu_eff=mu_eff_V, r_eff=r_eff_V, n=n_V) for l in p_]
disc_V = [mu_disc(l, mu0=mu0d_V, h=h_disc_V) for l in p_]
total_profile_V = total_mu_profile(bulge_V, disc_V)
fig = plt.figure(figsize=[8, 5])
plt.plot(r_photV, mu_photV, 'o')
plt.plot(p_, bulge_V, '-', label='bulge', color='red')
plt.plot(p_, disc_V, '-', label='disk', color='green')
plt.plot(p_, total_profile_V, '-', label='sum', color='m')
plt.xlim(-10, 175)
plt.ylim(30, 15)
plt.legend();
Вторая модель (2D):
TODO: какую модель брать?
In [85]:
# Pignatelli 2013
# у них там две модели, непонятно, какую лучше брать
mu_eff_V = 20.6
r_eff_V = 18.9
n_V = 4 #они использовали Вокулера
mu0d_V = 20.4
h_disc_V = 29.0
In [86]:
bulge_V = [mu_bulge(l, mu_eff=mu_eff_V, r_eff=r_eff_V, n=n_V) for l in p_]
disc_V = [mu_disc(l, mu0=mu0d_V, h=h_disc_V) for l in p_]
total_profile_V = total_mu_profile(bulge_V, disc_V)
fig = plt.figure(figsize=[8, 5])
plt.plot(r_photV, mu_photV, 'o')
plt.plot(p_, bulge_V, '-', label='bulge', color='red')
plt.plot(p_, disc_V, '-', label='disk', color='green')
plt.plot(p_, total_profile_V, '-', label='sum', color='m')
plt.xlim(-10, 175)
plt.ylim(30, 15)
plt.legend();
У них также есть $M/L$, для диска это 4.2
In [87]:
M_to_L_V = bell_mass_to_light(0.90, 'V', 'B-V')
print M_to_L_V
Достаточно близко.
In [88]:
surf_V = [surf_density(mu=m, M_to_L=M_to_L_V, band='V') for m in disc_V]
plt.plot(p_, surf_R, '-', label='Noord R [M/L={:2.2f}]'.format(M_to_L_R))
plt.plot(p_, surf_B, '-', label='Noord B [M/L={:2.2f}]'.format(M_to_L_B))
plt.plot(p_, surf_J, '-', label='Mendez-Abreu J [M/L={:2.2f}]'.format(M_to_L_J))
plt.plot(p_, surf_R2, '-', label='Gutierrez R [M/L={:2.2f}]'.format(M_to_L_R))
plt.plot(p_, surf_R2_, '-', label='Gutierrez R (new) [M/L={:2.2f}]'.format(M_to_L_R))
plt.plot(p_, surf_V, '-', label='Pignatelli V [M/L={:2.2f}]'.format(M_to_L_V))
plt.legend()
Out[88]:
Хм, занимательно - похоже имеем две конкурирующие модели фотометрии. Причем максимальный диск по ноордермееру даст что-то среднее.
In [89]:
all_photometry.append(('Pignatelli 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')))
Для $V$ можно вычислить калибровки McGaugh:
In [90]:
mcgaugh_mass_to_light(0.90, 'V')
Out[90]:
Не сильно отличается от используемой, среднее почти такое же.
In [91]:
show_all_photometry_table(all_photometry, scale)
Можно провести тест-сравнение с кривой вращения тонкого диска при заданной фотометрии, если она слишком массивная - то не брать ее (это ограничение сверху).
In [92]:
fig = plt.figure(figsize=[10,6])
plt.plot(test_points, spl_gas(test_points), '-', label='spline')
plt.errorbar(r_wsrt, vel_wsrt, yerr=e_vel_wsrt, fmt='.', marker='.', mew=0, label = 'WSRT')
for ind, photom in enumerate(all_photometry):
if type(photom[5]) == tuple:
plt.plot(test_points, map(lambda l: disc_vel(l, photom[7][0](0), photom[5][0], scale,
Sigma0_2=photom[7][1](0), h_2=photom[5][1]), test_points), next(linecycler), label=photom[0])
else:
plt.plot(test_points, map(lambda l: disc_vel(l, photom[7](0), photom[5], scale), test_points), next(linecycler), label=photom[0])
plt.ylim(0, 300)
plt.xlim(0, 300)
plt.legend(loc='best');
Тут Гутиеррез имеет массивнее хвост, чем Пигнателли, да и в целом стоит его проверить. По максимумам все похожи.
TODO: понять и добавить
TODO: $H_{\alpha}$
Изображение в $H_{\alpha}$ из Hameed 2005 http://iopscience.iop.org/article/10.1086/430211/pdf: (масштабная полоска - 1 кпк)
In [93]:
Image('hameed2005_halpha.png')
Out[93]:
Есть каринка в $\rm{H}_{\alpha}$ из Pignatelli http://mnras.oxfordjournals.org/content/323/1/188.full.pdf:
In [94]:
Image('pignatelli_halpha.png', width=700)
Out[94]:
Кривая вращения и карта из GHASP 2008 http://mnras.oxfordjournals.org/content/388/2/500.full.pdf:
In [95]:
Image('ghasp_2008_halpha_vel.png')
Out[95]:
In [96]:
Image('ghasp_2008_map.png')
Out[96]:
In [97]:
def plot_SF(ax):
ax.plot([70., 80.], [0., 0.], '-', lw=5., color='red')
ax.plot([190., 210.], [0., 0.], '-', lw=5., color='red')
ax.plot([0., 70.], [0., 0.], '-', lw=3., color='blue')
plot_SF(plt.gca())
plt.xlim(0, 350)
plt.ylim(0, 200)
Out[97]:
In [98]:
Image('diplom_results.png')
Out[98]:
Устойчиво, когда > 1: $$Q_g = \frac{\Sigma_g^{cr}}{\Sigma_g}=\frac{\kappa c_g}{\pi G \Sigma_g}$$ $$Q_s = \frac{\Sigma_s^{cr}}{\Sigma_s}=\frac{\sigma_R}{\sigma_R^{min}}=\frac{\kappa \sigma_R}{3.36 G \Sigma_s}$$
In [99]:
sound_vel = 6. #скорость звука в газе, км/с
data_lim = max(r_sig_ma)+10. #где заканчиваются данные
In [100]:
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=sound_vel, gas_density=l[1]*He_coeff) for l in
zip([epicyclicFreq_real(gas_approx, x, scale) for x in r_g_dens], gas_dens)], 's-', label='$Q_{gas*He_coeff}$')
plt.plot(r_g_dens, [1./Qg(epicycl=l[0], sound_vel=4., gas_density=l[1]*He_coeff) for l in
zip([epicyclicFreq_real(gas_approx, x, scale) for x in r_g_dens], gas_dens)], 's-', label='$Q_{gas*He_coeff}$ & c=4')
plt.plot(map(lambda l: l/scale, r_q_n), q_n, 'o', label='noordermeer 1F')
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(mu_disc(l_, mu0=mu0d_R, h=h_disc_R), M_to_L_R, 'R') for l_ 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(mu_disc(l_, mu0=mu0d_R, h=h_disc_R), M_to_L_R, 'R') for l_ in r_g_dens])], 's-', label='$Q_{star}^{max}$')
plt.axhline(y=1, ls='--')
plt.legend()
plt.ylabel('$Q^{-1}$', fontsize=15);
Одножидкостный уже не похож на диплом, но зато похоже на Ноордермеера (хоть и не совпадает с точностью).
НЕ ИСПРАВЛЕНО ЗА 1.6! Что верно и не надо.
Кинетическое приближение: $$\frac{1}{Q_{\mathrm{eff}}}=\frac{2}{Q_{\mathrm{s}}}\frac{1}{\bar{k}}\left[1-e^{-\bar{k}^{2}}I_{0}(\bar{k}^{2})\right]+\frac{2}{Q_{\mathrm{g}}}s\frac{\bar{k}}{1+\bar{k}^{2}s^{2}}>1\,$$
Гидродинамическое приближение: $$\frac{2\,\pi\, G\, k\,\Sigma_{\mathrm{s}}}{\kappa+k^{2}\sigma_{\mathrm{s}}}+\frac{2\,\pi\, G\, k\,\Sigma_{\mathrm{g}}}{\kappa+k^{2}c_{\mathrm{g}}}>1$$ или $$\frac{1}{Q_{\mathrm{eff}}}=\frac{2}{Q_{\mathrm{s}}}\frac{\bar{k}}{1+\bar{k}^{2}}+\frac{2}{Q_{\mathrm{g}}}s\frac{\bar{k}}{1+\bar{k}^{2}s^{2}}>1$$ для безразмерного волнового числа ${\displaystyle \bar{k}\equiv\frac{k\,\sigma_{\mathrm{s}}}{\kappa}},\, s=c/\sigma$
In [101]:
total_gas_data = zip(r_g_dens, map(lambda l: He_coeff*l, gas_dens))[1:7]
disk_scales = []
for l in all_photometry:
try:
disk_scales.append((l[5][0], l[0].split(' ')[1])) #внутренний диск только
except TypeError:
disk_scales.append((l[5], l[0].split(' ')[1]))
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_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='g', alpha=0.3, 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();
In [102]:
total_gas_data = zip(r_g_dens, map(lambda l: He_coeff*l, gas_dens))[1:7]
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=gutierrez_surf,
star_density_min=gutierrez_surf,
data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='Gutierrez R 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=gutierrez_surf_,
star_density_min=gutierrez_surf_,
data_lim=data_lim, color='y', alpha=0.3, disk_scales=disk_scales, label='Gutierrez R (new) maj max/min')
plt.ylim(0., 2.5)
plt.axhline(y=1., ls='-', color='grey')
plot_SF(ax)
plt.grid();
In [103]:
from matplotlib.animation import FuncAnimation
fig = plt.gcf()
plt.figure(figsize=(10,6))
ax = plt.gca()
def animate(i):
ax.cla()
if all_photometry[i][-1] is not None:
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=tot_dens(all_photometry[i][-1]),
star_density_min=tot_dens(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)
In [104]:
anim.save('..\\..\pics\\'+name+'.gif', writer='imagemagick', fps=1)
In [105]:
from IPython.display import HTML
HTML(anim.to_html5_video())
Out[105]:
Существует ограничение на максимальный диск в ~0.85 (изотермическое гало) и на субмаксимальный в 0.55-0.6 (NFW гало). Попробуем дотянуть фотметрию до максимальных дисков и посмотрим, как изменятся M/L (скорость зависит как корень из M/L):
In [106]:
fig = plt.figure(figsize=[12,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 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)
values = map(disc_v, test_points)
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)
if type(photom[5]) == tuple:
plt.plot(test_points, map(lambda l: disc_vel(l, max_coeff**2 *photom[7][0](0), photom[5][0], scale,
Sigma0_2=max_coeff**2 *photom[7][1](0), h_2=photom[5][1]), test_points), next(linecycler), label=photom[0] + '_MAX')
else:
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, 300)
plt.xlim(0, 300)
plt.legend(loc=(1.01, 0.));
В принципе все смотрится достаточно неплохо, хотя конечно скорее максимум правее. $B$ не проходит по M/L.
In [107]:
from matplotlib.animation import FuncAnimation
fig = plt.gcf()
plt.figure(figsize=(10,6))
ax = plt.gca()
def animate(i):
ax.cla()
if all_photometry[i][7] is not None:
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]*tot_dens(all_photometry[i][-1])(l),
star_density_min=lambda l: max_coeffs[all_photometry[i][0]][0]*tot_dens(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]*tot_dens(all_photometry[i][-1])(l),
star_density_min=lambda l: max_coeffs[all_photometry[i][0]][1]*tot_dens(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)
return ax
anim = FuncAnimation(plt.gcf(), animate, frames=len(all_photometry), interval=1000)
In [108]:
anim.save('..\\..\pics\\'+name+'_MAXDISCS.gif', writer='imagemagick', fps=1)
In [109]:
from IPython.display import HTML
HTML(anim.to_html5_video())
Out[109]:
In [110]:
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, 30.), '--', alpha=0.5)
plt.plot(r_g_dens, map(lambda l: y_interp_(l[0], 30.) + l[1], zip(r_g_dens, gas_dens)), '--', alpha=0.5);
In [111]:
fig = plt.figure(figsize=[10, 6])
total_gas_data_ = zip(r_g_dens, [He_coeff*(y_interp_(l[0], 29.) + l[1]) for l in zip(r_g_dens, gas_dens)])[:7]
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_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='g', alpha=0.3, 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();
Остальные фотометрии:
In [112]:
from matplotlib.animation import FuncAnimation
fig = plt.gcf()
plt.figure(figsize=(10,6))
ax = plt.gca()
def animate(i):
ax.cla()
if type(all_photometry[i][-1]) == tuple:
total_gas_data_ = zip(r_g_dens, [He_coeff*(y_interp_(l[0], all_photometry[i][-3][0]) + l[1]) for l in zip(r_g_dens, gas_dens)])[:7]
else:
total_gas_data_ = zip(r_g_dens, [He_coeff*(y_interp_(l[0], all_photometry[i][-3]) + l[1]) for l in zip(r_g_dens, gas_dens)])[:7]
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=tot_dens(all_photometry[i][-1]),
star_density_min=tot_dens(all_photometry[i][-1]),
data_lim=data_lim, color=np.random.rand(3), alpha=0.3, disk_scales=disk_scales)
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)
In [113]:
anim.save('..\\..\pics\\'+name+'_molec.gif', writer='imagemagick', fps=1)
In [114]:
from IPython.display import HTML
HTML(anim.to_html5_video())
Out[114]:
Приподнимается, но совсем чуть-чуть.
И оценки с максимальным диском:
In [115]:
from matplotlib.animation import FuncAnimation
fig = plt.gcf()
plt.figure(figsize=(10,6))
ax = plt.gca()
def animate(i):
ax.cla()
if type(all_photometry[i][5]) == tuple:
total_gas_data_ = zip(r_g_dens, [He_coeff*(y_interp_(l[0], all_photometry[i][-3][0]) + l[1]) for l in zip(r_g_dens, gas_dens)])[:7]
else:
total_gas_data_ = zip(r_g_dens, [He_coeff*(y_interp_(l[0], all_photometry[i][-3]) + l[1]) for l in zip(r_g_dens, gas_dens)])[:7]
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: max_coeffs[all_photometry[i][0]][0]*tot_dens(all_photometry[i][-1])(l),
star_density_min=lambda l: max_coeffs[all_photometry[i][0]][0]*tot_dens(all_photometry[i][-1])(l),
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., 1.5)
plot_Q_levels(ax, [1., 1.5, 2., 3.])
return ax
anim = FuncAnimation(plt.gcf(), animate, frames=len(all_photometry), interval=1000)
In [116]:
anim.save('..\\..\pics\\'+name+'_molec_MAX.gif', writer='imagemagick', fps=1)
In [117]:
from IPython.display import HTML
HTML(anim.to_html5_video())
Out[117]:
In [ ]:
summary_imgs_path = '..\\..\pics\\notebook_summary\\'
def save_model_plot(path):
fig, axes = plt.subplots(1, 5, figsize=[40,7])
fig.tight_layout()
axes[0].imshow(ImagePIL.open('ngc3898_SDSS_labeled.jpeg'), aspect='auto')
axes[0].set_title(name)
try:
axes[1].errorbar(r_sig_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='.', mew=0, color='red', label='$\sigma_{los}^{maj}$')
axes[1].plot(points, map(sig_R_maj_min, points))
axes[1].plot(points, map(sig_R_maj_max, points))
axes[1].plot(points, map(sig_R_maj_maxmaxtrue, points))
except Exception:
pass
try:
axes[1].errorbar(r_sig_mi, sig_mi, yerr=e_sig_mi, fmt='.', marker='.', mew=0, color='b', label='$\sigma_{los}^{min}$')
axes[1].plot(points, map(sig_R_minor_min, points), '--')
axes[1].plot(points, map(sig_R_minor_max, points), '--')
except Exception:
pass
axes[1].set_ylim(0,250)
axes[1].set_xlim(0, 105)
axes[1].grid()
axes[1].legend()
axes[1].set_title('Dispersions')
for photom in all_photometry:
axes[2].plot(r_g_dens, map(tot_dens(photom[-1]), r_g_dens), '-', label='{} [M/L={:2.2f}]'.format(photom[0], photom[-2]))
axes[2].set_xlim(0, 150)
axes[2].set_ylim(0, 300)
axes[2].set_title('Photometry')
axes[2].grid()
axes[2].legend()
axes[3].plot(r_g_dens, gas_dens, 'd-')
axes[3].plot(r_g_dens, [He_coeff*(y_interp_(l[0], h_disc_R) + l[1]) for l in zip(r_g_dens, gas_dens)], '*-')
axes[3].plot(r_g_dens, [y_interp_(l, h_disc_R) for l in r_g_dens], '--', label='H2 (R)')
axes[3].set_title('Gas')
axes[3].grid()
axes[3].set_xlim(0, 200)
axes[3].legend()
#change this
plot_2f_vs_1f(ax=axes[4], 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)])[1:7],
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_R, h=h_disc_R), 7.80, 'R'),
star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_R, h=h_disc_R), 7.80, 'R'),
data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='R Noord maxdisc')
plot_2f_vs_1f(ax=axes[4], total_gas_data=zip(r_g_dens, [He_coeff*(y_interp_(l[0], h_approx) + l[1]) for l in zip(r_g_dens, gas_dens)])[1:7],
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: tot_dens((lambda l1: surf_density(mu=mu_disc(l1, mu0=inner_approx[0], h=h_approx), M_to_L=2.93, band='R'),
lambda l2: surf_density(mu=mu_disc(l2, mu0=mu0_out, h=h_out), M_to_L=2.93, band='R')))(l),
star_density_min=lambda l: tot_dens((lambda l1: surf_density(mu=mu_disc(l1, mu0=inner_approx[0], h=h_approx), M_to_L=2.93, band='R'),
lambda l2: surf_density(mu=mu_disc(l2, mu0=mu0_out, h=h_out), M_to_L=2.93, band='R')))(l),
data_lim=data_lim, color='m', alpha=0.2, disk_scales=disk_scales, label='R Gutierrez 2d maxdisc')
axes[4].set_ylim(0., 2.5)
axes[4].set_xlim(0., 130.)
axes[4].axhline(y=1., ls='-', color='grey')
plot_SF(axes[4])
axes[4].grid()
axes[4].set_title('Instability')
# plt.savefig(path+name+'.png', format='png', bbox_inches='tight');
save_model_plot(summary_imgs_path)