Skip to content
GitLab
Commit a620f019 authored by NUNEZ Arturo's avatar NUNEZ Arturo
Browse files

Automatic commit samedi 14 juillet 2018, 16:30:02 (UTC+0200)

parent bf0d946c
Hide whitespace changes
Inline Side-by-side
LAM2LUPM/HaloA_2LUPM.ipynb
View file @ a620f019
This source diff could not be displayed because it is too large. You can view the blob instead.
LAM2LUPM/HaloB_2LUPM.ipynb
View file @ a620f019
This source diff could not be displayed because it is too large. You can view the blob instead.
LAM2LUPM/HaloC_2LUPM.ipynb
View file @ a620f019
This source diff could not be displayed because it is too large. You can view the blob instead.
LAM2LUPM/VDF_comparison_Plot_HALOA.ipynb
View file @ a620f019
This source diff could not be displayed because it is too large. You can view the blob instead.
LAM2LUPM/VDF_comparison_Plot_HALOB.ipynb
View file @ a620f019
This source diff could not be displayed because it is too large. You can view the blob instead.
LAM2LUPM/VDF_comparison_Plot_HaloC.ipynb
View file @ a620f019
This source diff could not be displayed because it is too large. You can view the blob instead.
Simulations/ALLPLOTS/AllPlots_HALOC.ipynb
View file @ a620f019
%% Cell type:code id: tags:
 
``` python
%matplotlib notebook
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/POL/HALOC_19/Hydro/output_00440"
#path = "/media/arturo/ARTUROTECA/OUTPUTS/HaloB/output_00417"
myhalo = wkbl.Galaxy_Hound(path)
 
ok,rho,_= CF.getDensity(np.array(myhalo.st.pos3d.reshape(len(myhalo.st.pos3d)*3),dtype=np.float32), myhalo.st.mass)
centro_rho = myhalo.st.pos3d[np.where(rho == rho.max())][0]
print "density",centro_rho
myhalo.center_shift(centro_rho)
myhalo.r_virial(600,n=3.)
```
 
%% Output
 
loading Dark matter..
loading Stars..
loading Gas..
density [9868.701 9744.951 9766.655]
 
/home/arturo/Documents/git/WKBL/wkbl/astro/galaxy_peeker.py:74: RuntimeWarning: divide by zero encountered in divide
rho_s = np.cumsum(mhist) / vol_bin
/home/arturo/Documents/git/WKBL/wkbl/astro/_dark_matter.py:68: RuntimeWarning: invalid value encountered in arccos
self.theta = np.arccos(np.copy(self.pos3d[:,0]),np.copy(self.r))
/home/arturo/Documents/git/WKBL/wkbl/astro/_stars.py:59: RuntimeWarning: invalid value encountered in arccos
self.theta = np.arccos(np.copy(self.pos3d[:,0]),np.copy(self.r))
/home/arturo/Documents/git/WKBL/wkbl/astro/_stars.py:63: RuntimeWarning: invalid value encountered in divide
self.vR = (vx*self.pos3d[:,0] + vy*self.pos3d[:,1])/ self.R
/home/arturo/Documents/git/WKBL/wkbl/astro/_stars.py:64: RuntimeWarning: invalid value encountered in divide
self.vr = (vx*self.pos3d[:,0] + vy*self.pos3d[:,1] + vz*self.pos3d[:,2])/ self.r
/home/arturo/Documents/git/WKBL/wkbl/astro/_stars.py:65: RuntimeWarning: invalid value encountered in divide
self.vphi = (-vx*self.pos3d[:,1] + vy*self.pos3d[:,0] )/ self.R
 
| r_200 = 182.23
| Diagonal matrix computed
| | 18, 0, 0|
| D =| 0, 13, 0|
| | 0, 0, 3|
 
%% Cell type:code id: tags:
 
``` python
myhalo.dm.m
```
%% Cell type:code id: tags:
``` python
disc = np.where((myhalo.st.r<8.5)&(myhalo.st.r>7.5)&(np.abs(myhalo.st.pos3d[:,2])<3))
hist40 , bins = np.histogram(myhalo.st.vphi[disc],bins=100)
```
 
%% Cell type:code id: tags:
 
``` python
fig, ax= plt.subplots()
ax.plot(bins[:-1],hist40)
```
 
%% Output
 
 
 
[<matplotlib.lines.Line2D at 0x7fd9cc5d8cd0>]
 
%% 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:markdown id: tags:
 
# print "{0:1.4e}".format(myhalo.dm.mass[0])
 
%% 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)
####
legend = ax.legend(loc='bottom left', ncol=1, shadow=False, fontsize=12)
```
 
%% Output
 
 
 
/usr/lib/python2.7/dist-packages/matplotlib/legend.py:325: UserWarning: Unrecognized location "bottom left". Falling back on "best"; valid locations are
right
center left
upper right
lower right
best
center
lower left
center right
upper left
upper center
lower center
six.iterkeys(self.codes))))
 
%% 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_gs(myhalo,[-length,length],150,thikness=0.3)#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,"SF0",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
 
 
 
143242288.25
346893217.1320877
 
<matplotlib.colorbar.Colorbar at 0x7fd96789a9d0>
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots()
length=1.2*myhalo.r200
SF1140_faceOn,edges= face_on_dm(myhalo,[-length,length],300)#H.T
 
mass_1 = ax.imshow(SF1140_faceOn+1e3, interpolation='nearest', origin='low',cmap="magma",
extent=[edges[0], edges[-1], edges[0], edges[-1]],
norm=LogNorm(vmin=5e6)
)
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
 
 
 
<matplotlib.text.Text at 0x7fd9674ee2d0>
 
%% 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.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+1e3, interpolation='nearest', origin='low',cmap="bone",
extent=[edges[0], edges[-1], edges[0], edges[-1]],
norm=LogNorm(vmin=1e5)
)
 
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 = ax2.imshow(SF1140_edgeOn+1e3, interpolation='nearest', origin='low',cmap="bone",
extent=[edges[0], edges[-1], edges[0], edges[-1]],
norm=LogNorm(vmin=1e5)
)
 
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',
extent=[edges[0], edges[-1], edges[0], edges[-1]],
norm=LogNorm(vmin=1e5)
)
 
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',
extent=[edges[0], edges[-1], edges[0], edges[-1]],
norm=LogNorm(vmin=1e5)
)
 
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
 
 
 
<matplotlib.text.Text at 0x7fd9662f5290>
 
%% 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 )
x = 0.03
y = -.02
#ax.scatter(np.log10(SF1_30.dm.total_m),np.log10(SF1_30.st.total_m),marker='s',c='r',s=30)
#ax.text(np.log10(SF1_30.dm.total_m)+x,np.log10(SF1_30.st.total_m)+y,'30 Myr',fontsize=13)
#ax.scatter(np.log10(SF1_90.dm.total_m),np.log10(SF1_90.st.total_m),marker='^',c='g',s=30)
#ax.text(np.log10(SF1_90.dm.total_m)+x,np.log10(SF1_90.st.total_m)+y,'90 Myr',fontsize=13)
ax.scatter(np.log10(myhalo.dm.total_m),np.log10(myhalo.st.total_m),marker='o',c='b',s=30,label='HALO B')
ax.scatter(12.158, 11.1386,marker='o',c='r',s=30,label='SF1 Mochima')
ax.scatter(12.1239, 11.1429,marker='o',c='k',s=30,label='SF1 Andorra, ep =4.5')
 
#ax.text(np.log10(myhalo.dm.total_m)+x,np.log10(myhalo.st.total_m)+y,'140 Myr',fontsize=13)
#ax.scatter(np.log10(SF0.dm.total_m),np.log10(SF0.st.total_m),marker='>',c='m',s=30)
#ax.text(np.log10(SF0.dm.total_m)+x,np.log10(SF0.st.total_m)+y,'Old Star Formation',fontsize=13)
 
legend = ax.legend(loc='upper right', ncol=2, shadow=False, fontsize=14)
frame = legend.get_frame()
```
 
%% Output
 
 
 
%% 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
 
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
```
 
%% Output
 
22376970 22376970
49257447 49257447
 
%% 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
 
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')
```
 
%% Output
 
/usr/local/lib/python2.7/dist-packages/ipykernel/__main__.py:18: RuntimeWarning: invalid value encountered in double_scalars
 
%% Cell type:code id: tags:
 
``` python
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.
```
 
%% Cell type:code id: tags:
 
``` python
file = open("vcdata.dat")
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
 
---------------------------------------------------------------------------
IOError Traceback (most recent call last)
<ipython-input-17-75f6c61d767c> in <module>()
----> 1 file = open("vcdata.dat")
2 R_mw = np.array([])
3 R_err = np.array([])
4 vc_mw = np.array([])
5 vc_err = np.array([])
IOError: [Errno 2] No such file or directory: 'vcdata.dat'
['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
 
 
 
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-18-0e1fe12d2f87> in <module>()
5
6
----> 7 print len(vc_mw),len(vc_err)
8 ax.scatter(R_mw,vc_mw,color='gray',marker='s',alpha=0.5,s=2,label="MW data\n(arXiv:1703.00020)")
9 ax.plot(r, vc_all,'b-',label='all')
NameError: name 'vc_mw' is not defined
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
Sf1_x, Sf1_y =11.0904, 2.5110828926
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="HALO B")
ax.scatter(Sf1_x,Sf1_y,color='r',marker='s',label="Mochima SF1")
 
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 B (Mollitor)aprox")
legend = ax.legend(loc='lower right', ncol=1, shadow=False, fontsize=15)
frame = legend.get_frame()
 
```
 
%% Output
 
 
 
%% 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,35])
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],age_hist[0]/(thikness*1e9))
```
 
%% Output
 
 
 
[<matplotlib.lines.Line2D at 0x7ff0a3edb2d0>]
 
%% Cell type:code id: tags:
 
``` python
print len(myhalo.dm.r[(myhalo.dm.r<7.4)&(myhalo.dm.r>7.5)])
```
 
%% Output
 
7989
 
%% Cell type:code id: tags:
 
``` python
```
 
%% Cell type:code id: tags:
 
``` python
```
 
%% 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