In [1]:
%run ../../../utils/load_notebook.py

In [2]:
return_back_path = '../../notebooks/2f/test_short/'

In [3]:
import datetime
start = datetime.datetime.now()
start


Out[3]:
datetime.datetime(2017, 7, 4, 22, 27, 29, 536000)

In [4]:
%%time
from n338 import *


importing Jupyter notebook from n338.ipynb
Using matplotlib backend: Qt5Agg
Populating the interactive namespace from numpy and matplotlib
importing Jupyter notebook from photometry.ipynb
Using matplotlib backend: Qt5Agg
Populating the interactive namespace from numpy and matplotlib
importing Jupyter notebook from instabilities.ipynb
Using matplotlib backend: Qt5Agg
Populating the interactive namespace from numpy and matplotlib
ERROR:root:File `u'../../utils/load_notebook.py'` not found.
importing Jupyter notebook from utils.ipynb
Wall time: 1.57 s
C:\Anaconda\lib\site-packages\scipy\integrate\quadpack.py:364: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
  If increasing the limit yields no improvement it is advised to analyze 
  the integrand in order to determine the difficulties.  If the position of a 
  local difficulty can be determined (singularity, discontinuity) one will 
  probably gain from splitting up the interval and calling the integrator 
  on the subranges.  Perhaps a special-purpose integrator should be used.
  warnings.warn(msg, IntegrationWarning)
9.82081483752
11.1404465179
I : 1.15; B : 1.06; R: 1.24.
+------+-----------+---------+----------+------+---------+----------+-------+-------------+-----------+
|      | Name      |   r_eff |   mu_eff |    n |   mu0_d |   h_disc |   M/L | M_d/M_sun   |   Sigma_0 |
|------+-----------+---------+----------+------+---------+----------+-------+-------------+-----------|
| 0.00 | Noorder R |   15.00 |    20.95 | 3.70 |   21.92 |    18.30 |  1.24 | 9.43E+09.   |        53 |
| 1.00 | Noorder B |   15.00 |    22.70 | 3.70 |   22.53 |    17.70 |  1.06 | 1.14E+10.   |        68 |
| 2.00 | Noorder I |   15.00 |    20.41 | 3.70 |   19.79 |    12.90 |  1.15 | 2.27E+10.   |       255 |
+------+-----------+---------+----------+------+---------+----------+-------+-------------+-----------+
Noorder R      : M/L was 1.24 and for max it equal 20.91, for submax equal 10.42
Noorder B      : M/L was 1.06 and for max it equal 14.02, for submax equal 6.99
Noorder I      : M/L was 1.15 and for max it equal 6.23, for submax equal 3.10
C:\Anaconda\lib\site-packages\matplotlib\pyplot.py:524: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  max_open_warning, RuntimeWarning)
Wall time: 2min 13s
C:\Anaconda\lib\site-packages\matplotlib\axes\_axes.py:545: UserWarning: No labelled objects found. Use label='...' kwarg on individual plots.
  warnings.warn("No labelled objects found. "

In [5]:
name


Out[5]:
'N0338'

In [6]:
n338dict = dict(globals(), **locals())

In [7]:
n338dict['name']


Out[7]:
'N0338'

In [8]:
os.chdir(return_back_path)

In [9]:
%%time
from n1167 import *


importing Jupyter notebook from n1167.ipynb
Using matplotlib backend: Qt5Agg
Populating the interactive namespace from numpy and matplotlib
14.7100823358
B : 3.80; R: 2.71.
+------+----------------+---------+----------+------+---------+----------+-------+-------------+-----------+
|      | Name           |   r_eff |   mu_eff |    n |   mu0_d |   h_disc |   M/L | M_d/M_sun   |   Sigma_0 |
|------+----------------+---------+----------+------+---------+----------+-------+-------------+-----------|
| 0.00 | Noorder R      |    6.70 |    19.40 | 1.70 |   20.12 |    24.20 |  2.71 | 2.30E+11.   |       605 |
| 1.00 | Noorder B      |    6.70 |    20.73 | 1.70 |   21.71 |    27.50 |  3.80 | 2.55E+11.   |       520 |
| 2.00 | Noorder R_max  |    6.70 |    19.40 | 1.70 |   20.12 |    24.20 |  4.00 | 3.39E+11.   |       893 |
| 3.00 | califa g (g-i) |    7.57 |    21.47 | 2.21 |   21.29 |    24.45 |  5.59 | 3.89E+11.   |      1006 |
| 4.00 | califa r (g-r) |    9.09 |    20.70 | 2.66 |   20.60 |    25.82 |  3.88 | 3.04E+11.   |       703 |
| 5.00 | califa i (r-i) |   11.43 |    20.61 | 3.21 |   20.32 |    27.51 |  2.95 | 2.85E+11.   |       581 |
+------+----------------+---------+----------+------+---------+----------+-------+-------------+-----------+
Noorder R      : M/L was 2.71 and for max it equal 5.53, for submax equal 2.76
Noorder B      : M/L was 3.80 and for max it equal 7.93, for submax equal 3.95
Noorder R_max  : M/L was 4.00 and for max it equal 5.53, for submax equal 2.76
califa g (g-i) : M/L was 5.59 and for max it equal 6.80, for submax equal 3.39
califa r (g-r) : M/L was 3.88 and for max it equal 6.39, for submax equal 3.18
califa i (r-i) : M/L was 2.95 and for max it equal 5.52, for submax equal 2.75
Wall time: 2min 25s

In [10]:
name


Out[10]:
'N1167'

In [11]:
n1167dict = dict(globals(), **locals())

In [12]:
n1167dict['name']


Out[12]:
'N1167'

In [13]:
del n1167dict['n338dict']

In [14]:
os.chdir(return_back_path)

In [15]:
%%time
from n2985 import *


importing Jupyter notebook from n2985.ipynb
Using matplotlib backend: Qt5Agg
Populating the interactive namespace from numpy and matplotlib
Abs B : 0.47; R: 0.75.
Rel B : 0.53; R: 0.81.
0.58 0.63
+------+------------------+---------+----------+--------+----------------+----------------+-------+-------------+-----------+
|      | Name             |   r_eff |   mu_eff |      n | mu0_d          | h_disc         |   M/L | M_d/M_sun   |   Sigma_0 |
|------+------------------+---------+----------+--------+----------------+----------------+-------+-------------+-----------|
| 0.00 | b:Noorder R      |   25.10 |    20.41 |   3.90 | 21.32          | 52.2           |  0.81 | 8.47E+09.   |        60 |
| 1.00 | b:Noorder B      |   25.10 |    21.86 |   3.90 | 22.06          | 57.6           |  0.53 | 9.08E+09.   |        53 |
| 2.00 | b:Mendez-Abreu J |   13.20 |    17.94 |   2.92 | 18.22          | 25.8           |  0.76 | 1.66E+10.   |       479 |
| 3.00 | Gutierrez R 2d   |  nan    |   nan    | nan    | (18.57, 21.76) | (18.1, 81.0)   |  0.75 | 2.44E+10.   |       732 |
| 4.00 | Heidt J          |   12.08 |    17.73 |   2.86 | 18.05          | 30.63          |  0.76 | 2.73E+10.   |       560 |
| 5.00 | Heidt K          |   14.06 |    17.23 |   3.03 | 17.32          | 31.08          |  0.65 | 3.39E+10.   |       674 |
| 6.00 | Heidt H          |   12.81 |    17.03 |   2.86 | 17.27          | 26.1           |  0.70 | 2.76E+10.   |       779 |
| 7.00 | S4G 2d           |    6.32 |   nan    |   2.82 | (18.55, 20.84) | (12.78, 48.89) |  0.67 | 3.42E+10.   |      1630 |
| 8.00 | S4G_AM 2d        |  nan    |   nan    | nan    | (18.28, 20.74) | (11.71, 47.23) |  0.67 | 3.58E+10.   |      2067 |
| 9.00 | Noorder R_max    |   25.10 |    20.41 |   3.90 | 21.32          | 52.2           |  6.00 | 6.29E+10.   |       444 |
+------+------------------+---------+----------+--------+----------------+----------------+-------+-------------+-----------+
b:Noorder R    : M/L was 0.81 and for max it equal 12.08, for submax equal 6.02
b:Noorder B    : M/L was 0.53 and for max it equal 8.18, for submax equal 4.08
b:Mendez-Abreu J: M/L was 0.76 and for max it equal 2.80, for submax equal 1.40
Gutierrez R 2d : M/L was 0.75 and for max it equal 2.48, for submax equal 1.24
Heidt J        : M/L was 0.76 and for max it equal 2.03, for submax equal 1.01
Heidt K        : M/L was 0.65 and for max it equal 1.42, for submax equal 0.71
Heidt H        : M/L was 0.70 and for max it equal 1.55, for submax equal 0.77
S4G 2d         : M/L was 0.67 and for max it equal 1.35, for submax equal 0.67
S4G_AM 2d      : M/L was 0.67 and for max it equal 1.18, for submax equal 0.59
Noorder R_max  : M/L was 6.00 and for max it equal 12.08, for submax equal 6.02
Wall time: 12min 57s

In [16]:
n2985dict = dict(globals(), **locals())

In [17]:
n2985dict['name']


Out[17]:
'N2985'

In [18]:
del n2985dict['n338dict']
del n2985dict['n1167dict']

In [19]:
os.chdir(return_back_path)

In [20]:
%%time
from n3898 import *


importing Jupyter notebook from n3898.ipynb
Using matplotlib backend: Qt5Agg
Populating the interactive namespace from numpy and matplotlib
Wall time: 179 ms
3.71203889959
0.189557760061 0.0590237999569
3.51965422525
+------+-------------------+---------+----------+--------+----------------+---------------+-------+-------------+-----------+
|      | Name              |   r_eff |   mu_eff |      n | mu0_d          | h_disc        |   M/L | M_d/M_sun   |   Sigma_0 |
|------+-------------------+---------+----------+--------+----------------+---------------+-------+-------------+-----------|
| 0.00 | Noorder R         |    8.80 |    18.37 |   2.30 | 20.49          | 36.2          |  1.95 | 2.16E+10.   |       310 |
| 1.00 | Noorder B         |    8.80 |    19.80 |   2.30 | 22.0           | 42.9          |  2.22 | 2.28E+10.   |       233 |
| 2.00 | Mendez-Abreu J    |   11.90 |    18.13 |   3.75 | 19.07          | 29.2          |  1.16 | 1.51E+10.   |       332 |
| 3.00 | Gutierrez R (new) |  nan    |   nan    | nan    | (19.03, 21.53) | (19.11, 59.9) |  1.95 | 4.58E+10.   |      1310 |
| 4.00 | Pignatelli V      |   18.90 |    20.60 |   4.00 | 20.4           | 29.0          |  3.52 | 3.96E+10.   |       886 |
+------+-------------------+---------+----------+--------+----------------+---------------+-------+-------------+-----------+
Noorder R      : M/L was 1.95 and for max it equal 7.80, for submax equal 3.89
Noorder B      : M/L was 2.22 and for max it equal 10.20, for submax equal 5.08
Mendez-Abreu J : M/L was 1.16 and for max it equal 4.96, for submax equal 2.47
Gutierrez R (new): M/L was 1.95 and for max it equal 2.93, for submax equal 1.46
Pignatelli V   : M/L was 3.52 and for max it equal 5.68, for submax equal 2.83
Wall time: 1min 3s

In [21]:
n3898dict = dict(globals(), **locals())

In [22]:
n3898dict['name']


Out[22]:
'N3898'

In [23]:
del n3898dict['n338dict']
del n3898dict['n1167dict']
del n3898dict['n2985dict']

In [24]:
os.chdir(return_back_path)

In [25]:
%%time
from n4258 import *


importing Jupyter notebook from n4258.ipynb
Using matplotlib backend: Qt5Agg
Populating the interactive namespace from numpy and matplotlib
19.9002390653
+------+--------------------+---------+----------+--------+---------+----------+-------+-------------+-----------+
|      | Name               |   r_eff |   mu_eff |      n |   mu0_d |   h_disc |   M/L | M_d/M_sun   |   Sigma_0 |
|------+--------------------+---------+----------+--------+---------+----------+-------+-------------+-----------|
| 0.00 | YOSHINO I          |   14.89 |    18.48 |   1.50 |   18.26 |    74.24 |  0.97 | 5.62E+10.   |       878 |
| 1.00 | YOSHINO V          |   10.34 |    18.57 |   1.50 |   19.54 |    78.55 |  1.74 | 6.93E+10.   |       967 |
| 2.00 | YOSHINO J          |    7.31 |    16.43 |   1.50 |   16.37 |    38.95 |  0.95 | 5.76E+10.   |      3271 |
| 3.00 | SPITZER 3.6        |   15.03 |    17.41 |   2.80 |   18.82 |    80.74 |  0.65 | 8.45E+10.   |      1116 |
| 4.00 | SPITZER 3.6 faceon |   15.03 |    17.41 |   2.80 |   19.90 |    80.74 |  0.65 | 3.12E+10.   |       413 |
| 5.00 | GALFIT K           |    6.27 |   nan    |   3.26 |   17.80 |   146.00 |  0.79 | 1.30E+11.   |       523 |
| 6.00 | S4G 3.6            |  nan    |   nan    | nan    |   20.50 |   178.76 |  0.65 | 8.80E+10.   |       237 |
+------+--------------------+---------+----------+--------+---------+----------+-------+-------------+-----------+
r = 9.88; gas_d = 170.26; epicycl = 828.67; sig = 119.96; star_d = 768.96
	Qs = 8.91; Qg = 2.15; Qeff = 2.10
r = 28.64; gas_d = 90.39; epicycl = 272.92; sig = 109.81; star_d = 597.24
	Qs = 3.46; Qg = 1.33; Qeff = 1.28
r = 45.43; gas_d = 35.24; epicycl = 212.11; sig = 93.27; star_d = 476.37
	Qs = 2.86; Qg = 2.66; Qeff = 2.36
r = 63.20; gas_d = 16.97; epicycl = 91.21; sig = 81.75; star_d = 374.95
	Qs = 1.37; Qg = 2.38; Qeff = 1.19
r = 80.98; gas_d = 14.67; epicycl = 82.19; sig = 81.75; star_d = 295.11
	Qs = 1.57; Qg = 2.48; Qeff = 1.35
r = 98.75; gas_d = 10.43; epicycl = 75.31; sig = 81.75; star_d = 232.28
	Qs = 1.83; Qg = 3.19; Qeff = 1.58
r = 116.53; gas_d = 10.96; epicycl = 63.68; sig = 81.75; star_d = 182.82
	Qs = 1.96; Qg = 2.57; Qeff = 1.66
YOSHINO I                : M/L was 0.97 and for max it equal 0.82, for submax equal 0.41
YOSHINO V                : M/L was 1.74 and for max it equal 1.22, for submax equal 0.61
YOSHINO J                : M/L was 0.95 and for max it equal 0.53, for submax equal 0.27
SPITZER 3.6              : M/L was 0.65 and for max it equal 0.38, for submax equal 0.19
SPITZER 3.6 faceon       : M/L was 0.65 and for max it equal 1.03, for submax equal 0.51
GALFIT K                 : M/L was 0.79 and for max it equal 0.51, for submax equal 0.26
S4G 3.6                  : M/L was 0.65 and for max it equal 0.80, for submax equal 0.40
Wall time: 27min 39s
Error in callback <function post_execute at 0x0000000009A21C18> (for post_execute):
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
C:\Anaconda\lib\site-packages\matplotlib\pyplot.pyc in post_execute()
    147             def post_execute():
    148                 if matplotlib.is_interactive():
--> 149                     draw_all()
    150 
    151             # IPython >= 2

C:\Anaconda\lib\site-packages\matplotlib\_pylab_helpers.pyc in draw_all(cls, force)
    148         for f_mgr in cls.get_all_fig_managers():
    149             if force or f_mgr.canvas.figure.stale:
--> 150                 f_mgr.canvas.draw_idle()
    151 
    152 atexit.register(Gcf.destroy_all)

C:\Anaconda\lib\site-packages\matplotlib\backend_bases.pyc in draw_idle(self, *args, **kwargs)
   2038         if not self._is_idle_drawing:
   2039             with self._idle_draw_cntx():
-> 2040                 self.draw(*args, **kwargs)
   2041 
   2042     def draw_cursor(self, event):

C:\Anaconda\lib\site-packages\matplotlib\backends\backend_agg.pyc in draw(self)
    462 
    463         try:
--> 464             self.figure.draw(self.renderer)
    465         finally:
    466             RendererAgg.lock.release()

C:\Anaconda\lib\site-packages\matplotlib\artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs)
     61     def draw_wrapper(artist, renderer, *args, **kwargs):
     62         before(artist, renderer)
---> 63         draw(artist, renderer, *args, **kwargs)
     64         after(artist, renderer)
     65 

C:\Anaconda\lib\site-packages\matplotlib\figure.pyc in draw(self, renderer)
   1142 
   1143             mimage._draw_list_compositing_images(
-> 1144                 renderer, self, dsu, self.suppressComposite)
   1145 
   1146             renderer.close_group('figure')

C:\Anaconda\lib\site-packages\matplotlib\image.pyc in _draw_list_compositing_images(renderer, parent, dsu, suppress_composite)
    137     if not_composite or not has_images:
    138         for zorder, a in dsu:
--> 139             a.draw(renderer)
    140     else:
    141         # Composite any adjacent images together

C:\Anaconda\lib\site-packages\matplotlib\artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs)
     61     def draw_wrapper(artist, renderer, *args, **kwargs):
     62         before(artist, renderer)
---> 63         draw(artist, renderer, *args, **kwargs)
     64         after(artist, renderer)
     65 

C:\Anaconda\lib\site-packages\matplotlib\axes\_base.pyc in draw(self, renderer, inframe)
   2424             renderer.stop_rasterizing()
   2425 
-> 2426         mimage._draw_list_compositing_images(renderer, self, dsu)
   2427 
   2428         renderer.close_group('axes')

C:\Anaconda\lib\site-packages\matplotlib\image.pyc in _draw_list_compositing_images(renderer, parent, dsu, suppress_composite)
    137     if not_composite or not has_images:
    138         for zorder, a in dsu:
--> 139             a.draw(renderer)
    140     else:
    141         # Composite any adjacent images together

C:\Anaconda\lib\site-packages\matplotlib\artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs)
     61     def draw_wrapper(artist, renderer, *args, **kwargs):
     62         before(artist, renderer)
---> 63         draw(artist, renderer, *args, **kwargs)
     64         after(artist, renderer)
     65 

C:\Anaconda\lib\site-packages\matplotlib\axis.pyc in draw(self, renderer, *args, **kwargs)
   1148         self._update_label_position(ticklabelBoxes, ticklabelBoxes2)
   1149 
-> 1150         self.label.draw(renderer)
   1151 
   1152         self._update_offset_text_position(ticklabelBoxes, ticklabelBoxes2)

C:\Anaconda\lib\site-packages\matplotlib\artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs)
     61     def draw_wrapper(artist, renderer, *args, **kwargs):
     62         before(artist, renderer)
---> 63         draw(artist, renderer, *args, **kwargs)
     64         after(artist, renderer)
     65 

C:\Anaconda\lib\site-packages\matplotlib\text.pyc in draw(self, renderer)
    760             posy = float(textobj.convert_yunits(textobj._y))
    761             if not np.isfinite(posx) or not np.isfinite(posy):
--> 762                 raise ValueError("posx and posy should be finite values")
    763             posx, posy = trans.transform_point((posx, posy))
    764             canvasw, canvash = renderer.get_canvas_width_height()

ValueError: posx and posy should be finite values
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
C:\Anaconda\lib\site-packages\IPython\core\formatters.pyc in __call__(self, obj)
    305                 pass
    306             else:
--> 307                 return printer(obj)
    308             # Finally look for special method names
    309             method = get_real_method(obj, self.print_method)

C:\Anaconda\lib\site-packages\IPython\core\pylabtools.pyc in <lambda>(fig)
    225 
    226     if 'png' in formats:
--> 227         png_formatter.for_type(Figure, lambda fig: print_figure(fig, 'png', **kwargs))
    228     if 'retina' in formats or 'png2x' in formats:
    229         png_formatter.for_type(Figure, lambda fig: retina_figure(fig, **kwargs))

C:\Anaconda\lib\site-packages\IPython\core\pylabtools.pyc in print_figure(fig, fmt, bbox_inches, **kwargs)
    117 
    118     bytes_io = BytesIO()
--> 119     fig.canvas.print_figure(bytes_io, **kw)
    120     data = bytes_io.getvalue()
    121     if fmt == 'svg':

C:\Anaconda\lib\site-packages\matplotlib\backend_bases.pyc in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, **kwargs)
   2198                     orientation=orientation,
   2199                     dryrun=True,
-> 2200                     **kwargs)
   2201                 renderer = self.figure._cachedRenderer
   2202                 bbox_inches = self.figure.get_tightbbox(renderer)

C:\Anaconda\lib\site-packages\matplotlib\backends\backend_agg.pyc in print_png(self, filename_or_obj, *args, **kwargs)
    543 
    544     def print_png(self, filename_or_obj, *args, **kwargs):
--> 545         FigureCanvasAgg.draw(self)
    546         renderer = self.get_renderer()
    547         original_dpi = renderer.dpi

C:\Anaconda\lib\site-packages\matplotlib\backends\backend_agg.pyc in draw(self)
    462 
    463         try:
--> 464             self.figure.draw(self.renderer)
    465         finally:
    466             RendererAgg.lock.release()

C:\Anaconda\lib\site-packages\matplotlib\artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs)
     61     def draw_wrapper(artist, renderer, *args, **kwargs):
     62         before(artist, renderer)
---> 63         draw(artist, renderer, *args, **kwargs)
     64         after(artist, renderer)
     65 

C:\Anaconda\lib\site-packages\matplotlib\figure.pyc in draw(self, renderer)
   1142 
   1143             mimage._draw_list_compositing_images(
-> 1144                 renderer, self, dsu, self.suppressComposite)
   1145 
   1146             renderer.close_group('figure')

C:\Anaconda\lib\site-packages\matplotlib\image.pyc in _draw_list_compositing_images(renderer, parent, dsu, suppress_composite)
    137     if not_composite or not has_images:
    138         for zorder, a in dsu:
--> 139             a.draw(renderer)
    140     else:
    141         # Composite any adjacent images together

C:\Anaconda\lib\site-packages\matplotlib\artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs)
     61     def draw_wrapper(artist, renderer, *args, **kwargs):
     62         before(artist, renderer)
---> 63         draw(artist, renderer, *args, **kwargs)
     64         after(artist, renderer)
     65 

C:\Anaconda\lib\site-packages\matplotlib\axes\_base.pyc in draw(self, renderer, inframe)
   2424             renderer.stop_rasterizing()
   2425 
-> 2426         mimage._draw_list_compositing_images(renderer, self, dsu)
   2427 
   2428         renderer.close_group('axes')

C:\Anaconda\lib\site-packages\matplotlib\image.pyc in _draw_list_compositing_images(renderer, parent, dsu, suppress_composite)
    137     if not_composite or not has_images:
    138         for zorder, a in dsu:
--> 139             a.draw(renderer)
    140     else:
    141         # Composite any adjacent images together

C:\Anaconda\lib\site-packages\matplotlib\artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs)
     61     def draw_wrapper(artist, renderer, *args, **kwargs):
     62         before(artist, renderer)
---> 63         draw(artist, renderer, *args, **kwargs)
     64         after(artist, renderer)
     65 

C:\Anaconda\lib\site-packages\matplotlib\axis.pyc in draw(self, renderer, *args, **kwargs)
   1148         self._update_label_position(ticklabelBoxes, ticklabelBoxes2)
   1149 
-> 1150         self.label.draw(renderer)
   1151 
   1152         self._update_offset_text_position(ticklabelBoxes, ticklabelBoxes2)

C:\Anaconda\lib\site-packages\matplotlib\artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs)
     61     def draw_wrapper(artist, renderer, *args, **kwargs):
     62         before(artist, renderer)
---> 63         draw(artist, renderer, *args, **kwargs)
     64         after(artist, renderer)
     65 

C:\Anaconda\lib\site-packages\matplotlib\text.pyc in draw(self, renderer)
    760             posy = float(textobj.convert_yunits(textobj._y))
    761             if not np.isfinite(posx) or not np.isfinite(posy):
--> 762                 raise ValueError("posx and posy should be finite values")
    763             posx, posy = trans.transform_point((posx, posy))
    764             canvasw, canvash = renderer.get_canvas_width_height()

ValueError: posx and posy should be finite values
<matplotlib.figure.Figure at 0x1b57afd0>

In [26]:
n4258dict = dict(globals(), **locals())

In [27]:
n4258dict['name']


Out[27]:
'N4258'

In [28]:
del n4258dict['n338dict']
del n4258dict['n1167dict']
del n4258dict['n2985dict']
del n4258dict['n3898dict']

In [29]:
os.chdir(return_back_path)

In [30]:
%%time
from n4725 import *


importing Jupyter notebook from n4725.ipynb
Using matplotlib backend: Qt5Agg
Populating the interactive namespace from numpy and matplotlib
+------+-------------------+---------+----------+--------+---------+----------+-------+-------------+-----------+
|      | Name              |   r_eff |   mu_eff |      n |   mu0_d |   h_disc |   M/L | M_d/M_sun   |   Sigma_0 |
|------+-------------------+---------+----------+--------+---------+----------+-------+-------------+-----------|
| 0.00 | S4G 3.6           |   10.14 |   nan    |   2.21 |   20.34 |    73.20 |  0.68 | 9.26E+10.   |       286 |
| 1.00 | Heidt J           |  nan    |   nan    | nan    |   17.78 |    49.99 |  0.75 | 1.06E+11.   |       703 |
| 2.00 | Heidt H           |  nan    |   nan    | nan    |   17.11 |    50.28 |  0.69 | 1.36E+11.   |       891 |
| 3.00 | Heidt K           |  nan    |   nan    | nan    |   17.01 |    54.66 |  0.65 | 1.60E+11.   |       889 |
| 4.00 | infra 3.6         |   15.72 |    17.49 |   3.61 |   19.65 |    71.86 |  0.68 | 1.68E+11.   |       539 |
| 5.00 | infra 3.6 face-on |   15.72 |    17.49 |   3.61 |   20.03 |    71.86 |  0.68 | 1.18E+11.   |       379 |
+------+-------------------+---------+----------+--------+---------+----------+-------+-------------+-----------+
316.386554622
41.1302521008
17.8888888889
2.16129032258
2.2905027933
S4G 3.6        : M/L was 0.68 and for max it equal 1.61, for submax equal 0.80
Heidt J        : M/L was 0.75 and for max it equal 0.94, for submax equal 0.47
Heidt H        : M/L was 0.69 and for max it equal 0.68, for submax equal 0.34
Heidt K        : M/L was 0.65 and for max it equal 0.61, for submax equal 0.31
infra 3.6      : M/L was 0.68 and for max it equal 0.87, for submax equal 0.43
infra 3.6 face-on: M/L was 0.68 and for max it equal 1.24, for submax equal 0.62
Wall time: 12min 2s

In [31]:
n4725dict = dict(globals(), **locals())

In [32]:
n4725dict['name']


Out[32]:
'N4725'

In [33]:
del n4725dict['n338dict']
del n4725dict['n1167dict']
del n4725dict['n2985dict']
del n4725dict['n3898dict']
del n4725dict['n4258dict']

In [34]:
os.chdir(return_back_path)

In [35]:
%%time
from n5533 import *


importing Jupyter notebook from n5533.ipynb
Using matplotlib backend: Qt5Agg
Populating the interactive namespace from numpy and matplotlib
Abs B : 2.34; R: 2.01.
Rel B : 2.52; R: 2.11.
1.21 1.24
+------+----------------+---------+----------+--------+---------+----------+-------+-------------+-----------+
|      | Name           |   r_eff |   mu_eff |      n |   mu0_d |   h_disc |   M/L | M_d/M_sun   |   Sigma_0 |
|------+----------------+---------+----------+--------+---------+----------+-------+-------------+-----------|
| 0.00 | Noorder R      |    9.90 |    19.75 |   2.70 |   21.27 |    34.40 |  2.11 | 8.52E+10.   |       163 |
| 1.00 | Noorder B      |    9.90 |    21.45 |   2.70 |   22.39 |    32.40 |  2.52 | 8.56E+10.   |       185 |
| 2.00 | Noorder R_max  |    9.90 |    19.75 |   2.70 |   21.27 |    34.40 |  5.00 | 2.02E+11.   |       387 |
| 3.00 | Mendez-Abreu J |    5.10 |    17.70 |   2.34 |   18.03 |    14.10 |  1.18 | 7.70E+10.   |       878 |
| 4.00 | 1991 R         |  nan    |    17.64 | nan    |   20.73 |    43.02 |  2.01 | 2.09E+11.   |       256 |
| 5.00 | 1998 V         |   19.60 |    21.62 |   4.00 |   22.43 |    40.70 |  2.26 | 6.41E+10.   |        88 |
| 6.00 | Silchenko V    |  nan    |   nan    | nan    |   20.20 |    18.00 |  2.26 | 9.78E+10.   |       684 |
| 7.00 | califa g (r-i) |    8.80 |    20.87 |   3.33 |   21.42 |    28.17 |  3.82 | 2.15E+11.   |       613 |
| 8.00 | califa r (g-i) |    8.88 |    20.01 |   3.17 |   20.72 |    28.01 |  2.91 | 1.63E+11.   |       472 |
| 9.00 | califa i (g-r) |    9.72 |    19.75 |   3.21 |   20.37 |    28.01 |  2.35 | 1.52E+11.   |       441 |
+------+----------------+---------+----------+--------+---------+----------+-------+-------------+-----------+
96.426
89.1
87.12
Noorder R      : M/L was 2.11 and for max it equal 7.15, for submax equal 3.56
Noorder B      : M/L was 2.52 and for max it equal 8.15, for submax equal 4.06
Noorder R_max  : M/L was 5.00 and for max it equal 7.15, for submax equal 3.56
Mendez-Abreu J : M/L was 1.18 and for max it equal 2.06, for submax equal 1.03
1991 R         : M/L was 2.01 and for max it equal 3.27, for submax equal 1.63
1998 V         : M/L was 2.26 and for max it equal 11.53, for submax equal 5.75
Silchenko V    : M/L was 2.26 and for max it equal 3.97, for submax equal 1.98
califa g (r-i) : M/L was 3.82 and for max it equal 4.42, for submax equal 2.20
califa r (g-i) : M/L was 2.91 and for max it equal 4.41, for submax equal 2.20
califa i (g-r) : M/L was 2.35 and for max it equal 3.81, for submax equal 1.90
[ 38.34400115]
[ 24.67098717]
Wall time: 1min 49s

In [36]:
n5533dict = dict(globals(), **locals())

In [37]:
n5533dict['name']


Out[37]:
'N5533'

In [38]:
del n5533dict['n338dict']
del n5533dict['n1167dict']
del n5533dict['n2985dict']
del n5533dict['n3898dict']
del n5533dict['n4258dict']
del n5533dict['n4725dict']

In [39]:
end = datetime.datetime.now()
end


Out[39]:
datetime.datetime(2017, 7, 4, 23, 31, 4, 41000)

In [40]:
finish = end-start

In [41]:
finish.total_seconds()/3600.


Out[41]:
1.0595847222222223

Конец вычислений


In [ ]:


In [42]:
# %%time
# for ind, name1 in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
#     for name2 in ['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533'][ind+1:]:
#         dict1 = locals()[name1+'dict']
#         dict2 = locals()[name2+'dict']
#         for key in dict1.keys():
#             if key in dict2.keys():
#                 try:
#                     if type(dict1[key]) == type(dict2[key]) == np.ndarray:
#                         if (dict1[key] == dict2[key]).all():
#                             del dict2[key]
#                     elif dict1[key] == dict2[key]:
#                         del dict2[key]
#                 except Exception:
#                     print key, name1, name2

In [ ]:


In [43]:
n338dict['star_approx'] == n2985dict['star_approx']


Out[43]:
False

In [ ]:


In [ ]:


In [ ]:


In [44]:
tex_imgs_dir = 'C:\\Users\\root\\Dropbox\\RotationCurves\\PhD\\paper2\\imgs'

In [45]:
os.chdir(tex_imgs_dir)

fig, axs = plt.subplots(nrows=7, ncols=3, sharex=False, sharey=False, figsize=[16,30])

names  = [r'$\rm{NGC\, 338}$', r'$\rm{NGC\, 1167}$', r'$\rm{NGC\, 2985}$', r'$\rm{NGC\, 3898}$', r'$\rm{NGC\, 4258}$', r'$\rm{NGC\, 4725}$', r'$\rm{NGC\, 5533}$',]

ax1 = axs[0,0]
ax2 = axs[1,0]
ax3 = axs[2,0]
ax4 = axs[3,0]
ax5 = axs[4,0]
ax6 = axs[5,0]
ax7 = axs[6,0]

axes = [ax1, ax2, ax3, ax4, ax5, ax6, ax7]

# ===================================================
# Дисперии
# ===================================================

for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
    ax = axes[ind]
    try:
        ax.errorbar(map(abs, locals()[name+'dict']['r_sig_ma']), locals()[name+'dict']['sig_ma'], yerr=locals()[name+'dict']['e_sig_ma'], 
                    fmt='.', marker='.', mew=0, color='red', label='$\sigma_{los}^{maj}$')
        ax.plot(points, map(locals()[name+'dict']['sig_R_maj_min'], locals()[name+'dict']['points']))
        ax.plot(points, map(locals()[name+'dict']['sig_R_maj_max'], locals()[name+'dict']['points']))
#         ax.plot(points, map(locals()[name+'dict']['sig_R_maj_maxmaxtrue'], locals()[name+'dict']['points']))
    except Exception:
        print 'WARNING S:{}'.format(name)

    ax.xaxis.set_major_locator(MultipleLocator(10))
    ax.xaxis.set_minor_locator(MultipleLocator(2))
#     ax.plot([reb, reb],[0, (50 - 30*(ind/2))], color='black', lw=2)
#     ax.axvline(x=reb, ls='--', color='black')
    ax.text(0.8, 0.9, names[ind], fontsize=20, ha='center', va='center', transform=ax.transAxes)
    ax.yaxis.set_label_coords(-0.055, 0.5)
    ax.set_ylabel(r'$\sigma_{\rm{los}},\, \rm{km/s}$', fontsize=20)
    ax.set_xlabel(r'$R,\, \rm{arcsec}$', fontsize=20)
    ax.yaxis.set_major_locator(MultipleLocator(100))

    
# ===================================================
# Кривая вращения
# ===================================================

ax1 = axs[0,1]
ax2 = axs[1,1]
ax3 = axs[2,1]
ax4 = axs[3,1]
ax5 = axs[4,1]
ax6 = axs[5,1]
ax7 = axs[6,1]

axes = [ax1, ax2, ax3, ax4, ax5, ax6, ax7]

for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
    ax = axes[ind]
    try:
        ax.plot(points, map(locals()[name+'dict']['spl_gas'], locals()[name+'dict']['points']), '--')
        ax.plot(locals()[name+'dict']['r_g_b'], locals()[name+'dict']['vel_g_b'], '.')
    except Exception:
        print 'WARNING V:{}'.format(name)
        
    try:
        ax.plot(locals()[name+'dict']['r_wsrt'], locals()[name+'dict']['vel_wsrt'], '.')
    except Exception:
        print 'WARNING V:{}'.format(name)
        
    try:
        ax.plot(locals()[name+'dict']['r_noord'], locals()[name+'dict']['vel_noord'], '.')
        ax.plot(locals()[name+'dict']['r_ma_n'], locals()[name+'dict']['vel_ma_n'], 's')
    except Exception:
        print 'WARNING V:{}'.format(name)    
        
    
    ax.text(0.8, 0.2, names[ind], fontsize=20, ha='center', va='center', transform=ax.transAxes)
    ax.yaxis.set_label_coords(-0.055, 0.5)
    ax.set_ylabel(r'$V,\, \rm{km/s}$', fontsize=20)
    ax.set_xlabel(r'$R,\, \rm{arcsec}$', fontsize=20)
    ax.set_ylim(0, 400)
    ax.xaxis.set_major_locator(MultipleLocator(10))
#     ax.xaxis.set_minor_locator(MultipleLocator(2))
    ax.yaxis.set_major_locator(MultipleLocator(100))
    
# ax2.set_ylim(0, 400)


# ===================================================
# Поверхностная плотность
# ===================================================
    
ax1 = axs[0,2]
ax2 = axs[1,2]
ax3 = axs[2,2]
ax4 = axs[3,2]
ax5 = axs[4,2]
ax6 = axs[5,2]
ax7 = axs[6,2]

axes = [ax1, ax2, ax3, ax4, ax5, ax6, ax7]



#     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-photom)')
#     axes[3].set_title('Gas')
#     axes[3].grid()
#     axes[3].set_xlim(0, 200)
#     axes[3].legend()
    
#     axes[3].plot(r_g_dens, gas_dens, 'd-')
#     axes[3].plot(r_g_dens, [He_coeff*(y_interp_(l[0], h_disc_I) + l[1]) for l in zip(r_g_dens, gas_dens)], '*-')
#     axes[3].plot(r_g_dens, [y_interp_(l, h_disc_I) for l in r_g_dens], '--', label='H2 (I-photom)')
#     axes[3].set_title('Gas')
#     axes[3].grid()
#     axes[3].set_xlim(0, 200)
#     axes[3].legend()
    
#     axes[3].plot(r_HI_dens, HI_dens, '--', label='HI')
#     axes[3].plot(zip(*total_gas_data)[0], zip(*total_gas_data)[1], '*-')
#     axes[3].plot(r_mol_dens, mol_dens, '--', label='mol')
#     axes[3].set_title('Gas')
#     axes[3].grid()
#     axes[3].set_xlim(0, 200)
#     axes[3].legend()


for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
    ax = axes[ind]
    try:
        ax.plot(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'], 'o-')
    except Exception:
        print 'WARNING G:{}'.format(name)
        
    try:
        ax.plot(locals()[name+'dict']['r_HI_dens'], locals()[name+'dict']['HI_dens'], 'o-')
    except Exception:
        print 'WARNING G:{}'.format(name)
        
    try:
        ax.plot(zip(*locals()[name+'dict']['total_gas_data'])[0], zip(*locals()[name+'dict']['total_gas_data'])[1], 'o-')
    except Exception:
        print 'WARNING G:{}'.format(name)
        
    try:
        ax.plot(zip(*locals()[name+'dict']['total_gas_data_'])[0], zip(*locals()[name+'dict']['total_gas_data_'])[1], 'o-')
    except Exception:
        print 'WARNING G:{}'.format(name)

    ax.xaxis.set_major_locator(MultipleLocator(100))
    ax.xaxis.set_minor_locator(MultipleLocator(50))
#     ax.plot([reb, reb],[0, (50 - 30*(ind/2))], color='black', lw=2)
#     ax.axvline(x=reb, ls='--', color='black')
    ax.text(0.8, 0.9, names[ind], fontsize=20, ha='center', va='center', transform=ax.transAxes)
    ax.yaxis.set_label_coords(-0.055, 0.5)
    ax.set_ylabel(r'$\Sigma_{\rm{los}},\, \rm{M_{sun}/{pc}^2}$', fontsize=20)
    ax.set_xlabel(r'$R,\, \rm{arcsec}$', fontsize=20)
    ax.yaxis.set_major_locator(MultipleLocator(10))
    ax.yaxis.set_minor_locator(MultipleLocator(1))
    if name == 'n4258':
        ax.set_ylim(0, 30)

# ax1.set_xlim(0, 45)
# ax1.set_ylim(0)    

# ax2.set_ylim(0, 240)
# ax2.set_xlim(0, 45)
# # ax2.set_ylabel(r'$\sigma_{\rm{los}},\, \rm{km/s}$', fontsize=15)

# ax3.set_ylim(0, 225)
# ax3.set_xlim(0, 60)
# # ax3.axes.yaxis.set_ticklabels([])

# # ax4.set_ylim(0, 150)
# # ax4.set_xlim(0, 72)
# # ax4.set_ylabel(r'$\sigma_{\rm{los}},\, \rm{km/s}$', fontsize=10)
# # ax4.set_xlabel(r'$R,\, \rm{arcsec}$', fontsize=15)

# ax4.set_ylim(0, 110)
# # ax5.set_xlim(0, 50)
# # ax5.axes.yaxis.set_ticklabels([])
# # ax5.set_ylabel(r'$\sigma_{\rm{gas}},\, \rm{km/s}$', fontsize=20)
# # ax5.set_xlabel(r'$R,\, \rm{arcsec}$', fontsize=20)


fig.subplots_adjust(wspace=0.15, hspace=0.25)
# plt.savefig('observ_data.eps', format='eps')
# plt.savefig('observ_data.png', format='png')
# plt.savefig('observ_data.pdf', format='pdf', dpi=150)
plt.show()


WARNING V:n338
WARNING G:n338
WARNING G:n1167
WARNING G:n2985
WARNING G:n3898

In [ ]:

Kennicutt-like plot


In [46]:
def plot_kennicutt(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, sfrange=None):

    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))

    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))
   
    rr = zip(*total_gas_data)[0]
    
#     ax.fill_between(rr, invQeff_min, invQeff_max, color=color, alpha=alpha, label=label)
    ax.plot(np.array(rr)/sfrange, invQeff_min, '.-', color=color, alpha=0.8)
    ax.plot(np.array(rr)/sfrange, invQeff_max, '--', color=color, alpha=0.8)
#     ax.plot(rr, invQg, 'v-', color='b')

#     ax.set_ylim(0., 1.5)
#     ax.set_xlim(0., data_lim+50.)
#     plot_SF(ax)
#     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()

In [47]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
            'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
            'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
            'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
            'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
            'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
            'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}

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

for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
    ax = plt.gca()
    
    color = cm.rainbow(np.linspace(0, 1, 7))[ind]
    
    phot = kenn_photom[name]
    mud = locals()[name+'dict'][phot[0]]
    h = locals()[name+'dict'][phot[1]]
#     print name
#     print mud
#     print h
    
    if name in ['n4258', 'n4725']:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
    elif name == 'n2985':
        total_gas_data=locals()[name+'dict']['total_gas_data_']
    else:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in 
                                                                   zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
            
            
    if total_gas_data[0][0] < 0.1:
        total_gas_data = total_gas_data[1:]
    
#     print total_gas_data

    plot_kennicutt(ax=ax, total_gas_data=total_gas_data,
                      epicycl=locals()[name+'dict']['epicyclicFreq_real'], 
                      gas_approx=locals()[name+'dict']['spl_gas'], 
                      sound_vel=locals()[name+'dict']['sound_vel'], 
                      scale=locals()[name+'dict']['scale'], 
                      sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
                      sigma_min=locals()[name+'dict']['sig_R_maj_min'],
                      star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      data_lim=locals()[name+'dict']['data_lim'], 
                  color=color, alpha=0.2, 
                  label=name+' '+phot[3],
                  sfrange=phot[4])
    plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
    
plt.axvline(x=1, alpha=0.5)
plot_Q_levels(ax, [1.5, 2., 3.])
plt.legend()
plt.xlim(0, 3)
plt.ylim(0, 1.3)
plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
plt.xlabel(r'$R/R_{SF}$', fontsize=20)
plt.show()


Среднее значение


In [54]:
def plot_kennicutt_half(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, sfrange=None):

    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))

    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))
   
    rr = zip(*total_gas_data)[0]
    ax.plot(np.array(rr)/sfrange, (np.array(invQeff_min)+np.array(invQeff_max))/2., '.-', color=color, alpha=0.8)

In [55]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
            'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
            'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
            'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
            'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
            'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
            'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}

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

for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
    ax = plt.gca()
    
    color = cm.rainbow(np.linspace(0, 1, 7))[ind]
    
    phot = kenn_photom[name]
    mud = locals()[name+'dict'][phot[0]]
    h = locals()[name+'dict'][phot[1]]
#     print name
#     print mud
#     print h
    
    if name in ['n4258', 'n4725']:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
    elif name == 'n2985':
        total_gas_data=locals()[name+'dict']['total_gas_data_']
    else:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in 
                                                                   zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
            
            
    if total_gas_data[0][0] < 0.1:
        total_gas_data = total_gas_data[1:]
    
#     print total_gas_data

    plot_kennicutt_half(ax=ax, total_gas_data=total_gas_data,
                      epicycl=locals()[name+'dict']['epicyclicFreq_real'], 
                      gas_approx=locals()[name+'dict']['spl_gas'], 
                      sound_vel=locals()[name+'dict']['sound_vel'], 
                      scale=locals()[name+'dict']['scale'], 
                      sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
                      sigma_min=locals()[name+'dict']['sig_R_maj_min'],
                      star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      data_lim=locals()[name+'dict']['data_lim'], 
                  color=color, alpha=0.2, 
                  label=name+' '+phot[3],
                  sfrange=phot[4])
    plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
    
plt.axvline(x=1, alpha=0.5)
plot_Q_levels(ax, [1.5, 2., 3.])
plt.legend()
plt.xlim(0, 3)
plt.ylim(0, 1.3)
plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
plt.xlabel(r'$R/R_{SF}$', fontsize=20)
plt.show()


Только для галактик с молек. газом:


In [58]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
            'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
            'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
            'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
            'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
            'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
            'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}

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

for ind, name in enumerate(['n2985', 'n4258', 'n4725']):
    ax = plt.gca()
    
    color = cm.rainbow(np.linspace(0, 1, 7))[ind]
    
    phot = kenn_photom[name]
    mud = locals()[name+'dict'][phot[0]]
    h = locals()[name+'dict'][phot[1]]
#     print name
#     print mud
#     print h
    
    if name in ['n4258', 'n4725']:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
    elif name == 'n2985':
        total_gas_data=locals()[name+'dict']['total_gas_data_']
    else:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in 
                                                                   zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
            
            
    if total_gas_data[0][0] < 0.1:
        total_gas_data = total_gas_data[1:]
    
#     print total_gas_data

    plot_kennicutt_half(ax=ax, total_gas_data=total_gas_data,
                      epicycl=locals()[name+'dict']['epicyclicFreq_real'], 
                      gas_approx=locals()[name+'dict']['spl_gas'], 
                      sound_vel=locals()[name+'dict']['sound_vel'], 
                      scale=locals()[name+'dict']['scale'], 
                      sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
                      sigma_min=locals()[name+'dict']['sig_R_maj_min'],
                      star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      data_lim=locals()[name+'dict']['data_lim'], 
                  color=color, alpha=0.2, 
                  label=name+' '+phot[3],
                  sfrange=phot[4])
    plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
    
plt.axvline(x=1, alpha=0.5)
plot_Q_levels(ax, [1.5, 2., 3.])
plt.legend()
plt.xlim(0, 3)
plt.ylim(0, 1.3)
plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
plt.xlabel(r'$R/R_{SF}$', fontsize=20)
plt.show()


ВАУ.

Если изменить масштаб для 4725 (там непонятно), то тоже хорошо, но хуже.


In [59]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
            'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
            'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
            'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
            'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
            'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 100.],
            'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}

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

for ind, name in enumerate(['n2985', 'n4258', 'n4725']):
    ax = plt.gca()
    
    color = cm.rainbow(np.linspace(0, 1, 7))[ind]
    
    phot = kenn_photom[name]
    mud = locals()[name+'dict'][phot[0]]
    h = locals()[name+'dict'][phot[1]]
#     print name
#     print mud
#     print h
    
    if name in ['n4258', 'n4725']:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
    elif name == 'n2985':
        total_gas_data=locals()[name+'dict']['total_gas_data_']
    else:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in 
                                                                   zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
            
            
    if total_gas_data[0][0] < 0.1:
        total_gas_data = total_gas_data[1:]
    
#     print total_gas_data

    plot_kennicutt_half(ax=ax, total_gas_data=total_gas_data,
                      epicycl=locals()[name+'dict']['epicyclicFreq_real'], 
                      gas_approx=locals()[name+'dict']['spl_gas'], 
                      sound_vel=locals()[name+'dict']['sound_vel'], 
                      scale=locals()[name+'dict']['scale'], 
                      sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
                      sigma_min=locals()[name+'dict']['sig_R_maj_min'],
                      star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      data_lim=locals()[name+'dict']['data_lim'], 
                  color=color, alpha=0.2, 
                  label=name+' '+phot[3],
                  sfrange=phot[4])
    plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
    
plt.axvline(x=1, alpha=0.5)
plot_Q_levels(ax, [1.5, 2., 3.])
plt.legend()
plt.xlim(0, 3)
plt.ylim(0, 1.3)
plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
plt.xlabel(r'$R/R_{SF}$', fontsize=20)
plt.show()


Разброс для молек.


In [60]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
            'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
            'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
            'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
            'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
            'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
            'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}

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

for ind, name in enumerate(['n2985', 'n4258', 'n4725']):
    ax = plt.gca()
    
    color = cm.rainbow(np.linspace(0, 1, 7))[ind]
    
    phot = kenn_photom[name]
    mud = locals()[name+'dict'][phot[0]]
    h = locals()[name+'dict'][phot[1]]
#     print name
#     print mud
#     print h
    
    if name in ['n4258', 'n4725']:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
    elif name == 'n2985':
        total_gas_data=locals()[name+'dict']['total_gas_data_']
    else:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in 
                                                                   zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
            
            
    if total_gas_data[0][0] < 0.1:
        total_gas_data = total_gas_data[1:]
    
#     print total_gas_data

    plot_kennicutt(ax=ax, total_gas_data=total_gas_data,
                      epicycl=locals()[name+'dict']['epicyclicFreq_real'], 
                      gas_approx=locals()[name+'dict']['spl_gas'], 
                      sound_vel=locals()[name+'dict']['sound_vel'], 
                      scale=locals()[name+'dict']['scale'], 
                      sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
                      sigma_min=locals()[name+'dict']['sig_R_maj_min'],
                      star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      data_lim=locals()[name+'dict']['data_lim'], 
                  color=color, alpha=0.2, 
                  label=name+' '+phot[3],
                  sfrange=phot[4])
    plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
    
plt.axvline(x=1, alpha=0.5)
plot_Q_levels(ax, [1.5, 2., 3.])
plt.legend()
plt.xlim(0, 3)
plt.ylim(0, 1.3)
plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
plt.xlabel(r'$R/R_{SF}$', fontsize=20)
plt.show()


Одножидкостный


In [64]:
def plot_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, sfrange=None):

    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))
   
    rr = zip(*total_gas_data)[0]
    ax.plot(np.array(rr)/sfrange, invQg, '--', color=color, alpha=0.8)

In [65]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
            'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
            'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
            'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
            'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
            'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
            'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}

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

for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
    ax = plt.gca()
    
    color = cm.rainbow(np.linspace(0, 1, 7))[ind]
    
    phot = kenn_photom[name]
    mud = locals()[name+'dict'][phot[0]]
    h = locals()[name+'dict'][phot[1]]
#     print name
#     print mud
#     print h
    
    if name in ['n4258', 'n4725']:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
    elif name == 'n2985':
        total_gas_data=locals()[name+'dict']['total_gas_data_']
    else:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in 
                                                                   zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
            
            
    if total_gas_data[0][0] < 0.1:
        total_gas_data = total_gas_data[1:]
    
#     print total_gas_data

    plot_1f(ax=ax, total_gas_data=total_gas_data,
                      epicycl=locals()[name+'dict']['epicyclicFreq_real'], 
                      gas_approx=locals()[name+'dict']['spl_gas'], 
                      sound_vel=locals()[name+'dict']['sound_vel'], 
                      scale=locals()[name+'dict']['scale'], 
                      sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
                      sigma_min=locals()[name+'dict']['sig_R_maj_min'],
                      star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      data_lim=locals()[name+'dict']['data_lim'], 
                  color=color, alpha=0.2, 
                  label=name+' '+phot[3],
                  sfrange=phot[4])
    plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
    
plt.axvline(x=1, alpha=0.5)
plot_Q_levels(ax, [1.5, 2., 3.])
plt.legend()
plt.xlim(0, 3)
plt.ylim(0, 1.3)
plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
plt.xlabel(r'$R/R_{SF}$', fontsize=20)
plt.show()


Только с молек.:


In [66]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
            'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
            'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
            'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
            'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
            'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
            'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}

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

for ind, name in enumerate(['n2985', 'n4258', 'n4725']):
    ax = plt.gca()
    
    color = cm.rainbow(np.linspace(0, 1, 7))[ind]
    
    phot = kenn_photom[name]
    mud = locals()[name+'dict'][phot[0]]
    h = locals()[name+'dict'][phot[1]]
#     print name
#     print mud
#     print h
    
    if name in ['n4258', 'n4725']:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
    elif name == 'n2985':
        total_gas_data=locals()[name+'dict']['total_gas_data_']
    else:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in 
                                                                   zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
            
            
    if total_gas_data[0][0] < 0.1:
        total_gas_data = total_gas_data[1:]
    
#     print total_gas_data

    plot_1f(ax=ax, total_gas_data=total_gas_data,
                      epicycl=locals()[name+'dict']['epicyclicFreq_real'], 
                      gas_approx=locals()[name+'dict']['spl_gas'], 
                      sound_vel=locals()[name+'dict']['sound_vel'], 
                      scale=locals()[name+'dict']['scale'], 
                      sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
                      sigma_min=locals()[name+'dict']['sig_R_maj_min'],
                      star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      data_lim=locals()[name+'dict']['data_lim'], 
                  color=color, alpha=0.2, 
                  label=name+' '+phot[3],
                  sfrange=phot[4])
    plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
    
plt.axvline(x=1, alpha=0.5)
plot_Q_levels(ax, [1.5, 2., 3.])
plt.legend()
plt.xlim(0, 3)
plt.ylim(0, 1.3)
plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
plt.xlabel(r'$R/R_{SF}$', fontsize=20)
plt.show()


Vs 2F:


In [67]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
            'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
            'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
            'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
            'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
            'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
            'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}

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

for ind, name in enumerate(['n2985', 'n4258', 'n4725']):
    ax = plt.gca()
    
    color = cm.rainbow(np.linspace(0, 1, 7))[ind]
    
    phot = kenn_photom[name]
    mud = locals()[name+'dict'][phot[0]]
    h = locals()[name+'dict'][phot[1]]
#     print name
#     print mud
#     print h
    
    if name in ['n4258', 'n4725']:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
    elif name == 'n2985':
        total_gas_data=locals()[name+'dict']['total_gas_data_']
    else:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in 
                                                                   zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
            
            
    if total_gas_data[0][0] < 0.1:
        total_gas_data = total_gas_data[1:]
    
#     print total_gas_data

    plot_1f(ax=ax, total_gas_data=total_gas_data,
                      epicycl=locals()[name+'dict']['epicyclicFreq_real'], 
                      gas_approx=locals()[name+'dict']['spl_gas'], 
                      sound_vel=locals()[name+'dict']['sound_vel'], 
                      scale=locals()[name+'dict']['scale'], 
                      sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
                      sigma_min=locals()[name+'dict']['sig_R_maj_min'],
                      star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      data_lim=locals()[name+'dict']['data_lim'], 
                  color=color, alpha=0.2, 
                  label=name+' '+phot[3],
                  sfrange=phot[4])
    
    plot_kennicutt_half(ax=ax, total_gas_data=total_gas_data,
                      epicycl=locals()[name+'dict']['epicyclicFreq_real'], 
                      gas_approx=locals()[name+'dict']['spl_gas'], 
                      sound_vel=locals()[name+'dict']['sound_vel'], 
                      scale=locals()[name+'dict']['scale'], 
                      sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
                      sigma_min=locals()[name+'dict']['sig_R_maj_min'],
                      star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      data_lim=locals()[name+'dict']['data_lim'], 
                  color=color, alpha=0.2, 
                  label=name+' '+phot[3],
                  sfrange=phot[4])
    
    plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
    
plt.axvline(x=1, alpha=0.5)
plot_Q_levels(ax, [1.5, 2., 3.])
plt.legend()
plt.xlim(0, 3)
plt.ylim(0, 1.3)
plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
plt.xlabel(r'$R/R_{SF}$', fontsize=20)
plt.show()


1F vs 2F

1F vs 2F (среднее):


In [68]:
def plot_1f_vs_2f(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, sfrange=None):

    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))

    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))
   
    rr = zip(*total_gas_data)[0]
    ax.plot(invQg, (np.array(invQeff_min)+np.array(invQeff_max))/2., '.', color=color, alpha=0.8)

In [81]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
            'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
            'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
            'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
            'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
            'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
            'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}

fig = plt.figure(figsize=[8, 8])

for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
    ax = plt.gca()
    
    color = cm.rainbow(np.linspace(0, 1, 7))[ind]
    
    phot = kenn_photom[name]
    mud = locals()[name+'dict'][phot[0]]
    h = locals()[name+'dict'][phot[1]]
#     print name
#     print mud
#     print h
    
    if name in ['n4258', 'n4725']:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
    elif name == 'n2985':
        total_gas_data=locals()[name+'dict']['total_gas_data_']
    else:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in 
                                                                   zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
            
            
    if total_gas_data[0][0] < 0.1:
        total_gas_data = total_gas_data[1:]
    
#     print total_gas_data

    plot_1f_vs_2f(ax=ax, total_gas_data=total_gas_data,
                      epicycl=locals()[name+'dict']['epicyclicFreq_real'], 
                      gas_approx=locals()[name+'dict']['spl_gas'], 
                      sound_vel=locals()[name+'dict']['sound_vel'], 
                      scale=locals()[name+'dict']['scale'], 
                      sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
                      sigma_min=locals()[name+'dict']['sig_R_maj_min'],
                      star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      data_lim=locals()[name+'dict']['data_lim'], 
                  color=color, alpha=0.2, 
                  label=name+' '+phot[3],
                  sfrange=phot[4])
    plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
    
plt.plot([0., 1.], [0., 1.], '--', alpha=0.1)
plt.legend(loc='lower right')
plt.xlim(0, 1.)
plt.ylim(0, 1.)
plt.ylabel(r'$Q_{eff}^{-1}$', fontsize=20)
plt.xlabel(r'$Q_g^{-1}$', fontsize=20)
plt.show()


Только молек


In [82]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
            'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
            'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
            'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
            'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
            'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
            'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}

fig = plt.figure(figsize=[8, 8])

for ind, name in enumerate(['n2985', 'n4258', 'n4725']):
    ax = plt.gca()
    
    color = cm.rainbow(np.linspace(0, 1, 7))[ind]
    
    phot = kenn_photom[name]
    mud = locals()[name+'dict'][phot[0]]
    h = locals()[name+'dict'][phot[1]]
#     print name
#     print mud
#     print h
    
    if name in ['n4258', 'n4725']:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
    elif name == 'n2985':
        total_gas_data=locals()[name+'dict']['total_gas_data_']
    else:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in 
                                                                   zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
            
            
    if total_gas_data[0][0] < 0.1:
        total_gas_data = total_gas_data[1:]
    
#     print total_gas_data

    plot_1f_vs_2f(ax=ax, total_gas_data=total_gas_data,
                      epicycl=locals()[name+'dict']['epicyclicFreq_real'], 
                      gas_approx=locals()[name+'dict']['spl_gas'], 
                      sound_vel=locals()[name+'dict']['sound_vel'], 
                      scale=locals()[name+'dict']['scale'], 
                      sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
                      sigma_min=locals()[name+'dict']['sig_R_maj_min'],
                      star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      data_lim=locals()[name+'dict']['data_lim'], 
                  color=color, alpha=0.2, 
                  label=name+' '+phot[3],
                  sfrange=phot[4])
    plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])

plt.plot([0., 1.], [0., 1.], '--', alpha=0.1)
plt.legend(loc='lower right')
plt.xlim(0, 1.)
plt.ylim(0, 1.)
plt.ylabel(r'$Q_{eff}^{-1}$', fontsize=20)
plt.xlabel(r'$Q_g^{-1}$', fontsize=20)
plt.show()


Верхнее и нижнее значения:


In [83]:
def plot_1f_vs_2f_full(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, sfrange=None):

    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))

    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))
   
    rr = zip(*total_gas_data)[0]
    ax.plot(invQg, invQeff_min, 'o', color=color, alpha=0.8)
    ax.plot(invQg, invQeff_max, '*', color=color, alpha=0.8)

In [84]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
            'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
            'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
            'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
            'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
            'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
            'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}

fig = plt.figure(figsize=[8, 8])

for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
    ax = plt.gca()
    
    color = cm.rainbow(np.linspace(0, 1, 7))[ind]
    
    phot = kenn_photom[name]
    mud = locals()[name+'dict'][phot[0]]
    h = locals()[name+'dict'][phot[1]]
#     print name
#     print mud
#     print h
    
    if name in ['n4258', 'n4725']:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
    elif name == 'n2985':
        total_gas_data=locals()[name+'dict']['total_gas_data_']
    else:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in 
                                                                   zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
            
            
    if total_gas_data[0][0] < 0.1:
        total_gas_data = total_gas_data[1:]
    
#     print total_gas_data

    plot_1f_vs_2f_full(ax=ax, total_gas_data=total_gas_data,
                      epicycl=locals()[name+'dict']['epicyclicFreq_real'], 
                      gas_approx=locals()[name+'dict']['spl_gas'], 
                      sound_vel=locals()[name+'dict']['sound_vel'], 
                      scale=locals()[name+'dict']['scale'], 
                      sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
                      sigma_min=locals()[name+'dict']['sig_R_maj_min'],
                      star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      data_lim=locals()[name+'dict']['data_lim'], 
                  color=color, alpha=0.2, 
                  label=name+' '+phot[3],
                  sfrange=phot[4])
    plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
    
plt.plot([0., 1.], [0., 1.], '--', alpha=0.1)
plt.legend(loc='lower right')
plt.xlim(0, 1.)
plt.ylim(0, 1.)
plt.ylabel(r'$Q_{eff}^{-1}$', fontsize=20)
plt.xlabel(r'$Q_g^{-1}$', fontsize=20)
plt.show()


Только молек


In [85]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
            'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
            'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
            'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
            'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
            'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
            'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}

fig = plt.figure(figsize=[8, 8])

for ind, name in enumerate(['n2985', 'n4258', 'n4725']):
    ax = plt.gca()
    
    color = cm.rainbow(np.linspace(0, 1, 7))[ind]
    
    phot = kenn_photom[name]
    mud = locals()[name+'dict'][phot[0]]
    h = locals()[name+'dict'][phot[1]]
#     print name
#     print mud
#     print h
    
    if name in ['n4258', 'n4725']:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
    elif name == 'n2985':
        total_gas_data=locals()[name+'dict']['total_gas_data_']
    else:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in 
                                                                   zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
            
            
    if total_gas_data[0][0] < 0.1:
        total_gas_data = total_gas_data[1:]
    
#     print total_gas_data

    plot_1f_vs_2f_full(ax=ax, total_gas_data=total_gas_data,
                      epicycl=locals()[name+'dict']['epicyclicFreq_real'], 
                      gas_approx=locals()[name+'dict']['spl_gas'], 
                      sound_vel=locals()[name+'dict']['sound_vel'], 
                      scale=locals()[name+'dict']['scale'], 
                      sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
                      sigma_min=locals()[name+'dict']['sig_R_maj_min'],
                      star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      data_lim=locals()[name+'dict']['data_lim'], 
                  color=color, alpha=0.2, 
                  label=name+' '+phot[3],
                  sfrange=phot[4])
    plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
    
plt.plot([0., 1.], [0., 1.], '--', alpha=0.1)
plt.legend(loc='lower right')
plt.xlim(0, 1.)
plt.ylim(0, 1.)
plt.ylabel(r'$Q_{eff}^{-1}$', fontsize=20)
plt.xlabel(r'$Q_g^{-1}$', fontsize=20)
plt.show()



In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [48]:
print locals()[name+'dict']['total_gas_data_']
print locals()[name+'dict']['epicyclicFreq_real']
print locals()[name+'dict']['spl_gas'](15.)
print locals()[name+'dict']['sound_vel']
print locals()[name+'dict']['scale']
print locals()[name+'dict']['sig_R_maj_max'](15.)
print locals()[name+'dict']['sig_R_maj_min'](15.)
print surf_density(mu=mu_disc(15., mu0=locals()[name+'dict']['mudK'], h=locals()[name+'dict']['hK']), M_to_L=1.42, band='K')
print surf_density(mu=mu_disc(15., mu0=locals()[name+'dict']['mudK'], h=locals()[name+'dict']['hK']), M_to_L=1.42, band='K')
print locals()[name+'dict']['data_lim']
print locals()[name+'dict']['disk_scales']


[(0.0, 22.780000000000001), (15.0, 17.254108218335578), (30.0, 13.356413456740752), (45.0, 10.187361920100734), (60.0, 7.9384773506614295), (75.0, 6.6393114577321892), (90.0, 5.5382294229423552)]
<function epicyclicFreq_real at 0x000000000A69A5F8>
272.083511499
6.0
0.265
185.392460246
121.731459498
902.803666556
902.803666556
77.1001025481
[(34.4, 'R'), (32.4, 'B'), (34.4, 'R_max'), (14.1, 'J'), (43.0188679245283, 'R'), (40.7, 'V'), (18.0, 'V'), (28.1701, 'g'), (28.0063, 'r'), (28.0063, 'i')]

In [ ]:


In [ ]:


In [49]:
# 338
plot_2f_vs_1f(ax=axes[4], total_gas_data=zip(r_g_dens, [He_coeff*(y_interp_(l[0], h_disc_I) + l[1]) for l in zip(r_g_dens, gas_dens)])[1:], 
                  epicycl=epicyclicFreq_real_, 
                  gas_approx=spl_gas, 
                  sound_vel=sound_vel, 
                  scale=scale, 
                  sigma_max=sig_R_maj_max, 
                  sigma_min=sig_R_maj_min, 
                  star_density_max=lambda l: surf_density(mu_disc(l, mu0=mu0d_Ic, h=h_disc_I), 6.23, 'I'), 
                  star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_Ic, h=h_disc_I), 6.23, 'I'), 
                  data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='I maxdisc')
    
    plot_2f_vs_1f(ax=axes[4], total_gas_data=zip(r_g_dens, [He_coeff*(y_interp_(l[0], h_disc_B) + l[1]) for l in zip(r_g_dens, gas_dens)])[1:], 
                  epicycl=epicyclicFreq_real_, 
                  gas_approx=spl_gas, 
                  sound_vel=sound_vel, 
                  scale=scale, 
                  sigma_max=sig_R_maj_max, 
                  sigma_min=sig_R_maj_min, 
                  star_density_max=lambda l: surf_density(mu_disc(l, mu0=mu0d_Bc, h=h_disc_B), 6.99, 'B'), 
                  star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_Bc, h=h_disc_B), 6.99, 'B'), 
                  data_lim=data_lim, color='y', alpha=0.2, disk_scales=disk_scales, label='B submaxdisc')
    
# 1167
    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:], 
                  epicycl=epicyclicFreq_real, 
                  gas_approx=spl_gas, 
                  sound_vel=sound_vel, 
                  scale=scale, 
                  sigma_max=sig_R_maj_max, 
                  sigma_min=sig_R_maj_min, 
                  star_density_max=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), 5.53, 'R'), 
                  star_density_min=lambda l: surf_density(mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), 5.53, 'R'), 
                  data_lim=data_lim, color='g', alpha=0.3, disk_scales=disk_scales, label='R maxdisc')
    
# 2985
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=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.42, band='K'), 
              star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mudK, h=hK), M_to_L=1.42, band='K'), 
              data_lim=data_lim, color='g', alpha=0.2, disk_scales=disk_scales2, label='K Heidt maxdisc')

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: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu0d_s4g, h=h_disc_s4g), 1.35),
                   lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu0d_s4g_2, h=h_disc_s4g_2), 1.35)))(l), 
              star_density_min=lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu0d_s4g, h=h_disc_s4g), 1.35),
                   lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu0d_s4g_2, h=h_disc_s4g_2), 1.35)))(l), 
              data_lim=data_lim, color='y', alpha=0.2, disk_scales=disk_scales2, label='S4G 2d maxdisc')
    

#3898
    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')

#4258
    plot_2f_vs_1f(ax=axes[4], total_gas_data=zip(r_g_dens, map(lambda l: l, gas_dens))[:15], 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=mu_disc(l, mu0=mu0d_I, h=h_disc_I), M_to_L=0.82, band='I'), 
              star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_I, h=h_disc_I), M_to_L=0.82, band='I'), 
              data_lim=data_lim, color='g', alpha=0.2, disk_scales=disk_scales, label='I Yoshino maxdisc')
    
    plot_2f_vs_1f(ax=axes[4], total_gas_data=zip(r_g_dens, map(lambda l: l, gas_dens))[:15], 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: s4g_surf_density(mu_disc(l, mu0=mu0d_36, h=h_disc_36), M_to_L=0.38), 
              star_density_min=lambda l: s4g_surf_density(mu_disc(l, mu0=mu0d_36, h=h_disc_36), M_to_L=0.38), 
              data_lim=data_lim, color='m', alpha=0.2, disk_scales=disk_scales, label='SPITZER 3.6 maxdisc')

# 4725
    plot_2f_vs_1f(ax=axes[4], total_gas_data=zip(r_g_dens, map(lambda l: l, gas_dens))[:15], 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=mu_disc(l, mu0=mu0d_H, h=h_disc_H), M_to_L=0.68, band='H'), 
              star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_H, h=h_disc_H), M_to_L=0.68, band='H'), 
              data_lim=data_lim, color='y', alpha=0.3, disk_scales=disk_scales, label='H Heidt maxdisc')
    
    plot_2f_vs_1f(ax=axes[4], total_gas_data=zip(r_g_dens, map(lambda l: l, gas_dens))[:15], 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: s4g_surf_density(mu_disc(l, mu0=mu0d_s4g, h=h_disc_s4g), 1.61), 
              star_density_min=lambda l: s4g_surf_density(mu_disc(l, mu0=mu0d_s4g, h=h_disc_s4g), 1.61), 
              data_lim=data_lim, color='g', alpha=0.2, disk_scales=disk_scales, label='S4G maxdisc')
#  SF   55/100/175

# 5533

#change this
    total_gas_data_ = zip(r_g_dens, [He_coeff*(y_interp_(l[0], hi_r) + l[1]) for l in zip(r_g_dens, gas_dens)])[1:10]
    total_gas_data_2 = zip(r_g_dens, [He_coeff*(y_interp_2(l[0], hi_r) + l[1]) for l in zip(r_g_dens, gas_dens)])[1:10]

    plot_2f_vs_1f(ax=axes[4], 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=mu_disc(l, mu0=MU0_r, h=hi_r), M_to_L=4.41, band='r'), 
                  star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=MU0_r, h=hi_r), M_to_L=4.41, band='r'), 
                  data_lim=data_lim, color='g', alpha=0.2, disk_scales=disk_scales, label='r(g-i) maxdisc')
    
    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:10]
    total_gas_data_2 = zip(r_g_dens, [He_coeff*(y_interp_2(l[0], h_disc_R) + l[1]) for l in zip(r_g_dens, gas_dens)])[1:10]

    plot_2f_vs_1f(ax=axes[4], total_gas_data=total_gas_data_2, 
                  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=mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L=7.15, band='R'), 
                  star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mu0d_Rc, h=h_disc_R), M_to_L=7.15, band='R'), 
                  data_lim=data_lim, color='r', alpha=0.2, disk_scales=disk_scales, label='R Noord maxdisc 2')


  File "<ipython-input-49-ef18fd200c32>", line 13
    plot_2f_vs_1f(ax=axes[4], total_gas_data=zip(r_g_dens, [He_coeff*(y_interp_(l[0], h_disc_B) + l[1]) for l in zip(r_g_dens, gas_dens)])[1:],
    ^
IndentationError: unexpected indent

In [ ]:

Two at the same plot


In [16]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
            'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
            'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
            'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
            'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
            'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
            'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}

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

for ind, name in enumerate(['n1167', 'n3898']):
    ax = plt.gca()
    
    color = cm.rainbow(np.linspace(0, 1, 7))[ind]
    
    phot = kenn_photom[name]
    mud = locals()[name+'dict'][phot[0]]
    h = locals()[name+'dict'][phot[1]]
#     print name
#     print mud
#     print h
    
    if name in ['n4258', 'n4725']:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
    elif name == 'n2985':
        total_gas_data=locals()[name+'dict']['total_gas_data_']
    else:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in 
                                                                   zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
            
            
    if total_gas_data[0][0] < 0.1:
        total_gas_data = total_gas_data[1:]
    
#     print total_gas_data

    plot_2f_vs_1f(ax=ax, total_gas_data=total_gas_data,
                      epicycl=locals()[name+'dict']['epicyclicFreq_real'], 
                      gas_approx=locals()[name+'dict']['spl_gas'], 
                      sound_vel=locals()[name+'dict']['sound_vel'], 
                      scale=locals()[name+'dict']['scale'], 
                      sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
                      sigma_min=locals()[name+'dict']['sig_R_maj_min'],
                      star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      data_lim=locals()[name+'dict']['data_lim'], 
                  color=color, alpha=0.2, 
                  label=name+' '+phot[3])
    
    plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
    
# plt.axvline(x=1, alpha=0.5)
plot_Q_levels(ax, [1.5, 2., 3.])
plt.legend()
plt.xlim(0, 130.)
plt.ylim(0, 1.3)
plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
plt.xlabel(r'$R$', fontsize=20)
plt.show()


Wang & Silk approximation

Wang & Silk 2011 approximation:


In [45]:
def plot_WS_vs_2f(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, sfrange=None):
    
    rr = zip(*total_gas_data)[0]

    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))
    
    invQwff_WS_min = [1./Qg(epicycl=epicycl(gas_approx, r, scale), sound_vel=sound_vel, gas_density=gd) + 1./Qs(
        epicycl=epicycl(gas_approx, r, scale), sigma=sigma_max(r), star_density=star_density_min(r)) for r, gd in total_gas_data]
    
    ax.plot(invQwff_WS_min, invQeff_min, 'd', color=color, alpha=0.8)

    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))
        
    invQwff_WS_max = [1./Qg(epicycl=epicycl(gas_approx, r, scale), sound_vel=sound_vel, gas_density=gd) + 1./Qs(
        epicycl=epicycl(gas_approx, r, scale), sigma=sigma_min(r), star_density=star_density_max(r)) for r, gd in total_gas_data]
    
    ax.plot(invQwff_WS_max, invQeff_max, 'o', color=color, alpha=0.8)

In [46]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
            'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
            'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
            'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
            'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
            'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
            'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}

fig = plt.figure(figsize=[8, 8])

for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
    ax = plt.gca()
    
    color = cm.rainbow(np.linspace(0, 1, 7))[ind]
    
    phot = kenn_photom[name]
    mud = locals()[name+'dict'][phot[0]]
    h = locals()[name+'dict'][phot[1]]
#     print name
#     print mud
#     print h
    
    if name in ['n4258', 'n4725']:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
    elif name == 'n2985':
        total_gas_data=locals()[name+'dict']['total_gas_data_']
    else:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in 
                                                                   zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
            
            
    if total_gas_data[0][0] < 0.1:
        total_gas_data = total_gas_data[1:]
    
#     print total_gas_data

    plot_WS_vs_2f(ax=ax, total_gas_data=total_gas_data,
                      epicycl=locals()[name+'dict']['epicyclicFreq_real'], 
                      gas_approx=locals()[name+'dict']['spl_gas'], 
                      sound_vel=locals()[name+'dict']['sound_vel'], 
                      scale=locals()[name+'dict']['scale'], 
                      sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
                      sigma_min=locals()[name+'dict']['sig_R_maj_min'],
                      star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      data_lim=locals()[name+'dict']['data_lim'], 
                  color=color, alpha=0.2, 
                  label=name+' '+phot[3],
                  sfrange=phot[4])
    plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
    
plt.plot([0., 1.], [0., 1.], '--', alpha=0.1)
plt.legend(loc='lower right')
plt.xlim(0, 1.)
plt.ylim(0, 1.)
plt.ylabel(r'$Q_{eff}^{-1}$', fontsize=20)
plt.xlabel(r'$Q_{WS}^{-1}$', fontsize=20)
plt.grid()
plt.show()



In [42]:
def plot_WS_vs_2f_dist(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, sfrange=None):
    
    rr = zip(*total_gas_data)[0]

    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))
    
    invQwff_WS_min = [1./Qg(epicycl=epicycl(gas_approx, r, scale), sound_vel=sound_vel, gas_density=gd) + 1./Qs(
        epicycl=epicycl(gas_approx, r, scale), sigma=sigma_max(r), star_density=star_density_min(r)) for r, gd in total_gas_data]
    
    ax.plot(rr, [l[0]/l[1] for l in zip(invQwff_WS_min, invQeff_min)], 'd', color=color, alpha=0.8)

    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))
        
    invQwff_WS_max = [1./Qg(epicycl=epicycl(gas_approx, r, scale), sound_vel=sound_vel, gas_density=gd) + 1./Qs(
        epicycl=epicycl(gas_approx, r, scale), sigma=sigma_min(r), star_density=star_density_max(r)) for r, gd in total_gas_data]
    
    ax.plot(rr, [l[0]/l[1] for l in zip(invQwff_WS_max, invQeff_max)], 'o', color=color, alpha=0.8)

In [43]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
            'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
            'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
            'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
            'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
            'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
            'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}

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

for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
    ax = plt.gca()
    
    color = cm.rainbow(np.linspace(0, 1, 7))[ind]
    
    phot = kenn_photom[name]
    mud = locals()[name+'dict'][phot[0]]
    h = locals()[name+'dict'][phot[1]]
#     print name
#     print mud
#     print h
    
    if name in ['n4258', 'n4725']:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
    elif name == 'n2985':
        total_gas_data=locals()[name+'dict']['total_gas_data_']
    else:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in 
                                                                   zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
            
            
    if total_gas_data[0][0] < 0.1:
        total_gas_data = total_gas_data[1:]
    
#     print total_gas_data

    plot_WS_vs_2f_dist(ax=ax, total_gas_data=total_gas_data,
                      epicycl=locals()[name+'dict']['epicyclicFreq_real'], 
                      gas_approx=locals()[name+'dict']['spl_gas'], 
                      sound_vel=locals()[name+'dict']['sound_vel'], 
                      scale=locals()[name+'dict']['scale'], 
                      sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
                      sigma_min=locals()[name+'dict']['sig_R_maj_min'],
                      star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      data_lim=locals()[name+'dict']['data_lim'], 
                  color=color, alpha=0.2, 
                  label=name+' '+phot[3],
                  sfrange=phot[4])
    plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
    
# plt.plot([0., 1.], [0., 1.], '--', alpha=0.1)
plt.legend(loc='upper right')
plt.xlim(0)
plt.ylim(1.0)
plt.ylabel(r'$Q_{WS}^{-1}/Q_{eff}^{-1}$', fontsize=20)
plt.xlabel(r'$R$', fontsize=20)
plt.grid()
plt.show()


Все 2F в одной


In [44]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
            'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
            'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
            'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
            'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
            'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
            'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}

fi, axes = plt.subplots(ncols=1, nrows=7, figsize=[20., 26.], sharex=True)nbv

for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
# for ind, name in enumerate(['n338']):
    ax = axes[ind]
    
    color = cm.rainbow(np.linspace(0, 1, 7))[ind]
    
    phot = kenn_photom[name]
    mud = locals()[name+'dict'][phot[0]]
    h = locals()[name+'dict'][phot[1]]
#     print name
#     print mud
#     print h
    
    if name in ['n4258', 'n4725']:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
    elif name == 'n2985':
        total_gas_data=locals()[name+'dict']['total_gas_data_']
    else:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in 
                                                                   zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
            
            
    if total_gas_data[0][0] < 0.1:
        total_gas_data = total_gas_data[1:]
    
#     print total_gas_data

    plot_2f_vs_1f(ax=ax, total_gas_data=total_gas_data,
                      epicycl=locals()[name+'dict']['epicyclicFreq_real'], 
                      gas_approx=locals()[name+'dict']['spl_gas'], 
                      sound_vel=locals()[name+'dict']['sound_vel'], 
                      scale=locals()[name+'dict']['scale'], 
                      sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
                      sigma_min=locals()[name+'dict']['sig_R_maj_min'],
                      star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      data_lim=locals()[name+'dict']['data_lim'], 
                  color=color, alpha=0.2, 
                  label='')
    
    ax.set_xlim(0, 130)
    ax.set_ylim(0, 1.2)
    ax.xaxis.grid(True)
    
    if ind == 6:
        ax.set_xlabel(r'$R, arcsec$', fontsize=20)
    ax.set_ylabel(r'$Q_{eff}^{-1}$', fontsize=20)
    for q_ in [1., 1.5, 2., 3.]:
        ax.axhline(y=1./q_, lw=10, alpha=0.15, color='grey')
        
    ax.text(95, 0.85, 'NGC '+name[1:]+' $'+phot[3]+'$', fontsize=30)
#     ax.set_yticks(['', '0.2', '0.4', '0.6', '0.8', '1.0', ''])
    plt.setp(ax.get_yticklabels()[0], visible=False)    
    plt.setp(ax.get_yticklabels()[-1], visible=False)
    
    
#     plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
    
# plt.axvline(x=1, alpha=0.5)
# plot_Q_levels(ax, [1.5, 2., 3.])
# plt.legend()
# plt.xlim(0, 130.)
# plt.ylim(0, 1.3)
# plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
# plt.xlabel(r'$R$', fontsize=20)
plt.subplots_adjust(wspace=0, hspace=0)
plt.show()



In [51]:
kenn_photom={'n338' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60.],
            'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40.],
            'n2985' : ['mudK', 'hK', 1.42, 'K', 70.],
            'n3898' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80.],
            'n4258' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150.],
            'n4725' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55.],
            'n5533' : ['MU0_r', 'hi_r', 4.41, 'r', 100.]}

fi, axes = plt.subplots(ncols=1, nrows=7, figsize=[20., 26.], sharex=True)

for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
# for ind, name in enumerate(['n338']):
    ax = axes[ind]
    
    color = cm.rainbow(np.linspace(0, 1, 7))[ind]
    
    phot = kenn_photom[name]
    mud = locals()[name+'dict'][phot[0]]
    h = locals()[name+'dict'][phot[1]]
#     print name
#     print mud
#     print h
    
    if name in ['n4258', 'n4725']:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
    elif name == 'n2985':
        total_gas_data=locals()[name+'dict']['total_gas_data_']
    else:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in 
                                                                   zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]
            
            
    if total_gas_data[0][0] < 0.1:
        total_gas_data = total_gas_data[1:]
    
#     print total_gas_data

    plot_2f_vs_1f(ax=ax, total_gas_data=total_gas_data,
                      epicycl=locals()[name+'dict']['epicyclicFreq_real'], 
                      gas_approx=locals()[name+'dict']['spl_gas'], 
                      sound_vel=locals()[name+'dict']['sound_vel'], 
                      scale=locals()[name+'dict']['scale'], 
                      sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
                      sigma_min=locals()[name+'dict']['sig_R_maj_min'],
                      star_density_max=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      star_density_min=lambda l: surf_density(mu=mu_disc(l, mu0=mud, h=h), M_to_L=phot[2], band=phot[3]), 
                      data_lim=locals()[name+'dict']['data_lim'], 
                  color=color, alpha=0.2, 
                  label='')
    
    ax.set_xlim(0., 130)
    ax.set_ylim(0.01, 1.5)
    ax.xaxis.grid(True)
    
    if ind == 6:
        ax.set_xlabel(r'$R, arcsec$', fontsize=20)
    ax.set_ylabel(r'$Q_{eff}^{-1}$', fontsize=20)
    for q_ in [1., 1.5, 2., 3.]:
        ax.axhline(y=1./q_, lw=10, alpha=0.15, color='grey')
        
    ax.text(95, 0.85, 'NGC '+name[1:]+' $'+phot[3]+'$', fontsize=30)
#     ax.set_yticks(['', '0.2', '0.4', '0.6', '0.8', '1.0', ''])
    plt.setp(ax.get_yticklabels()[0], visible=False)    
    plt.setp(ax.get_yticklabels()[-1], visible=False)
    ax.set_yscale('log')
    
    
#     plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
    
# plt.axvline(x=1, alpha=0.5)
# plot_Q_levels(ax, [1.5, 2., 3.])
# plt.legend()
# plt.xlim(0, 130.)
# plt.ylim(0, 1.3)
# plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
# plt.xlabel(r'$R$', fontsize=20)
plt.subplots_adjust(wspace=0, hspace=0)
plt.show()



In [42]:
from random import shuffle
random.seed(42)
kenn_photom={'n338_I' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60., 15.],
             'n338_B' : ['mu0d_Bc', 'h_disc_B', 6.99, 'B', 60., 15.],
            'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40., 6.7],
            'n2985_K' : ['mudK', 'hK', 1.42, 'K', 70., 14.],
            'n2985_S4G' : [('mu0d_s4g', 'mu0d_s4g_2'), ('h_disc_s4g', 'h_disc_s4g_2'), 1.35, 'S4G', 70., 6.3],
            'n3898_R' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80., 8.8],
            'n3898_R2' : [('mu_inner_approx', 'mu_out'), ('h_approx', 'h_out'), 2.93, 'R', 80., 8.8],
            'n4258_I' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150., 15.],
            'n4258_S4G' : ['mu0d_36', 'h_disc_36', 0.38, 'S4G', 150., 15.],
            'n4725_H' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55., 10.14],
            'n4725_S4G' : ['mu0d_s4g', 'h_disc_s4g', 1.61, 'S4G', 55., 10.14],
            'n5533_r' : ['MU0_r', 'hi_r', 4.41, 'r', 100., 8.9],
            'n5533_R' : ['mu0d_Rc', 'h_disc_R', 7.15, 'R', 100., 9.9]}
color_shuffle = range(13)
shuffle(color_shuffle)

In [43]:
%%time
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):
    '''Картинка сравнения 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))

    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))

#     invQg = map(lambda l: l*1.6, invQg)
#     invQeff_min = map(lambda l: l*1.6, invQeff_min)
#     invQeff_max = map(lambda l: l*1.6, invQeff_max)
    
    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_SF(ax)
    plot_data_lim(ax, data_lim)
    for h, annot in disk_scales:
        plot_disc_scale(h, ax, '')
    plot_Q_levels(ax, [1., 1.5, 2., 3.])
    ax.legend(fontsize=20)
    

fi, axes = plt.subplots(ncols=1, nrows=7, figsize=[20., 26.], sharex=True)

color_counter = 0


for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
# for ind, name in enumerate(['n3898']):
    ax = axes[ind]
    
    for key in kenn_photom.keys():
        
        if name in key:
            print key
            phot = kenn_photom[key]
            color = cm.rainbow(np.linspace(0, 1, 13))[color_shuffle[color_counter]]
            
            try:
                h = locals()[name+'dict'][phot[1][0]]
            except KeyError:
                h = locals()[name+'dict'][phot[1]]
        #     print name
        #     print mud
        #     print h

            if name in ['n4258', 'n4725']:
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
            elif name == 'n2985':
                total_gas_data=locals()[name+'dict']['total_gas_data_']
            else:
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in 
                                                                           zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]


            while total_gas_data[0][0] < phot[5]:
                total_gas_data = total_gas_data[1:]

#             print total_gas_data
        
            ML = phot[2]
            band = phot[3]
        
            if type(phot[0]) == tuple and phot[3] == 'S4G':
                h1 = globals()[name+'dict'][phot[1][0]]
                h2 = globals()[name+'dict'][phot[1][1]]
                mu01 = globals()[name+'dict'][phot[0][0]]
                mu02 = globals()[name+'dict'][phot[0][1]]
                star_density = lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu01, h=h1), ML),
                               lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu02, h=h2), ML)))(l) 
            elif type(phot[0]) == tuple:
                h1 = globals()[name+'dict'][phot[1][0]]
                h2 = globals()[name+'dict'][phot[1][1]]
                mu01 = globals()[name+'dict'][phot[0][0]]
                mu02 = globals()[name+'dict'][phot[0][1]]
                star_density = lambda l: tot_dens((lambda l1: surf_density(mu=mu_disc(l1, mu0=mu01, h=h1), M_to_L=ML, band=band),
                               lambda l2: surf_density(mu=mu_disc(l2, mu0=mu02, h=h2), M_to_L=ML, band=band)))(l)
            elif phot[3] == 'S4G':
                star_density = lambda l: s4g_surf_density(mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), ML)
            else:
                star_density = lambda l: surf_density(mu=mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), M_to_L=ML, band=band)
            
            
            if name == 'n338':
                epicycl = locals()[name+'dict']['epicyclicFreq_real_']
            else:
                epicycl = locals()[name+'dict']['epicyclicFreq_real']
                            
            plot_2f_vs_1f(ax=ax, total_gas_data=total_gas_data,
                              epicycl=epicycl,
                              gas_approx=locals()[name+'dict']['spl_gas'], 
                              sound_vel=locals()[name+'dict']['sound_vel'], 
                              scale=locals()[name+'dict']['scale'], 
                              sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
                              sigma_min=locals()[name+'dict']['sig_R_maj_min'],
                              star_density_max=star_density, 
                              star_density_min=star_density, 
                              data_lim=locals()[name+'dict']['data_lim'], 
                          color=color, alpha=0.3, 
                          label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$',
                          disk_scales=[(h, band)])            

            ax.set_xlim(0, 130)
            if name == 'n338':
                ax.set_ylim(0, 1.4)
            else:
                ax.set_ylim(0, 1.2)
            ax.xaxis.grid(True)
            locals()[name+'dict']['plot_SF'](ax)

            if ind == 6:
                ax.set_xlabel(r'$R, arcsec$', fontsize=20)
            ax.set_ylabel(r'$Q_{eff}^{-1}$', fontsize=20)
            for q_ in [1., 1.5, 2., 3.]:
                ax.axhline(y=1./q_, lw=10, alpha=0.15, color='grey')

#             ax.text(95, 0.85, 'NGC '+name[1:]+' $'+phot[3]+'$', fontsize=30)
        #     ax.set_yticks(['', '0.2', '0.4', '0.6', '0.8', '1.0', ''])
            plt.setp(ax.get_yticklabels()[0], visible=False)    
            plt.setp(ax.get_yticklabels()[-1], visible=False)
            color_counter+=1
    
    
#     plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
    
# plt.axvline(x=1, alpha=0.5)
# plot_Q_levels(ax, [1.5, 2., 3.])
# plt.legend()
# plt.xlim(0, 130.)
# plt.ylim(0, 1.3)
# plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
# plt.xlabel(r'$R$', fontsize=20)
plt.subplots_adjust(wspace=0, hspace=0)
plt.show()


n338_I
n338_B
n1167
n2985_S4G
n2985_K
n3898_R2
n3898_R
n4258_S4G
n4258_I
n4725_S4G
n4725_H
n5533_R
n5533_r
Wall time: 2min 48s

In [ ]:


In [49]:
%%time

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):
    '''Картинка сравнения 2F и 1F критерия для разных фотометрий и величин sig_R, 
    куда подается весь газ, результат НЕ исправляется за осесимметричные возмущения.'''

    invQg, invQs_mi, 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))

    invQg, invQs_ma, 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))

#     invQg = map(lambda l: l*1.6, invQg)
#     invQeff_min = map(lambda l: l*1.6, invQeff_min)
#     invQeff_max = map(lambda l: l*1.6, invQeff_max)
    
    rr = zip(*total_gas_data)[0]
    
#     ax.fill_between(rr, invQeff_min, invQeff_max, color=color, alpha=alpha, label=label)
    ax.plot(rr, invQs_mi, '--', color=color, alpha=0.75)
    ax.plot(rr, invQs_ma, '--', color=color, alpha=0.75)
    ax.plot(rr, invQg, '--', color='b')

fi, axes = plt.subplots(ncols=1, nrows=13, figsize=[20., 50.], sharex=True)  


for ind, key in enumerate(kenn_photom.keys()):

    ax = axes[ind]
    name = key.split('_')[0]
    print key, name
    phot = kenn_photom[key]
    color = cm.rainbow(np.linspace(0, 1, 13))[color_shuffle[ind]]

    try:
        h = locals()[name+'dict'][phot[1][0]]
    except KeyError:
        h = locals()[name+'dict'][phot[1]]
#     print name
#     print mud
#     print h

    if name in ['n4258', 'n4725']:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
    elif name == 'n2985':
        total_gas_data=locals()[name+'dict']['total_gas_data_']
    else:
        total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in 
                                                                   zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]


    while total_gas_data[0][0] < phot[5]:
        total_gas_data = total_gas_data[1:]

#             print total_gas_data

    ML = phot[2]
    band = phot[3]

    if type(phot[0]) == tuple and phot[3] == 'S4G':
        h1 = globals()[name+'dict'][phot[1][0]]
        h2 = globals()[name+'dict'][phot[1][1]]
        mu01 = globals()[name+'dict'][phot[0][0]]
        mu02 = globals()[name+'dict'][phot[0][1]]
        star_density = lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu01, h=h1), ML),
                       lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu02, h=h2), ML)))(l) 
    elif type(phot[0]) == tuple:
        h1 = globals()[name+'dict'][phot[1][0]]
        h2 = globals()[name+'dict'][phot[1][1]]
        mu01 = globals()[name+'dict'][phot[0][0]]
        mu02 = globals()[name+'dict'][phot[0][1]]
        star_density = lambda l: tot_dens((lambda l1: surf_density(mu=mu_disc(l1, mu0=mu01, h=h1), M_to_L=ML, band=band),
                       lambda l2: surf_density(mu=mu_disc(l2, mu0=mu02, h=h2), M_to_L=ML, band=band)))(l)
    elif phot[3] == 'S4G':
        star_density = lambda l: s4g_surf_density(mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), ML)
    else:
        star_density = lambda l: surf_density(mu=mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), M_to_L=ML, band=band)


    if name == 'n338':
        epicycl = locals()[name+'dict']['epicyclicFreq_real_']
    else:
        epicycl = locals()[name+'dict']['epicyclicFreq_real']

    plot_param_depend(ax=ax,
                      epicycl=epicycl,
                      gas_approx=locals()[name+'dict']['spl_gas'], 
                      sound_vel=list(np.linspace(4., 20., 50)),
                      N=50,
                      scale=locals()[name+'dict']['scale'], 
                      sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
                      sigma_min=locals()[name+'dict']['sig_R_maj_min'],
                      star_density_max=star_density, 
                      star_density_min=star_density, 
                      data_lim=locals()[name+'dict']['data_lim'], 
                      color=color, alpha=0.3, 
                      label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$',
                      total_gas_data=total_gas_data)
    
    plot_2f_vs_1f_(ax=ax, total_gas_data=total_gas_data,
                              epicycl=epicycl,
                              gas_approx=locals()[name+'dict']['spl_gas'], 
                              sound_vel=locals()[name+'dict']['sound_vel'], 
                              scale=locals()[name+'dict']['scale'], 
                              sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
                              sigma_min=locals()[name+'dict']['sig_R_maj_min'],
                              star_density_max=star_density, 
                              star_density_min=star_density, 
                              data_lim=locals()[name+'dict']['data_lim'], 
                          color=color, alpha=0.3, 
                          label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$',
                          disk_scales=[(h, band)])

    ax.set_xlim(0, 130)
    if name == 'n338' or name == 'n5533':
        ax.set_ylim(0, 1.4)
    else:
        ax.set_ylim(0, 1.25)
    ax.xaxis.grid(True)
    locals()[name+'dict']['plot_SF'](ax)

    if ind == 6:
        ax.set_xlabel(r'$R, arcsec$', fontsize=20)
    ax.set_ylabel(r'$Q_{eff}^{-1}$', fontsize=20)
    for q_ in [1., 1.5, 2., 3.]:
        ax.axhline(y=1./q_, lw=10, alpha=0.15, color='grey')

    ax.text(95, 0.95, 'NGC '+name[1:]+' $'+phot[3]+'$', fontsize=30)
#     ax.set_yticks(['', '0.2', '0.4', '0.6', '0.8', '1.0', ''])
    plt.setp(ax.get_yticklabels()[0], visible=False)    
    plt.setp(ax.get_yticklabels()[-1], visible=False)    
    
#     plt.plot([-1, -2], [-1, -2], '-', color=color, label=name+' '+phot[3])
    
# plt.axvline(x=1, alpha=0.5)
# plot_Q_levels(ax, [1.5, 2., 3.])
# plt.legend()
# plt.xlim(0, 130.)
# plt.ylim(0, 1.3)
# plt.ylabel(r'$\Sigma_g/\Sigma_{cr}$', fontsize=20)
# plt.xlabel(r'$R$', fontsize=20)
plt.subplots_adjust(wspace=0, hspace=0)
plt.show()


n2985_S4G n2985
n4725_S4G n4725
n4725_H n4725
n5533_R n5533
n4258_S4G n4258
n1167 n1167
n3898_R2 n3898
n338_I n338
n2985_K n2985
n4258_I n4258
n338_B n338
n5533_r n5533
n3898_R n3898
Wall time: 2h 37min 51s
Compiler : 146 ms
Parser   : 1.99 s

In [ ]:

WS RF


In [43]:
%%time

import scipy.interpolate
def plot_WS_vs_2f(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, sfrange=None):
    
    rr = zip(*total_gas_data)[0]

    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))
    
    invQwff_WS_min = [1./Qg(epicycl=epicycl(gas_approx, r, scale), sound_vel=sound_vel, gas_density=gd) + 1./Qs(
        epicycl=epicycl(gas_approx, r, scale), sigma=sigma_max(r), star_density=star_density_min(r)) for r, gd in total_gas_data]
    
    
    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))
        
    invQwff_WS_max = [1./Qg(epicycl=epicycl(gas_approx, r, scale), sound_vel=sound_vel, gas_density=gd) + 1./Qs(
        epicycl=epicycl(gas_approx, r, scale), sigma=sigma_min(r), star_density=star_density_max(r)) for r, gd in total_gas_data]
    
    ax.plot(invQwff_WS_min, invQeff_min, '^', color=color, alpha=0.4, ms=8, lw=1)
    ax.plot(invQwff_WS_max, invQeff_max, 'D', color=color, alpha=0.4, ms=8, lw=1)
    
    return invQwff_WS_min, invQwff_WS_max
    
    
def plot_RF13_vs_2F(ax=None, r_g_dens=None, HI_gas_dens=None, CO_gas_dens=None, epicycl=None, sound_vel=None, sigma_R_max=None, sigma_R_min=None,  
           star_density=None, alpha_max=None, alpha_min=None, scale=None, gas_approx=None, thin=True, show=False, color=None):
    '''Плотности газа передаются не скорр. за гелий.'''
    
    if thin:
        romeo_Q = romeo_Qinv_thin
    else:
        romeo_Q = romeo_Qinv
            
    totgas = zip(r_g_dens, [He_coeff*(l[0]+l[1]) for l in zip(HI_gas_dens, CO_gas_dens)])
    
    if show:
        print 'sig_R_max case:'
    romeo_min = []
    for r, g, co in zip(r_g_dens, HI_gas_dens, CO_gas_dens):
        rom, _ = romeo_Q(r=r, epicycl=epicycl, sound_vel=sound_vel, sigma_R=sigma_R_max, 
               star_density=star_density, HI_density=He_coeff*g, CO_density=He_coeff*co, 
                         alpha=alpha_min, scale=scale, gas_approx=gas_approx)
        romeo_min.append(rom)



    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 show:
        print 'sig_R_min case:'
    
    romeo_max = []
    for r, g, co in zip(r_g_dens, HI_gas_dens, CO_gas_dens):
        rom, _ = romeo_Q(r=r, epicycl=epicycl, sound_vel=sound_vel, sigma_R=sigma_R_min, 
               star_density=star_density, HI_density=He_coeff*g, CO_density=He_coeff*co, 
                         alpha=alpha_max, scale=scale, gas_approx=gas_approx)
        romeo_max.append(rom)



    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(romeo_min, invQeff_min, 's', color=color, alpha=0.4, ms=8)
    ax.plot(romeo_max, invQeff_max, 'o', color=color, alpha=0.4, ms=8)
    
    return invQeff_min, invQeff_max, romeo_min, romeo_max
    

markers = {'.': 'point', ',': 'pixel', 'o': 'circle', 'v': 'triangle_down', '^': 'triangle_up', '<': 'triangle_left', '>': 'triangle_right', '1': 'tri_down', '2': 
           'tri_up', '3': 'tri_left', '4': 'tri_right', '8': 'octagon', 's': 'square', 'p': 'pentagon', '*': 'star', 'h': 'hexagon1', 'H': 'hexagon2', '+': 'plus', 
           'x': 'x', 'D': 'diamond', 'd': 'thin_diamond', '|': 'vline', '_': 'hline', 'P': 'plus_filled', 'X': 'x_filled', 0: 'tickleft', 1: 'tickright', 2: 'tickup', 
           3: 'tickdown', 4: 'caretleft', 5: 'caretright', 6: 'caretup', 7: 'caretdown', 8: 'caretleftbase', 9: 'caretrightbase', 10: 'caretupbase', 11: 'caretdownbase'}

kenn_photom={
             'n338_I' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60., 15.],
             'n338_B' : ['mu0d_Bc', 'h_disc_B', 6.99, 'B', 60., 15.],
            'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40., 6.7],
            'n2985_K' : ['mudK', 'hK', 1.42, 'K', 70., 14.],
            'n2985_S4G' : [('mu0d_s4g', 'mu0d_s4g_2'), ('h_disc_s4g', 'h_disc_s4g_2'), 1.35, 'S4G', 70., 6.3],
            'n3898_R' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80., 8.8],
            'n3898_R2' : [('mu_inner_approx', 'mu_out'), ('h_approx', 'h_out'), 2.93, 'R', 80., 8.8],
            'n4258_I' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150., 15.],
            'n4258_S4G' : ['mu0d_36', 'h_disc_36', 0.38, 'S4G', 150., 15.],
            'n4725_H' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55., 10.14],
            'n4725_S4G' : ['mu0d_s4g', 'h_disc_s4g', 1.61, 'S4G', 55., 10.14],
            'n5533_r' : ['MU0_r', 'hi_r', 4.41, 'r', 100., 8.9],
            'n5533_R' : ['mu0d_Rc', 'h_disc_R', 7.15, 'R', 100., 9.9]
            }
    

fig, (ax, ax2) = plt.subplots(ncols=2, nrows=1, figsize=[16, 8])

color_counter = 0
color_shuffle = range(13)
shuffle(color_shuffle)

saved_data_ = []

for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
# for ind, name in enumerate(['n4258']): 
    for key in kenn_photom.keys():
        
        if name in key:
            print key
            phot = kenn_photom[key]
            color = cm.rainbow(np.linspace(0, 1, 13))[color_shuffle[color_counter]]
            try:
                h = locals()[name+'dict'][phot[1][0]]
            except KeyError:
                h = locals()[name+'dict'][phot[1]]
        #     print name
        #     print mud
        #     print h

            if name in ['n4258', 'n4725']:
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
            elif name == 'n2985':
                total_gas_data=locals()[name+'dict']['total_gas_data_']
            else:
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in 
                                                                           zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]


            while total_gas_data[0][0] < phot[5]:
                total_gas_data = total_gas_data[1:]
                
            while total_gas_data[0][-1] > 130.:
                total_gas_data = total_gas_data[:-1]

#             print total_gas_data
        
            ML = phot[2]
            band = phot[3]
        
            if type(phot[0]) == tuple and phot[3] == 'S4G':
                h1 = globals()[name+'dict'][phot[1][0]]
                h2 = globals()[name+'dict'][phot[1][1]]
                mu01 = globals()[name+'dict'][phot[0][0]]
                mu02 = globals()[name+'dict'][phot[0][1]]
                star_density = lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu01, h=h1), ML),
                               lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu02, h=h2), ML)))(l) 
            elif type(phot[0]) == tuple:
                h1 = globals()[name+'dict'][phot[1][0]]
                h2 = globals()[name+'dict'][phot[1][1]]
                mu01 = globals()[name+'dict'][phot[0][0]]
                mu02 = globals()[name+'dict'][phot[0][1]]
                star_density = lambda l: tot_dens((lambda l1: surf_density(mu=mu_disc(l1, mu0=mu01, h=h1), M_to_L=ML, band=band),
                               lambda l2: surf_density(mu=mu_disc(l2, mu0=mu02, h=h2), M_to_L=ML, band=band)))(l)
            elif phot[3] == 'S4G':
                star_density = lambda l: s4g_surf_density(mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), ML)
            else:
                star_density = lambda l: surf_density(mu=mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), M_to_L=ML, band=band)
            
            
            if name == 'n338':
                epicycl = locals()[name+'dict']['epicyclicFreq_real_']
            else:
                epicycl = locals()[name+'dict']['epicyclicFreq_real']   
            
            invQwff_WS_min, invQwff_WS_max = plot_WS_vs_2f(ax=ax, total_gas_data=total_gas_data,
                              epicycl=epicycl,
                              gas_approx=locals()[name+'dict']['spl_gas'], 
                              sound_vel=locals()[name+'dict']['sound_vel'], 
                              scale=locals()[name+'dict']['scale'], 
                              sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
                              sigma_min=locals()[name+'dict']['sig_R_maj_min'],
                              star_density_max=star_density, 
                              star_density_min=star_density, 
                              data_lim=locals()[name+'dict']['data_lim'], 
                          color=color, alpha=0.5, 
                          label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$',
                          sfrange=phot[4])
            ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$')
            
            if name in ['n4258', 'n4725']:
                
                r_mol_dens = locals()[name+'dict']['r_mol_dens']
                mol_dens = locals()[name+'dict']['mol_dens']
                y_interp = scipy.interpolate.interp1d(list(r_mol_dens), list(mol_dens))

                def y_interp_(r):
                    if r <= min(r_mol_dens):
                        return y_interp(min(r_mol_dens))
                    elif r < max(r_mol_dens):
                        return y_interp(r)
                    else:
                        return 0.
                
                # не совсем точно так, но все же
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [l[1]/He_coeff - y_interp_(l[0]) for l in zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])], 
                                   [y_interp_(l) for l in locals()[name+'dict']['r_g_dens']])[1:15]
                
            elif name == 'n2985':
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'], 
                                   [locals()[name+'dict']['y_interp_'](l) for l in locals()[name+'dict']['r_g_dens']])[1:15]
                
            else:
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'], 
                                   [locals()[name+'dict']['y_interp_'](l, h) for l in locals()[name+'dict']['r_g_dens']])[1:15]


            while total_gas_data[0][0] < phot[5]:
                total_gas_data = total_gas_data[1:]
                
            while total_gas_data[0][-1] > 130.:
                total_gas_data = total_gas_data[:-1]
                
            invQeff_min, invQeff_max, romeo_min, romeo_max = plot_RF13_vs_2F(ax=ax, r_g_dens=zip(*total_gas_data)[0], 
                            HI_gas_dens=zip(*total_gas_data)[1], 
                            CO_gas_dens=zip(*total_gas_data)[2], 
                            epicycl=epicycl, 
                            sound_vel=6., 
                            sigma_R_max=locals()[name+'dict']['sig_R_maj_max'], 
                            sigma_R_min=locals()[name+'dict']['sig_R_maj_min'],  
                            star_density=star_density, 
                            alpha_max=0.7, 
                            alpha_min=0.3, 
                            scale=locals()[name+'dict']['scale'], 
                            gas_approx=locals()[name+'dict']['spl_gas'], 
                            thin=True,
                            color=color)
            
            saved_data_.append(zip(zip(*total_gas_data)[0], invQeff_min, invQeff_max, romeo_min, romeo_max, invQwff_WS_min, invQwff_WS_max))
            
            ax2.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$')
            color_counter+=1

            
def label_line(line, label, x, y, color='0.5', size=12):
    """Add a label to a line, at the proper angle.

    Arguments
    ---------
    line : matplotlib.lines.Line2D object,
    label : str
    x : float
        x-position to place center of text (in data coordinated
    y : float
        y-position to place center of text (in data coordinates)
    color : str
    size : float
    """
    xdata, ydata = line.get_data()
    x1 = xdata[0]
    x2 = xdata[-1]
    y1 = ydata[0]
    y2 = ydata[-1]

    ax = line.get_axes()
    text = ax.annotate(label, xy=(x, y), xytext=(-10, 0),
                       textcoords='offset points',
                       size=size, color=color,
                       horizontalalignment='left',
                       verticalalignment='bottom')

    sp1 = ax.transData.transform_point((x1, y1))
    sp2 = ax.transData.transform_point((x2, y2))

    rise = (sp2[1] - sp1[1])
    run = (sp2[0] - sp1[0])

    slope_degrees = np.degrees(np.arctan2(rise, run))
    text.set_rotation(slope_degrees)
    return text

ax.legend(loc='lower right')
ax.set_xlim(0, 1.75)
ax.set_ylim(0, 1.05)
ax.set_ylabel(r'$Q_{2F}^{-1}$', fontsize=20)
ax.set_xlabel(r'$Q_{WS,RF}^{-1}$', fontsize=20)
ax.grid()

ax2.legend(loc='lower right')
ax2.set_xlim(0, 1.05)
ax2.set_ylim(0, 1.05)
ax2.set_ylabel(r'$Q_{eff}^{-1}$', fontsize=20)
ax2.set_xlabel(r'$Q_{RF}^{-1}$', fontsize=20)
ax2.grid()

for _ in np.linspace(1., 2., 5):
       
    line, = ax.plot([0., _*2], [0., 2.0], '--', alpha=0.3, color='gray')
    label_line(line, 'y={:2.2f}x'.format(1./_), 0.8*_, 0.8, color='black')
    
    line2, = ax.plot([0., 1.0], [0., _], '--', alpha=0.3, color='gray')
    label_line(line2, 'y={:2.2f}x'.format(_), 0.8/_, 0.8, color='black')
#     print np.arctan(1./_)*(180./np.pi)
#     ax.text(0.8*_, 0.8, 'y={}x'.format(_), rotation=np.arctan(1./_)*(180./np.pi))

ax.plot([0., 2.0], [0., 2.0], '--', alpha=0.9, color='gray', lw=3)

plt.show()


n338_I
n338_B
n1167
n2985_S4G
n2985_K
n3898_R2
n3898_R
n4258_S4G
n4258_I
n4725_S4G
n4725_H
n5533_R
n5533_r
C:\Anaconda\lib\site-packages\matplotlib\artist.py:233: MatplotlibDeprecationWarning: get_axes has been deprecated in mpl 1.5, please use the
axes property.  A removal date has not been set.
  stacklevel=1)
Wall time: 8min 18s

In [44]:
os.chdir(return_back_path)

In [45]:
np.save('2f_appr_data.npy', saved_data_)

In [46]:
approx_data = np.load('2f_appr_data.npy')

In [47]:
fig = plt.figure(figsize=[12,12])
ax = plt.gca()

color_counter = 0
color_shuffle = range(13)
shuffle(color_shuffle)

for ind2, key in enumerate(['n338_I', 'n338_B','n1167', 'n2985_S4G','n2985_K','n3898_R2','n3898_R','n4258_S4G','n4258_I','n4725_S4G','n4725_H','n5533_R','n5533_r']):
    name = key.split('_')[0]
    phot = kenn_photom[key]
    rr, invQeff_min, invQeff_max, romeo_min, romeo_max, invQwff_WS_min, invQwff_WS_max = zip(*approx_data[ind2])
    color = cm.rainbow(np.linspace(0, 1, 13))[color_shuffle[color_counter]]
    ax.plot(romeo_min, invQeff_min, 's', color=color, alpha=0.4, ms=8)
    ax.plot(romeo_max, invQeff_max, 'o', color=color, alpha=0.4, ms=8)
    ax.plot(invQwff_WS_min, invQeff_min, '^', color=color, alpha=0.4, ms=8, lw=1)
    ax.plot(invQwff_WS_max, invQeff_max, 'D', color=color, alpha=0.4, ms=8, lw=1)
    ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$', lw=5)
    color_counter += 1
        

ax.legend(loc='lower right')
ax.set_xlim(0, 1.75)
ax.set_ylim(0, 1.05)
ax.set_ylabel(r'$Q_{2F}^{-1}$', fontsize=20)
ax.set_xlabel(r'$Q_{WS,RF}^{-1}$', fontsize=20)
ax.grid()

def label_line(line, label, x, y, color='0.5', size=12):
    """Add a label to a line, at the proper angle.

    Arguments
    ---------
    line : matplotlib.lines.Line2D object,
    label : str
    x : float
        x-position to place center of text (in data coordinated
    y : float
        y-position to place center of text (in data coordinates)
    color : str
    size : float
    """
    xdata, ydata = line.get_data()
    x1 = xdata[0]
    x2 = xdata[-1]
    y1 = ydata[0]
    y2 = ydata[-1]

    ax = line.get_axes()
    text = ax.annotate(label, xy=(x, y), xytext=(-10, 0),
                       textcoords='offset points',
                       size=size, color=color,
                       horizontalalignment='left',
                       verticalalignment='bottom')

    sp1 = ax.transData.transform_point((x1, y1))
    sp2 = ax.transData.transform_point((x2, y2))

    rise = (sp2[1] - sp1[1])
    run = (sp2[0] - sp1[0])

    slope_degrees = np.degrees(np.arctan2(rise, run))
    text.set_rotation(slope_degrees)
    return text

for _ in np.linspace(1., 2., 5):
       
    line, = ax.plot([0., _*2], [0., 2.0], '--', alpha=0.3, color='gray')
    label_line(line, 'y={:2.2f}x'.format(1./_), 0.8*_, 0.8, color='black')
    
    line2, = ax.plot([0., 1.0], [0., _], '--', alpha=0.3, color='gray')
    label_line(line2, 'y={:2.2f}x'.format(_), 0.8/_, 0.8, color='black')

    
line, = ax.plot([0., 0.89*2], [0., 2.0], '--', alpha=0.3, color='gray')
label_line(line, 'y={:2.2f}x'.format(1./0.89), 0.8*0.89, 0.8, color='black')

ax.plot([0., 2.0], [0., 2.0], '--', alpha=0.9, color='gray', lw=3)

plt.show()



In [48]:
fig = plt.figure(figsize=[24,8])
ax = plt.gca()

color_counter = 0
color_shuffle = range(13)
shuffle(color_shuffle)

for ind2, key in enumerate(['n338_I', 'n338_B','n1167', 'n2985_S4G','n2985_K','n3898_R2','n3898_R','n4258_S4G','n4258_I','n4725_S4G','n4725_H','n5533_R','n5533_r']):
    name = key.split('_')[0]
    phot = kenn_photom[key]
    rr, invQeff_min, invQeff_max, romeo_min, romeo_max, invQwff_WS_min, invQwff_WS_max = zip(*approx_data[ind2])
    color = cm.rainbow(np.linspace(0, 1, 13))[color_shuffle[color_counter]]
    
    ax.plot(rr, [l[0]/l[1] for l in zip(invQeff_min, romeo_min)], 's-', color=color, alpha=0.4, ms=8)
    ax.plot(rr, [l[0]/l[1] for l in zip(invQeff_max, romeo_max)], 'o-', color=color, alpha=0.4, ms=8)
#     ax.plot(rr, romeo_min, '^-', color=color, alpha=0.4, ms=8)
#     ax.plot(rr, romeo_max, 'D-', color=color, alpha=0.4, ms=8)
#     ax.plot(invQwff_WS_min, invQeff_min, '^', color=color, alpha=0.4, ms=8, lw=1)
#     ax.plot(invQwff_WS_max, invQeff_max, 'D', color=color, alpha=0.4, ms=8, lw=1)
    ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$', lw=5)
    
    color_counter += 1
        

ax.legend(loc='lower right')
ax.set_xlim(0)
ax.set_ylim(0.9, 1.6)
# ax.set_ylabel(r'$Q_{2F}^{-1}$', fontsize=20)
# ax.set_xlabel(r'$Q_{WS,RF}^{-1}$', fontsize=20)
ax.grid()
plt.show()


Плохо видно - разнесем:


In [49]:
fig = plt.figure(figsize=[24,8])
ax = plt.gca()

color_counter = 0
color_shuffle = range(13)
shuffle(color_shuffle)

for ind2, key in enumerate(['n338_I', 'n338_B','n1167', 'n2985_S4G','n2985_K','n3898_R2','n3898_R','n4258_S4G','n4258_I','n4725_S4G','n4725_H','n5533_R','n5533_r']):
    name = key.split('_')[0]
    phot = kenn_photom[key]
    rr, invQeff_min, invQeff_max, romeo_min, romeo_max, invQwff_WS_min, invQwff_WS_max = zip(*approx_data[ind2])
    color = cm.rainbow(np.linspace(0, 1, 13))[color_shuffle[color_counter]]
    a_ = filter(lambda l: l[1] > 1.1 and l[1] < 1.2,  zip(rr, [l[0]/l[1] for l in zip(invQeff_min, romeo_min)]))
    try:
        ax.plot(zip(*a_)[0], map(lambda l: l+ind2*0.1, zip(*a_)[1]), 's-', color=color, alpha=0.4, ms=8)
    
#     ax.plot(rr, [ind2 + l[0]/l[1] for l in zip(invQeff_max, romeo_max)], 'o-', color=color, alpha=0.4, ms=8)
#     ax.plot(rr, romeo_min, '^-', color=color, alpha=0.4, ms=8)
#     ax.plot(rr, romeo_max, 'D-', color=color, alpha=0.4, ms=8)
#     ax.plot(invQwff_WS_min, invQeff_min, '^', color=color, alpha=0.4, ms=8, lw=1)
#     ax.plot(invQwff_WS_max, invQeff_max, 'D', color=color, alpha=0.4, ms=8, lw=1)
        ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$', lw=5)
    except Exception:
        print name
    
    color_counter += 1
        

ax.legend(loc='lower right')
ax.set_xlim(0)
ax.set_ylim(1.2, 2.5)
# ax.set_ylabel(r'$Q_{2F}^{-1}$', fontsize=20)
# ax.set_xlabel(r'$Q_{WS,RF}^{-1}$', fontsize=20)
ax.grid()
plt.show()


n338
n338
n5533

In [77]:
import scipy.interpolate

from math import pi

    
def plot_RF13_vs_2F(ax=None, r_g_dens=None, HI_gas_dens=None, CO_gas_dens=None, epicycl=None, sound_vel=6., sigma_R_max=None, sigma_R_min=None,  
           star_density=None, alpha_max=None, alpha_min=None, scale=None, gas_approx=None, thin=True, show=False, color=None, verbose=True):
    '''Плотности газа передаются не скорр. за гелий.'''
    
    if thin:
        romeo_Q = romeo_Qinv_thin
    else:
        romeo_Q = romeo_Qinv
            
    totgas = zip(r_g_dens, [He_coeff*(l[0]+l[1]) for l in zip(HI_gas_dens, CO_gas_dens)])
    
    if show:
        print 'sig_R_max case:'
    romeo_min = []
    for r, g, co in zip(r_g_dens, HI_gas_dens, CO_gas_dens):
        rom = romeo_Q(r=r, epicycl=epicycl, sound_vel=sound_vel, sigma_R=sigma_R_max, 
               star_density=star_density, HI_density=He_coeff*g, CO_density=He_coeff*co, 
                         alpha=alpha_min, scale=scale, gas_approx=gas_approx, verbose=verbose)
        romeo_min.append(rom)



    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 show:
        print 'sig_R_min case:'
    
    romeo_max = []
    for r, g, co in zip(r_g_dens, HI_gas_dens, CO_gas_dens):
        rom = romeo_Q(r=r, epicycl=epicycl, sound_vel=sound_vel, sigma_R=sigma_R_min, 
               star_density=star_density, HI_density=He_coeff*g, CO_density=He_coeff*co, 
                         alpha=alpha_max, scale=scale, gas_approx=gas_approx, verbose=verbose)
        romeo_max.append(rom)



    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(romeo_min, invQeff_min, 's', color=color, alpha=0.4, ms=8)
#     ax.plot(romeo_max, invQeff_max, 'o', color=color, alpha=0.4, ms=8)
    
    return invQeff_min, invQeff_max, romeo_min, romeo_max
    


kenn_photom={
             'n338_I' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60., 15.],
             'n338_B' : ['mu0d_Bc', 'h_disc_B', 6.99, 'B', 60., 15.],
            'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40., 6.7],
            'n2985_K' : ['mudK', 'hK', 1.42, 'K', 70., 14.],
            'n2985_S4G' : [('mu0d_s4g', 'mu0d_s4g_2'), ('h_disc_s4g', 'h_disc_s4g_2'), 1.35, 'S4G', 70., 6.3],
            'n3898_R' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80., 8.8],
            'n3898_R2' : [('mu_inner_approx', 'mu_out'), ('h_approx', 'h_out'), 2.93, 'R', 80., 8.8],
            'n4258_I' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150., 15.],
            'n4258_S4G' : ['mu0d_36', 'h_disc_36', 0.38, 'S4G', 150., 15.],
            'n4725_H' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55., 10.14],
            'n4725_S4G' : ['mu0d_s4g', 'h_disc_s4g', 1.61, 'S4G', 55., 10.14],
            'n5533_r' : ['MU0_r', 'hi_r', 4.41, 'r', 100., 8.9],
            'n5533_R' : ['mu0d_Rc', 'h_disc_R', 7.15, 'R', 100., 9.9]
            }
    

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

color_counter = 0
color_shuffle = range(13)
shuffle(color_shuffle)

ax = plt.gca()

for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
# for ind, name in enumerate(['n1167']): 
    for key in kenn_photom.keys():
        
        if name in key:
            print key
            phot = kenn_photom[key]
            color = cm.rainbow(np.linspace(0, 1, 13))[color_shuffle[color_counter]]
            try:
                h = locals()[name+'dict'][phot[1][0]]
            except KeyError:
                h = locals()[name+'dict'][phot[1]]
        #     print name
        #     print mud
        #     print h

            if name in ['n4258', 'n4725']:
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
            elif name == 'n2985':
                total_gas_data=locals()[name+'dict']['total_gas_data_']
            else:
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in 
                                                                           zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]


            while total_gas_data[0][0] < phot[5]:
                total_gas_data = total_gas_data[1:]
                
            while total_gas_data[0][-1] > 130.:
                total_gas_data = total_gas_data[:-1]

#             print total_gas_data
        
            ML = phot[2]
            band = phot[3]
        
            if type(phot[0]) == tuple and phot[3] == 'S4G':
                h1 = globals()[name+'dict'][phot[1][0]]
                h2 = globals()[name+'dict'][phot[1][1]]
                mu01 = globals()[name+'dict'][phot[0][0]]
                mu02 = globals()[name+'dict'][phot[0][1]]
                star_density = lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu01, h=h1), ML),
                               lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu02, h=h2), ML)))(l) 
            elif type(phot[0]) == tuple:
                h1 = globals()[name+'dict'][phot[1][0]]
                h2 = globals()[name+'dict'][phot[1][1]]
                mu01 = globals()[name+'dict'][phot[0][0]]
                mu02 = globals()[name+'dict'][phot[0][1]]
                star_density = lambda l: tot_dens((lambda l1: surf_density(mu=mu_disc(l1, mu0=mu01, h=h1), M_to_L=ML, band=band),
                               lambda l2: surf_density(mu=mu_disc(l2, mu0=mu02, h=h2), M_to_L=ML, band=band)))(l)
            elif phot[3] == 'S4G':
                star_density = lambda l: s4g_surf_density(mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), ML)
            else:
                star_density = lambda l: surf_density(mu=mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), M_to_L=ML, band=band)
            
            
            if name == 'n338':
                epicycl = locals()[name+'dict']['epicyclicFreq_real_']
            else:
                epicycl = locals()[name+'dict']['epicyclicFreq_real']   
            
#             invQwff_WS_min, invQwff_WS_max = plot_WS_vs_2f(ax=ax, total_gas_data=total_gas_data,
#                               epicycl=epicycl,
#                               gas_approx=locals()[name+'dict']['spl_gas'], 
#                               sound_vel=locals()[name+'dict']['sound_vel'], 
#                               scale=locals()[name+'dict']['scale'], 
#                               sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
#                               sigma_min=locals()[name+'dict']['sig_R_maj_min'],
#                               star_density_max=star_density, 
#                               star_density_min=star_density, 
#                               data_lim=locals()[name+'dict']['data_lim'], 
#                           color=color, alpha=0.5, 
#                           label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$',
#                           sfrange=phot[4])
#             ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$')
            
            if name in ['n4258', 'n4725']:
                
                r_mol_dens = locals()[name+'dict']['r_mol_dens']
                mol_dens = locals()[name+'dict']['mol_dens']
                y_interp = scipy.interpolate.interp1d(list(r_mol_dens), list(mol_dens))

                def y_interp_(r):
                    if r <= min(r_mol_dens):
                        return y_interp(min(r_mol_dens))
                    elif r < max(r_mol_dens):
                        return y_interp(r)
                    else:
                        return 0.
                
                # не совсем точно так, но все же
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [l[1]/He_coeff - y_interp_(l[0]) for l in zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])], 
                                   [y_interp_(l) for l in locals()[name+'dict']['r_g_dens']])[1:15]
                
            elif name == 'n2985':
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'], 
                                   [locals()[name+'dict']['y_interp_'](l) for l in locals()[name+'dict']['r_g_dens']])[1:15]
                
            else:
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'], 
                                   [locals()[name+'dict']['y_interp_'](l, h) for l in locals()[name+'dict']['r_g_dens']])[1:15]


            while total_gas_data[0][0] < phot[5]:
                total_gas_data = total_gas_data[1:]
                
            while total_gas_data[0][-1] > 130.:
                total_gas_data = total_gas_data[:-1]
                
            invQeff_min, invQeff_max, romeo_min, romeo_max = plot_RF13_vs_2F(ax=ax, r_g_dens=zip(*total_gas_data)[0], 
                            HI_gas_dens=zip(*total_gas_data)[1], 
                            CO_gas_dens=zip(*total_gas_data)[2], 
                            epicycl=epicycl, 
                            sound_vel=sound_vel,
                            sigma_R_max=locals()[name+'dict']['sig_R_maj_max'], 
                            sigma_R_min=locals()[name+'dict']['sig_R_maj_min'],  
                            star_density=star_density, 
                            alpha_max=0.7, 
                            alpha_min=0.3, 
                            scale=locals()[name+'dict']['scale'], 
                            gas_approx=locals()[name+'dict']['spl_gas'], 
                            thin=True,
                            color=color, verbose=False)
            
            romeo_min, components = zip(*romeo_min)
            a_ = filter(lambda l: l[1] > 1.1 and l[1] < 1.2,  zip(rr, [l[0]/l[1] for l in zip(invQeff_min, romeo_min)], components))
#             print a_
            try:
                for _ in a_:
                    if _[2] == 'star':
                        marker = '*' 
                    elif _[2] == 'HI':
                        marker = 'o'
                    else:
                        marker = 's'
#                     ax.plot([_[0], _[1]], [_[0]+0.01, _[1]+0.01])
                    ax.scatter(_[0], _[1]+color_counter*0.1, 30, marker=marker, color=color)
#                 ax.plot(zip(*a_)[0], zip(*a_)[1],'-', lw=10)
#                 print zip(*a_)[0]
            except Exception:
                print 'Exception for ' + name
            
            
            
            
#             ax.plot(romeo_min, romeo_min2, '.', ms=9, color=color)
#             ax.plot(romeo_max, romeo_max2, '.', ms=9, color=color)
            
            
            ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$')
            color_counter+=1

            
ax.legend(loc='lower right')
ax.set_xlim(0, 100)
ax.set_ylim(1.0, 2.5)
ax.set_ylabel(r'$Q_{RF, thick}$', fontsize=20)
ax.set_xlabel(r'$Q_{RF}^{-1}$', fontsize=20)
ax.grid()

plt.show()


n338_I
n338_B
n1167
n2985_S4G
n2985_K
n3898_R2
n3898_R
n4258_S4G
n4258_I
n4725_S4G
n4725_H
n5533_R
n5533_r

Влияние разных дисперсий H2 и HI:


In [50]:
import scipy.interpolate

from math import pi

def romeo_Qinv_thin(r=None, epicycl=None, sound_vel_CO=11., sound_vel_HI=6., sigma_R=None, star_density=None, 
               HI_density=None, CO_density=None, alpha=None, scale=None, gas_approx=None, verbose=False, show=False):
    G = 4.32
    kappa = epicycl(gas_approx, r, scale)
    Q_star = kappa*sigma_R(r)/(pi*G*star_density(r))
    Q_CO = kappa*sound_vel_CO/(pi*G*CO_density)
    Q_HI = kappa*sound_vel_HI/(pi*G*HI_density)
    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
        
    #     TODO: оставить только show или verbose
    if show:
        print 'r={:7.3f} Qg={:7.3f} Qs={:7.3f} Qg^-1={:7.3f} Qs^-1={:7.3f}'.format(r, Q_HI, Q_star, 1./Q_HI, 1./Q_star)
    
    dispersions = [sigma_R(r), sound_vel_HI, sound_vel_CO]
    QTs = [Q_star*T_star, Q_HI*T_HI, Q_CO*T_CO]
    components = ['star', 'HI', 'H2']
    index = QTs.index(min(QTs))
    
    if verbose:
        print 'QTs: {}'.format(QTs)
        print 'min index: {}'.format(index)
        print 'min component: {}'.format(components[index])
    
    sig_m = dispersions[index]
    def W_i(sig_m, sig_i):
        return 2*sig_m*sig_i/(sig_m**2 + sig_i**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[index]
    
def plot_RF13_vs_2F(ax=None, r_g_dens=None, HI_gas_dens=None, CO_gas_dens=None, epicycl=None, sound_vel_CO=11., sound_vel_HI=6., sigma_R_max=None, sigma_R_min=None,  
           star_density=None, alpha_max=None, alpha_min=None, scale=None, gas_approx=None, thin=True, show=False, color=None):
    '''Плотности газа передаются не скорр. за гелий.'''
    
    if thin:
        romeo_Q = romeo_Qinv_thin
    else:
        romeo_Q = romeo_Qinv
            
    totgas = zip(r_g_dens, [He_coeff*(l[0]+l[1]) for l in zip(HI_gas_dens, CO_gas_dens)])
    
    if show:
        print 'sig_R_max case:'
    romeo_min = []
    for r, g, co in zip(r_g_dens, HI_gas_dens, CO_gas_dens):
        rom, _ = romeo_Q(r=r, epicycl=epicycl, sound_vel_CO=sound_vel_CO, sound_vel_HI=sound_vel_HI, sigma_R=sigma_R_max, 
               star_density=star_density, HI_density=He_coeff*g, CO_density=He_coeff*co, 
                         alpha=alpha_min, scale=scale, gas_approx=gas_approx)
        romeo_min.append(rom)



    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 show:
        print 'sig_R_min case:'
    
    romeo_max = []
    for r, g, co in zip(r_g_dens, HI_gas_dens, CO_gas_dens):
        rom, _ = romeo_Q(r=r, epicycl=epicycl, sound_vel_CO=sound_vel_CO, sound_vel_HI=sound_vel_HI, sigma_R=sigma_R_min, 
               star_density=star_density, HI_density=He_coeff*g, CO_density=He_coeff*co, 
                         alpha=alpha_max, scale=scale, gas_approx=gas_approx)
        romeo_max.append(rom)



    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(romeo_min, invQeff_min, 's', color=color, alpha=0.4, ms=8)
#     ax.plot(romeo_max, invQeff_max, 'o', color=color, alpha=0.4, ms=8)
    
    return invQeff_min, invQeff_max, romeo_min, romeo_max
    


kenn_photom={
             'n338_I' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60., 15.],
             'n338_B' : ['mu0d_Bc', 'h_disc_B', 6.99, 'B', 60., 15.],
            'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40., 6.7],
            'n2985_K' : ['mudK', 'hK', 1.42, 'K', 70., 14.],
            'n2985_S4G' : [('mu0d_s4g', 'mu0d_s4g_2'), ('h_disc_s4g', 'h_disc_s4g_2'), 1.35, 'S4G', 70., 6.3],
            'n3898_R' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80., 8.8],
            'n3898_R2' : [('mu_inner_approx', 'mu_out'), ('h_approx', 'h_out'), 2.93, 'R', 80., 8.8],
            'n4258_I' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150., 15.],
            'n4258_S4G' : ['mu0d_36', 'h_disc_36', 0.38, 'S4G', 150., 15.],
            'n4725_H' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55., 10.14],
            'n4725_S4G' : ['mu0d_s4g', 'h_disc_s4g', 1.61, 'S4G', 55., 10.14],
            'n5533_r' : ['MU0_r', 'hi_r', 4.41, 'r', 100., 8.9],
            'n5533_R' : ['mu0d_Rc', 'h_disc_R', 7.15, 'R', 100., 9.9]
            }
    

fig  = plt.figure(figsize=[8, 8])

color_counter = 0
color_shuffle = range(13)
shuffle(color_shuffle)

ax = plt.gca()

for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
    for key in kenn_photom.keys():
        
        if name in key:
            print key
            phot = kenn_photom[key]
            color = cm.rainbow(np.linspace(0, 1, 13))[color_shuffle[color_counter]]
            try:
                h = locals()[name+'dict'][phot[1][0]]
            except KeyError:
                h = locals()[name+'dict'][phot[1]]
        #     print name
        #     print mud
        #     print h

            if name in ['n4258', 'n4725']:
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
            elif name == 'n2985':
                total_gas_data=locals()[name+'dict']['total_gas_data_']
            else:
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in 
                                                                           zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]


            while total_gas_data[0][0] < phot[5]:
                total_gas_data = total_gas_data[1:]
                
            while total_gas_data[0][-1] > 130.:
                total_gas_data = total_gas_data[:-1]

#             print total_gas_data
        
            ML = phot[2]
            band = phot[3]
        
            if type(phot[0]) == tuple and phot[3] == 'S4G':
                h1 = globals()[name+'dict'][phot[1][0]]
                h2 = globals()[name+'dict'][phot[1][1]]
                mu01 = globals()[name+'dict'][phot[0][0]]
                mu02 = globals()[name+'dict'][phot[0][1]]
                star_density = lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu01, h=h1), ML),
                               lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu02, h=h2), ML)))(l) 
            elif type(phot[0]) == tuple:
                h1 = globals()[name+'dict'][phot[1][0]]
                h2 = globals()[name+'dict'][phot[1][1]]
                mu01 = globals()[name+'dict'][phot[0][0]]
                mu02 = globals()[name+'dict'][phot[0][1]]
                star_density = lambda l: tot_dens((lambda l1: surf_density(mu=mu_disc(l1, mu0=mu01, h=h1), M_to_L=ML, band=band),
                               lambda l2: surf_density(mu=mu_disc(l2, mu0=mu02, h=h2), M_to_L=ML, band=band)))(l)
            elif phot[3] == 'S4G':
                star_density = lambda l: s4g_surf_density(mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), ML)
            else:
                star_density = lambda l: surf_density(mu=mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), M_to_L=ML, band=band)
            
            
            if name == 'n338':
                epicycl = locals()[name+'dict']['epicyclicFreq_real_']
            else:
                epicycl = locals()[name+'dict']['epicyclicFreq_real']   
            
#             invQwff_WS_min, invQwff_WS_max = plot_WS_vs_2f(ax=ax, total_gas_data=total_gas_data,
#                               epicycl=epicycl,
#                               gas_approx=locals()[name+'dict']['spl_gas'], 
#                               sound_vel=locals()[name+'dict']['sound_vel'], 
#                               scale=locals()[name+'dict']['scale'], 
#                               sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
#                               sigma_min=locals()[name+'dict']['sig_R_maj_min'],
#                               star_density_max=star_density, 
#                               star_density_min=star_density, 
#                               data_lim=locals()[name+'dict']['data_lim'], 
#                           color=color, alpha=0.5, 
#                           label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$',
#                           sfrange=phot[4])
#             ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$')
            
            if name in ['n4258', 'n4725']:
                
                r_mol_dens = locals()[name+'dict']['r_mol_dens']
                mol_dens = locals()[name+'dict']['mol_dens']
                y_interp = scipy.interpolate.interp1d(list(r_mol_dens), list(mol_dens))

                def y_interp_(r):
                    if r <= min(r_mol_dens):
                        return y_interp(min(r_mol_dens))
                    elif r < max(r_mol_dens):
                        return y_interp(r)
                    else:
                        return 0.
                
                # не совсем точно так, но все же
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [l[1]/He_coeff - y_interp_(l[0]) for l in zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])], 
                                   [y_interp_(l) for l in locals()[name+'dict']['r_g_dens']])[1:15]
                
            elif name == 'n2985':
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'], 
                                   [locals()[name+'dict']['y_interp_'](l) for l in locals()[name+'dict']['r_g_dens']])[1:15]
                
            else:
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'], 
                                   [locals()[name+'dict']['y_interp_'](l, h) for l in locals()[name+'dict']['r_g_dens']])[1:15]


            while total_gas_data[0][0] < phot[5]:
                total_gas_data = total_gas_data[1:]
                
            while total_gas_data[0][-1] > 130.:
                total_gas_data = total_gas_data[:-1]
                
            invQeff_min, invQeff_max, romeo_min, romeo_max = plot_RF13_vs_2F(ax=ax, r_g_dens=zip(*total_gas_data)[0], 
                            HI_gas_dens=zip(*total_gas_data)[1], 
                            CO_gas_dens=zip(*total_gas_data)[2], 
                            epicycl=epicycl, 
                            sound_vel_HI=6.,
                            sound_vel_CO=6., 
                            sigma_R_max=locals()[name+'dict']['sig_R_maj_max'], 
                            sigma_R_min=locals()[name+'dict']['sig_R_maj_min'],  
                            star_density=star_density, 
                            alpha_max=0.7, 
                            alpha_min=0.3, 
                            scale=locals()[name+'dict']['scale'], 
                            gas_approx=locals()[name+'dict']['spl_gas'], 
                            thin=True,
                            color=color)
            
            invQeff_min, invQeff_max, romeo_min2, romeo_max2 = plot_RF13_vs_2F(ax=ax, r_g_dens=zip(*total_gas_data)[0], 
                            HI_gas_dens=zip(*total_gas_data)[1], 
                            CO_gas_dens=zip(*total_gas_data)[2], 
                            epicycl=epicycl, 
                            sound_vel_HI=6.,
                            sound_vel_CO=11., 
                            sigma_R_max=locals()[name+'dict']['sig_R_maj_max'], 
                            sigma_R_min=locals()[name+'dict']['sig_R_maj_min'],  
                            star_density=star_density, 
                            alpha_max=0.7, 
                            alpha_min=0.3, 
                            scale=locals()[name+'dict']['scale'], 
                            gas_approx=locals()[name+'dict']['spl_gas'], 
                            thin=True,
                            color=color)
            
            ax.plot(romeo_min, romeo_min2, 's', ms=9, color=color)
            ax.plot(romeo_max, romeo_max2, 'o', ms=9, color=color)
            
            ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$')
            color_counter+=1

            
ax.legend(loc='lower right')
ax.set_xlim(0, 1.05)
ax.set_ylim(0, 1.05)
ax.set_ylabel(r'$Q_{RF, \sigma_{CO}=11}^{-1}$', fontsize=20)
ax.set_xlabel(r'$Q_{RF}^{-1}$', fontsize=20)
ax.grid()

line, = ax.plot([0., 1*2], [0., 2.0], '--', alpha=0.3, color='gray')
label_line(line, 'y=x', 0.8*1, 0.8, color='black')
plt.show()


n338_I
n338_B
n1167
n2985_S4G
n2985_K
n3898_R2
n3898_R
n4258_S4G
n4258_I
n4725_S4G
n4725_H
n5533_R
n5533_r

Толщины:


In [55]:
import scipy.interpolate

from math import pi

def romeo_Qinv_thin(r=None, epicycl=None, sound_vel=11., sigma_R=None, star_density=None, 
               HI_density=None, CO_density=None, alpha=None, scale=None, gas_approx=None, verbose=False, show=False):
    G = 4.32
    kappa = epicycl(gas_approx, r, scale)
    Q_star = kappa*sigma_R(r)/(pi*G*star_density(r))
    Q_CO = kappa*sound_vel/(pi*G*CO_density)
    Q_HI = kappa*sound_vel/(pi*G*HI_density)
    
#     TODO: оставить только show или verbose
    if show:
        print 'r={:7.3f} Qg={:7.3f} Qs={:7.3f} Qg^-1={:7.3f} Qs^-1={:7.3f}'.format(r, Q_HI, Q_star, 1./Q_HI, 1./Q_star)
        
    dispersions = [sigma_R(r), sound_vel, sound_vel]
    QTs = [Q_star, Q_HI, Q_CO]
    components = ['star', 'HI', 'H2']
    index = QTs.index(min(QTs))
    
    if verbose:
        print 'QTs: {}'.format(QTs)
        print 'min index: {}'.format(index)
        print 'min component: {}'.format(components[index])
    
    sig_m = dispersions[index]
    def W_i(sig_m, sig_i):
        return 2*sig_m*sig_i/(sig_m**2 + sig_i**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[index]


def romeo_Qinv(r=None, epicycl=None, sound_vel=11., sigma_R=None, star_density=None, 
               HI_density=None, CO_density=None, alpha=None, scale=None, gas_approx=None, verbose=False, show=False):
    G = 4.32
    kappa = epicycl(gas_approx, r, scale)
    Q_star = kappa*sigma_R(r)/(pi*G*star_density(r))
    Q_CO = kappa*sound_vel/(pi*G*CO_density)
    Q_HI = kappa*sound_vel/(pi*G*HI_density)
    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
        
    #     TODO: оставить только show или verbose
    if show:
        print 'r={:7.3f} Qg={:7.3f} Qs={:7.3f} Qg^-1={:7.3f} Qs^-1={:7.3f}'.format(r, Q_HI, Q_star, 1./Q_HI, 1./Q_star)
    
    dispersions = [sigma_R(r), sound_vel, sound_vel]
    QTs = [Q_star*T_star, Q_HI*T_HI, Q_CO*T_CO]
    components = ['star', 'HI', 'H2']
    index = QTs.index(min(QTs))
    
    if verbose:
        print 'QTs: {}'.format(QTs)
        print 'min index: {}'.format(index)
        print 'min component: {}'.format(components[index])
    
    sig_m = dispersions[index]
    def W_i(sig_m, sig_i):
        return 2*sig_m*sig_i/(sig_m**2 + sig_i**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[index]

    
def plot_RF13_vs_2F(ax=None, r_g_dens=None, HI_gas_dens=None, CO_gas_dens=None, epicycl=None, sound_vel=6., sigma_R_max=None, sigma_R_min=None,  
           star_density=None, alpha_max=None, alpha_min=None, scale=None, gas_approx=None, thin=True, show=False, color=None, verbose=True):
    '''Плотности газа передаются не скорр. за гелий.'''
    
    if thin:
        romeo_Q = romeo_Qinv_thin
    else:
        romeo_Q = romeo_Qinv
            
    totgas = zip(r_g_dens, [He_coeff*(l[0]+l[1]) for l in zip(HI_gas_dens, CO_gas_dens)])
    
    if show:
        print 'sig_R_max case:'
    romeo_min = []
    for r, g, co in zip(r_g_dens, HI_gas_dens, CO_gas_dens):
        rom, _ = romeo_Q(r=r, epicycl=epicycl, sound_vel=sound_vel, sigma_R=sigma_R_max, 
               star_density=star_density, HI_density=He_coeff*g, CO_density=He_coeff*co, 
                         alpha=alpha_min, scale=scale, gas_approx=gas_approx, verbose=verbose)
        romeo_min.append(rom)



    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 show:
        print 'sig_R_min case:'
    
    romeo_max = []
    for r, g, co in zip(r_g_dens, HI_gas_dens, CO_gas_dens):
        rom, _ = romeo_Q(r=r, epicycl=epicycl, sound_vel=sound_vel, sigma_R=sigma_R_min, 
               star_density=star_density, HI_density=He_coeff*g, CO_density=He_coeff*co, 
                         alpha=alpha_max, scale=scale, gas_approx=gas_approx, verbose=verbose)
        romeo_max.append(rom)



    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(romeo_min, invQeff_min, 's', color=color, alpha=0.4, ms=8)
#     ax.plot(romeo_max, invQeff_max, 'o', color=color, alpha=0.4, ms=8)
    
    return invQeff_min, invQeff_max, romeo_min, romeo_max
    


kenn_photom={
             'n338_I' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60., 15.],
             'n338_B' : ['mu0d_Bc', 'h_disc_B', 6.99, 'B', 60., 15.],
            'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40., 6.7],
            'n2985_K' : ['mudK', 'hK', 1.42, 'K', 70., 14.],
            'n2985_S4G' : [('mu0d_s4g', 'mu0d_s4g_2'), ('h_disc_s4g', 'h_disc_s4g_2'), 1.35, 'S4G', 70., 6.3],
            'n3898_R' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80., 8.8],
            'n3898_R2' : [('mu_inner_approx', 'mu_out'), ('h_approx', 'h_out'), 2.93, 'R', 80., 8.8],
            'n4258_I' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150., 15.],
            'n4258_S4G' : ['mu0d_36', 'h_disc_36', 0.38, 'S4G', 150., 15.],
            'n4725_H' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55., 10.14],
            'n4725_S4G' : ['mu0d_s4g', 'h_disc_s4g', 1.61, 'S4G', 55., 10.14],
            'n5533_r' : ['MU0_r', 'hi_r', 4.41, 'r', 100., 8.9],
            'n5533_R' : ['mu0d_Rc', 'h_disc_R', 7.15, 'R', 100., 9.9]
            }
    

fig  = plt.figure(figsize=[8, 8])

color_counter = 0
color_shuffle = range(13)
shuffle(color_shuffle)

ax = plt.gca()

for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
# for ind, name in enumerate(['n338']): 
    for key in kenn_photom.keys():
        
        if name in key:
            print key
            phot = kenn_photom[key]
            color = cm.rainbow(np.linspace(0, 1, 13))[color_shuffle[color_counter]]
            try:
                h = locals()[name+'dict'][phot[1][0]]
            except KeyError:
                h = locals()[name+'dict'][phot[1]]
        #     print name
        #     print mud
        #     print h

            if name in ['n4258', 'n4725']:
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
            elif name == 'n2985':
                total_gas_data=locals()[name+'dict']['total_gas_data_']
            else:
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in 
                                                                           zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]


            while total_gas_data[0][0] < phot[5]:
                total_gas_data = total_gas_data[1:]
                
            while total_gas_data[0][-1] > 130.:
                total_gas_data = total_gas_data[:-1]

#             print total_gas_data
        
            ML = phot[2]
            band = phot[3]
        
            if type(phot[0]) == tuple and phot[3] == 'S4G':
                h1 = globals()[name+'dict'][phot[1][0]]
                h2 = globals()[name+'dict'][phot[1][1]]
                mu01 = globals()[name+'dict'][phot[0][0]]
                mu02 = globals()[name+'dict'][phot[0][1]]
                star_density = lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu01, h=h1), ML),
                               lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu02, h=h2), ML)))(l) 
            elif type(phot[0]) == tuple:
                h1 = globals()[name+'dict'][phot[1][0]]
                h2 = globals()[name+'dict'][phot[1][1]]
                mu01 = globals()[name+'dict'][phot[0][0]]
                mu02 = globals()[name+'dict'][phot[0][1]]
                star_density = lambda l: tot_dens((lambda l1: surf_density(mu=mu_disc(l1, mu0=mu01, h=h1), M_to_L=ML, band=band),
                               lambda l2: surf_density(mu=mu_disc(l2, mu0=mu02, h=h2), M_to_L=ML, band=band)))(l)
            elif phot[3] == 'S4G':
                star_density = lambda l: s4g_surf_density(mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), ML)
            else:
                star_density = lambda l: surf_density(mu=mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), M_to_L=ML, band=band)
            
            
            if name == 'n338':
                epicycl = locals()[name+'dict']['epicyclicFreq_real_']
            else:
                epicycl = locals()[name+'dict']['epicyclicFreq_real']   
            
#             invQwff_WS_min, invQwff_WS_max = plot_WS_vs_2f(ax=ax, total_gas_data=total_gas_data,
#                               epicycl=epicycl,
#                               gas_approx=locals()[name+'dict']['spl_gas'], 
#                               sound_vel=locals()[name+'dict']['sound_vel'], 
#                               scale=locals()[name+'dict']['scale'], 
#                               sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
#                               sigma_min=locals()[name+'dict']['sig_R_maj_min'],
#                               star_density_max=star_density, 
#                               star_density_min=star_density, 
#                               data_lim=locals()[name+'dict']['data_lim'], 
#                           color=color, alpha=0.5, 
#                           label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$',
#                           sfrange=phot[4])
#             ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$')
            
            if name in ['n4258', 'n4725']:
                
                r_mol_dens = locals()[name+'dict']['r_mol_dens']
                mol_dens = locals()[name+'dict']['mol_dens']
                y_interp = scipy.interpolate.interp1d(list(r_mol_dens), list(mol_dens))

                def y_interp_(r):
                    if r <= min(r_mol_dens):
                        return y_interp(min(r_mol_dens))
                    elif r < max(r_mol_dens):
                        return y_interp(r)
                    else:
                        return 0.
                
                # не совсем точно так, но все же
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [l[1]/He_coeff - y_interp_(l[0]) for l in zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])], 
                                   [y_interp_(l) for l in locals()[name+'dict']['r_g_dens']])[1:15]
                
            elif name == 'n2985':
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'], 
                                   [locals()[name+'dict']['y_interp_'](l) for l in locals()[name+'dict']['r_g_dens']])[1:15]
                
            else:
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'], 
                                   [locals()[name+'dict']['y_interp_'](l, h) for l in locals()[name+'dict']['r_g_dens']])[1:15]


            while total_gas_data[0][0] < phot[5]:
                total_gas_data = total_gas_data[1:]
                
            while total_gas_data[0][-1] > 130.:
                total_gas_data = total_gas_data[:-1]
                
            invQeff_min, invQeff_max, romeo_min, romeo_max = plot_RF13_vs_2F(ax=ax, r_g_dens=zip(*total_gas_data)[0], 
                            HI_gas_dens=zip(*total_gas_data)[1], 
                            CO_gas_dens=zip(*total_gas_data)[2], 
                            epicycl=epicycl, 
                            sound_vel=sound_vel,
                            sigma_R_max=locals()[name+'dict']['sig_R_maj_max'], 
                            sigma_R_min=locals()[name+'dict']['sig_R_maj_min'],  
                            star_density=star_density, 
                            alpha_max=0.7, 
                            alpha_min=0.3, 
                            scale=locals()[name+'dict']['scale'], 
                            gas_approx=locals()[name+'dict']['spl_gas'], 
                            thin=True,
                            color=color, verbose=False)
            
            invQeff_min, invQeff_max, romeo_min2, romeo_max2 = plot_RF13_vs_2F(ax=ax, r_g_dens=zip(*total_gas_data)[0], 
                            HI_gas_dens=zip(*total_gas_data)[1], 
                            CO_gas_dens=zip(*total_gas_data)[2], 
                            epicycl=epicycl, 
                            sound_vel=sound_vel,
                            sigma_R_max=locals()[name+'dict']['sig_R_maj_max'], 
                            sigma_R_min=locals()[name+'dict']['sig_R_maj_min'],  
                            star_density=star_density, 
                            alpha_max=0.7, 
                            alpha_min=0.3, 
                            scale=locals()[name+'dict']['scale'], 
                            gas_approx=locals()[name+'dict']['spl_gas'], 
                            thin=False,
                            color=color, verbose=False)
            
            ax.plot(romeo_min, romeo_min2, '.', ms=9, color=color)
            ax.plot(romeo_max, romeo_max2, '.', ms=9, color=color)
            
            ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$')
            color_counter+=1

            
ax.legend(loc='lower right')
ax.set_xlim(0, 1.05)
ax.set_ylim(0, 1.05)
ax.set_ylabel(r'$Q_{RF, thick}$', fontsize=20)
ax.set_xlabel(r'$Q_{RF}^{-1}$', fontsize=20)
ax.grid()

line, = ax.plot([0., 1*2], [0., 2.0], '--', alpha=0.3, color='gray')
label_line(line, 'y=x', 0.8*1, 0.8, color='black')
plt.show()


n338_I
n338_B
n1167
n2985_S4G
n2985_K
n3898_R2
n3898_R
n4258_S4G
n4258_I
n4725_S4G
n4725_H
n5533_R
n5533_r

Хорошо видно проблему, которая была с газом:


In [69]:
%%time
import scipy.interpolate
def plot_WS_vs_2f(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, sfrange=None):
    
    rr = zip(*total_gas_data)[0]
    
    print total_gas_data

    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))
#     print invQeff_min
    
    invQwff_WS_min = [1./Qg(epicycl=epicycl(gas_approx, r, scale), sound_vel=sound_vel, gas_density=gd) + 1./Qs(
        epicycl=epicycl(gas_approx, r, scale), sigma=sigma_max(r), star_density=star_density_min(r)) for r, gd in total_gas_data]
    
    ax.plot(rr, invQeff_min, 'o', color=color, alpha=0.8, ms=10)
#     ax.plot(rr, invQwff_WS_min, 's', color=color, alpha=0.8, ms=10)
    

    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))
    
#     print invQeff_max
        
    invQwff_WS_max = [1./Qg(epicycl=epicycl(gas_approx, r, scale), sound_vel=sound_vel, gas_density=gd) + 1./Qs(
        epicycl=epicycl(gas_approx, r, scale), sigma=sigma_min(r), star_density=star_density_max(r)) for r, gd in total_gas_data]

#     ax.plot(invQwff_WS_max, invQeff_max, 'o', color=color, alpha=0.8, ms=10)
    
    
def plot_RF13_vs_2F(ax=None, r_g_dens=None, HI_gas_dens=None, CO_gas_dens=None, epicycl=None, sound_vel=None, sigma_R_max=None, sigma_R_min=None,  
           star_density=None, alpha_max=None, alpha_min=None, scale=None, gas_approx=None, thin=True, show=False, color=None):
    '''Плотности газа передаются не скорр. за гелий.'''
    
    if thin:
        romeo_Q = romeo_Qinv_thin
    else:
        romeo_Q = romeo_Qinv
            
    totgas = zip(r_g_dens, [He_coeff*(l[0]+l[1]) for l in zip(HI_gas_dens, CO_gas_dens)])
    print totgas
    
    if show:
        print 'sig_R_max case:'
    romeo_min = []
    for r, g, co in zip(r_g_dens, HI_gas_dens, CO_gas_dens):
        rom, _ = romeo_Q(r=r, epicycl=epicycl, sound_vel=sound_vel, sigma_R=sigma_R_max, 
               star_density=star_density, HI_density=He_coeff*g, CO_density=He_coeff*co, 
                         alpha=alpha_min, scale=scale, gas_approx=gas_approx)
        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))
    
#     print invQeff_min

    if show:
        print 'sig_R_min case:'
    
    romeo_max = []
    for r, g, co in zip(r_g_dens, HI_gas_dens, CO_gas_dens):
        rom, _ = romeo_Q(r=r, epicycl=epicycl, sound_vel=sound_vel, sigma_R=sigma_R_min, 
               star_density=star_density, HI_density=He_coeff*g, CO_density=He_coeff*co, 
                         alpha=alpha_max, scale=scale, gas_approx=gas_approx)
        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))
    
#     print invQeff_max

    
#     ax.plot(r_g_dens, romeo_min, 's', color=color, alpha=0.8, ms=10)
    ax.plot(r_g_dens, invQeff_min, 'o', color=color, alpha=0.8, ms=10)
    ax.plot(r_g_dens, invQg, '-', color='b', alpha=0.8)
    
#     ax.plot(romeo_min, invQeff_min, 's', color=color, alpha=0.8, ms=10)
#     ax.plot(romeo_max, invQeff_max, 'o', color=color, alpha=0.8, ms=10)
    

markers = {'.': 'point', ',': 'pixel', 'o': 'circle', 'v': 'triangle_down', '^': 'triangle_up', '<': 'triangle_left', '>': 'triangle_right', '1': 'tri_down', '2': 
           'tri_up', '3': 'tri_left', '4': 'tri_right', '8': 'octagon', 's': 'square', 'p': 'pentagon', '*': 'star', 'h': 'hexagon1', 'H': 'hexagon2', '+': 'plus', 
           'x': 'x', 'D': 'diamond', 'd': 'thin_diamond', '|': 'vline', '_': 'hline', 'P': 'plus_filled', 'X': 'x_filled', 0: 'tickleft', 1: 'tickright', 2: 'tickup', 
           3: 'tickdown', 4: 'caretleft', 5: 'caretright', 6: 'caretup', 7: 'caretdown', 8: 'caretleftbase', 9: 'caretrightbase', 10: 'caretupbase', 11: 'caretdownbase'}

kenn_photom={'n338_I' : ['mu0d_Ic', 'h_disc_I', 6.23, 'I', 60., 15.],
             'n338_B' : ['mu0d_Bc', 'h_disc_B', 6.99, 'B', 60., 15.],
            'n1167' : ['mu0d_Rc', 'h_disc_R', 5.53, 'R', 40., 6.7],
            'n2985_K' : ['mudK', 'hK', 1.42, 'K', 70., 14.],
            'n2985_S4G' : [('mu0d_s4g', 'mu0d_s4g_2'), ('h_disc_s4g', 'h_disc_s4g_2'), 1.35, 'S4G', 70., 6.3],
            'n3898_R' : ['mu0d_R', 'h_disc_R', 7.80, 'R', 80., 8.8],
            'n3898_R2' : [('mu_inner_approx', 'mu_out'), ('h_approx', 'h_out'), 2.93, 'R', 80., 8.8],
            'n4258_I' : ['mu0d_I', 'h_disc_I', 0.82, 'I', 150., 15.],
#             'n4258_S4G' : ['mu0d_36', 'h_disc_36', 0.38, 'S4G', 150., 15.],
            'n4725_H' : ['mu0d_H', 'h_disc_H', 0.68, 'H', 55., 10.14],
            'n4725_S4G' : ['mu0d_s4g', 'h_disc_s4g', 1.61, 'S4G', 55., 10.14],
            'n5533_r' : ['MU0_r', 'hi_r', 4.41, 'r', 100., 8.9],
            'n5533_R' : ['mu0d_Rc', 'h_disc_R', 7.15, 'R', 100., 9.9]}
    

fig, (ax, ax2) = plt.subplots(ncols=2, nrows=1, figsize=[16, 8])


color_counter = 0
color_shuffle = range(13)
shuffle(color_shuffle)

# for ind, name in enumerate(['n338', 'n1167', 'n2985', 'n3898', 'n4258', 'n4725', 'n5533']):
for ind, name in enumerate(['n4258']): 
    for key in kenn_photom.keys():
        
        if name in key:
            print key
            phot = kenn_photom[key]
            color = cm.rainbow(np.linspace(0, 1, 13))[color_shuffle[color_counter]]
            try:
                h = locals()[name+'dict'][phot[1][0]]
            except KeyError:
                h = locals()[name+'dict'][phot[1]]
        #     print name
        #     print mud
        #     print h

            if name in ['n4258', 'n4725']:
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], map(lambda l: l[0], zip(locals()[name+'dict']['gas_dens'])))[1:15]
            elif name == 'n2985':
                total_gas_data=locals()[name+'dict']['total_gas_data_']
            else:
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], [He_coeff*(locals()[name+'dict']['y_interp_'](l[0], h) + l[1]) for l in 
                                                                           zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'])])[1:15]


            while total_gas_data[0][0] < phot[5]:
                total_gas_data = total_gas_data[1:]

#             print total_gas_data
        
            ML = phot[2]
            band = phot[3]
        
            if type(phot[0]) == tuple and phot[3] == 'S4G':
                h1 = globals()[name+'dict'][phot[1][0]]
                h2 = globals()[name+'dict'][phot[1][1]]
                mu01 = globals()[name+'dict'][phot[0][0]]
                mu02 = globals()[name+'dict'][phot[0][1]]
                star_density = lambda l: tot_dens((lambda l1: s4g_surf_density(mu_disc(l1, mu0=mu01, h=h1), ML),
                               lambda l2: s4g_surf_density(mu_disc(l2, mu0=mu02, h=h2), ML)))(l) 
            elif type(phot[0]) == tuple:
                h1 = globals()[name+'dict'][phot[1][0]]
                h2 = globals()[name+'dict'][phot[1][1]]
                mu01 = globals()[name+'dict'][phot[0][0]]
                mu02 = globals()[name+'dict'][phot[0][1]]
                star_density = lambda l: tot_dens((lambda l1: surf_density(mu=mu_disc(l1, mu0=mu01, h=h1), M_to_L=ML, band=band),
                               lambda l2: surf_density(mu=mu_disc(l2, mu0=mu02, h=h2), M_to_L=ML, band=band)))(l)
            elif phot[3] == 'S4G':
                star_density = lambda l: s4g_surf_density(mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), ML)
            else:
                star_density = lambda l: surf_density(mu=mu_disc(l, mu0=globals()[name+'dict'][phot[0]], h=h), M_to_L=ML, band=band)
            
            
            if name == 'n338':
                epicycl = locals()[name+'dict']['epicyclicFreq_real_']
            else:
                epicycl = locals()[name+'dict']['epicyclicFreq_real']   
            
            plot_WS_vs_2f(ax=ax, total_gas_data=total_gas_data,
                              epicycl=epicycl,
                              gas_approx=locals()[name+'dict']['spl_gas'], 
                              sound_vel=locals()[name+'dict']['sound_vel'], 
                              scale=locals()[name+'dict']['scale'], 
                              sigma_max=locals()[name+'dict']['sig_R_maj_max'], 
                              sigma_min=locals()[name+'dict']['sig_R_maj_min'],
                              star_density_max=star_density, 
                              star_density_min=star_density, 
                              data_lim=locals()[name+'dict']['data_lim'], 
                          color=color, alpha=0.5, 
                          label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$',
                          sfrange=phot[4])
            ax.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$')
            
            if name in ['n4258', 'n4725']:
                
                r_mol_dens = locals()[name+'dict']['r_mol_dens']
                mol_dens = locals()[name+'dict']['mol_dens']
                y_interp = scipy.interpolate.interp1d(list(r_mol_dens), list(mol_dens))

                def y_interp_(r):
                    if r <= min(r_mol_dens):
                        return y_interp(min(r_mol_dens))
                    elif r < max(r_mol_dens):
                        return y_interp(r)
                    else:
                        return 0.
                
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['HI_dens'], 
                                   [y_interp_(l) for l in locals()[name+'dict']['r_g_dens']])[1:15]
                
            elif name == 'n2985':
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'], 
                                   [locals()[name+'dict']['y_interp_'](l) for l in locals()[name+'dict']['r_g_dens']])[1:15]
                
            else:
                total_gas_data=zip(locals()[name+'dict']['r_g_dens'], locals()[name+'dict']['gas_dens'], 
                                   [locals()[name+'dict']['y_interp_'](l, h) for l in locals()[name+'dict']['r_g_dens']])[1:15]


            while total_gas_data[0][0] < phot[5]:
                total_gas_data = total_gas_data[1:]
                
            plot_RF13_vs_2F(ax=ax, r_g_dens=zip(*total_gas_data)[0], 
                            HI_gas_dens=zip(*total_gas_data)[1], 
                            CO_gas_dens=zip(*total_gas_data)[2], 
                            epicycl=epicycl, 
                            sound_vel=6., 
                            sigma_R_max=locals()[name+'dict']['sig_R_maj_max'], 
                            sigma_R_min=locals()[name+'dict']['sig_R_maj_min'],  
                            star_density=star_density, 
                            alpha_max=0.7, 
                            alpha_min=0.3, 
                            scale=locals()[name+'dict']['scale'], 
                            gas_approx=locals()[name+'dict']['spl_gas'], 
                            thin=True,
                            color=color)    
            ax2.plot([-1, -2], [-1, -2], '-', color=color, label=r'$\rm{NGC\, '+name[1:]+'}\:\: '+phot[3]+'$')
            color_counter+=1
    
for _ in np.linspace(1., 2., 5):
    ax.plot([0., _], [0., 1.0], '--', alpha=0.3, color='gray')
    ax2.plot([0., 1.0], [0., _], '--', alpha=0.3, color='gray')

# ax.legend(loc='upper left')
# ax.set_xlim(0, 2.05)
ax.set_ylim(0, 1.05)
# ax.set_ylabel(r'$Q_{eff}^{-1}$', fontsize=20)
# ax.set_xlabel(r'$Q_{WS}^{-1}$', fontsize=20)
ax.grid()

# ax2.legend(loc='lower right')
# ax2.set_xlim(0, 1.05)
ax2.set_ylim(0, 1.05)
# ax2.set_ylabel(r'$Q_{eff}^{-1}$', fontsize=20)
# ax2.set_xlabel(r'$Q_{RF}^{-1}$', fontsize=20)
ax2.grid()

plt.show()


n4258_I
[(28.637732681578012, 90.39105833987557), (45.425369081123776, 35.243717313379562), (63.200513504172214, 16.972905545436063), (80.97565792722061, 14.668946801443612), (98.750802350269083, 10.432432612369912), (116.52594677331751, 10.956817095588416), (134.30109119636595, 9.7821991152168621), (153.0637436429171, 7.7949504380469508), (171.8263960894682, 7.0733900568403678), (188.61403248901397, 4.7911829205668797), (206.3891769120624, 3.0412119646485007), (223.17681331160813, 2.6283423968522013), (241.93946575815926, 2.271601615072282), (259.71461018120777, 3.2484021805019161)]
[(28.637732681578012, 85.134131416689144), (45.425369081123776, 33.116938249739718), (63.200513504172214, 16.530191458907307), (80.97565792722061, 14.452695431893424), (98.750802350269083, 10.09319199004219), (116.52594677331751, 10.875391264106273), (134.30109119636595, 7.9222111161279516), (153.0637436429171, 7.7969821602114839), (171.8263960894682, 6.7387100284340047), (188.61403248901397, 4.7925143650408817), (206.3891769120624, 3.1425102957409288), (223.17681331160813, 2.5450759879746965), (241.93946575815926, 2.1997506205034361), (259.71461018120777, 3.3025745319310271)]
Wall time: 30.8 s

In [ ]:


In [ ]: