Commit daab818a authored by NUNEZ Arturo's avatar NUNEZ Arturo
Browse files

Automatic commit vendredi 29 juin 2018, 16:30:02 (UTC+0200)

parent 6d18476e
......@@ -2726,44 +2726,44 @@
}
],
"source": [
" fig, ax = plt.subplots(2,1,gridspec_kw = {'height_ratios':[3.5, 1,3.5,1]},figsize=[6,11],sharex=True)\n",
" length=17\n",
" ax[0].set_ylim([-length,length]);ax[0].set_xlim([-length,length])\n",
" ax[1].set_ylim([-length,length]);ax[1].set_xlim([-length,length])\n",
"fig, ax = plt.subplots(2,1,gridspec_kw = {'height_ratios':[3.5, 1,3.5,1]},figsize=[6,11],sharex=True)\n",
"length=17\n",
"ax[0].set_ylim([-length,length]);ax[0].set_xlim([-length,length])\n",
"ax[1].set_ylim([-length,length]);ax[1].set_xlim([-length,length])\n",
"\n",
" SF1140_faceOn,edges= face_on_gs(myhalo,[-length,length],150,thikness=0.2)#H.T \n",
" print SF1140_faceOn.max()\n",
" mass_2 = ax[0].imshow(SF1140_faceOn+0.001, interpolation='nearest', origin='low',cmap='viridis',\n",
" extent=[edges[0], edges[-1], edges[0], edges[-1]],\n",
" norm=LogNorm(vmin=1e5,vmax=5e8)\n",
" )\n",
" \"\"\"\n",
" divider = make_axes_locatable(ax)\n",
" cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
" cbar = plt.colorbar(mass_2,cax=cax,label=r'mass [M$_{\\odot}$]')\n",
" \"\"\"\n",
" ax[0].text(-15,12,\"SF1\",fontsize=25,color='w')\n",
"SF1140_faceOn,edges= face_on_gs(myhalo,[-length,length],150,thikness=0.2)#H.T \n",
"print SF1140_faceOn.max()\n",
"mass_2 = ax[0].imshow(SF1140_faceOn+0.001, interpolation='nearest', origin='low',cmap='viridis',\n",
" extent=[edges[0], edges[-1], edges[0], edges[-1]],\n",
" norm=LogNorm(vmin=1e5,vmax=5e8)\n",
" )\n",
"\"\"\"\n",
"divider = make_axes_locatable(ax)\n",
"cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
"cbar = plt.colorbar(mass_2,cax=cax,label=r'mass [M$_{\\odot}$]')\n",
"\"\"\"\n",
"ax[0].text(-15,12,\"SF1\",fontsize=25,color='w')\n",
"\n",
"\n",
"\n",
"\n",
" SF1140_faceOn,edges= edge_on_gs(myhalo,[-length,length],150)#H.T \n",
" print SF1140_faceOn.max()\n",
" mass_2 = ax[1].imshow(SF1140_faceOn+0.001, interpolation='nearest', origin='low',cmap='viridis',\n",
" extent=[edges[0], edges[-1], edges[0], edges[-1]],\n",
" norm=LogNorm(vmin=1e5,vmax=5e8)\n",
" )\n",
"SF1140_faceOn,edges= edge_on_gs(myhalo,[-length,length],150)#H.T \n",
"print SF1140_faceOn.max()\n",
"mass_2 = ax[1].imshow(SF1140_faceOn+0.001, interpolation='nearest', origin='low',cmap='viridis',\n",
" extent=[edges[0], edges[-1], edges[0], edges[-1]],\n",
" norm=LogNorm(vmin=1e5,vmax=5e8)\n",
" )\n",
"\n",
" ax[0].set_ylabel(\"Kpc\",fontsize=15)\n",
" ax[1].set_xlabel(\"Kpc\",fontsize=15)\n",
" #divider = make_axes_locatable(fig)\n",
" #cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
" #cbar = plt.colorbar(mass_2,cax=cax,label=r'mass [M$_{\\odot}$]')\n",
" fig.tight_layout(h_pad=-1,w_pad=5,)\n",
" fig.colorbar(mass_2, ax=ax.ravel().tolist(),label=r'mass [M$_{\\odot}$]')\n",
"ax[0].set_ylabel(\"Kpc\",fontsize=15)\n",
"ax[1].set_xlabel(\"Kpc\",fontsize=15)\n",
"#divider = make_axes_locatable(fig)\n",
"#cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
"#cbar = plt.colorbar(mass_2,cax=cax,label=r'mass [M$_{\\odot}$]')\n",
"fig.tight_layout(h_pad=-1,w_pad=5,)\n",
"fig.colorbar(mass_2, ax=ax.ravel().tolist(),label=r'mass [M$_{\\odot}$]')\n",
"\n",
" #cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])\n",
" #fig.colorbar(mass_2, cax=cbar_ax)"
"#cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])\n",
"#fig.colorbar(mass_2, cax=cbar_ax)"
]
},
{
%% Cell type:code id: tags:
 
``` python
%matplotlib notebook
%load_ext autoreload
%autoreload 2
```
 
%%%% Output: stream
 
The autoreload extension is already loaded. To reload it, use:
%reload_ext autoreload
 
%% Cell type:code id: tags:
 
``` python
 
from scipy.stats import rv_continuous
from scipy.interpolate import interp1d
from matplotlib.patches import Circle
from scipy.special import gamma
import numpy as np
import emcee
from mpl_toolkits.axes_grid1 import make_axes_locatable
from numpy import exp, sqrt
from scipy.integrate import quad, dblquad
import matplotlib.patches as patches
from itertools import product
from scipy.integrate import quad
import scipy.optimize as optimize
import matplotlib.pyplot as plt
import matplotlib as mpl
from sklearn.neighbors import KDTree
import sys
import lmfit
from py_unsio import *
import pymc
import os
from pymodelfit import FunctionModel1DAuto
import wkbl
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm
import wkbl.astro.nbody_essentials as nbe
import cfalcon
CF =cfalcon.CFalcon()
import iminuit
from iminuit import Minuit, describe, Struct
import probfit
from matplotlib.colors import LogNorm
```
 
%% Cell type:code id: tags:
 
``` python
path = "/data/OWN/Mochima/SF1/Stable/output_00041"
myhalo= wkbl.Galaxy_Hound(path)
print "loaded"
myhalo.r_virial(600)
print "cutted"
nucenter = nbe.real_center(myhalo.dm.pos3d, myhalo.dm.mass)
myhalo.center_shift(nucenter)
myhalo.redefine(4.5)
```
 
%%%% Output: stream
 
loading Dark matter..
loading Stars..
loading Gas..
loaded
| r_200 = 228.890 kpc
---- taking particles inside 2.5 * r200
| number of praticles inside 2.5 * r200
| dm mass = 1.058e+12 M_sun
| p_dm_200 = 9.266e+05 particles
| stellar mass = 1.219e+11 M_sun
| p_st_200 = 5.090e+05 psrticles
| gas mass = 1.491e+11 M_sun
| p_gs_200 = 2.048e+06 particles
---- rotating galaxy
| Diagonal matrix computed
| |20, 0, 0|
| D =| 0,19, 0|
| | 0, 0, 1|
cutted
 
%% Cell type:code id: tags:
 
``` python
myh
```
 
%%%% Output: error
 
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-8-c533dd5972f2> in <module>()
----> 1 myh
 
NameError: name 'myh' is not defined
 
%% Cell type:code id: tags:
 
``` python
print np.unique(myhalo.dm.mass)
```
 
%%%% Output: stream
 
[ 1528510.6 12228085. 97824680. ]
 
%% Cell type:code id: tags:
 
``` python
ok,myhalo.dm.rho,_= CF.getDensity(np.array(myhalo.dm.pos3d.reshape(len(myhalo.dm.pos3d)*3),
dtype=np.float32), myhalo.dm.mass)
```
 
%% Cell type:code id: tags:
 
``` python
ok,myhalo.st.rho,_= CF.getDensity(np.array(myhalo.st.pos3d.reshape(len(myhalo.st.pos3d)*3),
dtype=np.float32), myhalo.st.mass)
```
 
%% Cell type:code id: tags:
 
``` python
fig, ax= plt.subplots(figsize=[7,6])
ax.set_xlim([1e-2,10**2.5])
ax.set_ylim([1,1e12])
ax.scatter(myhalo.st.r, myhalo.st.rho,s=0.5,alpha=0.5,lw=0,c='y')
ax.scatter(1e20,1e20,c='y',alpha=0.5, label="stars")
ax.scatter(1e20,1e20,c='b',alpha=0.5, label="DM")
 
ax.scatter(myhalo.dm.r, myhalo.dm.rho,s=0.5,alpha=0.5,lw=0)
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_ylim([1e0,1e12])
ax.set_xlabel(r"r [kpc]",fontsize=20)
ax.set_ylabel(r"$\rho$ [M$_{\odot}$ kpc$^{-3}$]",fontsize=20)
ax.set_title("",fontsize=15)
####
legend = ax.legend(loc='bottom left', ncol=1, shadow=False, fontsize=12)
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%% Cell type:code id: tags:
 
``` python
def face_on_dm(sim,lims,points):
edges = np.linspace(lims[0],lims[1],points)
H, xedges, yedges = np.histogram2d(sim.dm.pos3d[:,0],
sim.dm.pos3d[:,1],
bins=(edges, edges),
weights=sim.dm.mass)
result = H.T
return result, edges
 
def face_on_st(sim,lims,points,thikness=.5):
disk = (np.abs(sim.st.pos3d[:,2])<thikness)
edges = np.linspace(lims[0],lims[1],points)
H, xedges, yedges = np.histogram2d(sim.st.pos3d[disk,0],
sim.st.pos3d[disk,1],
bins=(edges, edges),
weights=sim.st.mass[disk])
result = H.T
return result, edges
 
def face_on_gs(sim,lims,points,thikness=.9):
disk = (np.abs(sim.gs.pos3d[:,2])<thikness)
edges = np.linspace(lims[0],lims[1],points)
H, xedges, yedges = np.histogram2d(sim.gs.pos3d[disk,0],
sim.gs.pos3d[disk,1],
bins=(edges, edges),
weights=sim.gs.mass[disk])
result = H.T
return result, edges
 
def edge_on_st(sim,lims,points):
#disk = sim.st.pos3d[:,2]
edges = np.linspace(lims[0],lims[1],points)
H, xedges, yedges = np.histogram2d(sim.st.pos3d[:,0],
sim.st.pos3d[:,2],
bins=(edges, edges),
weights=sim.st.mass)
result = H.T
return result, edges
 
def edge_on_gs(sim,lims,points):
edges = np.linspace(lims[0],lims[1],points)
H, xedges, yedges = np.histogram2d(sim.gs.pos3d[:,0],
sim.gs.pos3d[:,2],
bins=(edges, edges),
weights=sim.gs.mass)
result = H.T
return result, edges
 
```
 
%% Cell type:code id: tags:
 
``` python
# Moster et all
def M_1(z):
M10 ,M11 = 11.590, 1.195
log = M10 + M11*(z / (z+1))
return 10**(log)
 
def N(z):
N10 ,N11 = 0.0351, -0.0247
return N10 + N11*(z / (z+1))
 
 
def beta(z):
B10 ,B11 = 1.376, -0.826
return B10 + B11*(z / (z+1))
 
 
def gamma(z):
G10 ,G11 = 0.608, 0.329
return G10 + G11*(z / (z+1))
 
def mm(M,z=0):
one = ( M / M_1(z))**(-beta(z))
two = ( M / M_1(z))**gamma(z)
return 2*N(z) * M / (one +two)
 
def alpha(m):
return 0.15 / np.log10(m)
 
M = np.logspace(10,14,50)
m = mm(M)
al = np.sqrt(m)#alpha(m)
```
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots(2,1,gridspec_kw = {'height_ratios':[3.5, 1,3.5,1]},figsize=[6,11],sharex=True)
length=17
ax[0].set_ylim([-length,length]);ax[0].set_xlim([-length,length])
ax[1].set_ylim([-length,length]);ax[1].set_xlim([-length,length])
 
SF1140_faceOn,edges= face_on_st(myhalo,[-length,length],150,thikness=0.2)#H.T
print SF1140_faceOn.max()
mass_2 = ax[0].imshow(SF1140_faceOn+0.001, interpolation='nearest', origin='low',cmap='bone',
extent=[edges[0], edges[-1], edges[0], edges[-1]],
norm=LogNorm(vmin=1e5,vmax=3e9)
)
"""
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
cbar = plt.colorbar(mass_2,cax=cax,label=r'mass [M$_{\odot}$]')
"""
ax[0].text(-15,12,"SF1",fontsize=25,color='w')
 
 
 
 
SF1140_faceOn,edges= edge_on_st(myhalo,[-length,length],150)#H.T
print SF1140_faceOn.max()
mass_2 = ax[1].imshow(SF1140_faceOn+0.001, interpolation='nearest', origin='low',cmap='bone',
extent=[edges[0], edges[-1], edges[0], edges[-1]],
norm=LogNorm(vmin=1e5,vmax=3e9)
)
 
ax[0].set_ylabel("Kpc",fontsize=15)
ax[1].set_xlabel("Kpc",fontsize=15)
#divider = make_axes_locatable(fig)
#cax = divider.append_axes("right", size="5%", pad=0.05)
#cbar = plt.colorbar(mass_2,cax=cax,label=r'mass [M$_{\odot}$]')
fig.tight_layout(h_pad=-1,w_pad=5,)
fig.colorbar(mass_2, ax=ax.ravel().tolist(),label=r'mass [M$_{\odot}$]')
 
#cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
#fig.colorbar(mass_2, cax=cbar_ax)
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%%%% Output: stream
 
669135626.921875
1566738775.71875
 
%%%% Output: execute_result
 
<matplotlib.colorbar.Colorbar at 0x7f4caced8310>
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots(2,1,gridspec_kw = {'height_ratios':[3.5, 1,3.5,1]},figsize=[6,11],sharex=True)
length=17
ax[0].set_ylim([-length,length]);ax[0].set_xlim([-length,length])
ax[1].set_ylim([-length,length]);ax[1].set_xlim([-length,length])
SF1140_faceOn,edges= face_on_gs(myhalo,[-length,length],150,thikness=0.2)#H.T
print SF1140_faceOn.max()
mass_2 = ax[0].imshow(SF1140_faceOn+0.001, interpolation='nearest', origin='low',cmap='viridis',
extent=[edges[0], edges[-1], edges[0], edges[-1]],
norm=LogNorm(vmin=1e5,vmax=5e8)
)
"""
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
cbar = plt.colorbar(mass_2,cax=cax,label=r'mass [M$_{\odot}$]')
"""
ax[0].text(-15,12,"SF1",fontsize=25,color='w')
fig, ax = plt.subplots(2,1,gridspec_kw = {'height_ratios':[3.5, 1,3.5,1]},figsize=[6,11],sharex=True)
length=17
ax[0].set_ylim([-length,length]);ax[0].set_xlim([-length,length])
ax[1].set_ylim([-length,length]);ax[1].set_xlim([-length,length])
 
SF1140_faceOn,edges= face_on_gs(myhalo,[-length,length],150,thikness=0.2)#H.T
print SF1140_faceOn.max()
mass_2 = ax[0].imshow(SF1140_faceOn+0.001, interpolation='nearest', origin='low',cmap='viridis',
extent=[edges[0], edges[-1], edges[0], edges[-1]],
norm=LogNorm(vmin=1e5,vmax=5e8)
)
"""
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
cbar = plt.colorbar(mass_2,cax=cax,label=r'mass [M$_{\odot}$]')
"""
ax[0].text(-15,12,"SF1",fontsize=25,color='w')
 
 
 
SF1140_faceOn,edges= edge_on_gs(myhalo,[-length,length],150)#H.T
print SF1140_faceOn.max()
mass_2 = ax[1].imshow(SF1140_faceOn+0.001, interpolation='nearest', origin='low',cmap='viridis',
extent=[edges[0], edges[-1], edges[0], edges[-1]],
norm=LogNorm(vmin=1e5,vmax=5e8)
)
ax[0].set_ylabel("Kpc",fontsize=15)
ax[1].set_xlabel("Kpc",fontsize=15)
#divider = make_axes_locatable(fig)
#cax = divider.append_axes("right", size="5%", pad=0.05)
#cbar = plt.colorbar(mass_2,cax=cax,label=r'mass [M$_{\odot}$]')
fig.tight_layout(h_pad=-1,w_pad=5,)
fig.colorbar(mass_2, ax=ax.ravel().tolist(),label=r'mass [M$_{\odot}$]')
 
#cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
#fig.colorbar(mass_2, cax=cbar_ax)
SF1140_faceOn,edges= edge_on_gs(myhalo,[-length,length],150)#H.T
print SF1140_faceOn.max()
mass_2 = ax[1].imshow(SF1140_faceOn+0.001, interpolation='nearest', origin='low',cmap='viridis',
extent=[edges[0], edges[-1], edges[0], edges[-1]],
norm=LogNorm(vmin=1e5,vmax=5e8)
)
ax[0].set_ylabel("Kpc",fontsize=15)
ax[1].set_xlabel("Kpc",fontsize=15)
#divider = make_axes_locatable(fig)
#cax = divider.append_axes("right", size="5%", pad=0.05)
#cbar = plt.colorbar(mass_2,cax=cax,label=r'mass [M$_{\odot}$]')
fig.tight_layout(h_pad=-1,w_pad=5,)
fig.colorbar(mass_2, ax=ax.ravel().tolist(),label=r'mass [M$_{\odot}$]')
#cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
#fig.colorbar(mass_2, cax=cbar_ax)
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%%%% Output: stream
 
49808552.0
222334649.16015625
 
%%%% Output: execute_result
 
<matplotlib.colorbar.Colorbar at 0x7f4ca26e4ed0>
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots(2,1,gridspec_kw = {'height_ratios':[3.5, 1,3.5,1]},figsize=[6,11],sharex=True)
length=17
ax[0].set_ylim([-length,length]);ax[0].set_xlim([-length,length])
ax[1].set_ylim([-length,length]);ax[1].set_xlim([-length,length])
 
SF1140_faceOn,edges= face_on_gs(myhalo,[-length,length],150,thikness=0.2)#H.T
print SF1140_faceOn.max()
mass_2 = ax[0].imshow(SF1140_faceOn+0.001, interpolation='nearest', origin='low',cmap='viridis',
extent=[edges[0], edges[-1], edges[0], edges[-1]],
norm=LogNorm(vmin=1e5,vmax=1e8)
)
#####################################################################################
#####################################################################################
#####################################################################################
ax[0].text(-15,12,"SF1",fontsize=25,color='w')
 
 
 
 
SF1140_faceOn,edges= edge_on_gs(myhalo,[-length,length],150)#H.T
print SF1140_faceOn.max()
mass_2 = ax[1].imshow(SF1140_faceOn+0.001, interpolation='nearest', origin='low',cmap='viridis',
extent=[edges[0], edges[-1], edges[0], edges[-1]],
norm=LogNorm(vmin=1e5,vmax=1e9)
)
 
ax[0].set_ylabel("Kpc",fontsize=15)
ax[1].set_xlabel("Kpc",fontsize=15)
 
fig.tight_layout(h_pad=-1,w_pad=5,)
fig.colorbar(mass_2, ax=ax.ravel().tolist(),label=r'mass [M$_{\odot}$]')
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%%%% Output: stream
 
49808552.0
222334649.16015625
 
%%%% Output: execute_result
 
<matplotlib.colorbar.Colorbar at 0x7f8955d64dd0>
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots()
SF1140_faceOn,edges= face_on_gs(myhalo,[-length,length],150)#H.T
 
mass_2 = ax.imshow(SF1140_faceOn+1e5, interpolation='nearest', origin='low',cmap='viridis',
extent=[edges[0], edges[-1], edges[0], edges[-1]],
norm=LogNorm(vmin=1e5,vmax=1e9)
)
 
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
cbar = plt.colorbar(mass_2,cax=cax,label=r'mass [M$_{\odot}$]')
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%% Cell type:code id: tags:
 
``` python
 
SF130_faceon, edges = face_on_dm(myhalo,[-1.2*myhalo.r200,1.2*myhalo.r200],300)
length = 15.
fig,[[ax,ax1,ax2],[ax3,ax4,ax5]] = plt.subplots(2,3,figsize=[22,13])
fig.tight_layout(w_pad=3)
 
#######################################################################################################################3
mass_1 = ax[0].imshow(SF130_faceon+1e3, interpolation='nearest', origin='low',cmap="magma",
extent=[edges[0], edges[-1], edges[0], edges[-1]],
norm=LogNorm(vmin=1e7)
)
divider = make_axes_locatable(ax)
 
cax = divider.append_axes("right", size="5%", pad=0.05)
cbar = plt.colorbar(mass_1,cax=cax,label=r'mass [M$_{\odot}$]')
 
ax.add_artist(Circle(xy=(0, 0),radius=myhalo.r200,color='w',ls='--',lw=1.7,fill=False))
ax.text(-100,1.1*myhalo.r200,r"R$_{200}$ = "+str(int(myhalo.r200))+" Kpc ",color='w',fontsize=15)
#######################################################################################################################3
 
SF1140_faceOn,edges = face_on_st(myhalo,[-length,length],200)#H.T
 
mass_2 = ax1.imshow(SF1140_faceOn+1e2, interpolation='nearest', origin='low',cmap="bone",
extent=[edges[0], edges[-1], edges[0], edges[-1]],
norm=LogNorm(vmin=1e5,vmax=1e9)
)
 
divider = make_axes_locatable(ax1)
cax = divider.append_axes("right", size="5%", pad=0.05)
cbar = plt.colorbar(mass_2,cax=cax,label=r'mass [M$_{\odot}$]')
#######################################################################################################################3
 
SF1140_edgeOn,edges = edge_on_st(myhalo,[-length,length],150)#H.T
 
 
mass_2 = ax[1].imshow(SF1140_edgeOn+1e3, interpolation='nearest', origin='low',cmap="bone",
extent=[edges[0], edges[-1], edges[0], edges[-1]],
norm=LogNorm(vmin=1e5,vmax=1e9)
)
 
divider = make_axes_locatable(ax2)
cax = divider.append_axes("right", size="5%", pad=0.05)
cbar = plt.colorbar(mass_2,cax=cax,label=r'mass [M$_{\odot}$]')
 
SF1140_faceOn,edges= face_on_gs(myhalo,[-length,length],150)#H.T
 
mass_2 = ax4.imshow(SF1140_faceOn+1e5, interpolation='nearest', origin='low',cmap='viridis',
extent=[edges[0], edges[-1], edges[0], edges[-1]],
norm=LogNorm(vmin=1e5,vmax=1e9)
)
 
divider = make_axes_locatable(ax4)
cax = divider.append_axes("right", size="5%", pad=0.05)
cbar = plt.colorbar(mass_2,cax=cax,label=r'mass [M$_{\odot}$]')
 
#######################################################################################################################3
 
SF1140_edgeOn,edges = edge_on_gs(myhalo,[-length,length],150)#H.T
 
 
mass_2 = ax5.imshow(SF1140_edgeOn, interpolation='nearest', origin='low',cmap='viridis',
extent=[edges[0], edges[-1], edges[0], edges[-1]],
norm=LogNorm(vmin=1e5,vmax=1e9)
)
 
divider = make_axes_locatable(ax5)
cax = divider.append_axes("right", size="5%", pad=0.05)
cbar = plt.colorbar(mass_2,cax=cax,label=r'mass [M$_{\odot}$]')
 
 
ax3.scatter(myhalo.dm.r, myhalo.dm.rho,s=0.5,alpha=0.5,lw=0)
ax3.set_xscale("log")
ax3.set_yscale("log")
ax3.set_ylim([1e0,1e12])
ax3.set_xlabel(r"r [kpc]",fontsize=20)
ax3.set_ylabel(r"$\rho$ [M$_{\odot}$ kpc$^{-3}$]",fontsize=20)
ax3.set_title("sf_model = 0",fontsize=15)
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%%%% Output: error
 
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-17-b78eb0fd6c1c> in <module>()
6
7 #######################################################################################################################3
----> 8 mass_1 = ax[0].imshow(SF130_faceon+1e3, interpolation='nearest', origin='low',cmap="magma",
9 extent=[edges[0], edges[-1], edges[0], edges[-1]],
10 norm=LogNorm(vmin=1e7)
TypeError: 'AxesSubplot' object does not support indexing
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots(figsize = [6,6])
ax.set_xlim(11.75,12.75)
ax.set_ylim(10,11.5)
ax.set_xlabel(r'Log$_{10}$(M$_{Halo}$/M$_{\odot}$)',fontsize=18)
ax.set_ylabel(r'Log$_{10}$(M$_{Stars}$/M$_{\odot}$)',fontsize=18)
 
ax.fill_between(np.log10(M), np.log10(m)+.3,np.log10(m)-.3,color='blue',alpha=0.5 )
size = 70
ax.scatter(np.log10(myhalo.dm.total_m),np.log10(myhalo.st.total_m),marker='o',c='b',s=size)
ax.scatter(np.log10(myhalo.dm.total_m),np.log10(myhalo.st.fire_m),marker='o',c='gray',s=size,linewidths=0.2)
ax.scatter(np.log10(myhalo.dm.total_m),np.log10(myhalo.st.mass[(myhalo.st.r<17.)].sum()),alpha=0.6,marker='o',c='g',s=size)
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%%%% Output: execute_result
 
<matplotlib.collections.PathCollection at 0x7f89540cb950>
 
%% Cell type:code id: tags:
 
``` python
r_arr = np.arange(0,30,30./100.)
_= np.histogram(myhalo.dm.r,bins=r_arr,weights=myhalo.dm.vtheta)
 
```
 
%% Cell type:code id: tags:
 
``` python
 
def circular_speed(r,comp="all"):
if comp=='all':
m = mass
radii = np.sqrt(r2)
elif comp=="dm":
m = myhalo.dm.mass
radii = myhalo.dm.r
elif comp=="st":
m = myhalo.st.mass
radii = myhalo.st.r
elif comp=="gs":
m = myhalo.gs.mass
radii = myhalo.gs.r
else:
sys.exit("no valid component")
enclosed_m = np.sum(m[np.where(radii<r)])
return np.sqrt(myhalo.p.G * enclosed_m / r) * myhalo.p.kpctokm
 
 
 
pos_dm = np.array(myhalo.dm.pos3d.reshape(len(myhalo.dm.pos3d)*3),dtype=np.float32)
pos_gs = np.array(myhalo.gs.pos3d.reshape(len(myhalo.gs.pos3d)*3),dtype=np.float32)
pos_st = np.array(myhalo.st.pos3d.reshape(len(myhalo.st.pos3d)*3),dtype=np.float32)
pos = np.concatenate((pos_dm, pos_st, pos_gs))
mass = np.concatenate((myhalo.dm.mass,myhalo.st.mass,myhalo.gs.mass))
v = np.concatenate((myhalo.dm.v,myhalo.st.v,myhalo.gs.v))
print len(mass)*3, len(pos)
pos3d = pos.reshape(len(pos)/3,3)
r2 = pos3d[:,0]**2 + pos3d[:,1]**2 +pos3d[:,2]**2
 
get_vc= np.vectorize(circular_speed)
r = np.linspace(0,40,60)
vc_all = get_vc(r,comp='all')
vc_dm = get_vc(r,comp='dm')
vc_st = get_vc(r,comp='st')
vc_gs = get_vc(r,comp='gs')
 
vc_stars_vphi = np.array([])
std_stars_vphi = np.array([])
vc_gas_vphi = np.array([])
std_gas_vphi = np.array([])
for i in range(len(r)-1):
stars_condition = (myhalo.st.R>r[i])&(myhalo.st.R<r[i+1])&(np.abs(myhalo.st.pos3d[:,2])<0.5)
gas_condition = (myhalo.gs.R>r[i])&(myhalo.gs.R<r[i+1])&(np.abs(myhalo.gs.pos3d[:,2])<0.5)
vc_gas_vphi = np.append(vc_gas_vphi, np.average(myhalo.gs.vphi[gas_condition]))
std_gas_vphi = np.append(std_gas_vphi, np.std(myhalo.gs.vphi[gas_condition]))
#for stars
vc_stars_vphi = np.append(vc_stars_vphi, np.nanmean(myhalo.st.vphi[stars_condition]))
std_stars_vphi = np.append(std_stars_vphi, np.nanstd(myhalo.st.vphi[stars_condition]))
r_arrays = (r[1:] + r[:-1]) / 2.
```
 
%%%% Output: stream
 
10452216 10452216
 
%% Cell type:code id: tags:
 
``` python
file = open("../vcdata.dat")
R_mw = np.array([])
R_err = np.array([])
vc_mw = np.array([])
vc_err = np.array([])
 
for ln in file:
row = ln.split('\t')
if row[0][0]=='#':
continue
R_mw = np.append(R_mw, row[0])
R_err = np.append(R_err, row[1])
vc_mw = np.append(vc_mw, row[2])
vc_err = np.append(vc_err, row[3])
print vc_err
```
 
%%%% Output: stream
 
['4.5' '4.5' '4.5' ... '19.8466566348' '19.7816973913' '23.7338181969']
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots(figsize=[7,7])
ax.set_ylim(0,600)
ax.set_xlim(0,20)
ax1.set_xlim(0,20)
 
 
print len(vc_mw),len(vc_err)
ax.scatter(R_mw,vc_mw,color='gray',marker='s',alpha=0.5,s=2,label="MW data\n(arXiv:1703.00020)")
ax.plot(r, vc_all,'b-',label='all')
ax.plot(r, vc_dm, 'k--',label='dark matter')
ax.plot(r, vc_st, 'y-.',label='stars')
ax.plot(r, vc_gs ,'r--',label='gas')
texto = r"$v_c = \left(\frac{G M_{<r}}{r}\right)^{1/2}$"
ax.text(5,400,texto,fontsize=15)
ax.set_xlabel("r [kpc]",fontsize=15)
ax.set_ylabel(r"v$_{c}$ [km / s]",fontsize=15)
legend = ax.legend(loc='upper right', ncol=2, shadow=False, fontsize=14)
frame = legend.get_frame()
ax1.set_ylim(-200,700)
ax1.plot(r_arrays,-vc_stars_vphi+std_stars_vphi,'r--')
ax1.plot(r_arrays,-vc_stars_vphi-std_stars_vphi,'r--')
ax1.fill_between(r_arrays,-vc_stars_vphi+std_stars_vphi, -vc_stars_vphi-std_stars_vphi,color='r',alpha=0.2)
ax1.plot(r_arrays,-vc_stars_vphi,'ro-',markersize=4,markeredgewidth=0,
label=r"$\left<v_{\phi}\right>_{stars} \pm 1 \sigma$ ")
ax1.plot(r_arrays,-vc_gas_vphi-std_gas_vphi,'g--')
ax1.plot(r_arrays,-vc_gas_vphi+std_gas_vphi,'g--')
ax1.fill_between(r_arrays,-vc_gas_vphi+std_gas_vphi, -vc_gas_vphi-std_gas_vphi,color='g',alpha=0.2)
texto_dos = "cylindric rings of\n"
texto_dos += r"$|z| < 0.5$ kpc"
ax1.text(5,400,texto_dos,fontsize=15)
ax1.scatter(R_mw,vc_mw,color='gray',marker='s',alpha=0.5,s=2)
 
ax1.plot(r_arrays,-vc_gas_vphi,'go-',markersize=4,markeredgewidth=0,label=r"$\left<v_{\phi}\right>_{gas}\pm 1 \sigma$")
 
ax1.plot(r,vc_all,'b-',label=r"$v_c = \left(\frac{G M_{<r}}{r}\right)^{1/2}$")
 
ax1.set_xlabel("r [kpc]",fontsize=15)
ax1.set_ylabel(r"v$_{c}$ [km / s]",fontsize=15)
legend = ax1.legend(loc='upper right', ncol=2, shadow=False, fontsize=15)
frame = legend.get_frame()
plt.savefig("/home/arturo/Documents/LAM/LAM2LUPM/report07-17/HaloB/circularSpeedsCilyndric.pdf")
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%%%% Output: stream
 
2590 2590
 
%% Cell type:code id: tags:
 
``` python
# total stellar mass inside 0.1*r97
Mstar = np.sum(myhalo.st.mass[myhalo.st.r < (myhalo.r97/10.)])
# wheighted histogram of mass as the radius grows
hist = np.histogram(myhalo.st.r[myhalo.st.r < (myhalo.r97/10.)],bins=5012,
weights=myhalo.st.mass[myhalo.st.r < (myhalo.r97/10.)])
# find radius wher 80% of the stellar mass is contained
"""
note that this radius is decided like that by mollitor because
in his sim the rotatioon curve are already flat, in our case this does
not happens..
"""
r_80 = hist[1][np.argmin(np.abs((np.cumsum(hist[0])/hist[0].sum())-0.8))]
# circular speed
v_rot = circular_speed(r_80)
```
 
%% Cell type:code id: tags:
 
``` python
 
fig, ax = plt.subplots()
ax.set_xlim(9,12)
ax.set_ylim(1.6,2.6)
ax.set_xlabel(r"Log$_{10}$(M$_{stars}$/M$_{\odot}$)",fontsize=18)
ax.set_ylabel(r"Log$_{10}$(v/(km/s))",fontsize=18)
x = np.arange(9.,14,0.2)
ax.plot(x, 2.179+(x-10.3)*0.259,"k--",label="Dutton et al (2001)")
ax.scatter(np.log10(Mstar),np.log10(v_rot),label="Mochima")
print np.log10(Mstar),np.log10(v_rot)
ax.scatter(11.25,2.55,color="m",marker="^",edgecolor="k",label="Halo A (Mollitor)aprox")
ax.scatter(10.95,2.4,color="y",marker="^",edgecolor="k",label="Halo C (Mollitor)aprox")
ax.scatter(10.87,2.38,color="c",marker="^",edgecolor="k",label="Halo C (Mollitor)aprox")
legend = ax.legend(loc='lower right', ncol=1, shadow=False, fontsize=15)
frame = legend.get_frame()
 
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%%%% Output: stream
 
11.061818 2.4956806129034885
 
%% Cell type:code id: tags:
 
``` python
h=0.673
H0 = h/0.9777
def t_z(z):
return (2./3./H0)*(1-(1/(1+z)**(1.5)))
 
print t_z(3)
```
 
%% Cell type:code id: tags:
 
``` python
 
age_hist = np.histogram(myhalo.st.age,bins=512,weights=myhalo.st.mass)
fig,ax = plt.subplots()
ax.set_ylim([0,1.5])
ax.set_ylabel("SFR [Msun/year]", fontsize=18)
ax.set_xlabel("lookback time [Gyr]", fontsize=18)
thikness=age_hist[1][1]-age_hist[1][0]
 
ax.plot(-1*age_hist[1][:-1],np.log10(age_hist[0]/(thikness*1e9)/myhalo.p.aexp))
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%%%% Output: execute_result
 
[<matplotlib.lines.Line2D at 0x7f894b27b910>]
 
%% Cell type:code id: tags:
 
``` python
myhalo.st.metal.max()
```
 
%% Cell type:code id: tags:
 
``` python
myhalo.st.metal.min()
```
 
%%%% Output: execute_result
 
2e-05
 
%% Cell type:code id: tags:
 
``` python
myhalo.st.age.min()/myhalo.p.aexp
```
 
%%%% Output: execute_result
 
-13.509659750861164
 
%% Cell type:code id: tags:
 
``` python
 
def face_on_st_age(sim,lims,points,thikness=5):
disk = (np.abs(sim.st.pos3d[:,2])<thikness)
edges = np.linspace(lims[0],lims[1],points)
H, xedges, yedges = np.histogram2d(sim.st.pos3d[disk,0],
sim.st.pos3d[disk,1],
bins=(edges, edges),
weights=-sim.st.age[disk])
result = H.T
return result, edges
```
 
%% Cell type:code id: tags:
 
``` python
AGE,edges = face_on_st_age(myhalo,[-length,length],300)#H.T
print AGE.min()
```
 
%%%% Output: stream
 
-0.0005839324439875782
 
%% Cell type:code id: tags:
 
``` python
fig, ax1 = plt.subplots()
AGE,edges = face_on_st_age(myhalo,[-length,length],200)#H.T
 
mass_2 = ax1.imshow(AGE+1, interpolation='nearest', origin='low',cmap="bone",
extent=[edges[0], edges[-1], edges[0], edges[-1]],
norm=LogNorm(vmin=AGE.min()+1,vmax=AGE.max())
)
 
divider = make_axes_locatable(ax1)
cax = divider.append_axes("right", size="5%", pad=0.05)
cbar = plt.colorbar(mass_2,cax=cax,label=r'mass [M$_{\odot}$]')
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%% Cell type:code id: tags:
 
``` python
 
SF130_faceon, edges = face_on_dm(myhalo,[-1.2*myhalo.r200,1.2*myhalo.r200],300)
length = 15.
fig,ax = plt.subplots(figsize=[6,6])
fig.tight_layout(w_pad=3)
 
#######################################################################################################################3
mass_1 = ax.imshow(SF130_faceon+1e3, interpolation='nearest', origin='low',cmap="magma",
extent=[edges[0], edges[-1], edges[0], edges[-1]],
norm=LogNorm(vmin=1e7)
)
divider = make_axes_locatable(ax)
 
cax = divider.append_axes("right", size="5%", pad=0.05)
cbar = plt.colorbar(mass_1,cax=cax,label=r'mass [M$_{\odot}$]')
 
ax.add_artist(Circle(xy=(0, 0),radius=myhalo.r200,color='w',ls='--',lw=1.7,fill=False))
ax.text(-100,1.1*myhalo.r200,r"R$_{200}$ = "+str(int(myhalo.r200))+" Kpc ",color='w',fontsize=15)
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%%%% Output: execute_result
 
<matplotlib.text.Text at 0x7f2a86610b50>
 
%% Cell type:code id: tags:
 
``` python
```
 
%% Cell type:code id: tags:
 
``` python
```
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment