Commit 04c24cb6 authored by NUNEZ Arturo's avatar NUNEZ Arturo
Browse files

Automatic commit mercredi 30 mai 2018, 17:26:39 (UTC+0200)

parent 58264bc3
%% Cell type:code id: tags:
 
``` python
%matplotlib notebook
%load_ext autoreload
%autoreload 2
```
 
%% Cell type:code id: tags:
 
``` python
from scipy.stats import rv_continuous
from scipy.special import gamma
from scipy.interpolate import interp1d
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 scipy.integrate import simps
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
import warnings
from matplotlib.ticker import FormatStrFormatter
warnings.filterwarnings('ignore')
```
 
%% Cell type:code id: tags:
 
``` python
def abg_logprofile(x,p_s,r_s,al,be,ga):
x = 10**x
power = (be - ga) / (al)
denominator = ((x/(r_s))**ga) * ((1 + (x / (r_s))**al)**power)
return np.log10(10**p_s / denominator)
 
def abg_profile(x,po,r_s,al,be,ga):
power = (be - ga) / al
denominator = ((x/r_s)**ga) * ((1 + (x / r_s)**al)**power)
return (10**po) / denominator
```
 
%% Cell type:code id: tags:
 
``` python
"""
pathsim = "/data/POL/HALOB/hydro/output_00417"
#path = "/media/arturo/ARTUROTECA/OUTPUTS/HaloB/output_00417"
myhalo = wkbl.Galaxy_Hound(pathsim,"halo,stars,gas")
"""
 
path = "/data/OWN/SF1test/SF0/mstar1_T3600/output_00041"
myhalo= wkbl.Galaxy_Hound(path)
print "loaded"
```
 
%% Output
 
loading Dark matter..
loading Stars..
loading Gas..
loaded
 
%% Cell type:code id: tags:
 
``` python
myhalo.r_virial(600,n=24)
print "cutted"
nucenter = nbe.real_center(myhalo.dm.pos3d[myhalo.dm.r<myhalo.r200], myhalo.dm.mass[myhalo.dm.r<myhalo.r200])
myhalo.center_shift(nucenter)
myhalo.redefine(24)
```
 
%% Output
 
| r_200 = 228.997 kpc
---- taking particles inside 24 * r200
| number of praticles inside 24 * r200
| dm mass = 1.060e+12 M_sun
| p_dm_200 = 2.837e+06 particles
| stellar mass = 2.064e+11 M_sun
| p_st_200 = 9.015e+05 psrticles
| gas mass = 1.207e+13 M_sun
| p_gs_200 = 9.330e+06 particles
---- rotating galaxy
| Diagonal matrix computed
| |17, 0, 0|
| D =| 0,13, 0|
| | 0, 0, 2|
cutted
 
%% Cell type:code id: tags:
 
``` python
```
 
%% 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))
phi = np.concatenate((myhalo.dm.phi,myhalo.st.phi,myhalo.gs.phi))
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
 
39207024 39207024
 
%% Cell type:code id: tags:
 
``` python
ok, acc, Phy = CF.getGravity(pos,mass,0.15,G=myhalo.p.G)
Phy = Phy * myhalo.p.kpctokm**2
```
 
%% Cell type:code id: tags:
 
``` python
```
 
%% Cell type:markdown id: tags:
 
# Rmax from max of potential
# (spherical and/or line of sight of the neighbor galaxy)
 
 
### Spherical
 
 
%% Cell type:code id: tags:
 
``` python
 
bin_num = 512
 
pot_sph, bins_pot = np.histogram(r2,bins=bin_num,
weights=Phy)
n, _ = np.histogram(r2,bins=bin_num)
 
```
 
%% Cell type:code id: tags:
 
``` python
bins_pot = np.linspace(0.,myhalo.dm.r.max(),512)
print bins_pot.max()
```
 
%% Output
 
5495.915096235554
 
%% Cell type:code id: tags:
 
``` python
bins_pot = np.linspace(0.,myhalo.dm.r.max(),512)
mean_phi = np.array([])
std_phi = np.array([])
r_phi = np.array([])
for i in range(len(bins_pot)-1):
contidion = (r2>bins_pot[i]**2)&(r2<bins_pot[i+1]**2)
 
mean_phi = np.append(mean_phi,np.average(Phy[contidion]))
std_phi = np.append(std_phi,np.std((Phy[contidion])))
r_phi = np.append(r_phi, ((bins_pot[i])+(bins_pot[i+1]))/2)
 
 
```
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots(figsize=[7,6])
#ax.set_title("spherical shells",fontsize=30)
ax.set_ylim([-50e4,0])
ax.set_xlim([0,4500])
 
ax.set_ylabel(r'$\Phi$ [km $^2$ M$_{\odot}$ s$^{-2}$ ] ', fontsize=15)
#ax.set_xscale('log')
ax.set_xlabel(r'r [kpc] ', fontsize=15)
ax.plot(r_phi,mean_phi,'bo-',markersize=2)
ax.plot(r_phi,mean_phi+std_phi,'b--')
ax.plot(r_phi,mean_phi-std_phi,'b--',label=r"$1 \sigma$")
 
legend = ax.legend(loc='lower right', ncol=1, shadow=False, fontsize=14)
frame = legend.get_frame()
#fig.text(0.75,0.7,pathsim[10:14]+" "+pathsim[14],fontsize=18)
#ax.add_patch(
# patches.Rectangle((450, -1.5e4),
# 100,
# 1e4,
# fill=False # remove background
# )
#)
 
"""
left, bottom, width, height = [0.45, 0.23, 0.4, 0.4]
ax2 = fig.add_axes([left, bottom, width, height])
ax2.set_xlim([450,550])
ax2.set_ylim([-1.5e4,-5e3])
 
ax2.set_ylabel(r'$\Phi$ [kpc $^2$ M$_{\odot}$ s$^{-2}$ ] ', fontsize=15)
ax2.set_xlabel(r'r [kpc] ', fontsize=15)
ax2.plot(r_phi,mean_phi,'bo-',markersize=2)
ax2.plot(r_phi,mean_phi+std_phi,'b--')
ax2.plot(r_phi,mean_phi-std_phi,'b--')
"""
fig.tight_layout(w_pad=20)
ax.tick_params(axis='both', which='major', labelsize=15, size=5,width=1.2)
ax.tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
#fig.savefig("/home/arturo/Documents/LAM/LAM2LUPM/report07-17/HaloB/POTvsR_spherical.pdf",dpi=300)
```
 
%% Output
 
 
 
%% Cell type:markdown id: tags:
 
### L.O. S
 
%% Cell type:code id: tags:
 
``` python
ok,rho_dm,_= 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
# center of big neighbour
#neighbourg = rho_dm[np.where((myhalo.dm.r<4000)&(myhalo.dm.r>2000))].max()
#myneighbourg = myhalo.dm.pos3d[np.where(rho_dm==neighbourg)][0]
 
neighbourg = Phy[np.where((r2<4000**2)&(r2>2000**2))].min()
myneighbourg = pos3d[np.where(Phy==neighbourg)][0]
# its radius squared
r_nei2 = np.sum(myneighbourg**2)
phi_neigh = np.arctan2(myneighbourg[1],myneighbourg[0])
print np.sqrt(r_nei2)
the_neigh = np.arccos(myneighbourg[0]/np.sqrt(r_nei2))
# distance to every point proyected to the line conecting center to neigblourg
adyacent = (pos3d[:,0]*myneighbourg[0] + pos3d[:,1]*myneighbourg[1] + pos3d[:,2]*myneighbourg[2]) / np.sqrt(r_nei2)
# distance to each point
hipoteneuse = np.sqrt(pos3d[:,0]**2 + pos3d[:,1]**2 + pos3d[:,2]**2)
# angle of cone
alpha = np.radians(12)
cos_alpha = np.cos(alpha)
# cosine of all particles respective to their angle to the l.o.s
cos_all = adyacent / hipoteneuse
# final selection of cone
my_cone = pos3d[np.where((cos_all)>cos_alpha)]
```
 
%% Output
 
3871.106
 
%% Cell type:code id: tags:
 
``` python
# radii from center of neighbour
r_neig = np.sqrt((pos3d[:,0] - myneighbourg[0])**2 +(pos3d[:,1] - myneighbourg[1])**2 +(pos3d[:,2] - myneighbourg[2])**2 )
# calculating the R200 of neighbout
mhist, rhist = np.histogram(r_neig,range=(0.0,np.sqrt(r_nei2)),bins=512, weights=mass )
vol_bin = (4./3.)*np.pi*(rhist[:-1]**3)
r_bin = rhist[:-1]+ 0.5*(rhist[2]-rhist[1])
rho_s = np.cumsum(mhist) / vol_bin
r200_neigh = r_bin[np.argmin(np.abs(rho_s - (200 * myhalo.p.rho_crit)))]
# mass of neighbour
phi_neigh = np.arctan2(myneighbourg[1],myneighbourg[0])
mass_neigh = np.sum(mass[np.where(r_neig<r200_neigh)])
print "the biggest neighbourg has a mass of {0:.4e} Msun and a r200 of {1:.4f} kpc".format(mass_neigh,r200_neigh)
print "and it is located at {0:.4f} kpc from de center of our halo ".format(np.sqrt(r_nei2))
```
 
%% Output
 
the biggest neighbourg has a mass of 1.0808e+13 Msun and a r200 of 472.5471 kpc
and it is located at 3871.1060 kpc from de center of our halo
 
%% Cell type:code id: tags:
 
``` python
los_condition = np.where((cos_all)>cos_alpha)
cone_pos = my_cone
cone_mass = mass[los_condition]
cone_Phy = Phy[los_condition]
cone_r2 = my_cone[:,0]**2 + my_cone[:,1]**2 + my_cone[:,2]**2
cone_r = np.sqrt(cone_r2)
cone_phi = phi[los_condition]
 
#cone_prox = cone_r*np.cos(phi_neigh)
#cone_proy = cone_r*np.sin(phi_neigh)
 
# now the histogram
bin_num = 1000
pot_los, bins_pot_los = np.histogram(cone_r2,bins=bin_num,
weights=cone_Phy)
n_los, _ = np.histogram(cone_r2,bins=bin_num)
```
 
%% Cell type:code id: tags:
 
``` python
#a proyection Story
th,ph,ep = the_neigh, phi_neigh, -3*np.pi/4.
Ct, Cp, Ce = np.cos(th),np.cos(ph), np.cos(ep)
St, Sp, Se = np.sin(th),np.sin(ph), np.sin(ep)
#the rotation matrix
rot = np.array([[(Ce*Ct*Cp)+(Se*Sp),(Ce*Ct*Sp)-(Se*Cp),Ce*St],\
[(Se*Ct*Cp)-(Ce*Sp),(Se*Ct*Sp)+(Ce*Cp),Se*St],\
[-Cp*St,-St*Sp,Ct]])
cone_prox = rot[0,0]*cone_pos[:,0]+rot[0,1]*cone_pos[:,1]+rot[0,2]*cone_pos[:,2]
cone_proy = rot[1,0]*cone_pos[:,0]+rot[1,1]*cone_pos[:,1]+rot[1,2]*cone_pos[:,2]
```
 
%% Cell type:code id: tags:
 
``` python
```
 
%% Cell type:code id: tags:
 
``` python
bins_pot_los = np.linspace(0.,myhalo.dm.r.max(),512)
mean_phi_los = np.array([])
std_phi_los = np.array([])
r_phi_los = np.array([])
for i in range(len(bins_pot_los)-1):
contidion = (cone_r2>bins_pot_los[i]**2)&(cone_r2<bins_pot_los[i+1]**2)
mean_phi_los = np.append(mean_phi_los,np.average(cone_Phy[contidion]))
std_phi_los = np.append(std_phi_los,np.std((cone_Phy[contidion])))
r_phi_los = np.append(r_phi_los, ((bins_pot_los[i])+(bins_pot_los[i+1]))/2)
 
 
```
 
%% Cell type:code id: tags:
 
``` python
R_MAX = 1100.#r_phi_los[mean_phi_los==mean_phi_los.max()]
print R_MAX
```
 
%% Output
 
1100.0
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots(1,2,figsize=[12,5])
 
#fig, ax[0] = plt.subplots()
ax[0].set_xlim([0,5000])
ax[0].set_ylim([-7.7e5,0])
#ax.set_ylim([pot_sph.min(),np.abs(pot_sph.min())/6.])
ax[0].set_ylabel(r'$\Phi$ [kpc $^2$ M$_{\odot}$ s$^{-2}$ ] ', fontsize=15)
#ax.set_xscale('log')
ax[0].set_xlabel(r'r [kpc] ', fontsize=15)
ax[0].plot(r_phi_los, mean_phi_los+std_phi_los, 'b--')
ax[0].plot(r_phi_los, mean_phi_los-std_phi_los, 'b--')
ax[0].plot(r_phi_los, mean_phi_los, 'bo-',markersize=0.8,label="l.o.s")
ax[0].plot(r_phi,mean_phi,color='gray',ls='-', label='Shell')
ax[0].plot(r_phi,mean_phi+std_phi,color='gray',ls='--')
ax[0].plot(r_phi,mean_phi-std_phi,color='gray',ls='--')
ax[0].axvline(x=R_MAX,color='k',linestyle='--')
 
 
ax[1].set_ylim([0,4e3])
ax[1].set_xlim([0,4e3])
ax[1].set_xlabel(r' kpc ', fontsize=15)
ax[1].set_ylabel(r' kpc ', fontsize=15)
texto_neigh = "First massive neighbourg: \n"
texto_neigh += "M = {0:.4e}".format(mass_neigh)
texto_neigh += r" M$_{\odot}$"+"\n"
texto_neigh += r"r$_{200} = $ "+"{0:.2f} kpc\n".format(r200_neigh)
texto_neigh += r"r$_{max} = $ "+"{0:.2f} kpc".format(float(R_MAX))
fig.text(.6,.8,"L.O.S",fontsize=17)
fig.text(.6,.61,texto_neigh)
#ax[1].scatter(my_cone[:,0], my_cone[:,1],c="#902b08",lw=0,s=0.04,alpha=0.6)
ax[1].scatter(cone_prox, cone_proy,c="#902b08",lw=0,s=0.01,alpha=0.4)
 
#ax[1].scatter(myneighbourg[0],myneighbourg[1],s=20)
ax[1].scatter(np.sqrt(r_nei2)*np.cos(np.pi/4.),np.sqrt(r_nei2)*np.sin(np.pi/4.),s=30,lw=0)
circle = plt.Circle(((0),(0)),myhalo.r200, color='k',lw=1, fill=False)
ax[1].add_patch( circle )
circle = plt.Circle(((0),(0)),myhalo.r97, color='k',lw=1, fill=False)
ax[1].add_patch( circle )
circle = plt.Circle(((0),(0)),R_MAX, color='k',ls='--',lw=1, fill=False)
ax[1].add_patch( circle )
ax[1].text(-myhalo.r200-15,20,r'r$_{200}$')
ax[1].text(-myhalo.r97-100,-20,r'r$_{97}$')
 
circle = plt.Circle((np.sqrt(r_nei2)*np.cos(np.pi/4.),np.sqrt(r_nei2)*np.sin(np.pi/4.)),r200_neigh, color='k',lw=0.9, fill=False)
ax[1].add_patch( circle )
 
 
legend = ax[0].legend(loc='lower left', ncol=1, shadow=False, fontsize=16)
frame = legend.get_frame()
fig.tight_layout(w_pad=0)
ax[0].tick_params(axis='both', which='major', labelsize=15, size=5,width=1.2)
ax[0].tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
ax[1].tick_params(axis='both', which='major', labelsize=15, size=5,width=1.2)
ax[1].tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
#fig.savefig("/home/arturo/Documents/git/LAMtoLUPM_latex/Rmax.pdf")
```
 
%% Output
 
 
 
%% Cell type:code id: tags:
 
``` python
4500/166.
```
 
%% Output
 
27.10843373493976
 
%% Cell type:markdown id: tags:
 
# vesc(r)+-dvesc at different radii + comparison with the value calculated from the potential
 
$$v_{esc} = \sqrt(-2\left(\phi(r) - \phi(r_{max}\right))$$
 
%% Cell type:code id: tags:
 
``` python
bin_num = 60
times = 3.
bins = np.logspace(np.log10(3*myhalo.gs.hsml.min()),np.log10((times*myhalo.r200)**2),bin_num)
pot_sph_vesc, bins_pot_vesc = np.histogram(r2[(r2<(times*myhalo.r200)**2)], bins=bins, weights=Phy[(r2<(times*myhalo.r200)**2)])
```
 
%% Cell type:code id: tags:
 
``` python
diference = np.abs(bins_pot-R_MAX)
good_index = np.where(diference==diference.min())[0]
rmax = R_MAX
pot_max =np.sum((pot_sph/n)[good_index])/len(good_index)
print pot_max
```
 
%% Output
 
-99536.50121986923
 
%% Cell type:code id: tags:
 
``` python
tonei = (r_phi_los<np.sqrt(r_nei2))
r_relevant = r_phi_los[tonei]
R_MAXH = r_relevant[mean_phi_los[tonei]==(mean_phi_los[tonei]).max()][0]
pot_maxH = mean_phi_los[tonei].max()
print pot_max
```
 
%% Output
 
-99536.50121986923
 
%% Cell type:code id: tags:
 
``` python
def vesc_from_pot(limmin, limmax):
"""
a function calculating the vesc from potential inside interval
"""
contidion = (r2>limmin**2)&(r2<limmax**2)
mean = np.average(Phy[contidion])
sigma = np.std((Phy[contidion]))
v_esc = np.sqrt(2*np.abs(mean - pot_max))
sig_vesc = sigma / v_esc
return v_esc, sig_vesc
```
 
%% Cell type:code id: tags:
 
``` python
mean_phi_vesc = np.array([])
std_phi_vesc = np.array([])
r_phi_vesc = np.array([])
for i in range(len(bins_pot_vesc)-1):
contidion = (r2>bins_pot_vesc[i])&(r2<bins_pot_vesc[i+1])
mean_phi_vesc = np.append(mean_phi_vesc,np.average(Phy[contidion]))
std_phi_vesc = np.append(std_phi_vesc,np.std((Phy[contidion])))
r_phi_vesc = np.append(r_phi_vesc, (np.sqrt(bins_pot_vesc[i])+np.sqrt(bins_pot_vesc[i+1]))/2)
```
 
%% Cell type:code id: tags:
 
``` python
##################################
v_esc = np.sqrt(2*np.abs(mean_phi_vesc - pot_max))
sig_vesc = std_phi_vesc / v_esc
```
 
%% Cell type:code id: tags:
 
``` python
potdatHyd = np.loadtxt("/home/arturo/Downloads/Psi_Mochima2_DM_baryons_Rmax=1104.0kpc_dimensionful (1).txt")
```
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots()
ax.set_xscale('log')
ax.plot(r_phi_vesc, v_esc,'bo-',markersize=3,label=r"$v_{escape}$")
ax.plot(r_phi_vesc, v_esc+sig_vesc,'b--',markersize=3,label=r"1$\sigma$")
ax.plot(r_phi_vesc, v_esc-sig_vesc,'b--',markersize=3)
ax.set_xlabel("r [kpc]",fontsize=15)
ax.set_ylabel(r"v$_{esc}$ [km/s] ",fontsize=15)
ax.axvline(x=myhalo.r200,color='k',linestyle='--',label=r"r$_{200}$")
ax.axvline(x=8,color='y',linestyle='--',label=r"r$_{\odot}$")
ax.plot(potdatHyd[:,0],np.sqrt(2*potdatHyd[:,1]),"r--", lw=2,label="reconstructed")
 
legend = ax.legend(loc='lower left', ncol=1, shadow=False, fontsize=14)
frame = legend.get_frame()
```
 
%% Output
 
 
 
%% Cell type:code id: tags:
 
``` python
limmin, limmax = 7.5, 8.5
def slice_vmax(rmin, rmax,phimin,phimax):
in_the_shell = (myhalo.dm.r > rmin)&(myhalo.dm.r < rmax)
in_the_slice = (myhalo.dm.phi > phimin)&(myhalo.dm.phi < phimax)
result = myhalo.dm.v[in_the_shell&in_the_slice]
return result.max(),len(result)
slice_vesc = np.vectorize(slice_vmax)
 
def vesc_per_radius(rmin, rmax,n=20):
if not bool(int(rmax)%50):
print rmax
slices = np.linspace(-np.pi,np.pi,n)
v_maxs, N = slice_vesc(rmin,rmax,slices[:-1],slices[1:])
return np.mean(v_maxs), np.std(v_maxs), np.sum(N)
 
get_vesc_from_slices = np.vectorize(vesc_per_radius)
```
 
%% Cell type:code id: tags:
 
``` python
unbound_dm = dm_inside = st_inside = unbound_st = shell_vol = rho_max_bin = np.array([])
bins_lin = np.sqrt(bins_pot_vesc)
for i in range(len(bins_pot_vesc)-1):
inside_bin_dm = np.where((myhalo.dm.r>=bins_lin[i])&(myhalo.dm.r<bins_lin[i+1]))
inside_bin_st = np.where((myhalo.st.r>=bins_lin[i])&(myhalo.st.r<bins_lin[i+1]))
dm_inside = np.append(dm_inside,len(myhalo.dm.mass[inside_bin_dm]))
st_inside = np.append(st_inside,len(myhalo.st.mass[inside_bin_st]))
unbound_dm = np.append(unbound_dm, len(np.where(myhalo.dm.v[inside_bin_dm]>v_esc[i])[0]))
unbound_st = np.append(unbound_st,len(np.where(myhalo.st.v[inside_bin_st]>v_esc[i])[0]))
shell_vol = np.append(shell_vol, 4.* np.pi * ((bins_lin[i+1]**3)-(bins_lin[i]**3)) / 3.)
rho_max_bin = np.append(rho_max_bin, rho_dm[inside_bin_dm].mean())
```
 
%% Cell type:code id: tags:
 
``` python
slices = np.linspace(-np.pi,np.pi,20)
#slice_vesc(7.5,8.5,slices[:-1],slices[1:])
```
 
%% Cell type:code id: tags:
 
``` python
slices = np.linspace(-np.pi,np.pi,20)
an_array = np.linspace(0,myhalo.r200,60)
vesc_per_radius = get_vesc_from_slices(np.sqrt(bins_pot_vesc[:-1]),np.sqrt(bins_pot_vesc[1:]))
#vesc_per_radius = get_vesc_from_slices(an_array[:-1],an_array[1:])
```
 
%% Output
 
0.8483203065325621
0.848320306533
0.952148159965
 
%% Cell type:code id: tags:
 
``` python
fig, [ax,ax1,ax2] = plt.subplots(3,1,gridspec_kw = {'height_ratios':[3., 1,1]},figsize=[8,10.5],sharex=True)
ax.set_xlim([3.*myhalo.gs.hsml.min(),3*myhalo.r200])
ax.set_ylim([70,1000])
ax.set_xscale('log')
font = 17
##### from potential
ax.plot(r_phi_vesc, v_esc,'b-',markersize=3,label=r"$v_{esc}=\sqrt{-2(\Phi(r)-\Phi(r_{max}))}$")
ax.plot(r_phi_vesc, v_esc+sig_vesc,'b--',markersize=3)
#ax.plot(r_phi_vesc+2000, v_esc+sig_vesc+2000,color='gray', linestyle='--',label=r"1$\sigma$")
ax.plot(r_phi_vesc, v_esc-sig_vesc,'b--',markersize=3)
ax.fill_between(r_phi_vesc,v_esc+sig_vesc, v_esc-sig_vesc,color='b',alpha=0.2)
ax2.set_xlabel("r [kpc]",fontsize=font)
ax1.set_ylabel(r"$\frac{v_{max,shell}-v_{esc}}{v_{esc}}$",fontsize=1.5*font)
ax.set_ylabel(r"v$_{esc}$ [km/s] ",fontsize=font)
ax.axvline(x=myhalo.r200,color='k',linestyle='--')
ax1.axvline(x=myhalo.r200,color='k',linestyle='--')
ax2.axvline(x=myhalo.r200,color='k',linestyle='--')
ax.axvline(x=8,color='y',linestyle='--',linewidth=2)
ax1.axvline(x=8,color='y',linestyle='--',linewidth=2)
ax2.axvline(x=8,color='y',linestyle='--',linewidth=2)
rr = (an_array[:-1]+an_array[1:])/2
##### from shell
ax.plot(r_phi_vesc,vesc_per_radius[0],color="#902b08",linestyle="-",label=r"$v_{max,shell}$")
ax.plot(r_phi_vesc,vesc_per_radius[0]+vesc_per_radius[1],color="#902b08",linestyle="--")
ax.plot(r_phi_vesc,vesc_per_radius[0]-vesc_per_radius[1],color="#902b08",linestyle="--")
ax.fill_between(r_phi_vesc,vesc_per_radius[0]+vesc_per_radius[1],vesc_per_radius[0]-vesc_per_radius[1],color="#902b08",alpha=0.2)
#,yerr=vesc_per_radius[1], xerr=(an_array[1:]-an_array[:-1])/2,
######### second panel
ax1.set_ylim([-0.28,0.12])
ax1.plot(r_phi_vesc,[0 for i in r_phi_vesc],color='b',lw=2)
ax1.plot(r_phi_vesc,(vesc_per_radius[0]-v_esc)/v_esc,'r-',lw=2,label=r"$v_{esc}$ from mass")
ax.text(myhalo.r200+10,580,r"r$_{200}$",fontsize=font)
ax.text(8.6,v_esc[np.abs(r_phi_vesc-8)==np.abs(r_phi_vesc-8).min()],r"r$_{sun}$",fontsize=font)
legend = ax.legend(loc='lower left', ncol=1, shadow=False, fontsize=font)
frame = legend.get_frame()
fig.tight_layout(h_pad=-0.90)
ax.plot(potdatHyd[:,0],np.sqrt(2*potdatHyd[:,1]),"r--", lw=2,label="reconstructed")
 
 
 
ax2.set_ylim([-0.1,1.39])
ax2.set_xscale('log')
ax2.set_ylabel('Franction of\nunbound\nparticles\nper bin',fontsize=font)
ax2.plot(r_phi_vesc,unbound_dm / dm_inside,color='#212f3c',ls='-',lw=2,label="Dark Matter")
ax2.plot(r_phi_vesc,unbound_st / st_inside,color='#7b241c',ls='-',lw=2,label="Stars")
 
legend = ax2.legend(loc='center left', ncol=1, shadow=False, fontsize=font)
frame = legend.get_frame()
ax.grid()
ax1.grid()
ax2.grid()
 
 
ax.tick_params(axis='both', which='major', labelsize=15, size=5,width=1.2)
ax.tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
ax1.tick_params(axis='both', which='major', labelsize=15, size=5,width=1.2)
ax1.tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
ax2.tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
ax2.tick_params(axis='both', which='major', labelsize=15, size=5,width=1.2)
 
#plt.savefig("/home/arturo/Documents/git/LAMtoLUPM_latex/v_esc.png",dpi=300)
```
 
%% Output
 
 
 
%% Cell type:code id: tags:
 
``` python
for i in range(len(potdatHyd[:,1])):
print potdatHyd[i,0],potdatHyd[i,1],np.sqrt(2*potdatHyd[i,1])
```
 
%% Output
 
5.238999999999999e-06 395345.2978830042 889.2078473371726
5.340481945312715e-06 395345.2975131295 889.2078469212128
5.443929644628954e-06 395345.2969065218 889.2078462390239
5.549381175547536e-06 395345.29638144537 889.2078456485248
5.656875353248687e-06 395345.29580696294 889.207845002464
5.766451744781276e-06 395345.2955179537 889.2078446774451
5.878150683626916e-06 395345.2960499875 889.2078452757685
5.9920132845461125e-06 395345.2945314326 889.207843568007
6.108081458712054e-06 395345.29391366313 889.2078428732657
6.226397929137448e-06 395345.2936271269 889.207842551028
6.347006246400271e-06 395345.2925839087 889.2078413778285
6.4699508046741e-06 395345.2928799601 889.2078417107668
6.5952768580690034e-06 395345.2922529872 889.2078410056754
6.723030537288851e-06 395345.29194434785 889.2078406585806
6.853258866611406e-06 395345.2913192367 889.2078399555828
6.9860097811972565e-06 395345.289543719 889.2078379588419
7.121332144734107e-06 395345.2904370225 889.2078389634479
7.25927576742271e-06 395345.29049215006 889.2078390254441
7.399891424311337e-06 395345.28950279846 889.2078379128228
7.543230873985305e-06 395345.28893673356 889.2078372762282
7.689346877618626e-06 395345.2885216291 889.2078368094031
7.838293218394537e-06 395345.2879236819 889.2078361369539
7.99012472130237e-06 395345.28876822913 889.2078370867287
8.144897273317814e-06 395345.2757127089 889.2078224045365
8.302667843974146e-06 395345.29181446414 889.2078405125138
8.463494506331816e-06 395345.2847012853 889.207832513058
8.627436458354346e-06 395345.28595754574 889.207833925844
8.79455404469822e-06 395345.285311785 889.2078331996238
8.96490877892494e-06 395345.28482898127 889.2078326566644
9.138563366143178e-06 395345.284550564 889.2078323435574
9.31558172608968e-06 395345.2833239915 889.2078309641582
9.496029016657168e-06 395345.28354917833 889.2078312174026
9.679971657878082e-06 395345.2824981261 889.2078300353929
9.867477356372725e-06 395345.28260774014 889.2078301586645
1.0058615130271155e-05 395345.28224952926 889.2078297558218
1.0253455334617733e-05 395345.2818175328 889.20782927
1.0452069687267867e-05 395345.2823715759 889.207829893075
1.0654531295286193e-05 395345.28098867903 889.207828337874
1.0860914681856315e-05 395345.2804368275 889.2078277172637
1.107129581371167e-05 395345.27999926556 889.2078272251831
1.1285752129097787e-05 395345.27957744227 889.2078267508022
1.1504362566276252e-05 395345.27909805113 889.2078262116805
1.1727207592580604e-05 395345.27856233966 889.2078256092213
1.195436923403527e-05 395345.27806380036 889.2078250485657
1.2185931105548228e-05 395345.2777157411 889.2078246571396
1.2421978441688287e-05 395345.2773740657 889.2078242728926
1.2662598128058807e-05 395345.27682620694 889.2078236567726
1.2907878733278994e-05 395345.27590075345 889.2078226160107
1.315791054158471e-05 395345.2659771502 889.2078114559613
1.341278558606089e-05 395345.28111572895 889.207828480754
1.3672597682517431e-05 395345.2740849983 889.2078205740189
1.3937442464021535e-05 395345.2744851146 889.2078210239883
1.4207417416098803e-05 395345.2742441425 889.2078207529919
1.4482621912616357e-05 395345.27365862916 889.2078200945257
1.476315725236075e-05 395345.2731287702 889.2078194986481
1.5049126696324692e-05 395345.2727573767 889.2078190809804
1.5340635505715874e-05 395345.2723434889 889.2078186155235
1.56377909807022e-05 395345.2715994937 889.2078177788292
1.5940702499907258e-05 395345.27150174876 889.2078176689055
1.624948156067112e-05 395345.2708080611 889.2078168887869
1.6564241820090827e-05 395345.27083675214 889.2078169210527
1.6885099136856058e-05 395345.2707386728 889.207816810753
1.7212171613894845e-05 395345.26957231824 889.2078154990747
1.754557964184567e-05 395345.2697678462 889.2078157189648
1.788544594337144e-05 395345.2687644533 889.2078145905526
1.823189561833207e-05 395345.26816965913 889.2078139216492
1.858505618983173e-05 395345.2678392474 889.2078135500693
1.894505765115837e-05 395345.2673823859 889.2078130362845
1.9312032513632294e-05 395345.2669300693 889.2078125276108
1.968611585538183e-05 395345.26659013104 889.2078121453173
2.0067445371063386e-05 395345.2660064673 889.2078114889312
2.0456161422544994e-05 395345.265209517 889.2078105926837
2.085240709057145e-05 395345.2680741363 889.2078138142246
2.1256328227430364e-05 395345.2709228264 889.2078170178515
2.166807351063808e-05 395345.26463883533 889.2078099508971
2.2087794497665852e-05 395345.26394172333 889.2078091669273
2.2515645681725904e-05 395345.2635187146 889.2078086912131
2.2951784548638133e-05 395345.2627890842 889.2078078706734
2.3396371634798494e-05 395345.26302943303 889.2078081409688
2.384957058626988e-05 395345.26193135086 889.2078069060694
2.431154821901799e-05 395345.2626064673 889.2078076653031
2.4782474580314006e-05 395345.2612297965 889.2078061171039
2.526252301132626e-05 395345.2607713579 889.2078056015455
2.575187021092478e-05 395345.2604113112 889.2078051966381
2.6250696300721543e-05 395345.2618486995 889.20780681312
2.6759184891370567e-05 395345.2596162596 889.2078043025259
2.7277523150152576e-05 395345.25901428866 889.2078036255515
2.7805901869868193e-05 395345.25864068134 889.2078032053939
2.834451553906621e-05 395345.2582956863 889.2078028174137
2.8893562413632145e-05 395345.25773211295 889.2078021836211
2.9453244589763042e-05 395345.2572716603 889.2078016657977
3.002376807835636e-05 395345.2568463068 889.2078011874466
3.060534288083955e-05 395345.25638368266 889.2078006671811
3.1198183066468474e-05 395345.2558543886 889.2078000719389
3.1802506851123455e-05 395345.2544901271 889.2077985376951
3.241853667763099e-05 395345.2604586296 889.2078052498523
3.304649929764192e-05 395345.2538495836 889.2077978173421
3.3686625855095636e-05 395345.2542126653 889.2077982256626
3.4339151971300396e-05 395345.25394184445 889.2077979210984
3.50043178316624e-05 395345.25336773373 889.2077972754554
3.568236827409449e-05 395345.2529315295 889.2077967849017
3.6373552879137335e-05 395345.252430995 889.2077962220023
3.707812606182674e-05 395345.25145682285 889.2077951264517
3.779634716533968e-05 395345.25157407613 889.2077952583144
3.852848055645505e-05 395345.251234609 889.2077948765508
3.92747957228631e-05 395345.2506238581 889.2077941897024
4.0035567372360274e-05 395345.2503695131 889.2077939036669
4.081107553396467e-05 395345.24994729913 889.2077934288466
4.160160566099099e-05 395345.24990854145 889.2077933852598
4.240744873612163e-05 395345.249003475 889.2077923674252
4.322890137851348e-05 395345.2485671506 889.2077918767363
4.406626595297866e-05 395345.24798346055 889.2077912203205
4.491985068128097e-05 395345.2477335755 889.2077909393007
4.5789969755587645e-05 395345.24726764916 889.2077904153215
4.667694345411932e-05 395345.2468057465 889.2077898958673
4.758109825903928e-05 395345.24642726826 889.207789470232
4.850276697662713e-05 395345.2362254349 889.2077779972856
4.9442288859779824e-05 395345.23797146825 889.2077799608686
5.0400009732886086e-05 395345.2454142758 889.2077883310242
5.137628211911878e-05 395345.24452666927 889.207787332825
5.2371465370194e-05 395345.24410392717 889.2077868574107
5.3385925798643095e-05 395345.24363075383 889.2077863252816
5.442003681264692e-05 395345.2430570012 889.2077856800414
5.547417905348239e-05 395345.24286328285 889.2077854621864
5.6548740535630384e-05 395345.24232323485 889.2077848548503
5.764411678959854e-05 395345.241750301 889.2077842105308
5.8760711007510714e-05 395345.2414092656 889.2077838270036
5.9898934191515386e-05 395345.2405434015 889.2077828532558
6.105920530507013e-05 395345.2406334944 889.207782954574
6.224195142715584e-05 395345.239733334 889.2077819422567
6.344760790947817e-05 395345.23991262447 889.2077821438862
6.467661853671473e-05 395345.23902214144 889.2077811424521
6.592943568986499e-05 395345.2381454369 889.2077801565132
6.720652051276563e-05 395345.23776904866 889.2077797332282
6.850834308183118e-05 395345.23760876997 889.2077795529793
6.983538257908164e-05 395345.2371572015 889.207779045147
7.118812746852293e-05 395345.23857187526 889.2077806360843
7.256707567594331e-05 395345.258115156 889.20780261439
7.39727347721924e-05 395345.2374103437 889.2077793298299
7.540562216001138e-05 395345.23528354295 889.2077769380371
7.686626526448049e-05 395345.234910086 889.2077765180486
7.835520172715692e-05 395345.2348042293 889.2077763990026
7.987297960397352e-05 395345.23388110055 889.2077753608552
8.142015756696916e-05 395345.2333605405 889.2077747754352
8.299730510992847e-05 395345.2327856415 889.2077741289057
8.460500275800417e-05 395345.2323261061 889.2077736121138
8.624384228139966e-05 395345.231886878 889.2077731181594
8.79144269131919e-05 395345.2316958463 889.2077729033258
8.961737157137156e-05 395345.23115647474 889.2077722967504
9.13533030851861e-05 395345.23021871684 889.2077712421511
9.312286042586712e-05 395345.2294108157 889.2077703335882
9.492669494182446e-05 395345.229414251 889.2077703374515
9.676547059840073e-05 395345.2287044128 889.2077695391699
9.863986422226522e-05 395345.22834429145 889.2077691341788
0.00010055056575054479 395345.22786754026 889.207768598026
0.0001024982784847802 395345.227535592 889.207768224718
0.00010448371934980025 395345.2307152725 889.2077718005759
0.00010650761915761198 395345.22439111804 889.2077646884535
0.00010857072287640188 395345.22610101954 889.2077666114028
0.00011067378990474773 395345.22510613024 889.2077654925538
0.00011281759435114298 395345.22461339773 889.2077649384285
0.00011500292531893354 395345.2240444734 889.2077642986181
0.0001172305871967761 395345.22347837687 889.2077636619879
0.00011950139995472244 395345.2226734121 889.2077627567273
0.00012181619944603721 395345.22239020665 889.2077624382354
0.0001241758377148644 395345.2218179011 889.2077617946226
0.0001265811833098521 395345.2192592679 889.2077589171925
0.00012903312160385247 395345.2205876404 889.2077604110756
0.00013153255511981566 395345.2202475232 889.207760028581
0.00013408040386299393 395345.21972536633 889.2077594413652
0.00013667760565958293 395345.2188955335 889.2077585081379
0.0001393251165019228 395345.21823081543 889.2077577605983
0.0001420239109003836 395345.2175850323 889.2077570343528
0.00014477498224206917 395345.21694961865 889.2077563197688
0.00014757934315646869 395345.2165961224 889.2077559222281
0.00015043802588819033 395345.22076079843 889.2077606058085
0.00015335208267691684 395345.2147821722 889.2077538822659
0.00015632258614471797 395345.21427035314 889.207753306676
0.00015935062969086778 395345.2137202973 889.2077526880851
0.0001624373278943103 395345.21284246794 889.2077517008811
0.00016558381692391793 395345.2126250438 889.2077514563667
0.0001687912549566997 395345.2115477226 889.207750244815
0.0001720608226041096 395345.21073594794 889.2077493318959
0.00017539372334661186 395345.20994449494 889.2077484418306
0.00017879118397666615 395345.2098824606 889.2077483720669
0.0001822544550502901 395345.20856410667 889.2077468894506
0.00018578481134737167 395345.2076849117 889.207745900711
0.00018938355234089855 395345.2068029259 889.2077449088328
0.00019305200267527283 395345.20609579084 889.207744113591
0.00019679151265389514 395345.20538583776 889.2077433151803
0.00020060345873619125 395345.2044364456 889.2077422474971
0.0002044892440442664 395345.20351561427 889.2077412119332
0.00020845029887937594 395345.2025709733 889.2077401495933
0.000212488081248396 395345.1993384099 889.2077365142634
0.00021660407740049635 395345.19883287314 889.2077359457386
0.00022079980237420667 395345.19993398554 889.2077371840458
0.00022507680055508 395345.1988950222 889.2077360156312
0.00022943664624416028 395345.198010427 889.2077350208184
0.00023388094423745727 395345.19705348136 889.2077339446406
0.00023841133041664975 395345.1957126271 889.2077324367204
0.00024302947235123158 395345.19494180847 889.2077315698605
0.0002477370699123169 395345.1938179041 889.2077303059214
0.00025253585589833985 395345.1924765122 889.2077287973966
0.00025742759667287226 395345.1916236677 889.2077278382906
0.0002624140928147958 395345.19049567747 889.2077265697566
0.0002674971797810713 395345.1892054268 889.207725118745
0.0002726787285823403 395345.18784549 889.2077235893647
0.00027796064647161897 395345.18723317795 889.2077229007606
0.0002833448776463317 395345.18567833496 889.2077211521895
0.00028883340396393875 395345.184278783 889.207719578258
0.0002944282456714314 395345.17840354925 889.20771297099
0.0003001314621489552 395345.18219485664 889.2077172346815
0.0003059451526678371 395345.1801429299 889.2077149270916
0.0003118714571632996 395345.1787010937 889.2077133056075
0.00031791255702213583 395345.1772403218 889.2077116628283
0.0003240706758856484 395345.1757200609 889.207709953148
0.0003303480804681417 395345.1742621715 889.2077083136105
0.0003367470813912623 395345.17224367964 889.2077060436213
0.0003432700340345081 395345.1708889751 889.2077045201252
0.00034991933940220883 395345.1697722083 889.207703264213
0.0003566974450073 395345.1673420097 889.2077005312198
0.0003636068457722201 395345.16557009483 889.2076985385303
0.0003706500849472507 395345.16370722716 889.2076964435555
0.00037782975504665263 395345.1618179217 889.2076943188489
0.0003851484988029346 395345.15984758176 889.2076921030111
0.0003926090101396008 395345.15775165625 889.2076897459403
0.00040021403516274763 395345.1568038495 889.2076886800401
0.0004079663731718653 395345.15131891205 889.2076825116976
0.00041586887769021804 395345.15129492077 889.2076824847171
0.0004239244575151883 395345.14892060164 889.2076798145657
0.0004321360777889579 395345.1472783467 889.2076779676913
0.0004405067610899363 395345.14414429374 889.2076744431457
0.00044903958854533056 395345.1407263073 889.2076705992895
0.0004577377009652584 395345.1388878079 889.2076685317192
0.0004666042999988382 395345.136260005 889.2076655765007
0.00047564264931266914 395345.1343404276 889.2076634177504
0.0004848560757921382 395345.13056187995 889.2076591684081
0.0004942479707660026 395345.1276588173 889.2076559036334
0.000503821791254683 395345.1244934525 889.2076523438747
0.000513581061242746 395345.1213312455 889.207648787667
0.0005235293729760316 395345.1180589907 889.2076451077
0.0005336703882839049 395345.1145042028 889.2076411099972
0.0005440078399271275 395345.1082443213 889.207634070155
0.000554545532971825 395345.1131519718 889.2076395892826
0.0005652873461900838 395345.1038376256 889.2076291143993
0.0005762372334876654 395345.09984471695 889.2076246239873
0.0005873992253593863 395345.0959228478 889.207620213466
0.0005987774303726807 395345.0917624074 889.2076155346482
0.000610376036679909 395345.08738118364 889.2076106075382
0.0006221993135599509 395345.0830375215 889.2076057226698
0.0006342516129896704 395345.07789359626 889.207599937828
0.0006465373712458101 395345.0736367755 889.2075951506212
0.000659061110537926 395345.0686871811 889.2075895843233
0.0006718274406729409 395345.0635493625 889.2075838063489
0.0006848410607519548 395345.0582382132 889.2075778334473
0.0006981067608999084 395345.05292340275 889.2075718564286
0.0007116294240287653 395345.04726318095 889.2075654909611
0.0007254140276348364 395345.0397799676 889.2075570753632
0.0007394656456309257 395345.0331097759 889.2075495740867
0.0007537894502139664 395345.0289251719 889.2075448680941
0.0007683907137688265 395345.02256578545 889.2075377163483
0.0007832748108090033 395345.0157767993 889.2075300814757
0.0007984472199548981 395345.00889132527 889.2075223380932
0.0008139135259504226 395345.0016923201 889.2075142421145
0.0008296794217186564 395344.9941051725 889.2075057096318
0.0008457507104573353 395344.98587724566 889.2074964565309
0.0008621333077749181 395344.9785021686 889.207488162542
0.0008788332438680438 395344.9701926179 889.2074788176468
0.000895856665741151 395344.9614313696 889.2074689647739
0.0009132098394691097 395344.95243120164 889.2074588432124
0.0009308991525036647 395344.94326134276 889.207448530817
0.0009489311160245714 395344.93369996734 889.2074377781231
0.0009673123673362597 395344.9255503318 889.2074286130676
0.0009860496723109367 395344.9126398029 889.2074140939254
0.0010051499278789981 395344.9027029397 889.2074029189587
0.0010246201645676978 395344.89161275025 889.2073904469646
0.0010444675490889727 395344.88010764326 889.2073775083552
0.0010646993869774157 395344.8680096435 889.2073639029802
0.0010853231252793282 395344.85573046847 889.2073500938558
0.0011063463552938772 395344.84351746336 889.2073363591456
0.0011277768153673358 395344.8296700335 889.2073207863658
0.0011496223937414638 395344.81591470074 889.207305317158
0.0011718911314570468 395344.8016273661 889.207289249662
0.001194591225313695 395344.7868446798 889.2072726250948
0.001217731030886959 395344.77143903536 889.2072552999501
0.0012413190656039024 395344.7554559325 889.207237325397
0.0012653640118782373 395344.7393857216 889.2072192528822
0.0012898747203062017 395344.72170166404 889.2071993654392
0.0013148602129243276 395344.70397174626 889.2071794264217
0.0013403296865303314 395344.68558421155 889.2071587478495
0.0013662925160683106 395344.6663493451 889.2071371163696
0.0013927582580795374 395344.64652479225 889.2071148217295
0.0014197366542200707 395344.62587657 889.2070916007924
0.0014472376348465298 395344.60500023235 889.2070681233166
0.001475271322671304 395344.58245998324 889.2070427746097
0.0015038480364885872 395344.559433337 889.207016878901
0.0015329782949725654 395344.53544032475 889.2069898964186
0.0015626728205491997 395344.51095467294 889.2069623599142
0.0015929425433429834 395344.48491350334 889.2069330740774
0.0016237986052001759 395344.45898384164 889.2069039136411
0.0016552523637899427 395344.4318988332 889.2068734539035
0.0016873153967849545 395344.4020016568 889.2068398316072
0.0017199995061229643 395344.37222356163 889.2068063432282
0.00175331672235092 395344.34104040655 889.2067712747205
0.0017872793090532462 395344.3088242526 889.2067350445031
0.0018218997673658786 395344.2758361607 889.2066979461645
0.001857190840577765 395344.2422681425 889.2066601956403
0.001893165518821472 395344.2054818128 889.2066188258079
0.0019298370438546766 395344.16830619023 889.2065770181755
0.0019672189139342535 395344.12976602407 889.206533675978
0.0020053248887848007 395344.0897645202 889.2064886903606
0.0020441689946633797 395344.0480505543 889.2064417789092
0.0020837655295223947 395344.00397673465 889.2063922135677
0.002124129068272443 395343.9609460733 889.2063438213577
0.002165274468147147 395343.9135166973 889.2062904823574
0.0022072168741718687 395343.865363008 889.2062363287923
0.0022499717247383916 395343.8152483718 889.2061799699458
0.002293554757287551 395343.7631132961 889.2061213389122
0.0023379820141019763 395343.7095560605 889.2060611085155
0.0023832698482110044 395343.6541491234 889.2059987979427
0.0024294349294100106 395343.59515322896 889.205932451228
0.002476494250396301 395343.5346312673 889.2058643882948
0.0025244651330238966 395343.47299005854 889.2057950666522
0.002573365234679437 395343.40721800824 889.2057210994633
0.002623212554781629 395343.33971302246 889.2056451834103
0.0026740254414065506 395343.26439270197 889.2055604782305
0.002725822598041334 395343.1970705768 889.2054847678087
0.0027786230904686257 395343.12155415217 889.2053998420749
0.002832446353784447 395343.0433280503 889.205311869031
0.0028873121995519483 395342.9616483853 889.2052200121019
0.002943240823093782 395342.87732680136 889.2051251840617
0.0030002528109256857 395342.7887734812 889.2050255970006
0.0030583691483341 395342.69893788133 889.2049245678763
0.003117611227100529 395342.604472196 889.2048183317453
0.0031780008533755667 395342.5065462013 889.204708204136
0.003239560255705402 395342.40483783314 889.2045938228537
0.0033023120932138476 395342.2993053335 889.2044751409358
0.003366279463942809 395342.1910941497 889.2043534465514
0.003431485913354369 395342.0758229232 889.204223812419
0.0034979554429975003 395341.9576741944 889.2040909422251
0.003565712519342717 395341.83525626257 889.2039532708596
0.0036347820827877988 395341.70722915046 889.20380929138
0.0037051895568380204 395341.5755919636 889.2036612519806
0.003776960857464148 395341.4370190067 889.2035054125762
0.0038501224026417434 395341.2958738173 889.2033466804062
0.003924701122075252 395341.1480773424 889.2031804681565
0.0040007244671104066 395340.9945722109 889.2030078359057
0.004078220420838703 395340.83516855276 889.2028285701218
0.004157217508397537 395340.6700262626 889.2026428506189
0.00423774480746993 395340.4933016492 889.2024441055582
0.004319831958987564 395340.31971967977 889.2022488946817
0.004403509178041219 395340.13481310755 889.202040948071
0.0044888072650024705 395339.94246653345 889.2018246343554
0.004575757616860898 395339.74307342584 889.2016003960248
0.004664392238780825 395339.53584004537 889.2013673404302
0.004754743755881992 395339.32109384355 889.2011258358185
0.004846845425248354 395339.0977926529 889.2008747101556
0.004940731148169567 395338.8660994784 889.2006141467497
0.005036435482619529 395338.62559257715 889.2003436712979
0.0051339936559767 395338.37599676056 889.2000629743123
0.005233441577990754 395338.1117306464 889.1997657789237
0.005334815854000479 395337.8476007618 889.1994687366405
0.005438153798407629 395337.568195957 889.1991545159689
0.005543493448411856 395337.2781500289 889.1988283280955
0.005650873578011608 395336.97713930096 889.1984898089976
0.005760333712276324 395336.66453743994 889.1981382542813
0.0058719141418949895 395336.34004850127 889.1977733311091
0.005985655938006604 395336.003094789 889.1973943897823
0.006101600967317825 395335.65330000856 889.197001007098
0.006219791907513542 395335.2903620403 889.1965928432703
0.006340272262965867 395334.9137392788 889.1961692891831
0.006463086380747511 395334.5188916948 889.1957252390441
0.0065882794669552495 395334.1161513932 889.1952723124356
0.006715897603349684 395333.694564126 889.1947981900547
0.006845987764317215 395333.2568426564 889.1943059226778
0.006978597834160691 395332.80259103933 889.1937950649896
0.007113776624724862 395332.33070145064 889.1932643710823
0.007251573893363381 395331.84105066303 889.1927137023368
0.0073920403612537005 395331.33279510983 889.1921421100277
0.007535227732066859 395330.8050915013 889.1915486457361
0.007681188710998797 395330.25730957283 889.1909326006117
0.00782997702417044 395329.68385719176 889.1902876855906
0.007981647438403431 395329.09848708723 889.1896293671978
0.008136255781379083 395328.4856885655 889.188940201761
0.00829385896218766 395327.84948766243 889.1882247169746
0.008454514992275865 395327.18921199837 889.1874821566016
0.008618283006799949 395326.5037320485 889.1867112502846
0.008785223286392531 395325.79231518554 889.1859111740193
0.008955397279351051 395325.05347227404 889.1850802530079
0.009128867624255937 395324.2871272406 889.1842184016095
0.009305698173027031 395323.49134305667 889.1833234412987
0.009485954014426482 395322.66409579106 889.1823930958046
0.009669701498017045 395321.807993483 889.1814302980949
0.00985700825858433 395320.917839364 889.1804292036167
0.01004794324103225 395319.99407788034 889.1793903120791
0.010242576725760566 395319.03516200837 889.1783118835145
0.010440980354534174 395318.03993235744 889.1771926138878
0.010643227156853305 395317.00677333446 889.1760306860891
0.010849391576834712 395315.9344630168 889.1748247257306
0.011059549500613395 395314.82141066785 889.173572943627
0.011273778284275263 395313.6661757893 889.1722737195412
0.011492156782330704 395312.46974842023 889.1709281667055
0.011714765376739872 395311.22269215726 889.1695256723065
0.011941686006500045 395309.930987096 889.16807296157
0.012173002197806261 395308.5904255787 889.1665653021133
0.012408799094796008 395307.1992256636 889.1650006895949
0.012649163490889645 395305.75477296225 889.1633761834348
0.012894183860737707 395304.25631537114 889.1616909374483
0.013143950392787238 395302.7009592956 889.1599416969881
0.013398555022478746 395301.08666403417 889.1581261665825
0.013658091466086422 395299.4113346571 889.1562419897384
0.013922655255213617 395297.67223582504 889.1542860896809
0.01419234377195677 395295.8683697294 889.1522573437345
0.014467256284750226 395293.9959579562 889.1501515019341
0.01474749398490561 395292.0527006823 889.1479659771846
0.015033160023858794 395290.03630198963 889.1456981867366
0.015324359551138532 395287.9435888569 889.1433445613333
0.01562119975307038 395285.7721375088 889.1409023743186
0.015923789892230544 395283.5188103196 889.1383680961244
0.016232241347663722 395281.1804923425 889.1357382226208
0.016546667655880257 395278.7564583501 889.133011937303
0.016867184552647126 395276.23638067517 889.1301776238114
0.017193910015588763 395273.62391294073 889.1272393903369
0.017526964307612743 395270.91341062804 889.1241908874463
0.01786647002117699 395268.1007750332 889.1210275041674
0.01821255212341411 395265.1824185094 889.1177452042101
0.018565338002130144 395262.1546909335 889.1143398809104
0.01892495751269404 395259.0133082113 889.1108067144514
0.01929154302583559 395255.75387677556 889.1071407617594
0.01966522947636923 395252.3744340992 889.1033398138815
0.020046154412861464 395248.8639736766 889.0993914896991
0.020434458048260638 395245.2244832202 889.0952980229061
0.020830283311507148 395241.4484880044 889.0910510043439
0.021233775900143675 395237.5308750599 889.0866446810006
0.021645084333944156 395233.4679074423 889.0820748473589
0.022064360009581955 395229.25243985344 889.0773334641407
0.022491757256356597 395224.8796613937 889.0724151174568
0.022927433393000324 395220.34290649276 889.0673123071084
0.02337154878558465 395215.63867720903 889.0620210955016
0.023824266906548985 395210.75808482774 889.0565314813538
0.024285754394872294 395205.6958121861 889.0508374802715
0.02475618111741066 395200.44534859573 889.0449317650888
0.02523572023142261 395194.999597509 889.0388063493168
0.02572454824830603 395189.35088721616 889.0324525991344
0.026222845098569206 395183.49325089995 889.0258637980112
0.02673079419806087 395177.4180610325 889.0190302361727
0.02724858251548269 395171.11807867617 889.0119437652975
0.027776400641209977 395164.583430603 889.0045932734015
0.028314442857444982 395157.807607974 888.9969714323823
0.02886290720972963 395150.78123536584 888.989067689098
0.029421995579842937 395143.4957348856 888.9808723868986
0.029991913760111042 395135.9396891981 888.9723726744247
0.03057287152915611 395128.10580882843 888.9635603429742
0.03116508272911307 395119.9831889739 888.9544231162516
0.03176876534434154 395111.5612968362 888.9449491355875
0.032384141581662966 395102.8294174517 888.935126336508
0.0330114379521515 395093.7767645322 888.9249425733673
0.03365088535450964 395084.3917682169 888.9143848180396
0.03430271916005854 395074.662267504 888.9034393762959
0.034967179299374915 395064.5758112148 888.892092226289
0.03564451035060586 395054.1205914363 888.8803300686052
0.03633496162949472 395043.2827038806 888.8681372440803
0.03703878728115034 395032.0448515866 888.8554942751792
0.03775624637359452 395020.4044770796 888.8423982653838
0.038487602993120824 395008.33609368006 888.828820520217
0.03923312634150132 394995.8277947932 888.8147476215651
0.03999309083507546 394982.86421740963 888.8001622607972
0.0407677762057591 394969.43029014394 888.7850474553945
0.04155746760400947 394955.5083299714 888.7693832822679
0.0423624557037852 394941.0811299385 888.7531503515905
0.043183036809538944 394926.13156163716 888.7363293594306
0.044019512965282864 394910.64230476366 888.7189007833283
0.04487219206576678 394894.59363511315 888.7008423931117
0.045741387969809456 394877.9640243166 888.6821299253368
0.04662742061582586 394860.73651650164 888.6627442584747
0.047530616139591454 394842.8895662899 888.6426611032017
0.048451306994288484 394824.4008472118 888.6218552874016
0.04938983207287671 394805.2456675436 888.6002989731026
0.050346536832835465 394785.4096053337 888.5779758753125
0.05132177342332109 394764.86162811035 888.5548510115854
0.05231590081478837 394743.5784433602 888.5308981046863
0.05332928493112198 394721.53588072007 888.5060898842731
0.0543622987843279 394698.70885356 888.4803980432657
0.05541532261183404 394675.0697584621 888.4537914359554
0.05648874401644955 394650.5911910595 888.4262391341889
0.05758295810903618 394625.2422699913 888.39770628924
0.058698367653942644 394598.99824066175 888.3681649413849
0.05983538321725685 394571.826076297 888.3375778118327
0.060994423317928564 394543.69341880555 888.3059083658123
0.06217591458182065 394514.5710383592 888.2731235812093
0.06338029189874368 394484.42530851736 888.2391854771071
0.06460799858253301 394453.2181899986 888.2040510941149
0.06585948653422544 394420.91949994484 888.1676863069775
0.06713521640839795 394387.49137591047 888.1300483329122
0.06843565778272771 394352.89580107125 888.0910942026964
0.06976128933083739 394317.0923093107 888.0507781757873
0.07111259899848706 394280.0441153992 888.0090586423082
0.07249008418318062 394241.71136355953 887.9658905200802
0.07389425191725027 394202.05071534956 887.9212247889444
0.07532561905448837 394161.0130173827 887.8750058621796
0.07678471246039294 394118.56280502334 887.8271935517895
0.0782720692060998 394074.6515675223 887.777732957436
0.07978823676607034 394029.2292357626 887.7265674021056
0.08133377321960965 393982.246977686 887.6736415797036
0.08290924745628635 393933.6603707958 887.618905128542
0.08451523938533335 393883.4152488557 887.5622966855406
0.08615234014910367 393831.4544422257 887.5037514762691
0.08782115234066194 393777.7278402681 887.4432126511173
0.08952229022558936 393722.1833530802 887.3806211013177
0.09125637996808646 393664.75761328265 887.3159049777961
0.09302405986145482 393605.3895670788 887.2489950031826
0.09482598056304424 393544.0253030542 887.1798299139292
0.09666280533374925 393480.60422408086 887.1083408739667
0.09853521028214651 393415.05508560955 887.0344470037334
0.10044388461336043 393347.3113296619 886.9580726614555
0.1023895308827494 393277.31557145045 886.8791525021326
0.10437286525450692 393204.99479295796 886.7976035070889
0.1063946177652703 393130.2718695957 886.713337972984
0.1084555325928371 393053.07749567856 886.6262769574096
0.11055636833008595 392973.34516732185 886.536344621383
0.11269789826420434 392890.9919869442 886.4434465739416
0.11488091066132325 392805.93224734365 886.3474851855153
0.11710620905666728 392718.09715259896 886.2483818350237
0.11937461255032389 392627.4062287765 886.1460446549163
0.12168695610874317 392533.76554587774 886.0403665137134
0.12404409087207555 392437.0868669755 885.9312466179026
0.12644688446746505 392337.29574055114 885.8185996472993
0.12889622132840972 392234.29815088 885.7023181079295
0.13139300302030998 392127.99104055896 885.5822841956121
0.13393814857232 392018.28133880376 885.458391274038
0.13653259481563032 391905.09053680865 885.330549045732
0.13917729672830095 391788.30729723076 885.198630022924
0.14187322778677644 391667.8196715641 885.0625058961249
0.14462138032420652 391543.53889286495 884.9220744143124
0.14742276589571074 391415.36810667405 884.777224058886
0.1502784156507167 391283.18541682756 884.6278148654693
0.1531893807125124 391146.8689984193 884.4737067866057
0.15615673256514762 391006.335495757 884.314803105497
0.15918156344783302 390861.46565360646 884.1509663554143
0.1622649867569767 390712.1259569579 883.9820427553468
0.16540813745601093 390558.1943270519 883.8078912603711
0.1686121724931533 390399.58233745384 883.6284087074769
0.17187827122726448 390236.1518908438 883.4434355303613
0.17520763586195193 390067.7565341237 883.2528024683802
0.17860149188808536 389894.28467701515 883.0563794877597
0.18206108853487976 389715.63049184106 882.8540428540168
0.18558769922972046 389531.63925601036 882.6456132061275
0.1891826220668919 389342.1529763519 882.4309071835052
0.1928471802853894 389147.08187461697 882.2098184384676
0.1965827227559818 388946.28524745133 881.9821826402746
0.20039062447771253 388739.59676684305 881.7478060838519
0.2042722870840164 388526.8606938065 881.5065067188177
0.20822913935864204 388307.9890467382 881.2581790221731
0.21226263776156387 388082.8178096794 881.0026308810654
0.2163742669650859 387851.16548312525 880.739649934219
0.22056554040032733 387612.8937466937 880.4690724229825
0.22483800081429436 387367.9040077691 880.1907793288556
0.22919322083774477 387116.0082951937 879.9045497043344
0.23363280356404773 386857.01399731 879.6101568277961
0.23815838313926116 386590.820807052 879.3074784249842
0.2427716253636362 386317.2936886512 878.996352311716
0.24747422830477603 386036.23030000774 878.6765392338727
0.2522679229226672 385747.4396704265 878.347812282158
0.25715447370682265 385450.85275782977 878.010082809793
0.2621356793257639 385146.29295426595 877.6631391989366
0.2672133732890869 384833.5450542173 877.3067252155512
0.27238942462234694 384512.44593893195 876.9406433036753
0.27766573855502247 384182.93204408675 876.5648088351331
0.283044257221801 383844.79690938274 876.1789736228355
0.28852696037745307 383497.80086414184 875.782850784533
0.2941158661255479 383141.8585047074 875.3763287920315
0.2998130316612904 382776.8655343603 874.959273948634
0.3056205540287436 382402.5919635801 874.5314081993627
0.31154057089272275 382018.8155582077 874.0924614229409
0.31757526132563507 381625.5145910985 873.6423920473393
0.323726846609567 381222.53034495824 873.1810011045342
0.329997591053905 380809.6242412957 872.7079972605908
0.3363898028287972 380386.61388688645 872.2231525095931
0.34290583481475395 379953.516655891 871.726467024939
0.34954808546871124 379510.1219134399 871.2176787846306
0.3563189997068667 379056.18532466 870.696485952091
0.3632210698046208 378591.6156533445 870.1627613881722
0.370256836313943 378116.4135007217 869.6164827103057
0.37742888899851457 377630.3377616095 869.0573488114688
0.3847398677869823 377133.14721933386 868.4850571188128
0.39219246374468036 376624.8943149975 867.8996420266544
0.39978942006416823 376105.49515071214 867.3009802262559
0.40753353307496315 375574.70389592944 866.6887606239387
0.4154276532728266 375032.3352654512 866.0627405280188
0.42347468636899444 374478.543393439 865.4230680926398
0.4316775943597212 373913.1671852935 864.7695267356424
0.4400393966165523 373335.9640410098 864.1018042349059
0.44856317099770837 372746.86678922654 863.4197898927572
0.4572520549810044 372146.0544940639 862.7236573713091
0.4661092468187042 371533.3116945839 862.0131225156423
0.47513800671475437 370908.40235105803 861.2878756270263
0.4843416580248163 370271.47168270074 860.548048260759
0.493723588479549 369622.6065600653 859.7937038151248
0.5032872514315779 368961.5791790236 859.0245388567472
0.5130361671266278 368288.2110562679 858.2403055744561
0.5229739239992733 367602.8460617181 857.4413636648493
0.5331041799937886 366905.43796614127 856.6276180069625
0.5434306639105875 366195.7755698774 855.7987795853385
0.5539571767787355 365473.8255414902 854.9547655186094
0.5646875932550566 364740.0228054739 854.0960400393785
0.5756258630503367 363994.22975980927 853.2223974554457
0.5867760123831596 363236.2503747771 852.3335619049353
0.5981421454618925 362466.3105002558 851.429751066118
0.6097284459953902 361684.76735719194 850.5113372050862
0.6215391787329555 360891.4489761603 849.5780705457978
0.6335786910341367 360086.202901115 848.6297224362519
0.6458514144689197 359269.56532935775 847.6668748150511
0.6583618664489304 358441.68862449925 846.6896581682089
0.6711146518902246 357602.4173630454 845.6978389035239
0.6841144649082979 356751.7519253697 844.6913660330259
0.6973660905459123 355890.38571082236 843.67100899678
0.7108744065344067 355018.31503346143 842.636712983076
0.7246443850891169 354135.408532641 841.5882705131304
0.7386810947395809 353241.9289881096 840.525941286894
0.7529897021951829 352338.5416594668 839.450465077561
0.7675754742469504 351425.15639478463 838.3616837556266
0.7824437797061797 350501.66984343465 837.2594219755722
0.7976000913806248 349568.7180534428 836.1443871167739
0.813049988088948 348626.70244639454 835.0170087446058
0.82879915671421 347675.53549333266 833.8771318285837
0.8448533942971241 346715.2299365435 832.7247203446569
0.8612186101698714 345746.63592518837 831.5607445342623
0.877900828131229 344769.909337475 830.3853434851497
0.8949061886638501 343784.98189706885 829.1983862708234
0.9122409511944829 342792.08770836564 828.000105927971
0.9299114963979833 341792.0600932984 826.7914611234182
0.9479243285459364 340784.9064051859 825.572415243128
0.9662860779007927 339770.5760305829 824.3428607449487
0.9850035031563675 338749.61448635446 823.103413291859
1.0040834939256265 337722.6061771533 821.8547392053578
1.0235330732766403 336689.5152120147 820.5967526282501
1.0433594003176818 335650.3508051664 819.3294219118052
1.0635697728323859 334605.8765158863 818.0536370139629
1.0841716299659676 333556.3597282981 816.7696856866054
1.105172554963449 332501.7766785546 815.477500215125
1.1265802779609513 331442.26786464703 814.1772139585423
1.1484026788310444 330378.55971231486 812.8696816985056
1.1706477900832102 329310.7266986873 811.5549601828423
1.1933237998205029 328238.7530555449 810.232994953359
1.2164390547534567 327162.9395539386 808.9041223209813
1.2400020632723998 326083.79177606414 807.5689342416091
1.2640214985792622 325001.30754796934 806.2273966418771
1.288506201880067 323915.48313644814 804.8794731342676
1.3134651856392365 322826.72518123087 803.5256376510097
1.338907636896957 321735.25616801827 802.1661376149186
1.3648429206507917 320641.06915921636 800.8009355129605
1.3912805833028115 319544.2011012517 799.4300483485115
1.4182303561734697 318444.96471873013 798.053838683494
1.445702159083573 317343.41796854156 796.6723516836034
1.473706104005615 316239.55625714967 795.2855540711772
1.5022524987858619 315133.424978155 793.8934751944432
1.5313518509385011 314025.1437000663 792.4962380984106
1.5610148715133183 312914.7138155995 791.0938172120921
1.5912524790382725 311802.1314661831 789.6861800312616
1.6220758035384604 310687.36446850194 788.2732577837484
1.6534961906329002 309570.3815770915 786.8549822897374
1.6855252057106997 308451.1805428766 785.4313216862141
1.7181746381881038 307329.74173086317 784.0022215923411
1.75145650584802 306205.8720674508 782.5674054897135
1.7853830592635707 305079.5075732105 781.1267599733228
1.8199667863073643 303950.642268669 779.6802450603311
1.8552204167480941 302819.18648416846 778.2277127989834
1.8911569269361956 301684.81046075525 776.7687049061069
1.9277895445802313 300547.47662084823 775.3031363548689
1.965131753615833 299407.1765824862 773.8309590375487
2.003197299168935 298263.6650198793 772.3518175286173
2.0420001926151734 297116.60670481634 770.8652368667514
2.081554716737244 295965.990648246 769.3711596469496
2.1218754309822057 294811.78986534075 767.8695069676106
2.1629771768205925 293653.5729932795 766.3596714249511
2.2048750832093615 292491.12981620873 764.8413297099062
2.247584572160616 291324.4445295166 763.314410357248
2.291121364418243 290153.4212128414 761.7787358713045
2.3355014852444813 288977.56806180027 760.2336062840162
2.3807412703186075 287796.7873773503 758.6788350512361
2.4268573717498345 286611.0583064864 757.1143352314582
2.473866764206731 285420.17359470954 755.5397720765063
2.5217867511653385 284223.7222043956 753.9545373620291
2.5706349712783547 283021.6661844168 752.3585131895787
2.6204294048676324 281813.97686426283 750.7515925581015
2.671188380542496 280600.3322952061 749.133275586135
2.7229305819462324 279380.49175651843 747.5031662227505
2.775675054633264 278154.4240432057 745.8611453121897
2.829441213079564 276922.06783167605 744.207051608188
2.884248847828807 275683.0969212394 742.5403651266905
2.9401181327770085 274437.39966234623 740.860850176801
2.997069632598228 273184.93951845693 739.1683698839621
3.055124310314153 271925.61745122506 737.4627006855669
3.1143035350102504 270659.2000096664 735.7434335550217
3.174629089701433 269385.6328638505 734.0103989234083
3.2361231793500687 268104.8755996592 732.263443850175
3.298808439039326 266816.8212142531 730.5023219870736
3.36270794230478 265521.3529505952 728.7267704024537
3.4278452096274767 264218.4280630998 726.9366245596651
3.494244217091465 262908.002993888 725.131716302477
3.561929405209073 261590.02145051002 723.3118572932564
3.630925687917068 260264.43485143565 721.4768670601097
3.701258461747137 258931.19918667743 719.6265686961224
3.772953615173977 257590.28661334512 717.7608050225996
3.8460375381444862 256241.73155636157 715.8795032075741
3.920537131791486 254885.50087890003 713.9824940135437
3.9964798183356613 253521.55081534912 712.0695904409191
4.073893551179271 252149.92093874433 710.1407197714328
4.152806825195419 250770.7129902133 708.1958952016219
4.23324868721656 249383.88778086312 706.234929440428
4.315248746726244 247989.40684484883 704.2576330361621
4.398837186757928 246587.4473351484 702.2641203068093
4.484044775004943 245178.09926338022 700.254381297797
4.570902875145584 243761.3237198711 698.2282201685508
4.659443458387664 242337.12439441116 696.1854988354917
4.7496991152366315 240905.79951907985 694.1265007461966
4.841703067491699 239467.38575499176 692.0511335948981
4.935489180474254 238021.84868607984 689.9591997880452
5.031091975493225 236569.29835576983 687.8507081566025
5.128546642551861 235110.0862190687 685.726018492909
5.227889053300698 233644.20223896494 683.5849650759808
5.329155774241338 232171.6159446276 681.4273489442988
5.432384080186079 230692.6329237993 679.2534621535606
5.537611967978226 229207.5348788271 677.0635640452484
5.644878170478208 227716.2966033782 674.8574613996325
5.754222170820528 226218.92897016008 672.6350109385626
5.865684216946969 224715.9319297523 670.3967958302788
5.979305336421272 223207.45735296293 668.1428849474682
6.095127351530772 221693.48520653506 665.8730888187854
6.213192894680626 220174.14228618407 663.5874355142419
6.33354542408611 218649.99797541773 661.286621633037
6.45622923976901 217121.09864095843 658.9705587368201
6.581289499863784 215587.4294301762 656.6390628498676
6.708772237239682 214049.32352776825 654.2924782202043
6.838724376444688 212507.27471184646 651.9313993233436
6.9711937509778235 210961.2770203692 649.5556589244208
7.1062291208959385 209411.33757936707 647.1651065676626
7.243880190761637 207858.0440370843 644.7604889214045
7.384197627938719 206301.68608056562 642.3420990104348
7.5272330812421675 204742.25651667136 639.9097694467109
7.6730391999492875 203179.86378109604 637.4635107691985
7.821669653179213 201615.22018795123 635.004283746104
7.973179149647635 200048.4485667352 632.5321313051776
8.12762345780435 198479.54626701458 630.0468970910254
8.285059426360785 196908.80743802185 627.5488944106615
8.44554500521525 195336.88657093095 625.03901729561
8.60913926678336 193763.80667199165 622.5171590759433
8.775902427741801 192189.56960460145 619.9831765533665
8.945895871193171 190614.7382866516 617.4378321526008
9.119182169260299 189039.73947420984 614.8816788199333
9.295825106118013 187464.57340598523 612.3145815771256
9.475889701471257 185889.30732027412 609.7365124712053
9.65944223448785 184314.67717483029 607.1485438915756
9.846550268194974 182740.88463994383 604.5508822918777
10.037282674348033 181167.92958712458 601.9434019691961
10.23170965878141 179596.02933468146 599.326337373357
10.429902787250189 178025.9104699457 596.7007800731379
10.631935011772553 176457.62926794525 594.0667121930755
10.837880697482241 174891.18373845576 591.4240166554885
11.047815650001354 173327.018570085 588.7733325654025
11.261817143343267 171765.66839030644 586.1154636934713
11.479963948356165 170207.12934868413 583.4503052508999
11.702336361717308 168651.42522615052 580.7777978300315
11.929016235489174 167099.20760464313 578.0989666218807
12.160087007247961 165550.7419771889 575.4141847003581
12.395633730795852 164006.01899983955 572.7233520642222
12.635743107467942 162465.16480233535 570.0266042955107
12.880503518045836 160928.86651940647 567.3250682270376
13.130005055289287 159397.20999661225 564.6188271685814
13.384339557098134 157870.18112765028 561.9077880358134
13.643600640316354 156348.07577845477 559.1924101388622
13.907883735191163 154831.45107329809 556.4736311332247
14.177286120499481 153320.30003127464 553.7513883165886
14.451906959354814 151814.60934269583 551.0256061975629
14.73184733570788 150314.87437507804 548.2971354568215
15.017210291553955 148821.39470559277 545.5664848679633
15.30810086486127 147334.14435066227 542.8335736681405
15.604626128233905 145853.16652514657 540.0984475540483
15.906895228323771 144379.0225810154 537.3621173492144
16.215019426005746 142911.8153112244 534.6247568364645
16.52911213733129 141451.51193228812 531.8862884720533
16.849288975275215 139998.26676075408 529.1469866884892
17.175667792291257 138552.5669018287 526.4077638140013
17.508368723691653 137114.4038750617 523.6686048925632
17.84751423186736 135683.73851456976 520.9294357483934
18.19322915136461 134260.87610091077 518.1908453473696
18.545640734834834 132846.11295357667 515.4534177858882
18.904878699874306 131439.4026800781 512.7170812057623
19.271075276771455 130040.72706953257 509.98181745927485
19.644365257178848 128650.48005761475 507.24842051526343
20.024886043728202 127268.7640269009 504.5171236477527
20.41277770060604 125895.52661032462 501.78785678875215
20.80818300510948 124530.81319384705 499.06074418621034
21.21124750020039 123174.9853635552 496.3365498601835
21.62211954807781 121828.03224093428 493.61530008891395
22.04095038478769 120489.89568210636 490.89692539698467
22.46789417589082 119160.71436654555 488.1817578864363
22.903108073208895 117840.73144394753 485.4703522233825
23.346752272669978 116529.88680217743 482.7626472754027
23.798990073274027 115228.12404970352 480.0585881946151
24.25998793720109 113935.6697258292 477.3587115070368
24.729915551083494 112652.60914639568 474.6632683205973
25.208945888465273 111378.87406615014 471.9721899988391
25.697255273470894 110114.44011145433 469.28549969385233
26.195023445707825 108859.53448729882 466.6037601376543
26.702433626425965 107614.14072285747 463.92702168090506
27.219672585958968 106378.18669538641 461.255215028267
27.746930712471336 105151.7011805273 458.5884891283847
28.28440208203781 103934.84602446976 455.92728811614194
28.832284530079786 102727.5552927867 453.2715638395744
29.390779724185958 101529.75429540084 450.6212473805487
29.96009323834284 100341.53609084693 447.9766424510254
30.54043462860384 99162.9590851989 445.3379819534797
31.132017510223722 97993.94506996394 442.7051955194652
31.73505963628723 96834.43245607147 440.0782486241997
32.34978297786109 95684.53560815974 437.45750789798944
32.976413805697995 94544.23057376622 434.8430304690791
33.615182773523784 93413.43734515454 432.2347448902148
34.26632500293767 92292.12426576971 429.63269025010123
34.93008016995731 91180.38054581532 427.03718935431215
35.606692593239686 90078.13937318057 424.4482050219569
36.29641132401156 88985.31951159124 421.86566466492917
36.999490237741426 87901.92879503765 419.28970603876655
37.71618812758774 86827.99909731786 416.72052768568494
38.44676879965645 85763.44824593648 414.15805737891054
39.19150117010451 84708.1989778001 411.6022326902518
39.95065936412386 83662.29117606454 409.05327568927873
40.72452281684326 82625.69552499293 406.51124344842646
41.51337637618385 81598.32907012566 403.97606134553484
42.31751040770782 80580.13140571202 401.4477086887208
43.137220901497635 79571.14061803624 398.9264108028854
43.97280958110602 78571.29177392513 396.4121384971079
44.82458401461552 77580.5020581609 393.9048160613447
45.692857727850196 76598.73711791505 391.4044892893158
46.577950319779745 75626.0084471689 388.9113226615263
47.480187580159644 74662.2353560138 386.425245955835
48.39990160944906 73707.3361770862 383.94618419014455
49.337430941052595 72761.3073229396 381.4742647229026
50.293120665929344 71824.11862332898 379.00954769854803
51.26732255961625 70895.68746903112 376.55195516430695
52.260395211710986 69975.94171854446 374.10143468996336
53.27270415786403 69064.889584837 371.65814826218195
54.30462201432675 68162.47060509025 369.22207573516033
55.3565286151066 67268.60298061199 366.79313783278985
56.42881115177774 66383.23260697303 364.37132874849806
57.521864316001114 65506.35842340624 361.9567886458444
58.636090444804566 64637.90416722519 359.54945186225825
59.77189966867765 63777.788992912814 357.1492376945884
60.929710062534085 62925.98380169825 354.7562086890045
62.10994779959951 62082.46204751895 352.37043589812964
63.313047308279394 61247.143184497574 349.9918375748142
64.53945143206641 60419.95125606405 347.62034248894025
65.78961159254392 59600.87570674124 345.25606643979836
67.0639879555482 58789.86221621017 342.8990003374468
68.36304960054869 57986.8314286334 340.5490608668108
69.68727469330925 57191.71989729898 338.2062089829191
71.0371506618946 56404.51813934852 335.8705647696699
72.41317437608463 55625.15589123773 333.54206898452173
73.81585233026557 54853.55506810859 331.22063664001547
75.24570082986297 54089.670813612756 328.9062809178711
76.70324618138676 53333.47919023425 326.59907896451347
78.18902488615623 52584.90352775632 324.2989470465679
79.70358383777896 51843.8681307538 322.00580159603896
81.2474805234541 51110.34850038462 319.7197163153521
82.82128322917603 50384.29710406612 317.4406940014658
84.42557124891113 49665.637905489384 315.1686466179318
86.06093509782825 48954.30343830581 312.90351048943444
87.72797672965838 48250.276150273596 310.6453802980936
89.42730975826584 47553.49244049666 308.3942037084895
91.15955968350951 46863.877784592216 306.1498906894863
92.92536412148111 46181.37827577187 303.9124159219951
94.72537303920193 45505.969809033486 301.6818516551285
96.5602489938671 44837.5810469002 299.4581140891
98.4306673767222 44176.13907677911 297.24111114305543
100.33731666166564 43521.60761816952 295.03087166657497
102.28089865866559 42873.94542595237 292.8274079588602
104.26212877208687 42233.08040917397 290.63062608463673
106.28173626401974 41598.94465956873 288.4404432792625
108.34046452271143 40971.51275734812 286.2569222127148
110.43907133619622 40350.727842535634 284.0800163423525
112.57832917122715 39736.51949802133 281.9096291296958
114.7590254576084 39128.82946016195 279.7457040247873
116.9819628780378 38527.631532527215 277.58829778118246
119.24795966356213 37932.86049390614 275.4373267874423
121.55784989475713 37344.447641356644 273.29269160135493
123.91248380873901 36762.34799767534 271.15437668485214
126.31272811212537 36186.52526932531 269.02239783826667
128.7594663000569 35616.9119697948 266.89665404345106
131.25359898140007 35053.442255884 264.7770467993176
133.79604421024672 34496.08380538052 262.66360161004616
136.38773782383748 33944.786204073236 260.5562749352747
139.02963378702938 33399.48354915009 258.4549614503467
141.72270454343752 32860.116662786364 256.3595781818435
144.4679413733759 32326.656528546187 254.2701576219521
147.26635475873448 31799.04431273421 252.18661468338962
150.1189747549222 31277.215905506517 250.10883992976545
153.02685137001666 30761.121920011774 248.03677920829313
155.9910549512553 30250.728554559006 245.97043950263213
159.0126765790156 29745.973763771424 243.90971183522572
162.0928284684258 29246.795415121644 241.85448275821412
165.232644378755 28753.15583401408 239.80473654210454
168.4332800307344 28265.010607071657 237.76042819221055
171.69591353195992 27782.298581049705 235.72143975909236
175.02174581053734 27304.962075451265 233.68766366862957
178.41200105712622 26832.969373808715 231.6591002909608
181.86792717554798 26366.26793005743 229.63565894720023
185.39079624211905 25904.79842733674 227.61721563773133
188.98190497388535 25448.510315295585 225.60368044557956
192.64257520592403 24997.370914848132 223.5950398146083
196.37415437789355 24551.32399852366 221.59117310273737
200.1780160300043 24110.31208249707 219.5919492262732
204.05556030860038 23674.293807719958 217.59730608497873
208.0082144815324 23243.230114335387 215.6071896497674
212.03743346351573 22817.06472123636 213.62146297240994
216.14470035166102 22395.743007119723 211.63999152863204
220.33152697138243 21979.23116300971 209.66273470986545
224.59945443287788 21567.482614182 207.6895886373797
228.95005369839097 21160.44288635793 205.72040679698225
233.38492616045693 20758.062516801347 203.7550613692889
237.90570423135364 20360.309446072573 201.79350557474626
242.514051943968 19967.13328557941 199.83559885855877
247.21166556430452 19578.4813689338 197.88118338504952
252.00027421585293 19194.311129965143 195.93014637857618
256.88164051605577 18814.587481317023 193.98240889996714
261.85756122510105 18439.259265700315 192.03780495360968
266.929867907286 18068.275843404674 190.09616431377395
272.1004276051859 17701.602118691542 188.15739219436233
277.37114352688775 17339.197046095884 186.22135777668404
282.7439557465324 16981.01093200473 184.28787769142457
288.2208419184298 16626.99679477632 182.35677555153427
293.80381800500186 16277.122948140543 180.42795209246566
299.494939018831 15931.344599185637 178.5012302432991
305.29629977907933 15589.61380773492 176.57640730139983
311.2100356825638 15251.888774570418 174.65330672260643
317.238323489761 14918.137270786601 172.7317994509789
323.38338212604384 14588.313334677634 170.81167017904622
329.6474734984341 14262.370750347061 168.8926922655155
336.03290332818085 13940.274030238805 166.97469287433233
342.5420219994585 13621.987203844414 165.05748819029336
349.17722542451213 13307.465206067338 163.14082999707546
355.94095592555743 12996.664391850512 161.22446707525822
362.8357031337631 12689.553710839866 159.3082151732287
369.86400490565114 12386.093624417415 157.39182713481293
377.0284482572406 12086.24073834132 155.47501881872418
384.33167031629375 11789.955251598676 153.55751529377306
391.77635929300277 11497.207080291173 151.63909179556023
399.36525546948445 11207.95529179619 149.7194395647819
407.10115220843335 10922.158127592267 147.79822818689178
414.9868969813213 10639.780578937527 145.8751560680401
423.02539241650885 10360.790849263016 143.94992774755406
431.2195973676641 10085.148373834718 142.02216991607133
439.5725280028687 9812.813256042813 140.0914933608948
448.08725891482885 9543.755227128295 138.15755663103118
456.76692425258693 9277.939532349144 136.21996573446305
465.61471887515967 9015.327099372033 134.27827150639104
474.633899527514 8755.880908645226 132.33201357680028
483.8277860393316 8499.572853930396 130.3807720020893
493.1997625469878 8246.366754764766 128.42403789606342
502.7532787392073 7996.2250702521715 126.46125944535086
512.4918511268387 7749.114516634245 124.49188340317006
522.419064337235 7505.006792656467 122.51536060965145
532.5385724337011 7263.865718718876 120.53103931119881
542.8541002605067 7025.655269484624 118.53822395737693
553.3694448139356 6790.34641725924 116.53622970784014
564.0884766399196 6557.908828087281 114.52431032830786
575.0151412587231 6328.30746641122 112.50162191196374
586.1534606172257 6101.5084318349245 110.46726602785935
597.5075345693559 5877.485326301455 108.42034242983605
609.0815423851814 5656.20632459448 106.3598262935257
620.8797442892301 5437.637813631376 104.28458959627137
632.9064830286344 5221.748693549844 102.1934312326369
645.1661854716137 5008.513156771355 100.08509536161071
657.6633642369576 4797.899157365984 97.95814572934692
670.4026193550486 4589.87446590014 95.81100631869118
683.3886399610581 4384.411288433064 93.64199152552304
696.6262060209649 4181.482968517964 91.4492533432391
710.1201900909844 3981.0582445363666 89.23069252825921
723.8755591110762 3783.1064739776 86.98398098474914
737.8973762332234 3587.60283217829 84.7065857200996
752.1908026850824 3394.519415866244 82.3956238627543
766.7610996697755 3203.8262640223434 80.0478140116561
781.6136303024563 3015.494900802912 77.65944760044218
796.7538615843839 2829.5018618854433 75.22634992986757
812.1873664152696 2645.8189440740675 72.7436450018016
827.9198256445545 2464.417454580929 70.2056615178709
843.9570301624711 2285.2715347063836 67.60579168542269
860.3048830315759 2108.357746235549 64.93624174889626
876.9694016595705 1933.648382571898 62.18759333776952
893.9567200142435 1761.1160249863244 59.348395513043556
911.2730908812849 1590.737612210621 56.40456740744708
928.9248881658343 1422.4888238041647 53.3383318787561
946.9186092386491 1256.3430407833912 50.12670028604299
965.260877327662 1092.2745138733844 46.73915946769656
983.958443955923 930.2618882659725 43.13378926702296
1003.0181914267248 770.2804607351229 39.24997989133556
1022.4471353568672 612.30476557495 34.99442142899208
1042.2524272590274 456.3111078229925 30.20963779402171
1062.4413571741175 302.2786929779054 24.58774869636931
1083.021356354629 150.1831083836225 17.33107661881526
1103.9999999999995 1.3843301847344052e-05 0.005261806124771997
 
%% Cell type:code id: tags:
 
``` python
```
 
%% Cell type:markdown id: tags:
 
# histograms of f(v) at different radii with chi2 (Eddington vs. Maxwellian) - already done, except for the Maxwellian case.
 
%% Cell type:code id: tags:
 
``` python
```
 
%% Cell type:code id: tags:
 
``` python
def maxw(v,sigma):
N = np.sqrt(32 * np.pi) * v**2 / sigma**3
return N * np.exp(- v**2 / 2. / sigma**2)
 
get_maxw = np.vectorize(maxw)
```
 
%% Cell type:code id: tags:
 
``` python
def eddingtong_from_file(path):
v = np.array([])
fv = np.array([])
files = open(path)
for line in files:
row = line.split(' ')
if np.isnan(float(row[3][:-1])):
continue
v = np.append(v,float(row[0]))
fv = np.append(fv,float(row[3][:-1]))
return v,fv
 
def fdv_plot_chi2_max_edd(ax, path, limmin, limmax,save=False,outname="/home/arturo/Pictures/ploto",width=None):
# get eddington fdv from file
v, fv = eddingtong_from_file(path)
# interpolation of eddington points
f = interp1d(v,fv)
# get data points from simulations errors included
# particles in the shell between limmin and limmax
shell_wc = (myhalo.dm.r>limmin)&(myhalo.dm.r<limmax)
if width == None:
width = limmin /6.
disc = (np.abs(myhalo.dm.pos3d[:,2])<width)
local_v_wc = myhalo.dm.v[shell_wc]
v_disc_in_the_shell = myhalo.dm.v[shell_wc&disc]
Ntot=len(local_v_wc)
# speed distribution in the shell
bins = np.linspace(0,local_v_wc.max(),30)
hist_wc, bins_wc = np.histogram(myhalo.dm.v[shell_wc],bins=bins,normed=True)
x = (bins_wc[1:]+bins_wc[:-1])/2
# considering informations from theory inside the bins
fv_the = np.array([])
N = np.array([])
test = np.array([])
means = np.array([])
N_disc = np.array([])
for i in range(len(bins_wc)-1):
try:
fv_the = np.append(fv_the, quad(f,bins_wc[i],bins_wc[i+1])[0] / (bins_wc[i+1] - bins_wc[i]))
test = np.append(test,(bins_wc[i+1]+ bins_wc[i])/2)
pop = local_v_wc[(local_v_wc>bins_wc[i])&(local_v_wc<=bins_wc[i+1])]
disc_in_the_bin = v_disc_in_the_shell[(v_disc_in_the_shell>bins_wc[i])&(v_disc_in_the_shell<=bins_wc[i+1])]
means, N = np.append(means,np.average(pop)), np.append(N,len(pop))
N_disc = np.append(N_disc,len(disc_in_the_bin))
except:
test = np.append(test,(bins_wc[i+1]+ bins_wc[i])/2)
fv_the = np.append(fv_the,0.0)
pop = local_v_wc[(local_v_wc>bins_wc[i])&(local_v_wc<=bins_wc[i+1])]
disc_in_the_bin = v_disc_in_the_shell[(v_disc_in_the_shell>bins_wc[i])&(v_disc_in_the_shell<=bins_wc[i+1])]
means, N = np.append(means,np.average(pop)), np.append(N,len(pop))
N_disc = np.append(N_disc,len(disc_in_the_bin))
 
dof = len(test)
# normalizations factor for data histogram
Ntot = np.trapz(N,x=means)
## MAXWELLIAN Fdv
# circular velocity
enclosed_m = (np.sum(mass[np.where(r2<limmin**2)])+ np.sum(mass[np.where(r2<limmax**2)]))/2.
vc = np.sqrt(myhalo.p.G * enclosed_m / ((limmin + limmax) /2 )) * myhalo.p.kpctokm
# for gas
vc_gas = np.mean(myhalo.gs.vphi[(myhalo.gs.R>limmin)&(myhalo.gs.R<limmax)&(np.abs(myhalo.gs.pos3d[:,2])<2.)])
std_vc_gas = np.std(myhalo.gs.vphi[(myhalo.gs.R>limmin)&(myhalo.gs.R<limmax)&(np.abs(myhalo.gs.pos3d[:,2])<2.)])
#for stars
vc_stars = np.nanmean(myhalo.st.vphi[(myhalo.st.R>limmin)&(myhalo.st.R<limmax)&(np.abs(myhalo.st.pos3d[:,2])<2.)])
print np.nanmean(myhalo.st.vphi[(myhalo.st.r<limmin)&(myhalo.st.r<limmax)])
std_vc_stars = np.nanstd(myhalo.st.vphi[(myhalo.st.r>limmin)&(myhalo.st.r<limmax)])
leno = len(myhalo.st.vphi[(myhalo.st.r<limmin)&(myhalo.st.r<limmax)])
# this sigma is a parameter in the fdv not an error of any kind
#sigma_8 = vc / np.sqrt(2.)
#sigma_8 = vc_gas / np.sqrt(2.)
sigma_8 = vc_stars / np.sqrt(2.)
print "sigma = {0:0.1f}, vc = {1:.1f} km/s".format(sigma_8,vc)
maxwellian = get_maxw(test,sigma_8)
# normalizarion factor
N0 = np.trapz(maxwellian, x=test)
maxwellian /= N0
## CHI2 s
#error in data
sigma = np.sqrt(N) / Ntot
sigma_disc = np.sqrt(N_disc) / Ntot
# eddington
chi2 = np.sum(((N/Ntot) - fv_the)**2 / sigma**2)
# maxwellian
chi2_max = np.sum(((N/Ntot) - maxwellian)**2 / sigma**2)
# v escape from potential
v_esc , sig_vesc = vesc_from_pot(limmin,limmax)
# plotting
ax.set_xlim([0,1.2*v_esc])
ax.set_ylim([0,1.3*fv_the.max()])
ax.plot(means,fv_the,'k-',label="eddington")#,yerr=means/np.sqrt(N),xerr=sigma_x)
ax.errorbar(test,N/Ntot,yerr=np.sqrt(N)/Ntot,xerr=((bins[1:]-bins[:-1])/2.)[:dof],c='g',lw=1.5, label="data")
ax.errorbar(test,N_disc/Ntot,yerr=np.sqrt(N_disc)/Ntot,xerr=((bins[1:]-bins[:-1])/2.)[:dof],c='gray',
lw=1.2, ls='None', label=r"in disc of $\Delta$z "+"= {0} kpc".format(2*width))
ax.set_xlabel(r"$\vec{v}$ [km/s]",fontsize=16)
ax.set_ylabel(r"$f(v)_{normalized} $ ",fontsize=18)
ax.axvline(x=v_esc,color='k',linestyle='--')
ax.axvspan(v_esc-sig_vesc, v_esc+sig_vesc, alpha=0.5, color='red',label=r"$v_{esc} \pm 1\sigma$ ")
texto_max = "Maxwellian\n"
texto_max += r"$\chi_{red}^2 = $"
texto_max += "{0:.2f}".format(chi2_max/dof)
texto_edd = "Eddington\n"
texto_edd += r"$\chi_{red}^2 = $"
texto_edd += "{0:.2f}".format(chi2/dof)
#texto_vc = r"$v_c = $"
#texto_vc += " {0:.2f} km/s".format(vc)
#texto_vc += "\n"
#texto_vc = r"$v_c ^{gas} =$ "
#texto_vc += "{0:.2f}".format(vc_gas)
#texto_vc += r" $\pm$ "
#texto_vc +="{0:.2f}".format(std_vc_gas)
texto_vc = r"$v_c ^{stars} =$ "
texto_vc += "{0:.2f}".format(vc_stars,leno)
texto_vc += "\n"
texto_vc += "dof = {0}".format(dof)
 
fig.text(0.15,0.79,texto_max,fontsize=16,color='r')
fig.text(0.15,0.71,texto_edd,fontsize=16,color='k')
fig.text(0.15,0.60,texto_vc,fontsize=16,color='k')
 
ax.set_title(str(limmin)+" kpc < r < "+str(limmax)+" kpc ", fontsize=18)
max_v = np.linspace(0,1.5*test.max(),50)
ax.plot(max_v,get_maxw(max_v,sigma_8)/N0,'r-',label="maxwellian")
legend = ax.legend(loc='upper right', ncol=1, shadow=False, fontsize=14)
frame = legend.get_frame()
if (save):
plt.savefig(outname+".pdf",dpi=300)
```
 
%% Cell type:code id: tags:
 
``` python
"""
fig, ax = plt.subplots(figsize=[9,8])
path = '/home/arturo/Downloads/fv_df0_model13_bar3_r3_rmax500_div0 (1).dat'
fdv_plot_chi2_max_edd(ax, path, 2.5, 3.5,width=0.8,save=True,
outname='/home/arturo/Documents/LAM/LAM2LUPM/report07-17/HaloB/fdv_MAX_EDD_3kpc_vcstars')
fig, ax = plt.subplots(figsize=[9,8])
path = '/home/arturo/Downloads/fv_df0_model13_bar3_r8_rmax500_div0 (1).dat'
fdv_plot_chi2_max_edd(ax, path, 7.5, 8.5,save=True,
outname='/home/arturo/Documents/LAM/LAM2LUPM/report07-17/HaloB/fdv_MAX_EDD_8kpc_vcstars')
fig, ax = plt.subplots(figsize=[9,8])
path = '/home/arturo/Downloads/fv_df0_model13_bar3_r20_rmax500_div0 (1).dat'
fdv_plot_chi2_max_edd(ax, path, 19.5, 20.5,save=True,
outname='/home/arturo/Documents/LAM/LAM2LUPM/report07-17/HaloB/fdv_MAX_EDD_20kpc_vcstars')
fig, ax = plt.subplots(figsize=[9,8])
path = '/home/arturo/Downloads/fv_df0_model13_bar3_r50_rmax500_div0 (1).dat'
fdv_plot_chi2_max_edd(ax, path, 45, 55,save=True,
outname='/home/arturo/Documents/LAM/LAM2LUPM/report07-17/HaloB/fdv_MAX_EDD_50kpc_vcstars')
"""
```
 
%% Output
 
"\nfig, ax = plt.subplots(figsize=[9,8])\npath = '/home/arturo/Downloads/fv_df0_model13_bar3_r3_rmax500_div0 (1).dat'\nfdv_plot_chi2_max_edd(ax, path, 2.5, 3.5,width=0.8,save=True,\n outname='/home/arturo/Documents/LAM/LAM2LUPM/report07-17/HaloB/fdv_MAX_EDD_3kpc_vcstars')\nfig, ax = plt.subplots(figsize=[9,8])\npath = '/home/arturo/Downloads/fv_df0_model13_bar3_r8_rmax500_div0 (1).dat'\nfdv_plot_chi2_max_edd(ax, path, 7.5, 8.5,save=True,\n outname='/home/arturo/Documents/LAM/LAM2LUPM/report07-17/HaloB/fdv_MAX_EDD_8kpc_vcstars')\nfig, ax = plt.subplots(figsize=[9,8])\npath = '/home/arturo/Downloads/fv_df0_model13_bar3_r20_rmax500_div0 (1).dat'\nfdv_plot_chi2_max_edd(ax, path, 19.5, 20.5,save=True,\n outname='/home/arturo/Documents/LAM/LAM2LUPM/report07-17/HaloB/fdv_MAX_EDD_20kpc_vcstars')\nfig, ax = plt.subplots(figsize=[9,8])\npath = '/home/arturo/Downloads/fv_df0_model13_bar3_r50_rmax500_div0 (1).dat'\nfdv_plot_chi2_max_edd(ax, path, 45, 55,save=True,\n outname='/home/arturo/Documents/LAM/LAM2LUPM/report07-17/HaloB/fdv_MAX_EDD_50kpc_vcstars')\n"
 
%% Cell type:code id: tags:
 
``` python
```
 
%% 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,35,50)
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')
```
 
%% Cell type:code id: tags:
 
``` python
vc_stars_vphi = gas_mass_r = stars_mass_r = std_stars_vphi = vc_gas_vphi = 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)
# for gas
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]))
gas_mass_r= np.append(gas_mass_r, np.sum(myhalo.gs.mass[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]))
stars_mass_r= np.append(stars_mass_r, np.sum(myhalo.st.mass[stars_condition]))
 
r_arrays = (r[1:] + r[:-1]) / 2.
```
 
%% Cell type:code id: tags:
 
``` python
fig, [ax,ax1,ax2] = plt.subplots(3,1,figsize=[8,10],sharex=True)
ax.set_xlim(1,32)
ax.set_ylim(10,350)
#ax.set_title("Rotation curve\nHalo B",fontsize=18)
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,300,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=1, shadow=False, fontsize=14)
frame = legend.get_frame()
 
ax1.set_ylim(0,400)
ax1.set_xlim(1,32)
ax1.plot(r_arrays,np.abs(vc_stars_vphi)+std_stars_vphi,'r--')
ax1.plot(r_arrays,np.abs(vc_stars_vphi)-std_stars_vphi,'r--')
ax1.fill_between(r_arrays,np.abs(vc_stars_vphi)+std_stars_vphi, np.abs(vc_stars_vphi)-std_stars_vphi,color='r',alpha=0.2)
ax1.plot(r_arrays,np.abs(vc_stars_vphi),'ro-',markersize=4,markeredgewidth=0,label=r"$\left<v_{\phi}\right>_{stars} \pm 1 \sigma$ ")
ax1.plot(r_arrays,np.abs(vc_gas_vphi)+std_gas_vphi,'g--')
ax1.plot(r_arrays,np.abs(vc_gas_vphi)-std_gas_vphi,'g--')
ax1.fill_between(r_arrays,np.abs(vc_gas_vphi)+std_gas_vphi, np.abs(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,300,texto_dos,fontsize=15)
 
ax1.plot(r_arrays,np.abs(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}$")
 
ax2.set_xlabel("r [kpc]",fontsize=15)
ax1.set_ylabel(r"v$_{c}$ [km / s]",fontsize=15)
legend = ax1.legend(loc='upper right', ncol=1, shadow=False, fontsize=15)
frame = legend.get_frame()
fig.tight_layout(h_pad=-2.2)
ax.tick_params(axis='both', which='major', labelsize=15, size=5,width=1.2)
ax.tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
ax1.tick_params(axis='both', which='major', labelsize=15, size=5,width=1.2)
ax1.tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
ax2.set_ylabel("M [M$_{\odot}$]",fontsize=15)
ax2.plot(r_arrays,stars_mass_r,label="stars")
ax2.plot(r_arrays,gas_mass_r,label="gas")
ax2.set_yscale('log')
legend = ax2.legend(loc='upper right', ncol=1, shadow=False, fontsize=15)
frame = legend.get_frame()
ax2.tick_params(axis='both', which='major', labelsize=15, size=5,width=1.2)
ax2.tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
ax2.tick_params(axis='both', which='major', labelsize=15, size=5,width=1.2)
ax2.tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
#plt.savefig("/home/arturo/Documents/git/LAMtoLUPM_latex/v_c.png",dpi=300)
```
 
%% Output
 
 
 
%% Cell type:code id: tags:
 
``` python
fig, ax1 = plt.subplots(figsize=[8,7])
 
 
ax1.set_ylim(0,400)
ax1.plot(r_arrays,np.abs(vc_stars_vphi)+std_stars_vphi,'r--')
ax1.plot(r_arrays,np.abs(vc_stars_vphi)-std_stars_vphi,'r--')
ax1.fill_between(r_arrays,np.abs(vc_stars_vphi)+std_stars_vphi, np.abs(vc_stars_vphi)-std_stars_vphi,color='r',alpha=0.2)
ax1.plot(r_arrays,np.abs(vc_stars_vphi),'ro-',markersize=4,markeredgewidth=0,label=r"$\left<v_{\phi}\right>_{stars} \pm 1 \sigma$ ")
ax1.plot(r_arrays,np.abs(vc_gas_vphi)+std_gas_vphi,'g--')
ax1.plot(r_arrays,np.abs(vc_gas_vphi)-std_gas_vphi,'g--')
ax1.fill_between(r_arrays,np.abs(vc_gas_vphi)+std_gas_vphi, np.abs(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,300,texto_dos,fontsize=15)
 
ax1.plot(r_arrays,np.abs(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=1, shadow=False, fontsize=15)
frame = legend.get_frame()
ax1.tick_params(axis='both', which='major', labelsize=15, size=5,width=1.2)
ax1.tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
```
 
%% Output
 
 
 
%% Cell type:code id: tags:
 
``` python
myhalo.gs.vtheta
```
 
%% Output
 
array([ -7.62777163, -17.315813 , -3.02011556, ..., 215.65508344,
215.27539702, 69.70578333])
 
%% Cell type:code id: tags:
 
``` python
myhalo.gs.vphi
```
 
%% Output
 
array([ -9.60223707, -32.76665095, 0.24738798, ..., -61.61756659,
-49.45941557, 110.85658131])
 
%% Cell type:code id: tags:
 
``` python
myhalo.dm.phi = Phy[:len(myhalo.dm.mass)]
myhalo.gs.phi = Phy[len(myhalo.dm.mass):len(myhalo.gs.mass)]
myhalo.st.phi = Phy[len(myhalo.dm.mass)+len(myhalo.gs.mass):]
```
 
%% Cell type:markdown id: tags:
 
# $\mathcal{E}$ for DM particles inside R$_{200}$
 
%% Cell type:code id: tags:
 
``` python
print 2**5
```
 
%% Output
 
32
 
%% Cell type:code id: tags:
 
``` python
inside = np.where(myhalo.dm.r < myhalo.r200)
epsilon = -myhalo.dm.phi[inside] - 0.5*(myhalo.dm.v[inside]**2)
e_max = epsilon[np.where(myhalo.dm.r[inside]==myhalo.dm.r[inside].min())]
print e_max
binnum = 2**8 # 128
bins = np.linspace(0,e_max,binnum)
hist_ep, bins = np.histogram(epsilon,bins=bins,normed=1)
bins_c = (bins[:-1]+bins[1:])/2.
bins_width = (bins[1:]-bins[:-1])
```
 
%% Output
 
[542120.25]
 
%% Cell type:code id: tags:
 
``` python
fig,ax = plt.subplots()
ax.set_xlim([0,1])
#ax.set_xlabel('r[kpc]',fontsize=18)
ax.set_xlabel(r'$\mathcal{E}/\mathcal{E}_{max}$',fontsize=20)
ax.bar(bins_c/e_max,hist_ep,width=bins_width/e_max,lw=0)
fig.tight_layout()
ax.tick_params(axis='both', which='major', labelsize=15, size=5,width=1.2)
ax.tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
plt.savefig("/home/arturo/Documents/git/LAMtoLUPM_latex/f_E.pdf")
```
 
%% Output
 
 
 
%% Cell type:code id: tags:
 
``` python
```
 
%% Cell type:markdown id: tags:
 
# $\frac{\vec{L}}{m}$ for particles inside R$_{200}$
 
%% Cell type:code id: tags:
 
``` python
L = np.sqrt(myhalo.dm.vphi**2 + myhalo.dm.vtheta**2)*myhalo.dm.r
L = L[inside]
#bins = np.logspace(np.log10(myhalo.gs.hsml.min()),np.log10(myhalo.r200),512)
bins = np.linspace(0,L.max(),binnum)
hist_L, Lbins = np.histogram(L,bins=bins,normed=1)
Lbins_c = (Lbins[:-1]+Lbins[1:])/2.
Lbins_width = (Lbins[1:]-Lbins[:-1])
```
 
%% Cell type:code id: tags:
 
``` python
fig,ax = plt.subplots()
ax.set_xlim([0,1])
#ax.set_xlabel('r[kpc]',fontsize=18)
ax.set_xlabel(r'|$\vec{\mathcal{L}}$|/|$\vec{\mathcal{L}}_{max}$|',fontsize=20)
#ax.axvline(x=1,color="k",linestyle='--')
ax.bar(Lbins_c/L.max(),hist_L,width=Lbins_width/L.max(),lw=0)
fig.tight_layout()
ax.tick_params(axis='both', which='major', labelsize=15, size=5,width=1.2)
ax.tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
plt.savefig("/home/arturo/Documents/git/LAMtoLUPM_latex/f_L.png",dpi=300)
```
 
%% Output
 
 
 
%% Cell type:markdown id: tags:
 
# Correlation(?)
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots(figsize=[6,6])
ax.set_xlim([0,1])
ax.set_ylim([0,1])
ax.set_xlabel(r"$\mathcal{E}/\mathcal{E}_{max}$",fontsize=20)
ax.set_ylabel(r"$|\vec{\mathcal{L}}|/|\vec{\mathcal{L}}_{max}|$",fontsize=20)
 
out = [np.where(myhalo.dm.r[inside]>30)]
ax.scatter(epsilon[out] / e_max, L[out]/L.max(),s=0.2,alpha=0.1,lw=0,c='#800b35')
ax.scatter(100,100,s=20,lw=0,c='#800b35',label='r>30kpc')
 
indisc = [np.where(myhalo.dm.r[inside]<30)]
ax.scatter(epsilon[indisc] / e_max, L[indisc]/L.max(),s=0.2,alpha=0.1,lw=0,c='#252489')
ax.scatter(100,100,s=20,lw=0,c='#252489',label='r<30kpc')
 
legend = ax.legend(loc='upper right', ncol=1, shadow=False, fontsize=15)
frame = legend.get_frame()
fig.tight_layout()
ax.tick_params(axis='both', which='major', labelsize=15, size=5,width=1.2)
ax.tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
plt.savefig("/home/arturo/Documents/git/LAMtoLUPM_latex/corrLE.png",dpi=300)
```
 
%% Output
 
 
 
%% Cell type:code id: tags:
 
``` python
x_edges = np.linspace(0,1,binnum)
y_edges = np.linspace(0,1,binnum)
 
hist2d,xbins2d,ybins2d = np.histogram2d(epsilon/e_max,L/L.max(),bins=[x_edges,y_edges])
```
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots(figsize=[6,6])
ax.set_xlabel(r"$f(\mathcal{E})$",fontsize=20)
ax.set_ylabel(r"$f(\vec{\mathcal{L}})$",fontsize=20)
 
mass_1 = ax.imshow(hist2d, interpolation='nearest', origin='low',cmap="Blues",
extent=[xbins2d[0], xbins2d[-1], ybins2d[0], ybins2d[-1]],
)
divider = make_axes_locatable(ax)
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
cbar = plt.colorbar(mass_1,cax=cax)
 
fig.tight_layout()
ax.tick_params(axis='both', which='major', labelsize=15, size=5,width=1.2)
ax.tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
plt.savefig("/home/arturo/Documents/git/LAMtoLUPM_latex/corrfLfE.png",dpi=300)
```
 
%% Output
 
 
 
%% Cell type:code id: tags:
 
``` python
fig = plt.figure()
 
ax1 = plt.subplot2grid((3, 3), (0, 0),colspan=2)
ax3 = plt.subplot2grid((3, 3), (1, 0), colspan=2, rowspan=2)
ax4 = plt.subplot2grid((3, 3), (1, 2), rowspan=2)
ax1.set_xlim([0,1])
ax4.set_ylim([0,1])
 
ax1.set_xlabel(r'$\mathcal{E}/\mathcal{E}_{max}$',fontsize=18)
ax1.bar(bins_c/e_max,hist_ep,width=bins_width/e_max)
 
mass_1 = ax3.imshow(hist2d, interpolation='nearest', origin='low',cmap="Blues",
extent=[xbins2d[0], xbins2d[-1], ybins2d[0], ybins2d[-1]],
)
divider = make_axes_locatable(ax3)
 
ax4.set_xlabel(r'$\vec{\mathcal{L}}$',fontsize=18)
#ax.axvline(x=1,color="k",linestyle='--')
ax4.barh(Lbins_c/L.max(),hist_L,height=Lbins_width/L.max())
fig.tight_layout(w_pad=-1)
```
 
%% Output
 
 
 
%% Cell type:markdown id: tags:
 
# Moments of f(v)
 
%% Cell type:markdown id: tags:
 
### First moment
 
%% Cell type:code id: tags:
 
``` python
 
bis = np.logspace(np.log10(myhalo.gs.hsml.min()),np.log(myhalo.r200),)
```
 
%% Cell type:code id: tags:
 
``` python
```
 
%% Cell type:code id: tags:
 
``` python
```
 
%% Cell type:code id: tags:
 
``` python
```
 
%% Cell type:markdown id: tags:
 
# Energy Jacobian K($\mathcal{E}$)
 
%% Cell type:code id: tags:
 
``` python
r_sorted_dm = np.argsort(myhalo.dm.r)
r_K_DM = myhalo.dm.r[r_sorted_dm]
pot_K_DM = myhalo.dm.phi[r_sorted_dm]
pot_of_r = interp1d(r_K_DM,pot_K_DM)
```
 
%% Cell type:code id: tags:
 
``` python
#fig, ax = plt.subplots()
#ax.set_yscale('log')
def K_of_E(E,rmin,rmax):
"""
Energy Jacobian
rmin and rmax should be in kpc
and E in km^2 Msun s^-2
retunrs:
integral and error
integral = (4\pi)^2 \int_{rmin}^{rmax} r^2 \sqrt(2(|\Phi|-\mathcal{E}))
"""
def integrand(r):
#if (np.abs(pot_of_r(r))>E):
return (4.* np.pi * r)**2 * np.sqrt(2*(np.abs(pot_of_r(r)) - E))
r_array = np.linspace(rmin,rmax,1e6)
inte_array= integrand(r_array)
# print inte_array
inte_array = np.nan_to_num(inte_array)
# ax.plot(r_array,-pot_of_r(r_array))
return simps(inte_array,r_array)
 
"""
E = np.linspace(0.,e_max,50)
 
r = 10**np.linspace(np.log10(1.),np.log10(myhalo.r200),100)
r_E = np.reshape(np.repeat(r,len(E)),(len(r),len(E)))
print integrand(r_E)
#integral = simps(integrand(r_E),r_E,axis=0)
"""
K_of_E_vec = np.vectorize(K_of_E)
E = np.linspace(-1.1,e_max,100)
KE = K_of_E_vec(E,0.3,4*myhalo.r200)
```
 
%% Cell type:code id: tags:
 
``` python
```
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots()
ax.set_yscale('log')
fK_of_E_ = interp1d(E,KE)
ax.plot(E,KE)
```
 
%% Output
 
 
 
[<matplotlib.lines.Line2D at 0x7fc47803e190>]
 
%% Cell type:code id: tags:
 
``` python
fig,ax = plt.subplots()
#ax.set_xlim([0,1])
#ax.set_xlabel('r[kpc]',fontsize=18)
ax.set_xlabel(r'$\mathcal{E}/\mathcal{E}_{max}$',fontsize=18)
ax.bar(bins_c/e_max,hist_ep,width=bins_width/e_max,lw=0)
```
 
%% Output
 
 
 
<Container object of 255 artists>
 
%% Cell type:code id: tags:
 
``` python
k_corr = fK_of_E_(bins_c)
```
 
%% Cell type:code id: tags:
 
``` python
fig,ax = plt.subplots()
ax.set_xlim([0,2])
#ax.set_xlabel('r[kpc]',fontsize=18)
ax.set_xlabel(r'$\mathcal{E}/\mathcal{E}_{max}$',fontsize=18)
k_corr = fK_of_E_(bins_c)
ax.bar(bins_c/e_max,hist_ep/k_corr,width=bins_width/e_max)
```
 
%% Output
 
 
 
<Container object of 255 artists>
 
%% Cell type:code id: tags:
 
``` python
qfile = open("../Simulations/Q_param.txt")
qzj17 = qsh04 = K_halo = U_halo = mod = M_200 = Simname = np.array([])
for l in qfile:
if l[0]=="#":
continue
row = l.split(',')
qzj17 = np.append(qzj17, float(row[0]))
qsh04 = np.append(qsh04, float(row[1]))
K_halo = np.append(K_halo, float(row[2]))
U_halo = np.append(U_halo, float(row[3]))
mod = np.append(mod, float(row[4]))
M_200 = np.append(M_200, float(row[5]))
Simname = np.append(Simname, row[6][:-1])
 
 
```
 
%% Cell type:code id: tags:
 
``` python
fig, [ax,ax1] = plt.subplots(2,1,gridspec_kw = {'height_ratios':[1, 1]},figsize=[6,4.6],sharex=True)
ax.set_ylabel(r"q",fontsize=0.8*font)
ax.set_xscale('log')
ax.set_xlim([1e11,1e13])
ax.set_ylim([-0.025,.35])
ax.set_yticks(np.arange(0.05,0.4,0.05))
ax.scatter(M_200,qzj17)
ax.arrow( 7e12, 0.08, 0.0, -0.05, fc="k", ec="k",head_width=5e11, head_length=0.01 )
ax.text(2e12, 0.03,"U dominated\n q<0")
ax.arrow( 7e12, 0.1, 0.0, 0.05, fc="k", ec="k",head_width=5e11, head_length=0.01 )
ax.text(2e12, 0.1,"K dominated\n q>0")
ax1.set_ylim([-0.2,1.23])
ax1.set_xlabel(r"M$_{200}$ [M$_{\odot}$]",fontsize=font)
ax1.set_ylabel(r"$\Delta_{center}$",fontsize=1.2*font)
ax1.axhline(y=0.190,color='k',linestyle='--')
ax1.axhline(y=0.090,color='k',linestyle='--')
ax1.scatter(M_200,mod)
ax1.text(3e12,-0.01,'hsml level 18')
ax1.text(3e12,0.23,'hsml level 17')
ax.text(1.2e11,0.28,"(DMO)", fontsize=font)
fig.tight_layout(h_pad=-0.19)
ax.tick_params(axis='both', which='major', labelsize=15, size=5,width=1.2)
ax.tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
ax1.tick_params(axis='both', which='major', labelsize=15, size=5,width=1.2)
ax1.tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
plt.savefig("/home/arturo/Documents/git/LAMtoLUPM_latex/Qparam.png",dpi=300)
```
 
%% Output
 
 
 
%% Cell type:code id: tags:
 
``` python
get_mat = np.vectorize(nbe.m_matrix_for_r)
```
 
%% Cell type:code id: tags:
 
``` python
```
 
%% Cell type:code id: tags:
 
``` python
#r = np.logspace(np.log10(myhalo.gs.hsml.min()),np.log10(2*myhalo.r200),15)
r = np.linspace(myhalo.gs.hsml.min(),2*myhalo.r200,300)
M_dm = np.array([nbe.m_matrix_for_r(myhalo,'halo',i) for i in r])
M_st = np.array([nbe.m_matrix_for_r(myhalo,'stars',i) for i in r])
a_dm, b_dm, c_dm = np.sqrt(M_dm[:,0,0]), np.sqrt(M_dm[:,1,1]), np.sqrt(M_dm[:,2,2])
a_st, b_st, c_st = np.sqrt(M_st[:,0,0]), np.sqrt(M_st[:,1,1]), np.sqrt(M_st[:,2,2])
S_dm = c_dm/a_dm
T_dm = ((a_dm**2) - (b_dm**2))/((a_dm**2) -(c_dm**2))
S_st = c_st/a_st
T_st = ((a_st**2) - (b_st**2))/((a_st**2) -(c_st**2))
```
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots(figsize=[10,4])
#ax.set_xlim(4,250)
ax.set_ylim(0,1.18)
ax.set_xlim(0,410)
ax.set_xlabel('r [kpc]',fontsize=font)
ax.plot([4,5],[1e3,2e3],color='gray',linestyle='--',lw=2,label="S")
ax.plot([4,5],[1e3,2e3],color='gray',linestyle='-',lw=2,label="T")
ax.plot(r,T_dm,'-',color='#155c9e',lw=2,label="Dark matter")
ax.plot(r,S_dm,'--',color='#155c9e',lw=2)
#ax.plot(r,T_st,'-',color='#a91e4f',lw=2,label='stars')
#ax.plot(r,np.sqrt(S_st),'--',color='#a91e4f',lw=2)
ax.axvline(x=myhalo.r200,color='k',linestyle='-.')
ax.text(myhalo.r200+2,0.2,r"R$_{200}$",fontsize=font)
ax.axvline(x=myhalo.r97,color='k',linestyle='-.')
ax.text(myhalo.r97+2,0.2,r"R$_{97}$",fontsize=font)
ax.axvline(x=8,color='k',linestyle='-.')
ax.text(10,0.2,r"R$_{\odot}$",fontsize=font)
 
 
legend = ax.legend(loc='upper right', ncol=2, shadow=False, fontsize=font)
frame = legend.get_frame()
ax.tick_params(axis='both', which='major', labelsize=15, size=5,width=1.2)
ax.tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
fig.tight_layout()
```
 
%% Output
 
 
 
%% Cell type:code id: tags:
 
``` python
myhalo.p.
```
 
%% Output
 
File "<ipython-input-64-6f8a25829487>", line 1
myhalo.p.
^
SyntaxError: invalid syntax
 
%% Cell type:code id: tags:
 
``` python
m_bin = np.linspace(myhalo.gs.hsml.min(),2*myhalo.r200,60)
 
mass_r, bins = np.histogram(myhalo.dm.r,bins=m_bin,weights=myhalo.dm.r)
b_center = (bins[:-1]+bins[1:])/2
```
 
%% 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
Pcrit = myhalo.dm._p.rho_crit
Mdm = myhalo.dm.mass.min()
myradiuses = myhalo.dm.r[np.argsort(myhalo.dm.r)]
 
tabN = np.cumsum(np.ones(len(myradiuses)))[1:]
myradiuses = myradiuses[1:]
Rp03 = np.sqrt(200/64.) * np.sqrt(4 * np.pi * Pcrit * tabN / 3. / Mdm ) * (myradiuses**1.5)/ np.log(tabN)
val =0.6
R_P03 = myradiuses[ np.where(Rp03 > val) ][0]
print R_P03
hsml= myhalo.gs.hsml.min()
 
# R array logarithmic Bining
r_p = np.logspace(np.log10(0.2*hsml),np.log10(hsml),15)
# histogram of dm particles per logarithmic bin
n_hydro,r = np.histogram(myhalo.dm.r,bins=r_p)
# edges of bins
r1,r2 =r[:-1],r[1:]
# shell's volume
vol = 4.* np.pi * ((r2**3)-(r1**3)) / 3.
# density per shell
profilehydro_in = n_hydro*myhalo.dm.mass.min()/vol
# center of bins
r_in = (r_p[:-1]+r_p[1:])/2.
 
 
 
# R array logarithmic Bining
r_p = np.logspace(np.log10(2.*hsml),np.log10(2.5*myhalo.r200),100)
# histogram of dm particles per logarithmic bin
n_hydro,r = np.histogram(myhalo.dm.r,bins=r_p)
# edges of bins
r1,r2 =r[:-1],r[1:]
# shell's volume
vol = 4.* np.pi * ((r2**3)-(r1**3)) / 3.
# density per shell
profilehydro = n_hydro*myhalo.dm.mass.min()/vol
# center of bins
r = (r_p[:-1]+r_p[1:])/2.
bin_size= (r_p[:-1]-r_p[1:])/2.
rr = r
# extra estatistics from Cfalcon density
mean_hydro = std_hydro = stdlog_hydro = n_hydro=np.array([])
for i in range(len(r_p)-1):
shell_hydro = np.where((myhalo.dm.r > r_p[i])&(myhalo.dm.r < r_p[i+1])&(myhalo.dm.r > hsml))
n_hydro = np.append(n_hydro,len(shell_hydro[0]))
mean_hydro = np.append(mean_hydro,np.mean(myhalo.dm.rho[shell_hydro]))
std_hydro = np.append(std_hydro,np.std(myhalo.dm.rho[shell_hydro]))
stdlog_hydro = np.append(stdlog_hydro,np.std(np.log10(myhalo.dm.rho[shell_hydro])))
 
m_obs_hydro = n_hydro*myhalo.dm.mass.min()
n_hydro = np.array([len(myhalo.dm.mass[myhalo.dm.r<i]) for i in r])
```
 
%% Cell type:code id: tags:
 
``` python
fig , ax = plt.subplots()
r = np.linspace(myhalo.gs.hsml.min(),2*myhalo.r200,300)
 
ax.set_ylim(0,1.18)
ax.set_xlim(0,410)
ax.set_xlabel('r [kpc]',fontsize=font)
ax.plot([4,5],[1e3,2e3],color='gray',linestyle='--',lw=2,label="S")
ax.plot([4,5],[1e3,2e3],color='gray',linestyle='-',lw=2,label="T")
ax.plot(r,T_dm,'-',color='#155c9e',lw=2)
ax.plot(r,S_dm,'--',color='#155c9e',lw=2)
#ax.plot(r,T_st,'-',color='#a91e4f',lw=2,label='stars')
#ax.plot(r,np.sqrt(S_st),'--',color='#a91e4f',lw=2)
ax.axvline(x=myhalo.r200,color='k',linestyle='-.')
ax.text(myhalo.r200+2,0.2,r"R$_{200}$",fontsize=font)
ax.axvline(x=myhalo.r97,color='k',linestyle='-.')
ax.text(myhalo.r97+2,0.2,r"R$_{97}$",fontsize=font)
ax.axvline(x=8,color='k',linestyle='-.')
ax.text(10,0.2,r"R$_{\odot}$",fontsize=font)
 
 
legend = ax.legend(loc='upper right', ncol=2, shadow=False, fontsize=font)
frame = legend.get_frame()
ax.tick_params(axis='both', which='major', labelsize=15, size=5,width=1.2)
ax.tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
```
 
%% Cell type:markdown id: tags:
 
# $\beta$ parameter
 
$$\beta = 1 - \frac{ \sigma_{\phi}^2 + \sigma_{\theta}^2}{2 \sigma_{r}^2}$$
 
%% Cell type:code id: tags:
 
``` python
point_num = 150
r_beta = np.logspace(-1,np.log10(4*myhalo.r200),point_num)
vphi = np.concatenate((myhalo.dm.vphi,myhalo.st.vphi,myhalo.gs.vphi))
vtheta = np.concatenate((myhalo.dm.vtheta,myhalo.st.vtheta,myhalo.gs.vtheta))
vr = np.concatenate((myhalo.dm.vr,myhalo.st.vr,myhalo.gs.vr))
```
 
%% Cell type:code id: tags:
 
``` python
def beta_param(i):
condition = np.where((r2>r_beta[i]**2)&(r2<=r_beta[i+1]**2))
v_r = vr[condition]
v_phi = vphi[condition]
v_theta = vtheta[condition]
#print (np.std(v_phi))**2 ,(np.std(v_theta))**2 , (np.std(v_r))**2
return 1 - ((np.std(v_phi))**2 +(np.std(v_theta))**2) / 2. / (np.std(v_r))**2
 
get_beta = np.vectorize(beta_param)
```
 
%% Cell type:code id: tags:
 
``` python
beta_r = get_beta(range(point_num-1))
```
 
%% Cell type:code id: tags:
 
``` python
def beta_param_2(i):
condition = np.where((r2>r_beta[i]**2)&(r2<=r_beta[i+1]**2))
v_r = vr[condition]
v_phi = vphi[condition]
v_theta = vtheta[condition]
#print (np.std(v_phi))**2 ,(np.std(v_theta))**2 , (np.std(v_r))**2
return 1 - (np.mean(v_phi**2) + np.mean(v_theta**2)) / 2. / (np.mean(v_r**2 ))
 
get_beta_2 = np.vectorize(beta_param_2)
```
 
%% Cell type:code id: tags:
 
``` python
beta_r_2 = get_beta_2(range(point_num-1))
```
 
%% Cell type:code id: tags:
 
``` python
def beta_param_3(i):
condition = np.where((r2>r_beta[i]**2)&(r2<=r_beta[i+1]**2))
v_r = vr[condition]
v_phi = vphi[condition]
v_theta = vtheta[condition]
#print (np.std(v_phi))**2 ,(np.std(v_theta))**2 , (np.std(v_r))**2
return 1 - ((np.mean(v_phi)**2) + (np.mean(v_theta)**2)) / 2. / (np.mean(v_r)**2 )
get_beta_3 = np.vectorize(beta_param_3)
```
 
%% Cell type:code id: tags:
 
``` python
beta_r_3 = get_beta_3(range(point_num-1))
```
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots(figsize=[9,5])
r_2 = 6.95
ax.set_xscale('log')
#ax.set_title(r"$\beta(r)$ parameter", fontsize=1.5*font)
ax.set_xlabel(r"r [kpc]",fontsize=font)
ax.set_ylabel(r"$\beta(r)$",fontsize=1.5*font)
ax.set_ylim([-1.5,1])
ax.plot((r_beta[1:]+r_beta[:-1])/2.,beta_r,lw=1.5)
#ax.plot((r_beta[1:]+r_beta[:-1])/2/r_2,beta_r_2,lw=1.5)
#ax.plot((r_beta[1:]+r_beta[:-1])/2/r_2,beta_r_3,lw=1.5)
 
 
ax.axvline(x=myhalo.r200, color='k',lw=2,linestyle='--')
ax.axvline(x=myhalo.r97, color='gray',lw=2,linestyle='--')
 
ax.axvline(x=8/r_2, color='y',linestyle='--',lw=2,label=r'r$_{\odot}$')
fig.tight_layout()
ax.tick_params(axis='both', which='major', labelsize=15, size=5,width=1.2)
ax.tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
#plt.savefig("/home/arturo/Documents/git/LAMtoLUPM_latex/beta_r.png",dpi=300)
```
 
%% Cell type:code id: tags:
 
``` python
def delta_beta(i):
condition = np.where((r2>r_beta[i]**2)&(r2<=r_beta[i+1]**2))
v_r = vr[condition]
v_phi = vphi[condition]
v_theta = vtheta[condition]
first = (np.std(v_phi**2))**2 / 4./ np.mean(v_phi**2)**2
secon = (np.std(v_theta**2))**2 / 4./ np.mean(v_theta**2)**2
third = ((np.mean(v_phi**2) + np.mean(v_theta**2))**2) * np.std(v_r**2)**2 / 2. / np.mean(v_r**2)**2
return first + secon + third
get_deltabeta = np.vectorize(delta_beta)
```
 
%% Cell type:code id: tags:
 
``` python
deltabeta_r = get_deltabeta(range(point_num-1))
```
 
%% Cell type:code id: tags:
 
``` python
```
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots()
r_2 = 6.95
ax.set_xscale('log')
ax.set_title(r"$\beta(r)$ parameter", fontsize=22)
ax.set_xlabel(r"r / r$_2$ ",fontsize=15)
ax.set_ylabel(r"$\beta(r)$",fontsize=18)
ax.set_ylim([-4,2])
ax.plot((r_beta[1:]+r_beta[:-1])/2/r_2,beta_r_2,'b',lw=1.3)
ax.plot((r_beta[1:]+r_beta[:-1])/2/r_2,beta_r_2+deltabeta_r,'b--',lw=1.1)
ax.plot((r_beta[1:]+r_beta[:-1])/2/r_2,beta_r_2-deltabeta_r,'b--',lw=1.1)
ax.fill_between((r_beta[1:]+r_beta[:-1])/2/r_2,
beta_r_2+deltabeta_r, beta_r_2-deltabeta_r,color='b',alpha=0.2,label=r"$1 \sigma$")
ax.axvline(x=myhalo.r200/r_2, color='k',lw=2,linestyle='--',label=r'r$_{200}$')
#ax.axvline(x=myhalo.r97/r_2, color='gray',lw=2,linestyle='--',label=r'r$_{97}$')
#ax.axvline(x=np.sqrt(r_nei2)/r_2, color='g',lw=2,linestyle='--',label=r'first neighbourg')
#ax.axvline(x=8/r_2, color='y',linestyle='--',lw=2,label=r'r$_{\odot}$')
 
legend = ax.legend(loc='lower center', ncol=4, shadow=False, fontsize=14)
frame = legend.get_frame()
plt.savefig("/home/arturo/Documents/LAM/LAM2LUPM/report07-17/HaloB/Beta_param.pdf",dpi=300)
```
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots()
r_2 = 6.95
ax.set_xscale('log')
ax.set_title(r"$\beta(r)$ parameter", fontsize=22)
ax.set_xlabel(r"r / r$_2$ ",fontsize=15)
ax.set_ylabel(r"$\beta(r)$",fontsize=18)
ax.set_ylim([-4,2])
ax.plot((r_beta[1:]+r_beta[:-1])/2/r_2,beta_r,'k',lw=1.2,
label=r"$\beta_1$")
ax.plot((r_beta[1:]+r_beta[:-1])/2/r_2,beta_r_2,'b',lw=1.2,label=r"$\beta_2$")
ax.plot((r_beta[1:]+r_beta[:-1])/2/r_2,beta_r_2+deltabeta_r,'b--',lw=1.1)
ax.plot((r_beta[1:]+r_beta[:-1])/2/r_2,beta_r_2-deltabeta_r,'b--',lw=1.1)
ax.fill_between((r_beta[1:]+r_beta[:-1])/2/r_2,
beta_r_2+deltabeta_r, beta_r_2-deltabeta_r,color='b',alpha=0.2,label=r"$\beta_2 \pm 1 \sigma$")
 
beta1 = r"$\beta_1 = 1 - \frac{ \sigma_{\phi}^2 + \sigma_{\theta}^2}{2 \sigma_{r}^2}$"
ax.text(0.08,-2.5,beta1,fontsize=15)
beta2 = r"$\beta_2 = 1 - \frac{ <v_{\phi}^2> + <v_{\theta}^2>}{2 <v_{r}^2>}$"
ax.text(1.,-2.5,beta2,fontsize=15)
 
#ax.axvline(x=myhalo.r200/r_2, color='k',lw=2,linestyle='--',label=r'r$_{200}$')
#ax.axvline(x=myhalo.r97/r_2, color='gray',lw=2,linestyle='--',label=r'r$_{97}$')
#ax.axvline(x=np.sqrt(r_nei2)/r_2, color='g',lw=2,linestyle='--',label=r'first neighbourg')
#ax.axvline(x=8/r_2, color='y',linestyle='--',lw=2,label=r'r$_{\odot}$')
 
legend = ax.legend(loc='lower center', ncol=4, shadow=False, fontsize=14)
frame = legend.get_frame()
plt.savefig("/home/arturo/Documents/LAM/LAM2LUPM/report07-17/HaloB/Beta_paramB1vsB2.pdf",dpi=300)
```
 
%% Cell type:markdown id: tags:
 
# F(v) moments
 
%% Cell type:code id: tags:
 
``` python
def moments(rmin,rmax):
selection= np.where((myhalo.dm.r>rmin)&(myhalo.dm.r<=rmax))
fdv, vs = np.histogram(myhalo.dm.v[selection],bins=50,normed=1)
vcenter = (vs[:-1]+vs[1:])/2.
m1 = simps(vcenter*fdv,x=vcenter)
m2 = simps((vcenter**2)*fdv,x=vcenter)
m3 = simps((vcenter**3)*fdv,x=vcenter)
 
return m1,m2,m3
 
fdv_moments = np.vectorize(moments)
```
 
%% Cell type:code id: tags:
 
``` python
r_v = np.logspace(np.log10(4*myhalo.gs.hsml.min()),np.log10(3*myhalo.r200),150)
m1,m2,m3 = fdv_moments(r_v[:-1],r_v[1:])
```
 
%% Cell type:code id: tags:
 
``` python
 
fig, ax = plt.subplots(3,1)
ax[2].set_xlabel('r [kpc]',fontsize=18)
ax[0].set_ylabel('[km/s]',fontsize=18)
ax[1].set_ylabel(r'[km/s]$^2$',fontsize=18)
ax[2].set_ylabel(r'[km/s]$^3$',fontsize=18)
ax[0].set_xscale('log')
ax[1].set_xscale('log')
ax[2].set_xscale('log')
ax[0].set_ylim([52,400])
ax[1].set_ylim([1e4,1.5e5])
ax[2].set_ylim([1e4,6.5e7])
ax[0].set_xlim([0.5,3*myhalo.r200])
ax[1].set_xlim([0.5,3*myhalo.r200])
ax[2].set_xlim([0.5,3*myhalo.r200])
ax[0].axvline(x=myhalo.r200,color='k',linestyle='--')
ax[1].axvline(x=myhalo.r200,color='k',linestyle='--')
ax[2].axvline(x=myhalo.r200,color='k',linestyle='--')
r_v_c = (r_v[:-1]+r_v[1:])/2.
ax[0].plot(r_v_c,m1)
ax[1].plot(r_v_c,m2)
ax[2].plot(r_v_c,m3)
 
fig.tight_layout(h_pad=-2.7)
ax[0].text(22,250,r'$\int_0 ^{v_{max}} v f(v) dv$',fontsize=17)
ax[1].text(22,70000,r'$\int_0 ^{v_{max}} v^2 f(v) dv$',fontsize=17)
ax[2].text(22,3e7,r'$\int_0 ^{v_{max}} v^3 f(v) dv$',fontsize=17)
ax[2].text(170,1e7,r'R$_{200}$',fontsize=17)
 
 
ax[0].tick_params(axis='y', which='major', labelsize=15, size=5,width=1.2)
ax[0].tick_params(axis='x', which='major', labelsize=0, size=5,width=1.2)
ax[0].tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
ax[1].tick_params(axis='y', which='major', labelsize=15, size=5,width=1.2)
ax[1].tick_params(axis='x', which='major', labelsize=0, size=5,width=1.2)
ax[1].tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
ax[2].tick_params(axis='both', which='major', labelsize=15, size=5,width=1.2)
ax[2].tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
plt.savefig("/home/arturo/Documents/git/LAMtoLUPM_latex/moments.png",dpi=300)
```
 
%% Cell type:code id: tags:
 
``` python
```
 
%% Cell type:code id: tags:
 
``` python
```
 
%% Cell type:code id: tags:
 
``` python
```
 
%% Cell type:markdown id: tags:
 
* histograms of f(v) at different radii with chi2 (Eddington vs. Maxwellian) - already done, except for the Maxwellian case.
* vesc(r)+-dvesc at different radii + comparison with the value calculated from the potential, $\sqrt(-2\,\phi(r))$
* beta(r)+-dbeta at different radii
* compute the "equilibrium" test variable (ask Manu)
 
%% Cell type:code id: tags:
 
``` python
```
 
%% Cell type:code id: tags:
 
``` python
```
 
%% Cell type:code id: tags:
 
``` python
limmin, limmax = 7.5,8.5
shell = (myhalo.dm.r>limmin)&(myhalo.dm.r<limmax)
disc = (np.abs(myhalo.dm.pos3d[:,2])<1.5)
nodisc = (np.abs(myhalo.dm.pos3d[:,2])>1.5)
v_disc_in_the_shell = myhalo.dm.v[shell&disc]
v_nodisc_in_the_shell = myhalo.dm.v[shell&nodisc]
bins = np.linspace(0,myhalo.dm.v[shell].max(),30)
```
 
%% Cell type:code id: tags:
 
``` python
plt.figure()
plt.hist(v_disc_in_the_shell, bins, stacked=True)
plt.hist(v_nodisc_in_the_shell, bins, stacked=True)
print len(v_disc_in_the_shell)
print len(v_nodisc_in_the_shell)
```
 
%% Cell type:code id: tags:
 
``` python
plt.figure()
plt.hist([v_disc_in_the_shell,v_nodisc_in_the_shell], bins, stacked=True, normed=True)
plt.show()
```
 
%% Cell type:code id: tags:
 
``` python
```
 
%% Cell type:code id: tags:
 
``` python
nbe.real_center(myhalo.dm.pos3d,myhalo.dm.mass)
```
 
%% Cell type:code id: tags:
 
``` python
```
 
%% Cell type:code id: tags:
 
``` python
```
 
%% Cell type:code id: tags:
 
``` python
```