TODO: сделать так, чтобы можно было импортировать
In [1]:
%run ../../utils/load_notebook.py
from instabilities import *
import numpy as np
Коэффициент для учета вклада гелия в массу газа (см. Notes):
In [1]:
He_coeff = 1.36
Коэффициент, с которым пересчитывается масса молекулярного газа:
In [1]:
X_CO = 1.9
Путь для картинок в статью:
In [1]:
paper_imgs_dir = r'C:\Users\root\Dropbox\RotationCurves\PhD\paper2\imgs\\'
Путь до сохраненных моделей:
In [1]:
models_path = '..\\..\\notebooks\\2f\\test_short\\models\\'
In [ ]:
Врапер, который позволяет сохранить модель и не пересчитывать ее бесконечно. Используется для plot_2f_vs_1f.
In [1]:
def save_model(path):
'''враппер для plot_2f_vs_1f, сохраняющий параметры модели по заданному пути до .npy файла'''
def real_decorator(function):
def wrapper(*args, **kwargs):
dictionary = kwargs.copy()
rr = zip(*kwargs['total_gas_data'])[0]
dictionary.pop('ax', None)
for i in dictionary.items():
if callable(i[1]):
if i[0] == 'epicycl':
dictionary[i[0]] = [i[1](kwargs['gas_approx'], r, kwargs['scale']) for r in rr]
else:
try:
dictionary[i[0]] = [i[1](r) for r in rr]
except ValueError:
dictionary[i[0]] = i[1](rr)
np.save(path, dictionary)
return function(*args, **kwargs)
return wrapper
return real_decorator
In [ ]:
In [1]:
def calc_scale(D):
'''Масштаб в кпк/arcsec из расстояния в Мпк'''
return (D*1000*2*3.141596)/(360*60*60.)
In [ ]:
In [2]:
def flat_end(argument):
'''декоратор для того, чтобы продолжать функцию на уровне последнего значения'''
def real_decorator(function):
def wrapper(*args, **kwargs):
if args[0] < argument:
return function(*args, **kwargs)
else:
return function(argument, *args[1:], **kwargs)
return wrapper
return real_decorator
Для большой оси: $\sigma^2_{maj} = \sigma^2_{\varphi}\sin^2 i + \sigma^2_{z}\cos^2 i$, следовательно примерные ограничения $$\sigma_{maj} < \frac{\sigma_{maj}}{\sqrt{\sin^2 i + 0.49\cos^2 i}}< \sigma_R = \frac{\sigma_{maj}}{\sqrt{f\sin^2 i + \alpha^2\cos^2 i}} ~< \frac{\sigma_{maj}}{\sqrt{0.5\sin^2 i + 0.09\cos^2 i}} < \frac{\sqrt{2}\sigma_{maj}}{\sin i} (или \frac{\sigma_{maj}}{\sqrt{f}\sin i}),$$ или можно более точную оценку дать, если построить $f$ (сейчас $0.5 < f < 1$).
Для малой оси: $\sigma^2_{min} = \sigma^2_{R}\sin^2 i + \sigma^2_{z}\cos^2 i$ и ограничения $$\sigma_{min} < \frac{\sigma_{min}}{\sqrt{\sin^2 i + 0.49\cos^2 i}} < \sigma_R = \frac{\sigma_{min}}{\sqrt{\sin^2 i + \alpha^2\cos^2 i}} ~< \frac{\sigma_{min}}{\sqrt{\sin^2 i + 0.09\cos^2 i}} < \frac{\sigma_{min}}{\sin i}$$
Соответственно имеем 5 оценок из maj и 4 оценки из min.
In [3]:
# sig_maj_lim=None
# spl_maj=None
# @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 [4]:
# sig_min_lim=None
# spl_min=None
# @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 [1]:
# TODO: move to proper place
def plot_data_lim(ax, data_lim):
'''Вертикальная линия, обозначающая конец данных'''
ax.axvline(x=data_lim, ls='-.', color='black', alpha=0.5)
def plot_disc_scale(scale, ax, text=None):
'''Обозначает масштаб диска'''
ax.plot([scale, scale], [0., 0.05], '-', lw=6., color='black')
if text:
ax.annotate(text, xy=(scale, 0.025), xytext=(scale, 0.065), textcoords='data', arrowprops=dict(arrowstyle="->"))
def plot_Q_levels(ax, Qs, style='--', color='grey', alpha=0.4):
'''Функция, чтобы рисовать горизонтальные линии различных уровней $Q^{-1}$:'''
for Q in Qs:
ax.axhline(y=1./Q, ls=style, color=color, alpha=alpha)
def plot_2f_vs_1f(ax=None, total_gas_data=None, epicycl=None, gas_approx=None, sound_vel=None, scale=None, sigma_max=None, sigma_min=None, star_density_max=None,
star_density_min=None, data_lim=None, color=None, alpha=0.3, disk_scales=[], label=None, verbose=False, **kwargs):
'''Картинка сравнения 2F и 1F критерия для разных фотометрий и величин sig_R,
куда подается весь газ, результат НЕ исправляется за осесимметричные возмущения.'''
invQg, invQs, invQeff_min = zip(*get_invQeff_from_data(gas_data=total_gas_data,
epicycl=epicycl,
gas_approx=gas_approx,
sound_vel=sound_vel,
scale=scale,
sigma=sigma_max,
star_density=star_density_min, verbose=verbose))
invQg, invQs, invQeff_max = zip(*get_invQeff_from_data(gas_data=total_gas_data,
epicycl=epicycl,
gas_approx=gas_approx,
sound_vel=sound_vel,
scale=scale,
sigma=sigma_min,
star_density=star_density_max, verbose=verbose))
rr = zip(*total_gas_data)[0]
ax.fill_between(rr, invQeff_min, invQeff_max, color=color, alpha=alpha, label=label)
ax.plot(rr, invQeff_min, 'd-', color=color, alpha=0.6)
ax.plot(rr, invQeff_max, 'd-', color=color, alpha=0.6)
ax.plot(rr, invQg, 'v-', color='b')
ax.set_ylim(0., 1.5)
ax.set_xlim(0., data_lim+50.)
plot_data_lim(ax, data_lim)
for h, annot in disk_scales:
plot_disc_scale(h, ax, annot)
plot_Q_levels(ax, [1., 1.5, 2., 3.])
ax.legend()
Для случая бесконечного тонкого диска: $$\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 [1]:
def epicyclicFreq_real(poly_gas, R, resolution):
'''Честное вычисление эпициклической частоты на расстоянии R для сплайна или полинома'''
try:
return sqrt(2.0) * poly_gas(R) * sqrt(1 + R * poly_gas.deriv()(R) / poly_gas(R)) / (R * resolution )
except:
return sqrt(2.0) * poly_gas(R) * sqrt(1 + R * poly_gas.derivative()(R) / poly_gas(R)) / (R * resolution )
In [ ]:
Два других механизма из http://iopscience.iop.org/article/10.1088/0004-6256/148/4/69/pdf:
Schaye (2004), 'cold gas phase': $$\Sigma_g > 6.1 f_g^{0.3} Z^{-0.3} I^{0.23}$$ или при constant metallicity of 0.1 $Z_{sun}$ and interstellar flux of ionizing photons 10^6 cm−2 s−1: $$\Sigma_g > 6.1 \frac{\Sigma_g}{\Sigma_g + \Sigma_s}$$
In [59]:
def Sigma_crit_S04(gas_dens, r_gas, star_surf_dens):
return 6.1 * gas_dens / (gas_dens + star_surf_dens(r_gas))
Hunter et al (1998), 'competition with shear' according to Leroy: $$\Sigma_A = \alpha_A\frac{\sigma_g A}{\pi G}$$
In [2]:
def oort_a(r, gas_vel):
try:
return 0.5 * (gas_vel(r)/r - gas_vel.deriv()(r))
except:
return 0.5 * (gas_vel(r)/r - gas_vel.derivative()(r))
In [3]:
def Sigma_crit_A(r, gas_vel, alpha, sound_vel):
G = 4.32
return alpha * (sound_vel*oort_a(r, gas_vel)) / (np.pi*G)
In [ ]:
Кривая вращения тонкого диска: $$\frac{v^2}{r} = 2\pi G \Sigma_0 \frac{r}{2h} \left[I_0(\frac{r}{2h})K_0(\frac{r}{2h}) - I_1(\frac{r}{2h})K_1(\frac{r}{2h})\right]$$
In [4]:
from scipy.special import i0, i1, k0, k1
def disc_vel(r, Sigma0, h, scale, Sigma0_2=None, h_2=None):
G = 4.3
bessels = i0(0.5*r/h)*k0(0.5*r/h) - i1(0.5*r/h)*k1(0.5*r/h)
if h_2 is None:
return np.sqrt(2*np.pi*G*Sigma0*r*scale * 0.5*r/h * bessels)
else: #двухдисковая модель
bessels2 = i0(0.5*r/h_2)*k0(0.5*r/h_2) - i1(0.5*r/h_2)*k1(0.5*r/h_2)
return np.sqrt(2*np.pi*G*Sigma0*r*scale * 0.5*r/h * bessels + 2*np.pi*G*Sigma0_2*r*scale * 0.5*r/h_2 * bessels2)
In [ ]:
Функция для печати статистики по фотометриям, куда добавлена также информация о полной массе диска $M_d = 2\pi h^2 \Sigma(0)$ (только нужно учесть, что там пк в arcsec надо перевести):
In [3]:
from tabulate import tabulate
import pandas as pd
def show_all_photometry_table(all_photometry, scale):
'''scale in kpc/arcsec'''
copy = [list(l) for l in all_photometry]
#все это дальше из-за того, что бывают двухдисковые модели и их нужно по другому использовать
for entry in copy:
if type(entry[5]) == tuple:
entry[5] = (round(entry[5][0], 2), round(entry[5][1], 2))
else:
entry[5] = round(entry[5], 2)
for entry in copy:
if type(entry[4]) == tuple:
entry[4] = (round(entry[4][0], 2), round(entry[4][1], 2))
else:
entry[4] = round(entry[4], 2)
for entry in copy:
if type(entry[5]) == tuple:
entry.append(2*math.pi*entry[5][0]**2 * entry[-1][0](0) * (scale * 1000.)**2 +
2*math.pi*entry[5][1]**2 * entry[-1][1](0) * (scale * 1000.)**2)
else:
entry.append(2*math.pi*entry[5]**2 * entry[-1](0) * (scale * 1000.)**2)
for entry in copy:
if type(entry[5]) == tuple:
entry.append(entry[7][0](0) + entry[7][1](0))
else:
entry.append(entry[7](0))
df = pd.DataFrame(data=copy, columns=['Name', 'r_eff', 'mu_eff', 'n', 'mu0_d', 'h_disc', 'M/L', 'surf', 'M_d/M_sun', 'Sigma_0'])
df['M/L'] = df['M/L'].apply(lambda l: '%2.2f'%l)
# df['Sigma_0'] = df['surf'].map(lambda l:l(0))
df['Sigma_0'] = df['Sigma_0'].apply(lambda l: '%2.0f' % l)
# df['M_d/M_sun'] = 2*math.pi*df['h_disc']**2 * df['surf'].map(lambda l:l(0)) * (scale * 1000.)**2
df['M_d/M_sun'] = df['M_d/M_sun'].apply(lambda l: '%.2E.' % l)
df.drop('surf', axis=1, inplace=True)
print tabulate(df, headers='keys', tablefmt='psql', floatfmt=".2f")
In [ ]:
Префикс для плохих фотометрий, чтобы их не брать в конце.
In [1]:
BAD_MODEL_PREFIX = 'b:'
In [ ]:
Функция, которая возвращает суммарный профиль плотности для двухдисковой модели, если это необходимо:
In [1]:
def tot_dens(dens):
if type(dens) == tuple:
star_density = lambda l: dens[0](l) + dens[1](l)
else:
star_density = lambda l: dens(l)
return star_density
In [ ]:
Цикл по стилям линий для того, что бы они отличались там, где их много на одной картинке и цвета сливаются.
In [1]:
from itertools import cycle
lines = ["-","--","-.",":"]
linecycler = cycle(lines)
In [ ]:
Раскрашиваем задний фон в "зебру" с некоей периодичностью.
In [1]:
def foreground_zebra(ax, step, alpha):
for i in range(int(ax.get_xlim()[1])+1):
if i%2 == 0:
ax.axvspan(i*step, (i+1)*step, color='grey', alpha=alpha)
In [ ]:
Сравним с оценкой Romeo & Falstad (2013) https://ui.adsabs.harvard.edu/#abs/2013MNRAS.433.1389R/abstract:
$$Q_N^{-1} = \sum_{i=1}^{N}\frac{W_i}{Q_iT_i}$$где
$$Q_i = \frac{\kappa\sigma_{R,i}}{\pi G\Sigma_i}$$$$T_i= \begin{cases} 1 + 0.6(\frac{\sigma_z}{\sigma_R})^2_i, & \mbox{if } 0.0 \le \frac{\sigma_z}{\sigma_R} \le 0.5, \\ 0.8 + 0.7(\frac{\sigma_z}{\sigma_R})_i, & \mbox{if } 0.5 \le \frac{\sigma_z}{\sigma_R} \le 1.0 \end{cases}$$$$W_i = \frac{2\sigma_{R,m}\sigma_{R,i}}{\sigma_{R,m}^2 + \sigma_{R,i}^2},$$$$m:\ index\ of\ min(T_iQ_i)$$В самой развитой модели 3 компонента - HI, CO, stars. В их же модели полагается $(\sigma_z/\sigma_R)_{CO} = (\sigma_z/\sigma_R)_{HI} = 1$, т.е. $T_{CO} = T_{HI} = 1.5$. Скорость звука я буду полагать в обеих средах равной 11 км/c. Звезды у меня в двух крайних случаях бывают с $\alpha$ равным 0.3 и 0.7, соответственно $T_{0.3} = 1.05;\ T_{0.7}=1.29$.
In [ ]:
from math import pi
def romeo_Qinv(r=None, epicycl=None, sound_vel=None, sound_vel_CO=6., sound_vel_HI=11., sigma_R=None, star_density=None,
HI_density=None, CO_density=None, alpha=None, verbose=False, thin=True, He_corr=False):
'''Возвращает приближение Q из Romeo&Falstad(2013), оно же трехкомпонентное Romeo&Wiegert(2011) и наиболее неустойчивую компоненту.'''
G = 4.32
kappa = epicycl
if not He_corr:
CO_density = He_coeff*CO_density
HI_density = He_coeff*HI_density
if sound_vel is not None:
sound_vel_CO=sound_vel
sound_vel_HI=sound_vel
Q_star = kappa*sigma_R/(pi*G*star_density) # не 3.36
Q_CO = kappa*sound_vel_CO/(pi*G*CO_density)
Q_HI = kappa*sound_vel_HI/(pi*G*HI_density)
if not thin:
T_CO, T_HI = 1.5, 1.5
if alpha > 0 and alpha <= 0.5:
T_star = 1. + 0.6*alpha**2
else:
T_star = 0.8 + 0.7*alpha
else:
T_CO, T_HI, T_star = 1., 1., 1.
dispersions = [sigma_R, sound_vel_HI, sound_vel_CO]
QTs = [Q_star*T_star, Q_HI*T_HI, Q_CO*T_CO]
components = ['star', 'HI', 'H2']
mindex = QTs.index(min(QTs))
if verbose:
print 'QTs: {}'.format(QTs)
print 'min index: {}'.format(mindex)
print 'min component: {}'.format(components[mindex])
sig_m = dispersions[mindex]
def W_i(sig_m, sig_i):
return 2*sig_m*sig_i/(sig_m**2 + sig_i**2)
if verbose:
print 'Ws/TQs={:5.3f} WHI/TQHI={:5.3f} WCO/TQCO={:5.3f}'.format(W_i(sig_m, dispersions[0])/QTs[0], W_i(sig_m, dispersions[1])/QTs[1], W_i(sig_m, dispersions[2])/QTs[2])
print 'Ws={:5.3f} WHI={:5.3f} WCO={:5.3f}'.format(W_i(sig_m, dispersions[0]), W_i(sig_m, dispersions[1]), W_i(sig_m, dispersions[2]))
return W_i(sig_m, dispersions[0])/QTs[0] + W_i(sig_m, dispersions[1])/QTs[1] + W_i(sig_m, dispersions[2])/QTs[2], components[mindex]
Двухкомпонентная версия:
In [ ]:
def romeo_Qinv2(r=None, epicycl=None, sound_vel=None, sigma_R=None, star_density=None,
HI_density=None, CO_density=None, alpha=None, verbose=False, thin=True, He_corr=False, **kwargs):
'''Возвращает двухкомпонентное приближение Q из Romeo&Falstad(2011).'''
G = 4.32
kappa = epicycl
if not He_corr:
CO_density_ = He_coeff*CO_density
HI_density_ = He_coeff*HI_density
else:
CO_density_ = CO_density
HI_density_ = HI_density
gas = CO_density_ + HI_density_
Q_star = kappa*sigma_R/(pi*G*star_density) # не 3.36
Q_g = kappa*sound_vel/(pi*G*gas)
if not thin:
T_g = 1.5
if alpha > 0 and alpha <= 0.5:
T_star = 1. + 0.6*alpha**2
else:
T_star = 0.8 + 0.7*alpha
else:
T_g, T_star = 1., 1.
dispersions = [sigma_R, sound_vel]
QTs = [Q_star*T_star, Q_g*T_g]
components = ['star', 'gas']
mindex = QTs.index(min(QTs))
if verbose:
print 'QTs: {}'.format(QTs)
print 'min index: {}'.format(mindex)
print 'min component: {}'.format(components[mindex])
sig_m = dispersions[mindex]
def W_i(sig_m, sig_i):
return 2*sig_m*sig_i/(sig_m**2 + sig_i**2)
if verbose:
print 'Ws/TQs={:5.3f} Wg/TQg={:5.3f}'.format(W_i(sig_m, dispersions[0])/QTs[0], W_i(sig_m, dispersions[1])/QTs[1])
print 'Ws={:5.3f} Wg={:5.3f}'.format(W_i(sig_m, dispersions[0]), W_i(sig_m, dispersions[1]))
return W_i(sig_m, dispersions[0])/QTs[0] + W_i(sig_m, dispersions[1])/QTs[1], components[mindex]
Функция-сравнение с наблюдениями:
In [ ]:
def plot_RF13_vs_2F(r_g_dens=None, HI_gas_dens=None, CO_gas_dens=None, epicycl=None, sound_vel_CO=6., sound_vel_HI=11., sound_vel=None, sigma_R_max=None, sigma_R_min=None,
star_density=None, alpha_max=None, alpha_min=None, thin=True, verbose=False, scale=None, gas_approx=None):
'''Плотности газа передаются не скорр. за гелий.'''
fig = plt.figure(figsize=[20, 5])
ax = plt.subplot(131)
totgas = zip(r_g_dens, [He_coeff*(l[0]+l[1]) for l in zip(HI_gas_dens, CO_gas_dens)])[1:]
if sound_vel is not None:
sound_vel_CO=sound_vel
sound_vel_HI=sound_vel
if verbose:
print 'sig_R_max case:'
romeo_min = []
for r, g, co in zip(r_g_dens, HI_gas_dens, CO_gas_dens):
rom, _ = romeo_Qinv(r=r, epicycl=epicycl(gas_approx, r, scale), sound_vel_CO=sound_vel_CO, sound_vel_HI=sound_vel_HI, sigma_R=sigma_R_max(r),
star_density=star_density(r), HI_density=He_coeff*g, CO_density=He_coeff*co,
alpha=alpha_min, verbose=verbose, He_corr=True, thin=thin)
romeo_min.append(rom)
if _ == 'star':
color = 'g'
elif _ == 'HI':
color = 'b'
else:
color = 'm'
ax.scatter(r, rom, 10, marker='o', color=color)
invQg, invQs, invQeff_min = zip(*get_invQeff_from_data(gas_data=totgas,
epicycl=epicycl,
gas_approx=gas_approx,
sound_vel=sound_vel,
scale=scale,
sigma=sigma_R_max,
star_density=star_density))
if verbose:
print 'sig_R_min case:'
romeo_max = []
for r, g, co in zip(r_g_dens, HI_gas_dens, CO_gas_dens):
rom, _ = romeo_Qinv(r=r, epicycl=epicycl(gas_approx, r, scale), sound_vel_CO=sound_vel_CO, sound_vel_HI=sound_vel_HI, sigma_R=sigma_R_min(r),
star_density=star_density(r), HI_density=He_coeff*g, CO_density=He_coeff*co,
alpha=alpha_max, verbose=verbose, He_corr=True, thin=thin)
romeo_max.append(rom)
if _ == 'star':
color = 'g'
elif _ == 'HI':
color = 'b'
else:
color = 'm'
ax.scatter(r, rom, 10, marker = 's', color=color)
invQg, invQs, invQeff_max = zip(*get_invQeff_from_data(gas_data=totgas,
epicycl=epicycl,
gas_approx=gas_approx,
sound_vel=sound_vel,
scale=scale,
sigma=sigma_R_min,
star_density=star_density))
ax.plot(r_g_dens[1:], invQeff_min, '-', alpha=0.5, color='r')
ax.plot(r_g_dens[1:], invQeff_max, '-', alpha=0.5, color='r')
plot_Q_levels(ax, [1., 1.5, 2., 3.])
ax.set_xlim(0)
ax.set_ylim(0)
ax.legend([matplotlib.lines.Line2D([0], [0], linestyle='none', mfc='g', mec='none', marker='o'),
matplotlib.lines.Line2D([0], [0], linestyle='none', mfc='b', mec='none', marker='o'),
matplotlib.lines.Line2D([0], [0], linestyle='none', mfc='m', mec='none', marker='o')],
['star', 'HI', 'H2'], numpoints=1, markerscale=1, loc='upper right') #add custom legend
ax.set_title('RF13: major component')
ax = plt.subplot(132)
ax.plot(romeo_min[1:], invQeff_min, 'o')
ax.plot(romeo_max[1:], invQeff_max, 'o', color='m', alpha=0.5)
ax.set_xlabel('Romeo')
ax.set_ylabel('2F')
ax.set_xlim(0., 1.)
ax.set_ylim(0., 1.)
ax.plot(ax.get_xlim(), ax.get_ylim(), '--')
ax = plt.subplot(133)
ax.plot(r_g_dens[1:], [l[1]/l[0] for l in zip(romeo_min[1:], invQeff_min)], 'o-')
ax.plot(r_g_dens[1:], [l[1]/l[0] for l in zip(romeo_max[1:], invQeff_max)], 'o-', color='m', alpha=0.5)
ax.set_xlabel('R')
ax.set_ylabel('[2F]/[Romeo]');
In [ ]:
Функция для исправления центральной поверхностной яркости за наклон (приведение диска к виду "плашмя"). Взята из http://www.astronet.ru/db/msg/1166765/node20.html, (61).
In [2]:
def mu_face_on(mu0d, cos_i):
return mu0d + 2.5*np.log10(1./cos_i)
In [ ]:
Функция для анализа влияния параметров. Берет стандартные параметры, делает из них списки и прогоняет для всех в списке, после чего измеряет среднее и std. Можно варьировать несколько параметров одновременно.
In [1]:
def plot_param_depend(ax=None, N=None, data_lim=None, color=None, alpha=0.3, disk_scales=[], label=None, max_range=False, **kwargs):
params = kwargs.copy()
for p in params.keys():
if p == 'total_gas_data':
depth = lambda L: isinstance(L, list) and max(map(depth, L))+1 #depth of nested lists
if depth(params[p]) == 1:
params[p] = [params[p]]*N
elif type(params[p]) is not list:
params[p] = [params[p]]*N
result = []
for i in range(N):
invQg, invQs, invQeff_min = zip(*get_invQeff_from_data(gas_data=params['total_gas_data'][i],
epicycl=params['epicycl'][i],
gas_approx=params['gas_approx'][i],
sound_vel=params['sound_vel'][i],
scale=params['scale'][i],
sigma=params['sigma_max'][i],
star_density=params['star_density_min'][i]))
invQg, invQs, invQeff_max = zip(*get_invQeff_from_data(gas_data=params['total_gas_data'][i],
epicycl=params['epicycl'][i],
gas_approx=params['gas_approx'][i],
sound_vel=params['sound_vel'][i],
scale=params['scale'][i],
sigma=params['sigma_min'][i],
star_density=params['star_density_max'][i]))
result.append((invQeff_min, invQeff_max))
rr = zip(*params['total_gas_data'][0])[0]
qmins = []
qmaxs = []
for ind, rrr in enumerate(rr):
qmin = [result[l][0][ind] for l in range(len(result))]
qmax = [result[l][1][ind] for l in range(len(result))]
if max_range:
qmins.append(((np.min(qmin)+np.max(qmin))/2., (np.max(qmin)-np.min(qmin))/2.))
qmaxs.append(((np.min(qmax)+np.max(qmax))/2., (np.max(qmax)-np.min(qmax))/2.))
else:
qmins.append((np.mean(qmin), np.std(qmin)))
qmaxs.append((np.mean(qmax), np.std(qmax)))
ax.errorbar(rr, zip(*qmins)[0], fmt='o-', yerr=zip(*qmins)[1], elinewidth=6, alpha=0.3);
ax.errorbar(rr, zip(*qmaxs)[0], fmt='o-', yerr=zip(*qmaxs)[1])
ax.axhline(y=1., ls='-', color='grey')
ax.set_ylim(0.)
ax.set_xlim(0.)
plot_data_lim(ax, data_lim)
plot_Q_levels(ax, [1., 1.5, 2., 3.]);
In [ ]:
In [ ]:
In [ ]: