NGC5440-checkpoint


Эксперименты по восстановлению профилей дисперсий в трех направлениях для NGC5440 (UGC8963)

Сначала всякие настройки и импорты


In [26]:
import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize
from math import *
from IPython.display import HTML
from IPython.display import Image
import os
import pandas as pd

#Размер изображений
import matplotlib.pylab as pylab
pylab.rcParams['figure.figsize'] = 12, 12

#Наклон галактики по данным NED
incl=65.7

# Масштаб пк/секунда из NED
scale=284

Всякие картинки и БД для большего удобства:


In [27]:
# Данные из SDSS DR9
HTML('<iframe src=http://skyserver.sdss3.org/dr9/en/tools/explore/obj.asp?ra=14:03:01.0&dec=+34:45:28 width=1000 height=350></iframe>')


Out[27]:

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


Out[28]:

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

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

#Выводим данные за header-а файла
for line in file("v_stars_ma.dat"):
    if line[0] == '#':
        print(line)


# UGC 8963 = NGC 5440

# Zasov+2012

# PA = 230 deg; PA_maj = 46.1 deg (from Mendez-Abreu+2008)

# not corrected for V_sys

# inclination = 65.7 deg (NED) (q_b = 0.36 - Mendez-Abreu+2008)

# data ARE NOT corrected for inclination

# r        vel       d_vel   sig     d_sig

# "        km/s      km/s    km/s    km/s

#(1)       (2)       (3)     (4)     (5)

Выведем также для удобства данные из обоих файлов:


In [31]:
ma = pd.read_csv('v_stars_ma.dat', skiprows=9, sep='   ', header=False)
# display(ma)

In [32]:
mi = pd.read_csv('v_stars_mi.dat', skiprows=9, sep='   ', header=False)
# display(mi)

Посмотрим теперь на все доступные даные по кривым вращения. Т.к. некоторые из них скорректированы, а некоторые нет, разобьем этот этап на несколько шагов.

Данные Коткова по газу вдоль главной оси:


In [33]:
n5440pa230_gas = pd.read_csv('.\\gas\\n5440pa230_gas_corr.txt', sep=' ')
display(n5440pa230_gas)


r," v_Hbeta e_V max_Hb e_max sig_Hb e_sig v_OIII e_v max_OIII_1 e_max.1 sig_OIII e_sig.1 [OIII]/Hbeta e_oiii_hb S/N(HB) S/N(OIII)
0 -87.70 3981.76 15.50 359.05 -0.0 0.00 4.04 3951.49 2.110000e+01 193.96 1.020000e+03 0.00 62.900 -1.00 -1.00E+00 3.19 2.20
1 -86.74 4002.55 12.50 368.47 -0.0 0.00 4.05 4035.75 5.760000e+00 384.73 -0.000000e+00 0.00 1.360 -1.00 -1.00E+00 3.30 2.71
2 -85.70 3992.79 9.37 493.95 -0.0 0.00 3.10 3147.44 -0.000000e+00 71.17 -0.000000e+00 0.00 -0.000 -1.00 -1.00E+00 4.79 0.00
3 -84.77 4005.56 9.22 408.37 -0.0 4.34 3.16 4220.34 -0.000000e+00 62.20 -0.000000e+00 0.00 -0.000 -Infinity NaN 6.31 0.00
4 -84.08 4018.82 9.46 417.86 -0.0 0.00 3.05 4005.34 1.120000e+01 297.59 -0.000000e+00 0.00 1.720 -1.00 -1.00E+00 5.47 3.55
5 -83.31 4016.21 8.95 476.42 -0.0 0.00 2.62 3990.03 1.680000e+01 93.35 1.440000e+02 0.00 39.000 -1.00 -1.00E+00 5.45 2.71
6 -82.43 4010.11 7.02 686.52 -0.0 0.00 2.27 3983.04 1.220000e+01 312.66 9.830000e+02 0.00 39.500 -1.00 -1.00E+00 7.51 3.80
7 -81.55 4004.57 7.17 561.60 -0.0 0.00 2.34 3993.67 1.060000e+01 256.50 7.820000e+02 0.00 40.300 -1.00 -1.00E+00 7.27 4.03
8 -80.85 3992.30 7.22 506.11 -0.0 13.48 2.69 3979.04 1.220000e+01 121.76 1.120000e+02 0.00 26.000 -Infinity Infinity 7.86 4.14
9 -80.15 4001.32 8.08 442.60 178.0 15.76 14.90 4010.71 7.590000e+00 356.30 -0.000000e+00 0.00 0.916 -Infinity Infinity 7.11 3.86
10 -79.43 3988.34 6.85 596.78 -0.0 0.00 2.24 3995.83 2.650000e+00 493.29 -0.000000e+00 0.00 1.730 -1.00 -1.00E+00 7.90 4.52
11 -78.52 3982.14 6.02 493.83 -0.0 0.00 2.09 3975.25 -0.000000e+00 107.92 -0.000000e+00 0.00 -0.000 -1.00 -1.00E+00 9.17 0.00
12 -78.17 3981.28 6.16 489.30 -0.0 0.00 2.08 3976.09 1.290000e+01 78.32 6.240000e+01 0.00 25.500 -1.00 -1.00E+00 9.75 4.17
13 -77.81 3974.14 6.12 472.32 -0.0 3.94 2.17 3975.79 1.150000e+01 256.14 -0.000000e+00 0.00 1.640 -Infinity Infinity 12.10 4.90
14 -77.31 3987.57 5.39 721.54 -0.0 7.99 1.94 3997.38 9.930000e+00 367.13 -0.000000e+00 0.00 1.340 -Infinity Infinity 10.94 5.14
15 -76.74 3982.30 10.90 304.18 -0.0 0.00 2.90 3992.46 1.280000e+01 156.85 -0.000000e+00 0.00 1.780 -1.00 -1.00E+00 5.05 2.64
16 -76.03 3978.26 11.30 733.50 -0.0 0.00 3.63 3915.38 1.880000e+01 205.92 4.750000e+02 0.00 49.200 -1.00 -1.00E+00 4.13 2.11
17 -71.86 3010.12 28.40 379.65 -0.0 0.00 3.82 3877.30 2.680000e+01 136.61 9.290000e+02 0.00 97.600 -1.00 -1.00E+00 1.14 1.51
18 -71.66 4285.72 1040000.00 396.80 285000000.0 0.00 1260000.00 4505.68 -0.000000e+00 163.97 -0.000000e+00 0.00 -0.000 -1.00 -1.00E+00 0.23 0.00
19 -62.81 2991.46 -0.00 271.20 -0.0 0.00 -0.00 3788.57 4.450000e+04 342.04 -0.000000e+00 0.00 28200.000 -1.00 -1.00E+00 0.00 0.63
20 -57.74 3070.66 25.30 353.95 1370.0 0.00 73.90 4775.18 5.280000e+01 64.65 3.690000e+01 96.99 55.500 -1.00 -1.00E+00 1.72 1.31
21 -55.78 3977.96 29.30 87.76 38.4 64.22 37.30 3950.74 2.740000e+01 34.97 1.390000e+01 66.48 34.700 0.05 4.25E-01 1.85 1.95
22 -55.11 3962.73 11.60 359.13 -0.0 0.00 3.50 4274.96 4.870000e+00 222.99 -0.000000e+00 0.00 2.770 -1.00 -1.00E+00 4.68 1.45
23 -54.60 3960.21 9.09 289.74 -0.0 10.31 3.40 3975.42 1.860000e+01 43.00 2.530000e+01 31.29 29.100 0.09 4.99E-01 6.61 2.87
24 -54.25 3962.05 6.30 424.24 120.0 19.15 11.10 3950.61 6.720000e+00 290.98 -0.000000e+00 0.00 1.320 -Infinity Infinity 8.71 4.69
25 -53.89 3954.81 5.07 614.53 -0.0 15.92 1.95 3955.54 9.180000e+00 83.08 1.820000e+01 42.92 13.200 -0.01 1.72E-01 13.30 5.43
26 -53.53 3961.16 5.40 525.09 105.0 26.17 9.06 3951.55 6.200000e+00 186.76 9.340000e+01 0.00 13.600 -Infinity Infinity 13.51 8.38
27 -53.18 3952.24 6.13 356.81 54.3 42.35 8.85 3963.15 7.320000e+00 145.32 7.010000e+01 0.00 14.700 -Infinity Infinity 10.18 6.84
28 -52.82 3964.49 8.56 281.20 96.8 23.36 14.50 3980.27 9.210000e+00 87.20 2.360000e+01 34.40 14.000 0.09 3.74E-01 7.81 5.20
29 -52.25 3974.65 9.16 548.41 -0.0 0.00 2.86 3954.93 1.160000e+01 114.91 3.890000e+01 34.47 17.800 -1.00 -1.00E+00 4.87 3.82
30 -50.95 3977.48 22.20 162.84 55.2 63.94 28.60 3973.95 1.920000e+01 231.61 1.110000e+03 0.00 59.800 -Infinity Infinity 2.42 2.32
31 -48.47 3680.55 61.30 139.66 214.0 40.62 91.60 2834.98 9.620000e+08 400.10 1.100000e+10 0.00 15500000.000 -Infinity Infinity 0.75 1.05
32 -42.27 3028.66 68.20 73.78 122.0 41.67 100.00 4567.32 -0.000000e+00 139.53 -0.000000e+00 0.00 -0.000 -Infinity NaN 0.68 0.00
33 -39.82 2997.72 50.10 419.09 -0.0 0.00 6.90 3999.50 1.480000e+01 156.53 3.430000e+01 64.62 18.800 -1.00 -1.00E+00 0.82 3.16
34 -37.07 3053.81 38.50 347.09 -0.0 0.00 4.66 3973.26 2.470000e+01 70.46 1.740000e+01 91.19 28.600 -1.00 -1.00E+00 0.94 2.17
35 -35.28 3089.00 -0.00 298.27 -0.0 0.00 -0.00 3926.53 2.420000e+01 72.25 1.480000e+01 105.21 27.100 -1.00 -1.00E+00 0.00 2.30
36 -33.72 3969.43 23.90 327.74 -0.0 0.00 4.88 3206.46 -0.000000e+00 81.44 -0.000000e+00 0.00 -0.000 -1.00 -1.00E+00 2.13 0.00
37 -32.65 3996.09 21.60 248.39 605.0 0.00 55.20 3940.43 2.900000e+01 39.31 9.780000e+00 103.70 32.400 -1.00 -1.00E+00 2.24 2.19
38 -31.93 3973.42 17.40 357.25 -0.0 0.00 3.88 3963.85 1.100000e+01 327.38 -0.000000e+00 0.00 1.560 -1.00 -1.00E+00 2.88 3.54
39 -31.01 3933.13 19.20 399.05 -0.0 0.00 4.28 3958.86 1.240000e+01 204.45 2.690000e+02 0.00 30.400 -1.00 -1.00E+00 2.26 3.72
40 -30.16 4281.87 42.50 212.06 -0.0 0.00 5.40 3953.55 1.260000e+01 234.91 6.510000e+02 0.00 42.000 -1.00 -1.00E+00 0.68 3.92
41 -29.29 2829.64 47.10 199.53 -0.0 0.00 8.09 3948.80 1.300000e+01 117.72 5.650000e+01 24.58 20.900 -1.00 -1.00E+00 1.19 4.26
42 -28.30 3946.72 25.80 446.96 -0.0 0.00 4.97 3923.48 1.920000e+01 98.95 2.380000e+01 75.49 23.500 -1.00 -1.00E+00 1.67 2.69
43 -26.70 2996.52 -0.00 150.43 -0.0 0.00 -0.00 3940.97 2.380000e+01 61.73 1.320000e+01 99.61 26.900 -1.00 -1.00E+00 0.00 2.45
44 -25.51 3912.59 32.40 133.06 46.2 85.18 38.10 3903.50 2.500000e+01 67.44 1.500000e+01 100.64 28.300 0.21 2.91E-01 1.85 2.22
45 -23.93 3901.36 23.80 207.85 50.8 88.91 27.80 3960.45 1.360000e+01 141.74 2.590000e+01 69.50 16.700 0.16 2.16E-01 2.89 3.89
46 -22.61 3978.73 50.20 66.13 25.8 113.87 55.60 3957.97 1.920000e+01 65.40 1.530000e+01 76.42 23.100 0.25 3.18E-01 1.37 3.07
47 -21.94 4011.88 27.30 240.56 -0.0 0.00 6.32 4029.70 1.470000e+01 86.90 1.590000e+01 75.50 17.700 -1.00 -1.00E+00 1.91 3.99
48 -21.19 2725.71 -0.00 193.97 -0.0 0.00 -0.00 4019.80 8.910000e+00 163.18 3.830000e+01 37.68 13.100 -1.00 -1.00E+00 0.00 5.85
49 -20.69 3410.87 -0.00 104.13 -0.0 0.00 -0.00 4046.70 1.330000e+01 71.13 1.390000e+01 65.57 16.600 -1.00 -1.00E+00 0.00 4.80
50 -20.33 4039.30 26.90 137.49 264.0 0.00 57.00 4038.79 9.920000e+00 99.27 1.760000e+01 55.07 13.000 -1.00 -1.00E+00 2.36 7.15
51 -19.85 4002.58 18.40 204.16 83.7 44.68 26.00 4040.40 7.890000e+00 187.62 3.500000e+01 42.25 11.300 0.37 3.40E-01 2.80 6.93
52 -19.26 4255.62 -0.00 133.25 -0.0 0.00 -0.00 4028.80 1.080000e+01 92.09 1.480000e+01 64.56 13.500 -1.00 -1.00E+00 0.00 5.84
53 -18.73 2777.35 -0.00 197.93 -0.0 0.00 -0.00 4026.65 8.710000e+00 190.08 5.710000e+01 27.93 14.000 -1.00 -1.00E+00 0.00 5.48
54 -18.08 3854.23 26.30 152.86 112.0 36.39 40.20 4006.70 1.570000e+01 90.15 2.130000e+01 63.35 19.900 0.44 6.00E-01 2.24 3.46
55 -17.19 4730.15 -0.00 272.11 -0.0 0.00 -0.00 3948.14 1.160000e+01 151.50 2.510000e+01 66.80 14.500 -1.00 -1.00E+00 0.00 4.98
56 -16.28 3894.67 30.40 341.85 -0.0 0.00 4.69 3942.50 9.750000e+00 156.26 2.790000e+01 53.52 13.000 -1.00 -1.00E+00 1.58 5.39
57 -15.69 3898.68 24.80 160.40 -0.0 0.00 8.33 3951.13 1.250000e+01 89.30 2.260000e+01 48.82 17.000 -1.00 -1.00E+00 2.41 5.29
58 -15.20 3929.62 19.60 235.32 160.0 28.23 32.00 3940.43 7.560000e+00 210.71 2.920000e+01 53.28 10.000 0.66 5.82E-01 2.79 7.36
59 -14.42 3907.44 17.90 479.03 -0.0 0.00 3.76 3925.23 6.780000e+00 244.70 3.390000e+01 48.52 9.270 -1.00 -1.00E+00 2.50 6.94
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...

194 rows × 17 columns


In [34]:
plt.plot(n5440pa230_gas['r,"'], n5440pa230_gas['v_Hbeta'], 's')
plt.plot(n5440pa230_gas['r,"'], n5440pa230_gas['v_OIII'], 'x')
plt.show()



In [35]:
try:
    from PIL import Image
except:
    import Image

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

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

# Данные по кинематике газа Засова 2012 вдоль большой полуоси, не исправленные за наклон (они же Катков)
# Водород
r_g_H = n5440pa230_gas['r,"']
vel_g_H = n5440pa230_gas['v_Hbeta']
e_vel_g_H = n5440pa230_gas['e_V']

# Кислород
r_g_O = n5440pa230_gas['r,"']
vel_g_O = n5440pa230_gas['v_OIII']
e_vel_g_O = n5440pa230_gas['e_v']

fig, subplot = subplots(2, 1)
subplot[0].plot(r_g_H, vel_g_H, '^', label="gas Hbeta Zasov 2012", mfc='none')
subplot[0].plot(r_g_O, vel_g_O, 's', label="gas OIII Zasov 2012", mfc='none')
subplot[0].plot(r_ma, vel_ma, '.-', label="Zasov 2012, maj")
subplot[0].plot(r_mi, vel_mi, '.-', label="Zasov 2012, min")
subplot[0].legend(loc = 'lower left')
subplot[0].set_ylim(3300, 4100)

subplot[1].imshow(np.asarray(Image.open("zasov2012_fig4_part.png")))
subplot[1].set_title("From Z2012 fig.4 part")
plt.plot()


Out[35]:
[]

Попробуем как-нибудь автоматически вычистить значения по газу:


In [36]:
def bad_gas_data(data, diff=100):   
    pos = []
    neg = []
    level = None
    for d in filter(lambda l: l[0] >= 0, data):
        if level is None:
            level = d[1]
        elif abs(level - d[1]) < diff:
            level = d[1]
        else:
            pos.append(d)
    level = None
    for d in filter(lambda l: l[0] < 0, data):
        if level is None:
            level = d[1]
        elif abs(level - d[1]) < diff:
            level = d[1]
        else:
            neg.append(d)
    return neg + pos

H_data = zip(r_g_H, vel_g_H)
bad_H_r, bad_H_v = zip(*bad_gas_data(H_data))

O_data = filter(lambda l: l[0] < -32. or l[0] > 0., zip(r_g_O, vel_g_O))
bad_O_r, bad_O_v = zip(*bad_gas_data(O_data, diff=220.))

fig, subplot = subplots(2, 1)
subplot[0].plot(r_g_H, vel_g_H, '^', label="gas Hbeta Zasov 2012", mfc='none')
subplot[0].plot(r_g_O, vel_g_O, 's', label="gas OIII Zasov 2012", mfc='none')
subplot[0].plot(r_ma, vel_ma, '.', label="Zasov 2012, maj")
subplot[0].plot(r_mi, vel_mi, '.-', label="Zasov 2012, min")
subplot[0].plot(bad_H_r, bad_H_v, '^m')
subplot[0].plot(bad_O_r, bad_O_v, 'sm')
subplot[0].legend(loc = 'lower left')

# move legend outside of plot
box = subplot[0].get_position()
subplot[0].set_position([box.x0, box.y0, box.width * 0.8, box.height])
subplot[0].legend(loc='center left', bbox_to_anchor=(1, 0.5))


subplot[1].imshow(np.asarray(Image.open("zasov2012_fig4_part.png")))
subplot[1].set_title("From Z2012 fig.4 part")
plt.plot()


Out[36]:
[]

In [37]:
# Вычищаем и дочищаем:
good_H_data = zip(r_g_H, vel_g_H)
good_H_data = [q for q in good_H_data if q not in zip(bad_H_r, bad_H_v)]
good_H_data = filter(lambda l: l[0] < 51., good_H_data)
good_H_r, good_H_v = zip(*good_H_data)

good_O_data = zip(r_g_O, vel_g_O)
good_O_data = [q for q in good_O_data if q not in zip(bad_O_r, bad_O_v)]
good_O_data = filter(lambda l: l[0] < 51., good_O_data)
good_O_data = filter(lambda l: (l[0] < -61. and abs(l[1]-4000) < 20.) or l[0] > -61., good_O_data)
good_O_r, good_O_v = zip(*good_O_data)

fig, subplot = subplots(2, 1)
subplot[0].plot(good_H_r, good_H_v, '^', label="gas Hbeta Zasov 2012", mfc='none')
subplot[0].plot(good_O_r, good_O_v, 's', label="gas OIII Zasov 2012", mfc='none')
subplot[0].plot(r_ma, vel_ma, '.', label="Zasov 2012, maj")
subplot[0].plot(r_mi, vel_mi, '.-', label="Zasov 2012, min")
subplot[0].legend(loc = 'lower left')
subplot[0].set_ylim(3300, 4100)
subplot[0].set_xlim(-100, 60)

subplot[1].imshow(np.asarray(Image.open("zasov2012_fig4_part.png")))
subplot[1].set_title("From Z2012 fig.4 part")
plt.plot()


Out[37]:
[]

Теперь построим график дисперсий скоростей на луче зрения вдоль большой и малой осей по данным Засова:


In [38]:
plt.plot(r_ma, sig_ma, 's-', label='$\sigma_{los}^{maj}$')
plt.errorbar(r_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='.', mew=0, color='blue')
plt.plot(r_mi, sig_mi, 's-', label='$\sigma_{los}^{min}$')
plt.errorbar(r_mi, sig_mi, yerr=e_sig_mi, fmt='.', marker='.', mew=0, color='black')
plt.xlabel('$R$')
plt.ylabel('$\sigma$')
plt.legend()
plt.show()


Что мы видим? Даже без исправления значения вдоль малой оси на косинус угла видно, что центры двух кривых совсем не совпали, а это значит, что щели спектрометра не были точно соорентированы на центр и/или по осям галактики. Однако рассматриваемый случай оказывается достаточно простым - кривая вращения вдоль большой - симметрична, и вдоль малой тоже, но смещена - значит ${PA}_{min}$ и ${PA}_{maj}$ были определены правильно, просто щель для малой оси оказалась немного смещена. Данные и по скоростям, и по дисперсиям необходимо исправить, методика исправления описана ниже.


In [39]:
from mpl_toolkits.axes_grid.axislines import SubplotZero
import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure(1)
ax = SubplotZero(fig, 111)
fig.add_subplot(ax)

for direction in ["xzero", "yzero"]:
    ax.axis[direction].set_axisline_style("-|>")
    ax.axis[direction].set_visible(True)

for direction in ["left", "right", "bottom", "top"]:
    ax.axis[direction].set_visible(False)

x = np.linspace(-1., 1., 100)
y = np.zeros(100)
ax.plot(x, y)

from pylab import figure, show, rand
from matplotlib.patches import Ellipse
e = Ellipse(xy=(0,0), width=1.9, height=1.0, angle=0.0)

ax.add_artist(e)
e.set_clip_box(ax.bbox)
e.set_alpha(0.1)
e.set_facecolor('blue')

ax.set_ylim(-0.75, 0.75)
ax.axes.get_xaxis().set_ticks([])
ax.axes.get_yaxis().set_ticks([])

ax.text(0.98, 0.02, "$x$", fontsize=20)
ax.text(0.02, 0.73, "$y$", fontsize=20)
ax.plot((0.0, 0.2), (0.0, 0.2), '--', lw=0.3)
ax.text(0.05, 0.015, r"$\varphi$", fontsize=20)
ax.plot((0.2, 0.2), (-0.49, 0.49), '-', lw=1.5, color='r')
ax.plot((-0.87, 0.87), (0.2, 0.2), '-', lw=1.5, color='r')
ax.plot(0.2, 0.2, 's')
ax.text(0.01, 0.11, "$y_0$", fontsize=20)
ax.text(0.11, -0.03, "$x_0$", fontsize=20)
ax.plot((0.0, 0.2), (0.0, 0.0), '-', lw=2.0, color='k')
ax.plot((0.0, 0.0), (0.0, 0.2), '-', lw=2.0, color='k')
ax.text(0.22, -0.3, "$P_{obs}$", fontsize=20)
ax.arrow(0.2, 0., 0., -0.28, head_width=0.02, head_length=0.02, fc='k', ec='k', width=0.005)
ax.plot((0., 0.2), (0., -0.3), '--', lw=0.3)
ax.text(0.06, -0.06, r"$\theta$", fontsize=20)
ax.text(0.2, -0.16, "$R_{obs}$", fontsize=20)
ax.plot(0.2, -0.3, 's')

plt.show()


Пусть центр координат находится в истинном центре галактики, углы осей ${PA}_{maj}$ и ${PA}_{min}$ определены правильно, но щели спектрометра пересекаются в точке $(x_0,y_0)$. Угол $\varphi$ очевидно равен $$\varphi=\arctan \frac{y_0}{x_0}$$ Тогда расстояния и скорости для точки $P_{obs}$ для соответствующей оси (на рисунке выше - для большой) пересчитываются по следующему правилу: $$R_{real} = R_{obs}\times\frac{\sqrt{\sin^2\theta\cos^2i+\cos^2\theta}}{\cos i\sin \theta},$$ $$V_{real} = V_{obs}\times\frac{\sqrt{\sin^2\theta\cos^2i+\cos^2\theta}}{\sin i\cos i\sin \theta},$$ где $\theta$ есть угол на соответствующую точку.

В нашем случае из-за того, что $\sigma_{los, maj}(0) > \sigma_{los, min}(0)$, а также по причинам, которые были указаны выше, смещена щель только вдоль малой оси, угол $\varphi=0^{\circ}$. Угол $\theta$ и его тригонометрические функции в свою очередь считаются как $$\theta = \arctan \frac{R_{obs}}{x_0}, \sin \theta = \frac{R_{obs}}{\sqrt{R_{obs}^2 + x_0^2}}, \cos \theta = \frac{x_0}{\sqrt{R_{obs}^2 + x_0^2}}$$ После восстановления скоростей исправлять за наклон галактики профиль уже не нужно (случай $\theta=0^{\circ}$).

Остается лишь вопрос, как найти $y_0$. В рассматриваемом случае это сделать просто - достаточно взять максимальное значение скорости вдоль малой оси и найти схожее значениее на большой.


In [40]:
r_ma_b, vel_ma_b, e_vel_b = r_ma, map(lambda l: abs(l-3715.0), vel_ma), e_vel_ma
r_mi_b, vel_mi_b, e_vel_mi_b = r_mi, map(lambda l: abs(l-3715.0), vel_mi), e_vel_mi

plt.plot(r_ma_b, vel_ma_b, 'd', label = 'Zasov star maj')
plt.errorbar(r_ma_b, vel_ma_b, yerr=e_vel_b, fmt='.', marker='.', mew=0, color='blue')
plt.plot(r_mi_b, vel_mi_b, '.', label = 'Zasov star min', color='green')
plt.errorbar(r_mi_b, vel_mi_b, yerr=e_vel_mi_b, fmt='.', marker='.', mew=0, color='green')
plt.axvline(x = 9.)
plt.axvline(x = -7.)
plt.axhline(y =182.)
plt.xlim(-50., 50.)
plt.text(-40, 150, r"$x_0 = \frac{9+7}{2}=8$", fontsize=25)
plt.legend()
plt.plot()


Out[40]:
[]

In [41]:
x_0 = 8.*cos(incl * pi / 180)

def sin_theta(R_obs):
    return R_obs/sqrt(R_obs**2 + x_0**2)

def cos_theta(R_obs):
    return x_0/sqrt(R_obs**2 + x_0**2)

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

def correct_radii_min(R_obs):
    return R_obs * sqrt(sin_theta(R_obs)**2 * cos_i**2 + cos_theta(R_obs)**2) / cos_i /  abs(sin_theta(R_obs))

In [42]:
r_mi1 = map(correct_radii_min, r_mi)

In [43]:
def incline_velocity(v, angle):
    return v / sin(angle * pi / 180)

# Переносит центр в (r0,v0) и перегибает кривую вращения, 
# а также исправляет за наклон если необходимо
def correct_rotation_curve(rdata, vdata, dvdata, r0, v0, incl):
    rdata_tmp = [abs(r-r0) for r in rdata]
    vdata_tmp = [incline_velocity(abs(v-v0), incl) for v in vdata]
    data = zip(rdata_tmp, vdata_tmp, dvdata)
    data.sort()
    return zip(*data)

r_ma_b, vel_ma_b, e_vel_b = correct_rotation_curve(r_ma, vel_ma, e_vel_ma,  0., 3715., incl)
r_mi_b, vel_mi_b, e_vel_mi_b = correct_rotation_curve(r_mi1, vel_mi, e_vel_mi,  0., 3715., 90.)

plt.plot(r_ma_b, vel_ma_b, 'd', label = 'Zasov star maj')
plt.errorbar(r_ma_b, vel_ma_b, yerr=e_vel_b, fmt='.', marker='.', mew=0, color='blue')
plt.plot(r_mi_b, vel_mi_b, '.', label = 'Zasov star min', color='green')
plt.errorbar(r_mi_b, vel_mi_b, yerr=e_vel_mi_b, fmt='.', marker='.', mew=0, color='green')
plt.legend()
plt.ylim(0)
plt.plot()


Out[43]:
[]

В дальнейшем используем только засовские данные по звездам по большой полуоси, приблизим их полиномом.


In [44]:
poly_star = poly1d(polyfit(r_ma_b, vel_ma_b, deg=7))

plt.plot(r_ma_b, vel_ma_b, 'x-', color='blue', markersize=6)
test_points = np.arange(0.0, max(r_ma_b), 0.1)
plt.plot(test_points, poly_star(test_points), '-', color='red')
plt.xlabel('$R$'); plt.ylim(0)
plt.ylabel('$V^{maj}_{\phi}(R)$')
plt.show()


Кривая вращения нам нужна для нахождения соотношения $\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) и приближается гладко функцией $f=0.5(1+e^{-R/R_{0}}),$ где $R_{0}$ --- характерный масштаб.

${\bf Примечание:}$ Такое приближение оправдано следующими соображениями. Для равновесного диска верно уравнение, описанное выше. Для твердотельного участка вращения в центральных областях выражение в скобках равно 2, а $\sigma_{\varphi}^{2}/\sigma_{R}^{2}=1$. На плоском участке кривой вращения на периферии диска $\sigma_{\varphi}^{2}/\sigma_{R}^{2}\thickapprox0.5$. Функция $f$ как раз аппроксимирует такое поведение отношения $\sigma_{\varphi}^{2}/\sigma_{R}^{2}$.

Изобразим получившийся профиль $\sigma_{\varphi}^{2}/\sigma_{R}^{2}$, вычисляемый через производную полинома:


In [45]:
def sigPhi_to_sigR_real(R):
        return 0.5 * (1 + R*poly_star.deriv()(R) / poly_star(R))

plt.plot(test_points, [sigPhi_to_sigR_real(R) for R in test_points], 'd-', color='blue')
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)
plt.show()


Найдем теперь характерный масштаб $f=0.5(1+e^{-R/R_{0}})$:


In [46]:
def f(R, Ro):
    return 0.5*(1 + np.exp( -R/Ro ))

xdata = test_points
ydata = sigPhi_to_sigR_real(xdata)

from scipy.optimize import curve_fit
popt, pcov = curve_fit(f, xdata, ydata, p0=[1.0])
Ro = popt[0]

plt.plot(xdata, ydata, 'x-')
plt.plot(xdata, [f(p, Ro) for p in xdata], 's')
plt.axhline(y=0.5)
plt.axhline(y=0.0)
plt.title('$R_{0} = %s $' % Ro)
plt.ylim(0)
plt.show()


Теперь знаем значение отношения $\sigma_{\varphi}^{2}/\sigma_{R}^{2}$ в любой точке, заведем соответствующую функцию:


In [47]:
def sigPhi_to_sigR(R):
    return sqrt(f(R, Ro))

Построим графики дисперсий скоростей на луче зрения вдоль большой и малой оси ($\sigma_{los}^{maj}$ и $\sigma_{los}^{min}$), пересчитав расстояния для малой оси как было описано ранее:


In [48]:
# Исправляем значения вдоль малой оси на синус угла:    
def correct_min(R):    
    return R / cos(incl * pi / 180) 

# # r_mi_extend = map(correct_min, r_mi)
r_mi_extend = map(correct_radii_min, map(correct_min, r_mi))

plt.plot(r_ma, sig_ma, 's-', label='$\sigma_{los}^{maj}$')
plt.errorbar(r_ma, sig_ma, yerr=e_sig_ma, fmt='.', marker='.', mew=0, color='blue')
plt.plot(r_mi_extend, sig_mi, 's-', label='$\sigma_{los}^{min}$')
plt.errorbar(r_mi_extend, sig_mi, yerr=e_sig_mi, fmt='.', marker='.', mew=0, color='black')
plt.xlabel('$R$')
plt.ylabel('$\sigma$')
plt.legend()
plt.xlim(-50,50)
plt.show()



In [49]:
bind_curve = lambda p: (abs(p[0]), abs(p[1]), p[2])
sig_maj_data = zip(r_ma, sig_ma, e_sig_ma)
sig_maj_data = map(bind_curve, sig_maj_data)
sig_maj_data.sort()
radii_maj, sig_maj_p, e_sig_maj_p = zip(*sig_maj_data) 

poly_sig_maj = poly1d(polyfit(radii_maj, sig_maj_p, deg=9))

sig_min_data = zip(r_mi_extend, sig_mi, e_sig_mi)
sig_min_data = map(bind_curve, sig_min_data)
sig_min_data.sort()
radii_min, sig_min_p, e_sig_min_p = zip(*sig_min_data) 

# Добавляем лишние точки чтобы протянуть дальше
# num_fake_points = 10; expscale = 200.0
# fake_radii, fake_sig = zip(*[(31.0 + i, 115*exp(- i / expscale )) for i in range(1, num_fake_points+1)])
fake_radii, fake_sig = (),()

poly_sig_min = poly1d(polyfit(radii_min + fake_radii, sig_min_p + fake_sig, deg=3))

points = np.arange(0, max(radii_min), 0.1)
plt.plot(radii_maj, sig_maj_p, 's', label='$\sigma_{los}^{maj}$', color='blue')
plt.errorbar(radii_maj, sig_maj_p, yerr=e_sig_maj_p, fmt='.', marker='.', mew=0, color='blue')
plt.plot(points, poly_sig_maj(points), label = '$\sigma_{los}^{maj} polyfit$', color='blue')
plt.plot(radii_min, sig_min_p, 's', label='$\sigma_{los}^{min}$', color='red')
plt.errorbar(radii_min, sig_min_p, yerr=e_sig_min_p, fmt='.', marker='.', mew=0, color='red')
plt.plot(points, poly_sig_min(points), label = '$\sigma_{los}^{min} polyfit$', color='red')
plt.plot(fake_radii, fake_sig, 'bs', color='green', label='$fake points$')
plt.legend()
plt.ylim(70,250)
plt.xlim(0,50)
plt.show()


Посчитаем величину невязок для полученного приближения:


In [51]:
sqerr_maj = sum(power([poly_sig_maj(p[0]) - p[1] for p in sig_maj_data], 2))
sqerr_min = sum(power([poly_sig_min(p[0]) - p[1] for p in sig_min_data], 2))

print "Poly chi^2 for maj = %s, mean = %s" % (sqerr_maj, sqerr_maj / sig_maj_p.__len__())
print "Poly chi^2 for min = %s, mean = %s" % (sqerr_min, sqerr_min / sig_min_p.__len__())


Poly chi^2 for maj = 4322.98401725, mean = 61.7569145321
Poly chi^2 for min = 2867.43605769, mean = 77.4982718295

Методика восстановления профилей $\sigma_{R}(R)$, $\sigma_{\varphi}(R)$ и $\sigma_{z}(R)$ следующая. Представим, что $\sigma_{Z}/\sigma_{R} \equiv \alpha \equiv const$. Тогда, зная значения $\sigma_{\varphi}^{2}/\sigma_{R}^{2}=f(R)$ в каждой точке, получаем из уравнений, описанных выше: $$\sigma_{los,maj}^2=\sigma_R^2[f\sin^2i+\alpha^2\cos^2i]$$ $$\sigma_{los,min}^2=\sigma_R^2[\sin^2i+\alpha^2\cos^2i]$$ Представим теперь $\sigma_R(R)=\sigma_{R,0}\times F(R)$, где $F(0)=1$. Значение в квадратных скобках для $\sigma_{los,min}$ равно константе и, следуя предположению, получаем представление для дисперсии вдоль луча зрения для малой оси как $\sigma_{los,min}(R)=\sigma_{min,0}\times F(R)$. Очевидно $\sigma_{min,0} = \sigma_{los,min}(0)$, а значит мы знаем в каждой точке значение $F(R)=\sigma_{los,min(R)}/\sigma_{min,0}$. Описанная выше система уравнений вырождается в следующую: $$\sigma_{los,maj}^2(R)=\frac{\sigma_{R,0}^2\sigma_{los,min}^2(R)[f\sin^2i+\alpha^2\cos^2i]}{\sigma_{min,0}^2}$$ $$\sigma_{min,0}^2=\sigma_{R,0}^2\sin^2i+\sigma_{R,0}^2\alpha^2\cos^2i$$ Сделаем замену: $\sigma_{R,0}^2\sin^2 i \equiv A,\ \sigma_{R,0}^2\cos^2 i\times \alpha^2 \equiv B$. Окончательно, имеем $N+1$ линейное уравнение для $N$ точек, которые можем решить МНК: $$\left\{ \begin{array}{lr} \sigma_{los,maj}^2(R_j)\times \sigma_{min,0}^2 =\sigma_{los,min}^2(R_j)[Af(R_j)+B]\\ \sigma_{min,0}^2=A+B \end{array} \right. $$

Однако это все было верно для случая несмещенной щели. Для нашего случая исходное уравнение для дисперсий вдоль малой оси меняется на $$\sigma_{los,min}^2(R_{obs})=\sigma_R^2(R)[\cos^2\theta+f(R)\sin^2\theta]\sin^2i + \alpha^2\sigma_R^2(R)\cos^2i,$$ Соответственно неправомерно делать одно из ключевых предположений - $\sigma_{los,min}(R)=\sigma_{R}(R)\times const$ и систему нельзя записать как линейную. Но можно поступить следующим образом: сосчитаем разность $$\sigma_{los,min}^2(R_{obs}) - \sigma_{los,maj}^2\sin^2\theta = \sigma_R^2\cos^2\theta[\sin^2i + \alpha^2\cos^2i]$$ Исходя из этого выражения можем легко построить график невязок для большой оси для разных $\alpha$. (НЕ УВЕРЕH ЧТО ТАМ ДЕЛАТЬ С Robs!!)


In [74]:
alpha = 0.4
# Возьмем набор каких-то точек
points_obs = np.arange(8., 30., 0.5)
# Скорректируем
points_real = map(correct_radii_min, points_obs)
# Углы
cos2_theta = map(np.square, map(cos_theta, points_obs))
sin2_theta = map(np.square, map(sin_theta, points_obs))

# Вот не уверен я тут в obs или real - вроде бы мы уже исправили малую
sigR2_difr = [(poly_sig_min(p[0])**2 - poly_sig_maj(p[0])**2 * p[2]) / (p[1] * (sin_i**2 + alpha**2 * cos_i**2))  
              for p in zip(points_real, cos2_theta, sin2_theta)]

# plt.plot(points_real, map(sqrt, sigR2_difr), 'o')
# plt.plot(points_real, sigR2_difr, 'o')
# plt.plot(points_real, [(poly_sig_min(p[0])**2 - poly_sig_maj(p[0])**2 * p[1]) for p in zip(points_real, sin2_theta)], 'o')
# plt.plot(points_real, [(poly_sig_min(p[0])**2 - poly_sig_maj(p[1])**2 * p[2]) for p in zip(points_obs, points_real, sin2_theta)], 'o')

# plt.plot(radii_maj, sig_maj_p, 's', label='$\sigma_{los}^{maj}$', color='blue')
# plt.errorbar(radii_maj, sig_maj_p, yerr=e_sig_maj_p, fmt='.', marker='.', mew=0, color='blue')
# plt.plot(points, poly_sig_maj(points), label = '$\sigma_{los}^{maj} polyfit$', color='blue')
plt.plot(points_real, map(poly_sig_maj, points_real), 'x')
plt.plot(points_real, [sqrt(p[1]*(sigPhi_to_sigR(p[0])**2 * sin_i**2 + alpha**2 + cos_i**2)) for p in zip(points_real, sigR2_difr)], 'o')
plt.show()


Выклядит все это удручающе, но посчитать график мы все же можем.


In [93]:
alphas = np.arange(0.1, 0.9, 0.01)

sig_maj_difr_errs = []
for al in alphas: 
    sigR2_difr = [(poly_sig_min(p[0])**2 - poly_sig_maj(p[0])**2 * p[2]) / (p[1] * (sin_i**2 + al**2 * cos_i**2))  
              for p in zip(points_real, cos2_theta, sin2_theta)]
    sig2_maj_difr = [sqrt(p[1]*(sigPhi_to_sigR(p[0])**2 * sin_i**2 + al**2 + cos_i**2)) for p in zip(points_real, sigR2_difr)]
    sqerr_maj = sum(power([poly_sig_maj(p[0]) - p[1] for p in zip(points_real, sig2_maj_difr)], 2))/len(points_real)
    sig_maj_difr_errs.append(sqerr_maj)

plt.plot(alphas, sig_maj_difr_errs, '^')
plt.xlabel(r'$\alpha$')
plt.ylabel('$err$')
plt.show()


Проверить, что математика правильная!


In [26]:
#Новая функция f'
def f_prime(R_real, R_obs):
    return f(R_real, Ro)*sin_theta(R_obs)**2 + cos_theta(R_obs)**2

#ВОТ НЕ ПОНЯТНО - ВРОДЕ poly_sig_min(R_real)
def sigR_exp(R_real, R_obs, alp):
    return sqrt(poly_sig_min(R_real)/(f_prime(R_real, R_obs)*sin_i**2 + alp**2*cos_i**2))

def sigZ_exp(R_real, R_obs, alp):
    return alp * sigR_exp(R_real, R_obs, alp)

def sigPhi_exp(R_real, R_obs, alp):
    return sigPhi_to_sigR(R_real) * sigR_exp(R_real, R_obs, alp)

def sig_min_exp(R_real, R_obs, alp):
    return sqrt(sigR_exp(R_real, R_obs, alp)**2*f_prime(R_real, R_obs)*sin_i**2 + sigZ_exp(R_real, R_obs, alp)**2 * cos_i**2)

def sig_maj_exp(R_real, R_obs, alp):
    return sqrt(sigPhi_exp(R_real, R_obs, alp)**2 * sin_i**2 + sigZ_exp(R_real, R_obs, alp)**2 * cos_i**2)

Проверим, что все правильно. Построим профили для $\alpha=0.5$:


In [27]:
alpha = 0.5
points_obs = np.arange(0, max(radii_min), 0.1)
points_real = map(correct_radii_min, points_obs)

plt.plot(points_real, [sigR_exp(Rr, Ro, alpha) for (Rr, Ro) in zip(points_real, points_obs)], '-', color='red', label='$\sigma_{R, exp}$')
plt.plot(points_real, [sigPhi_exp(Rr, Ro, alpha) for (Rr, Ro) in zip(points_real, points_obs)], 
         '-', color='blue', label=r'$\sigma_{\varphi, exp}$')
plt.plot(points_real, [sigZ_exp(Rr, Ro, alpha) for (Rr, Ro) in zip(points_real, points_obs)], 
         '-', color='black', label='$\sigma_{Z, exp}$')
plt.legend()
plt.show()



In [28]:
plt.plot(points, poly_sig_maj(points), '-', label = '$\sigma_{los}^{maj} polyfit$', color='blue')
plt.plot(points, [sig_maj_exp(Rr, Robs, alpha) for (Rr, Robs) in zip(points_real, points_obs)], 
         '--', color='blue', label='$\sigma_{maj, exp}$')
plt.plot(points, poly_sig_min(points), '-', label = '$\sigma_{los}^{min} polyfit$', color='red')
plt.plot(points_real, [sig_min_exp(Rr, Robs, alpha) for (Rr, Robs) in zip(points_real, points_obs)], 
         '--', color='red', label='$\sigma_{min, exp}$')
plt.legend()
# plt.ylim(0, 300)
plt.show()



In [28]:


In [28]:


In [28]:

Теперь восстановим исходные неизвестные - $\alpha$ и $\sigma_{R, 0}$:


In [29]:
sig_R_0 = round( sqrt(A) / sin(incl*pi/180), 3)
alpha = round( sqrt(B)/ (cos(incl*pi/180) * sig_R_0), 3)

# sig_R_0 = 160
# alpha = 0.55

print "sig_R_0 = %s, alpha = %s" % (sig_R_0, alpha)


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-29-6acd7ae6cb97> in <module>()
----> 1 sig_R_0 = round( sqrt(A) / sin(incl*pi/180), 3)
      2 alpha = round( sqrt(B)/ (cos(incl*pi/180) * sig_R_0), 3)
      3 
      4 # sig_R_0 = 160
      5 # alpha = 0.55

NameError: name 'A' is not defined

Построим полученные профили дисперсий скоростей:


In [ ]:
def sigR_exp(R):
    return sig_R_0*poly_sig_min(R)/sig_min_0

def sigZ_exp(R):
    return alpha * sigR_exp(R)

def sigPhi_exp(R):
    return sigPhi_to_sigR(R) * sigR_exp(R)


plt.plot(points, [sigR_exp(R) for R in points], '-', color='red', label='$\sigma_{R, exp}$')
plt.plot(points, [sigPhi_exp(R) for R in points], '-', color='blue', label=r'$\sigma_{\varphi, exp}$')
plt.plot(points, [sigZ_exp(R) for R in points], '-', color='black', label='$\sigma_{Z, exp}$')
plt.legend()
plt.ylim(0, 350)
plt.xlim(0, 50)
plt.show()

И восстановим профили $\sigma_{los}^{maj}$ и $\sigma_{los}^{min}$. Связь профилей описывается следующими уравнениями: $$\sigma_{los,maj}^2=\sigma_{\varphi}^2\sin^2i+\sigma_Z^2\cos^2i$$ $$\sigma_{los,min}^2=\sigma_R^2\sin^2i+\sigma_Z^2\cos^2i$$


In [ ]:
def sig_maj_exp(R):
    return sqrt(sigPhi_exp(R)**2 * sin(incl*pi/180)**2 + sigZ_exp(R)**2 * cos(incl*pi/180)**2)

def sig_min_exp(R):
    return sqrt(sigR_exp(R)**2 * sin(incl*pi/180)**2 + sigZ_exp(R)**2 * cos(incl*pi/180)**2)

plt.plot(points, poly_sig_maj(points), '-', label = '$\sigma_{los}^{maj} polyfit$', color='blue')
plt.plot(points, [sig_maj_exp(R) for R in points], '--', color='blue', label='$\sigma_{maj, exp}$')
plt.plot(points, poly_sig_min(points), '-', label = '$\sigma_{los}^{min} polyfit$', color='red')
plt.plot(points, [sig_min_exp(R) for R in points], '--', color='red', label='$\sigma_{min, exp}$')
plt.legend()
plt.ylim(0, 350)
plt.xlim(0, 50)
plt.show()