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

Automatic commit samedi 21 juillet 2018, 16:30:01 (UTC+0200)

parent da823473
......@@ -98,6 +98,25 @@
"myhydro.r_virial(600,n=3)\n"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.152\n"
]
}
],
"source": [
"print \"{0:.3f}\".format(myhydro.gs.hsml.min())"
]
},
{
"cell_type": "code",
"execution_count": 4,
......
%% 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.interpolate import interp1d
from matplotlib.patches import Circle
from scipy.special import gamma
import numpy as np
import emcee
from mpl_toolkits.axes_grid1 import make_axes_locatable
from numpy import exp, sqrt
from scipy.integrate import quad, dblquad, nquad
import matplotlib.patches as patches
from itertools import product
from scipy.integrate import quad
import scipy.optimize as optimize
import matplotlib.pyplot as plt
import matplotlib as mpl
from sklearn.neighbors import KDTree
import sys
import lmfit
from py_unsio import *
import pymc
import os
from scipy.integrate import simps
from pymodelfit import FunctionModel1DAuto
import wkbl
from wkbl.astro.halo_info import *
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm
import wkbl.astro.nbody_essentials as nbe
import cfalcon
CF =cfalcon.CFalcon()
import iminuit
from iminuit import Minuit, describe, Struct
import probfit
from matplotlib.colors import LogNorm
from matplotlib.ticker import FormatStrFormatter
import warnings
warnings.filterwarnings('ignore')
```
%% Cell type:code id: tags:
``` python
#halo = HALOBHydro(where="home")
simname="HALO B"
pathsim = "/data/POL/HALOB/hydro/output_00417"
myhydro = wkbl.Galaxy_Hound(pathsim)
print myhydro.dm.pos3d[:,0].max()
zoom_reg = np.where(myhydro.dm.mass==myhydro.dm.mass.min())
nucenter = nbe.real_center(myhydro.dm.pos3d[zoom_reg], myhydro.dm.mass[zoom_reg])
print nucenter
myhydro.center_shift(nucenter)
myhydro.r_virial(600,n=3)
```
%%%% Output: stream
loading Dark matter..
loading Stars..
loading Gas..
19879.156
[9667.65173408 9866.02624312 9801.55717996]
| r_200 = 177.54
| Diagonal matrix computed
| | 20, 0, 0|
| D =| 0, 14, 0|
| | 0, 0, 2|
%% Cell type:code id: tags:
``` python
print "{0:.3f}".format(myhydro.gs.hsml.min())
```
%%%% Output: stream
0.152
%% Cell type:code id: tags:
``` python
simname_nospace = list(simname)
for i in range(len(simname)):
if simname[i]==" ":
simname_nospace[i]="_"
simname_nospace = "".join(simname_nospace)
```
%% Cell type:code id: tags:
``` python
pos_dm = np.array(myhydro.dm.pos3d.reshape(len(myhydro.dm.pos3d)*3),dtype=np.float32)
pos_gs = np.array(myhydro.gs.pos3d.reshape(len(myhydro.gs.pos3d)*3),dtype=np.float32)
pos_st = np.array(myhydro.st.pos3d.reshape(len(myhydro.st.pos3d)*3),dtype=np.float32)
pos = np.concatenate((pos_dm, pos_st, pos_gs))
phi_cord = np.concatenate((myhydro.dm.phi,myhydro.st.phi, myhydro.gs.phi))
mass = np.concatenate((myhydro.dm.mass,myhydro.st.mass,myhydro.gs.mass))
v = np.concatenate((myhydro.dm.v,myhydro.st.v,myhydro.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: stream
37736919 37736919
%% Cell type:markdown id: tags:
# Ellipticity T and S
%% Cell type:code id: tags:
``` python
get_mat = np.vectorize(nbe.m_matrix_for_r)
```
%% Cell type:code id: tags:
``` python
r = np.linspace(myhydro.gs.hsml.min(),2*myhydro.r200,300)
M_dm = np.array([nbe.m_matrix_for_r(myhydro,'halo',i) for i in r])
M_st = np.array([nbe.m_matrix_for_r(myhydro,'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
outputing = open("../../datafiles/"+simname_nospace+"_S_and_T.txt","w")
outputing.write("# "+simname+" ellipticity parameters\n")
outputing.write("# r200 = {0:.2f} kpc\n".format(myhydro.r200))
outputing.write("# format:\n")
outputing.write("# r (kpc), S , T\n")
for i in range(len(T_dm)):
outputing.write("{0:.2f} {1:.6f} {2:.6f} \n".format(r[i],S_dm[i],T_dm[i]))
outputing.close()
```
%% Cell type:code id: tags:
``` python
fig, ax = plt.subplots(figsize=[10,4])
#ax.set_xlim(4,250)
font=15
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=myhydro.r200,color='k',linestyle='-.')
ax.text(myhydro.r200+2,0.2,r"R$_{200}$",fontsize=font)
ax.axvline(x=myhydro.r97,color='k',linestyle='-.')
ax.text(myhydro.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)
ax.text(50,0.9,r"$\rm "+simname+"$",fontsize=1.5*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: display_data
%%%% Output: display_data
%% Cell type:markdown id: tags:
# Beta of r
%% Cell type:code id: tags:
``` python
point_num = 150
r_beta = np.logspace(-1,np.log10(4*myhydro.r200),point_num)
vphi = np.concatenate((myhydro.dm.vphi,myhydro.st.vphi,myhydro.gs.vphi))
vtheta = np.concatenate((myhydro.dm.vtheta,myhydro.st.vtheta,myhydro.gs.vtheta))
vr = np.concatenate((myhydro.dm.vr,myhydro.st.vr,myhydro.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
outputing = open("../../datafiles/"+simname_nospace+"_Beta.txt","w")
outputing.write("# "+simname+" Beta parameters\n")
outputing.write("# r200 = {0:.2f} kpc\n".format(myhydro.r200))
outputing.write("# format:\n")
outputing.write("# r (kpc), Beta(r) \n")
for i in range(len(beta_r)):
outputing.write("{0:.2f} {1:.6f} \n".format(((r_beta[1:]+r_beta[:-1])/2.)[i],beta_r[i]))
outputing.close()
```
%% 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=myhydro.r200, color='k',lw=2,linestyle='--')
ax.axvline(x=myhydro.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)
```
%%%% Output: display_data
%%%% Output: display_data
%% Cell type:markdown id: tags:
# F(v) moments
%% Cell type:code id: tags:
``` python
def moments(rmin,rmax):
selection= np.where((myhydro.dm.r>rmin)&(myhydro.dm.r<=rmax))
fdv, vs = np.histogram(myhydro.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*myhydro.gs.hsml.min()),np.log10(3*myhydro.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[0].yaxis.get_major_formatter().set_powerlimits((0, 2))
ax[1].yaxis.get_major_formatter().set_powerlimits((0, 2))
ax[2].yaxis.get_major_formatter().set_powerlimits((0, 2))
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,1.1*m1.max()])
ax[1].set_ylim([1e4,1.1*m2.max()])
ax[2].set_ylim([1e4,1.1*m3.max()])
ax[0].set_xlim([0.5,3*myhydro.r200])
ax[1].set_xlim([0.5,3*myhydro.r200])
ax[2].set_xlim([0.5,3*myhydro.r200])
ax[0].axvline(x=myhydro.r200,color='k',linestyle='--')
ax[1].axvline(x=myhydro.r200,color='k',linestyle='--')
ax[2].axvline(x=myhydro.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)
```
%%%% Output: display_data
%%%% Output: display_data
%% Cell type:code id: tags:
``` python
outputing = open("../../datafiles/"+simname_nospace+"_Moments.txt","w")
outputing.write("# "+simname+" Statistical moments of dm velocity\n")
outputing.write("# r200 = {0:.2f} kpc\n".format(myhydro.r200))
outputing.write("# format:\n")
outputing.write("# r (kpc), 1stmoment, 2ndmoment, 3rdmoment\n")
for i in range(len(r_v_c)):
outputing.write("{0:.2f} {1:.6e} {2:.6e} {3:.6e}\n".format(r_v_c[i],m1[i],m2[i],m3[i]))
outputing.close()
```
%% Cell type:code id: tags:
``` python
```
......
......@@ -101,6 +101,29 @@
"myhydro.r_virial(600,n=1)\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.281102\n",
"1536979.0\n",
"183658.06\n"
]
}
],
"source": [
"print myhydro.gs.hsml.min()\n",
"print myhydro.dm.mass.min()\n",
"print myhydro.st.mass.min()"
]
},
{
"cell_type": "code",
"execution_count": 4,
......@@ -139,7 +162,7 @@
},
{
"cell_type": "code",
"execution_count": 14,
"execution_count": 20,
"metadata": {
"collapsed": false
},
......@@ -161,12 +184,26 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 22,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": []
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"205.6640625\n",
"227993.1\n",
"/data/OWN/DMO/Adicora/output_00041\n"
]
}
],
"source": [
"print myDMO.r200\n",
"print myDMO.dm.mass.min()\n",
"print dmo.path"
]
},
{
"cell_type": "code",
......
%% 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.interpolate import interp1d
from matplotlib.patches import Circle
from scipy.special import gamma
import numpy as np
import emcee
from mpl_toolkits.axes_grid1 import make_axes_locatable
from numpy import exp, sqrt
from scipy.integrate import quad, dblquad, nquad
import matplotlib.patches as patches
from itertools import product
from scipy.integrate import quad
import scipy.optimize as optimize
import matplotlib.pyplot as plt
import matplotlib as mpl
from sklearn.neighbors import KDTree
import sys
import lmfit
from py_unsio import *
import pymc
import os
from pymodelfit import FunctionModel1DAuto
import wkbl
from wkbl.astro.halo_info import *
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm
import wkbl.astro.nbody_essentials as nbe
import cfalcon
CF =cfalcon.CFalcon()
import iminuit
from iminuit import Minuit, describe, Struct
import probfit
from matplotlib.colors import LogNorm
from matplotlib.ticker import FormatStrFormatter
import warnings
warnings.filterwarnings('ignore')
```
 
%% Cell type:markdown id: tags:
 
# Hydro
 
%% Cell type:code id: tags:
 
``` python
hydro = wkbl.astro.halo_info.AdicoraHydro()
simname=hydro.name
myhydro = wkbl.Galaxy_Hound(hydro.path)
print myhydro.dm.pos3d[:,0].max()
#zoom_reg = np.where(myhydro.dm.mass==myhydro.dm.mass.min())
#nucenter = nbe.real_center(myhydro.dm.pos3d[zoom_reg], myhydro.dm.mass[zoom_reg])
myhydro.center_shift(hydro.c_dm_com)
myhydro.r_virial(600,n=1)
```
 
%%%% Output: stream
 
loading Dark matter..
loading Stars..
loading Gas..
36844.594
| r_200 = 212.70
| Diagonal matrix computed
| | 20, 0, 0|
| D =| 0, 14, 0|
| | 0, 0, 4|
 
%% Cell type:code id: tags:
 
``` python
print myhydro.gs.hsml.min()
print myhydro.dm.mass.min()
print myhydro.st.mass.min()
```
%%%% Output: stream
0.281102
1536979.0
183658.06
%% Cell type:code id: tags:
``` python
pos_dm = np.array(myhydro.dm.pos3d.reshape(len(myhydro.dm.pos3d)*3),dtype=np.float32)
pos_gs = np.array(myhydro.gs.pos3d.reshape(len(myhydro.gs.pos3d)*3),dtype=np.float32)
pos_st = np.array(myhydro.st.pos3d.reshape(len(myhydro.st.pos3d)*3),dtype=np.float32)
pos = np.concatenate((pos_dm, pos_st, pos_gs))
phi_cord = np.concatenate((myhydro.dm.phi,myhydro.st.phi, myhydro.gs.phi))
 
mass = np.concatenate((myhydro.dm.mass,myhydro.st.mass,myhydro.gs.mass))
v = np.concatenate((myhydro.dm.v,myhydro.st.v,myhydro.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: stream
 
5973141 5973141
 
%% Cell type:markdown id: tags:
 
# DMO
 
%% Cell type:code id: tags:
 
``` python
dmo = wkbl.astro.halo_info.Adicoradmo()
myDMO = wkbl.Galaxy_Hound(dmo.path)
myDMO.center_shift(dmo.c_dm_com)
myDMO.r_virial(600,n=20)
```
 
%%%% Output: stream
 
loading Dark matter..
 
%% Cell type:code id: tags:
 
``` python
print myDMO.r200
print myDMO.dm.mass.min()
print dmo.path
```
 
%%%% Output: stream
205.6640625
227993.1
/data/OWN/DMO/Adicora/output_00041
%% Cell type:code id: tags:
 
``` python
def eddingtong_from_file(path):
files = np.loadtxt(path)
return files[:,0], files[:,1]
 
```
 
%% 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)
 
def eddingtong_from_file(path):
files = np.loadtxt(path)
return files[:,0], files[:,1]
 
def fdv_plot_chi2_max_edd(ax, path, limmin, limmax,index=0,spherical=False, save=False,outname="/home/arturo/Pictures/ploto.png",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 = (myhydro.dm.r>limmin)&(myhydro.dm.r<limmax)
if width == None:
width = limmin /6.
disc = (np.abs(myhydro.dm.pos3d[:,2])<width)
local_v_wc = myhydro.dm.v[shell_wc]
v_disc_in_the_shell = myhydro.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(myhydro.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 = N = test = means= binsize = np.array([])
for i in range(len(bins_wc)-1):
pop = local_v_wc[(local_v_wc>bins_wc[i])&(local_v_wc<=bins_wc[i+1])]
if len(pop)==0:
continue
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)
means, N = np.append(means,np.average(pop)), np.append(N,len(pop))
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])]
means, N = np.append(means,np.average(pop)), np.append(N,len(pop))
dof = len(test)
# normalizations factor for data histogram
ptot = np.sum(N)
Ntot = np.trapz(np.nan_to_num(N),x=np.nan_to_num(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(myhydro.p.G * enclosed_m / ((limmin + limmax) /2. )) * myhydro.p.kpctokm
# for gas
vc_gas = np.mean(myhydro.gs.vphi[(myhydro.gs.R>limmin)&(myhydro.gs.R<limmax)&(np.abs(myhydro.gs.pos3d[:,2])<2.)])
std_vc_gas = np.std(myhydro.gs.vphi[(myhydro.gs.R>limmin)&(myhydro.gs.R<limmax)&(np.abs(myhydro.gs.pos3d[:,2])<2.)])
vc_dm = np.nanmean(myhydro.dm.v[(myhydro.dm.R>limmin)&(myhydro.dm.R<limmax)&(np.abs(myhydro.dm.pos3d[:,2])<2.)])
#for stars
vc_stars = np.nanmean(myhydro.st.v[(myhydro.st.R>limmin)&(myhydro.st.R<limmax)&(np.abs(myhydro.st.pos3d[:,2])<2.)])
print myhydro.st.v[(myhydro.st.R>limmin)&(myhydro.st.R<limmax)&(np.abs(myhydro.st.pos3d[:,2])<2.)]
std_vc_stars = np.nanstd(myhydro.st.v[(myhydro.st.r>limmin)&(myhydro.st.r<limmax)])
leno = len(myhydro.st.v[(myhydro.st.r<limmin)&(myhydro.st.r<limmax)])
# this sigma is a parameter in the fdv not an error of any kind
sigma_8 = vc / np.sqrt(2.)
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
# eddington
chi2 = np.nansum((((N/Ntot) - fv_the)**2 / (sigma**2))[(fv_the>0)])
# maxwellian
chi2_max = np.nansum((((N/Ntot) - maxwellian)**2 / (sigma**2))[(fv_the>0)])
#for i in range(len(N)):
# print test[i],((((N/Ntot) - fv_the)**2 / (fv_the**2)))[i],((((N/Ntot) - maxwellian)**2 / (maxwellian**2)))[i]
# v escape from potential
if (spherical):
v_esc , sig_vesc = hydro.v_esc_sph[index], hydro.v_esc_sigma_sph[index]
else:
v_esc , sig_vesc = hydro.v_esc_iso[index], hydro.v_esc_sigma_iso[index]
# plotting
ax.set_xlim([0,1.2*v_esc])
ax.set_ylim([0,1.3*fv_the.max()])
x=np.linspace(v[0],v[-2],100)
ax.plot(x,f(x),'k-',lw=1.6,label=r"$\rm 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,fmt=".", label="data")
#ax.set_xlabel(r"$|\vec{v}|$ [km/s]",fontsize=16)
ax.set_ylabel(r"$f(v) $ ",fontsize=22)
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 = r"$\rm Maxwellian$"+"\n"
texto_max += r" $\chi^2_{red} = $"
texto_max += r"${0:.2f}$".format(chi2_max/dof)
texto_edd = r"$\rm Eddington$"+"\n"
texto_edd += r" $\chi^2_{red} = $"
texto_edd += r"${0:.2f}$".format(chi2/dof)
 
texto_vc = r"$v_c ^{stars} =$ "
texto_vc += r"${0:.2f}$".format(vc_stars,leno)
texto_vc += "\n"
texto_vc += "$dof = {0}$".format(dof)
 
 
 
 
fig.text(0.15,0.88,texto_max,fontsize=16,color='r')
fig.text(0.45,0.88,texto_edd,fontsize=16,color='k')
print "---->>>",(vc_stars>0), vc_stars
if (vc_stars>=0):
fig.text(0.15,0.80,texto_vc,fontsize=16,color='k')
fig.text(0.7,0.70,r"$\rm Hydro$",fontsize=26,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-',lw=1.6,label=r"$\rm 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)
 
 
def fdv_plot_chi2_max_eddDMO(ax, path, limmin, limmax,index=0, spherical=False, save=False,outname="/home/arturo/Pictures/ploto.png",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 = (myDMO.dm.r>limmin)&(myDMO.dm.r<limmax)
if width == None:
width = limmin /6.
disc = (np.abs(myDMO.dm.pos3d[:,2])<width)
local_v_wc = myDMO.dm.v[shell_wc]
print local_v_wc
v_disc_in_the_shell = myDMO.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(myDMO.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 = N = s = test = means = 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])]
means, N, s= np.append(means,np.average(pop)), np.append(N,len(pop)), np.append(s,np.std(pop))
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])]
means, N, s= np.append(means,np.average(pop)), np.append(N,len(pop)), np.append(s,np.std(pop))
 
dof = len(test)
# normalizations factor for data histogram
ptot=np.sum(N)
Ntot = np.trapz(N,x=means)
## MAXWELLIAN Fdv
# circular velocity
enclosed_m = (np.sum(myDMO.dm.mass[np.where(myDMO.dm.r<limmin)])+ np.sum(myDMO.dm.mass[np.where(myDMO.dm.r<limmax)]))/2.
vc = np.sqrt(myDMO.p.G * enclosed_m / ((limmin + limmax) /2. )) * myDMO.p.kpctokm
# this sigma is a parameter in the fdv not an error of any kind
sigma_8 = vc / np.sqrt(2.)
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
# eddington
chi2 = np.sum((((N/Ntot) - fv_the)**2 / (sigma)**2)[:-1])
# maxwellian
chi2_max = np.sum((((N/Ntot) - maxwellian)**2 / (sigma**2))[:-1])
# v escape from potential
if (spherical):
v_esc , sig_vesc = dmo.v_esc_sph[index], dmo.v_esc_sigma_sph[index]
else:
v_esc , sig_vesc = dmo.v_esc_iso[index], dmo.v_esc_sigma_iso[index]
# plotting
#ax.set_xlim([0,1.2*v_esc])
ax.set_ylim([0,1.3*fv_the.max()])
x=np.linspace(v[0],v[-2],100)
ax.plot(x,f(x),'k-',lw=2,label=r"$\rm 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='#2a0972',lw=1.5,fmt=".", label="data")
ax.set_xlabel(r"$v$ [km/s]",fontsize=22)
ax.set_ylabel(r"$f(v) $ ",fontsize=22)
ax.axvline(x=v_esc,color='k',linestyle='--')
ax.axvspan(v_esc-sig_vesc, v_esc+sig_vesc, alpha=0.5, color='#721c09',label=r"$v_{esc} \pm 1\sigma$ ")
texto_max = r"$\rm Maxwellian$"+"\n"
texto_max += r" $\chi^2_{red} = $"
texto_max += r"${0:.2f}$".format(chi2_max/dof)
texto_edd = r"$\rm Eddington$"+"\n"
texto_edd += r" $\chi^2_{red} = $"
texto_edd += r"${0:.2f}$".format(chi2/dof)
fig.text(0.15,0.45,texto_max,fontsize=16,color='r')
fig.text(0.45,0.45,texto_edd,fontsize=16,color='k')
fig.text(0.7,0.25,r"$\rm DMO$",fontsize=26,color='k')
fig.text(0.45,0.4,r"$dof = {0}$".format(dof),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(),100)
ax.plot(max_v,get_maxw(max_v,sigma_8)/N0,'#721c09',lw=1.6,label=r"$\rm 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)
```
 
%%%% Output: stream
 
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots(2,1,sharey=True,sharex=True,figsize=[10,10])
path = "/home/arturo/Documents/LAM/LAM2LUPM/Adicora/fdvs/fv_Eddington_Adicora_DM_baryons_Rmax=1453.96kpc_3kpc_no_divergence.txt"
fdv_plot_chi2_max_edd(ax[0],path,2,4,index=0)
path = "/home/arturo/Documents/LAM/LAM2LUPM/Adicora/fdvs/fv_Eddington_Adicora_DMO_Rmax=1418.67kpc_3kpc_no_divergence.txt"
fdv_plot_chi2_max_eddDMO(ax[1],path,2,4,index=0)
fig.tight_layout(h_pad=-0.4)
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)
plt.savefig("/home/arturo/Documents/LAM/LAM2LUPM/"+hydro.namenospace+"/vdf_comparison_"+hydro.namenospace+"_3kpc.pdf")
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%%%% Output: stream
 
[ 67.539696 119.1637 291.34207 ... 303.4103 349.57852 339.4761 ]
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots(2,1,sharey=True,sharex=True,figsize=[10,10])
path = "/home/arturo/Documents/LAM/LAM2LUPM/Adicora/fdvs/fv_Eddington_Adicora_DM_baryons_Rmax=1453.96kpc_8kpc_no_divergence.txt"
fdv_plot_chi2_max_edd(ax[0],path,7,9,index=1)
path = "/home/arturo/Documents/LAM/LAM2LUPM/Adicora/fdvs/fv_Eddington_Adicora_DMO_Rmax=1418.67kpc_8kpc_no_divergence.txt"
fdv_plot_chi2_max_eddDMO(ax[1],path,7,9,index=1)
fig.tight_layout(h_pad=-0.4)
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)
plt.savefig("/home/arturo/Documents/LAM/LAM2LUPM/"+hydro.namenospace+"/vdf_comparison_"+hydro.namenospace+"_8kpc.pdf")
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%%%% Output: stream
 
[212.12239 243.36891 232.73643 ... 143.06313 134.6469 276.95 ]
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots(2,1,sharey=True,sharex=True,figsize=[10,10])
path = "/home/arturo/Documents/LAM/LAM2LUPM/Adicora/fdvs/fv_Eddington_Adicora_DM_baryons_Rmax=1453.96kpc_20kpc_no_divergence.txt"
fdv_plot_chi2_max_edd(ax[0],path,19,21,index=2)
path = "/home/arturo/Documents/LAM/LAM2LUPM/Adicora/fdvs/fv_Eddington_Adicora_DMO_Rmax=1418.67kpc_20kpc_no_divergence.txt"
fdv_plot_chi2_max_eddDMO(ax[1],path,19,21,index=2)
fig.tight_layout(h_pad=-0.4)
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)
plt.savefig("/home/arturo/Documents/LAM/LAM2LUPM/"+hydro.namenospace+"/vdf_comparison_"+hydro.namenospace+"_20kpc.pdf")
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%%%% Output: stream
 
[140.53859 337.40323 225.82458 ... 151.16046 65.96902 225.92229]
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots(2,1,sharey=True,sharex=True,figsize=[10,10])
path = "/home/arturo/Documents/LAM/LAM2LUPM/Adicora/fdvs/fv_Eddington_Adicora_DM_baryons_Rmax=1453.96kpc_50kpc_no_divergence.txt"
fdv_plot_chi2_max_edd(ax[0],path,49,51,index=3)
path = "/home/arturo/Documents/LAM/LAM2LUPM/Adicora/fdvs/fv_Eddington_Adicora_DMO_Rmax=1418.67kpc_50kpc_no_divergence.txt"
fdv_plot_chi2_max_eddDMO(ax[1],path,49,51,index=3)
fig.tight_layout(h_pad=-0.4)
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)
plt.savefig("/home/arturo/Documents/LAM/LAM2LUPM/"+hydro.namenospace+"/vdf_comparison_"+hydro.namenospace+"_50kpc.pdf")
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%%%% Output: stream
 
[129.005 328.8708 97.40869 ... 90.98444 237.75586 251.83408]
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots(2,1,sharey=True,sharex=True,figsize=[10,10])
path = "/home/arturo/Documents/LAM/LAM2LUPM/Adicora/fdvs/fv_Eddington_Adicora_DM_baryons_Rmax=1453.96kpc_150kpc_no_divergence.txt"
fdv_plot_chi2_max_edd(ax[0],path,149,151,index=4)
path = "/home/arturo/Documents/LAM/LAM2LUPM/Adicora/fdvs/fv_Eddington_Adicora_DMO_Rmax=1418.67kpc_150kpc_no_divergence.txt"
fdv_plot_chi2_max_eddDMO(ax[1],path,149,151,index=4)
fig.tight_layout(h_pad=-0.4)
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)
plt.savefig("/home/arturo/Documents/LAM/LAM2LUPM/"+hydro.namenospace+"/vdf_comparison_"+hydro.namenospace+"_150kpc.pdf")
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%%%% Output: stream
 
[]
---->>> False nan
[159.34361 129.08572 175.99571 ... 127.824905 69.82763 173.92186 ]
 
%% Cell type:code id: tags:
 
``` python
```
......
......@@ -4,7 +4,7 @@
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
"collapsed": false
},
"outputs": [],
"source": [
......@@ -147,7 +147,7 @@
},
{
"cell_type": "code",
"execution_count": 20,
"execution_count": 3,
"metadata": {
"collapsed": true
},
......@@ -158,7 +158,7 @@
},
{
"cell_type": "code",
"execution_count": 5,
"execution_count": 4,
"metadata": {
"collapsed": false
},
......@@ -177,6 +177,27 @@
"myDMO.r_virial(600,n=20)\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"163.4765625\n",
"275776.28\n"
]
}
],
"source": [
"print myDMO.r200\n",
"print myDMO.dm.mass.min()"
]
},
{
"cell_type": "code",
"execution_count": 6,
......
%% 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.interpolate import interp1d
from matplotlib.patches import Circle
from scipy.special import gamma
import numpy as np
import emcee
from mpl_toolkits.axes_grid1 import make_axes_locatable
from numpy import exp, sqrt
from scipy.integrate import quad, dblquad, nquad
import matplotlib.patches as patches
from itertools import product
from scipy.integrate import quad
import scipy.optimize as optimize
import matplotlib.pyplot as plt
import matplotlib as mpl
from sklearn.neighbors import KDTree
import sys
import lmfit
from py_unsio import *
import pymc
import os
from pymodelfit import FunctionModel1DAuto
import wkbl
from wkbl.astro.halo_info import *
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm
import wkbl.astro.nbody_essentials as nbe
import cfalcon
CF =cfalcon.CFalcon()
import iminuit
from iminuit import Minuit, describe, Struct
import probfit
from matplotlib.colors import LogNorm
from matplotlib.ticker import FormatStrFormatter
import warnings
warnings.filterwarnings('ignore')
```
 
%% Cell type:markdown id: tags:
 
# Hydro
 
%% Cell type:code id: tags:
 
``` python
hydro = wkbl.astro.halo_info.HALOBHydro()
```
 
%% Cell type:code id: tags:
 
``` python
simname=hydro.name
myhydro = wkbl.Galaxy_Hound(hydro.path)
print myhydro.dm.pos3d[:,0].max()
myhydro.center_shift(hydro.c_dm_com)
myhydro.r_virial(600,n=1)
```
 
%%%% Output: stream
 
loading Dark matter..
loading Stars..
loading Gas..
19879.156
| r_200 = 177.54
| Diagonal matrix computed
| | 20, 0, 0|
| D =| 0, 14, 0|
| | 0, 0, 2|
 
%% Cell type:code id: tags:
 
``` python
pos_dm = np.array(myhydro.dm.pos3d.reshape(len(myhydro.dm.pos3d)*3),dtype=np.float32)
pos_gs = np.array(myhydro.gs.pos3d.reshape(len(myhydro.gs.pos3d)*3),dtype=np.float32)
pos_st = np.array(myhydro.st.pos3d.reshape(len(myhydro.st.pos3d)*3),dtype=np.float32)
pos = np.concatenate((pos_dm, pos_st, pos_gs))
phi_cord = np.concatenate((myhydro.dm.phi,myhydro.st.phi, myhydro.gs.phi))
 
mass = np.concatenate((myhydro.dm.mass,myhydro.st.mass,myhydro.gs.mass))
v = np.concatenate((myhydro.dm.v,myhydro.st.v,myhydro.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: stream
 
22910217 22910217
 
%% Cell type:markdown id: tags:
 
# DMO
 
%% Cell type:code id: tags:
 
``` python
dmo = wkbl.astro.halo_info.HALOBdmo()
```
 
%% Cell type:code id: tags:
 
``` python
myDMO = wkbl.Galaxy_Hound(dmo.path)
myDMO.center_shift(dmo.c_dm_com)
myDMO.r_virial(600,n=20)
```
 
%%%% Output: stream
 
loading Dark matter..
 
%% Cell type:code id: tags:
 
``` python
print myDMO.r200
print myDMO.dm.mass.min()
```
%%%% Output: stream
163.4765625
275776.28
%% Cell type:code id: tags:
``` python
def eddingtong_from_file(path):
files = np.loadtxt(path)
return files[:,0], files[:,1]
 
```
 
%% Cell type:code id: tags:
 
``` python
print hydro.name
```
 
%%%% Output: stream
 
Halo B
 
%% 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)
 
def eddingtong_from_file(path):
files = np.loadtxt(path)
return files[:,0], files[:,1]
 
def fdv_plot_chi2_max_edd(ax, path, limmin, limmax,index=0,spherical=False, save=False,outname="/home/arturo/Pictures/ploto.png",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 = (myhydro.dm.r>limmin)&(myhydro.dm.r<limmax)
if width == None:
width = limmin /6.
disc = (np.abs(myhydro.dm.pos3d[:,2])<width)
local_v_wc = myhydro.dm.v[shell_wc]
v_disc_in_the_shell = myhydro.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(myhydro.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 = N = test = means= binsize = np.array([])
for i in range(len(bins_wc)-1):
pop = local_v_wc[(local_v_wc>bins_wc[i])&(local_v_wc<=bins_wc[i+1])]
if len(pop)==0:
continue
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)
means, N = np.append(means,np.average(pop)), np.append(N,len(pop))
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])]
means, N = np.append(means,np.average(pop)), np.append(N,len(pop))
dof = len(test)
# normalizations factor for data histogram
ptot = np.sum(N)
Ntot = np.trapz(np.nan_to_num(N),x=np.nan_to_num(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(myhydro.p.G * enclosed_m / ((limmin + limmax) /2. )) * myhydro.p.kpctokm
# for gas
vc_gas = np.mean(myhydro.gs.vphi[(myhydro.gs.R>limmin)&(myhydro.gs.R<limmax)&(np.abs(myhydro.gs.pos3d[:,2])<2.)])
std_vc_gas = np.std(myhydro.gs.vphi[(myhydro.gs.R>limmin)&(myhydro.gs.R<limmax)&(np.abs(myhydro.gs.pos3d[:,2])<2.)])
vc_dm = np.nanmean(myhydro.dm.v[(myhydro.dm.R>limmin)&(myhydro.dm.R<limmax)&(np.abs(myhydro.dm.pos3d[:,2])<2.)])
#for stars
vc_stars = np.nanmean(myhydro.st.v[(myhydro.st.R>limmin)&(myhydro.st.R<limmax)&(np.abs(myhydro.st.pos3d[:,2])<2.)])
std_vc_stars = np.nanstd(myhydro.st.v[(myhydro.st.r>limmin)&(myhydro.st.r<limmax)])
leno = len(myhydro.st.v[(myhydro.st.r<limmin)&(myhydro.st.r<limmax)])
# this sigma is a parameter in the fdv not an error of any kind
sigma_8 = vc / np.sqrt(2.)
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
# eddington
chi2 = np.nansum((((N/Ntot) - fv_the)**2 / (sigma**2))[(fv_the>0)])
# maxwellian
chi2_max = np.nansum((((N/Ntot) - maxwellian)**2 / (sigma**2))[(fv_the>0)])
#for i in range(len(N)):
# print test[i],((((N/Ntot) - fv_the)**2 / (fv_the**2)))[i],((((N/Ntot) - maxwellian)**2 / (maxwellian**2)))[i]
# v escape from potential
if (spherical):
v_esc , sig_vesc = hydro.v_esc_sph[index], hydro.v_esc_sigma_sph[index]
else:
v_esc , sig_vesc = hydro.v_esc_iso[index], hydro.v_esc_sigma_iso[index]
# plotting
ax.set_xlim([0,1.2*v_esc])
ax.set_ylim([0,1.3*fv_the.max()])
x=np.linspace(v[0],v[-2],100)
ax.plot(x,f(x),'k-',lw=1.6,label=r"$\rm 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,fmt=".", label="data")
#ax.set_xlabel(r"$|\vec{v}|$ [km/s]",fontsize=16)
ax.set_ylabel(r"$f(v) $ ",fontsize=22)
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 = r"$\rm Maxwellian$"+"\n"
texto_max += r" $\chi^2_{red} = $"
texto_max += r"${0:.2f}$".format(chi2_max/dof)
texto_edd = r"$\rm Eddington$"+"\n"
texto_edd += r" $\chi^2_{red} = $"
texto_edd += r"${0:.2f}$".format(chi2/dof)
texto_vc = r"$v_c ^{stars} =$ "
texto_vc += r"${0:.2f}$".format(vc_stars,leno)
texto_vc += "\n"
texto_vc += "$dof = {0}$".format(dof)
 
 
 
fig.text(0.15,0.88,texto_max,fontsize=16,color='r')
fig.text(0.45,0.88,texto_edd,fontsize=16,color='k')
if (vc_stars>=0):
fig.text(0.15,0.80,texto_vc,fontsize=16,color='k')
fig.text(0.7,0.70,r"$\rm Hydro$",fontsize=26,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-',lw=1.6,label=r"$\rm 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)
 
 
def fdv_plot_chi2_max_eddDMO(ax, path, limmin, limmax,index=0, spherical=False, save=False,outname="/home/arturo/Pictures/ploto.png",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 = (myDMO.dm.r>limmin)&(myDMO.dm.r<limmax)
if width == None:
width = limmin /6.
disc = (np.abs(myDMO.dm.pos3d[:,2])<width)
local_v_wc = myDMO.dm.v[shell_wc]
v_disc_in_the_shell = myDMO.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(myDMO.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 = N = s = test = means = 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])]
means, N, s= np.append(means,np.average(pop)), np.append(N,len(pop)), np.append(s,np.std(pop))
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])]
means, N, s= np.append(means,np.average(pop)), np.append(N,len(pop)), np.append(s,np.std(pop))
 
dof = len(test)
# normalizations factor for data histogram
ptot=np.sum(N)
Ntot = np.trapz(N,x=means)
## MAXWELLIAN Fdv
# circular velocity
enclosed_m = (np.sum(myDMO.dm.mass[np.where(myDMO.dm.r<limmin)])+ np.sum(myDMO.dm.mass[np.where(myDMO.dm.r<limmax)]))/2.
vc = np.sqrt(myDMO.p.G * enclosed_m / ((limmin + limmax) /2. )) * myDMO.p.kpctokm
# this sigma is a parameter in the fdv not an error of any kind
sigma_8 = vc / np.sqrt(2.)
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
# eddington
chi2 = np.sum((((N/Ntot) - fv_the)**2 / (sigma)**2)[:-1])
# maxwellian
chi2_max = np.sum((((N/Ntot) - maxwellian)**2 / (sigma**2))[:-1])
# v escape from potential
if (spherical):
v_esc , sig_vesc = dmo.v_esc_sph[index], dmo.v_esc_sigma_sph[index]
else:
v_esc , sig_vesc = dmo.v_esc_iso[index], dmo.v_esc_sigma_iso[index]
# plotting
#ax.set_xlim([0,1.2*v_esc])
ax.set_ylim([0,1.3*fv_the.max()])
x=np.linspace(v[0],v[-2],100)
ax.plot(x,f(x),'k-',lw=2,label=r"$\rm 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='#2a0972',lw=1.5,fmt=".", label="data")
ax.set_xlabel(r"$v$ [km/s]",fontsize=22)
ax.set_ylabel(r"$f(v) $ ",fontsize=22)
ax.axvline(x=v_esc,color='k',linestyle='--')
ax.axvspan(v_esc-sig_vesc, v_esc+sig_vesc, alpha=0.5, color='#721c09',label=r"$v_{esc} \pm 1\sigma$ ")
texto_max = r"$\rm Maxwellian$"+"\n"
texto_max += r" $\chi^2_{red} = $"
texto_max += r"${0:.2f}$".format(chi2_max/dof)
texto_edd = r"$\rm Eddington$"+"\n"
texto_edd += r" $\chi^2_{red} = $"
texto_edd += r"${0:.2f}$".format(chi2/dof)
fig.text(0.15,0.45,texto_max,fontsize=16,color='r')
fig.text(0.45,0.45,texto_edd,fontsize=16,color='k')
fig.text(0.7,0.25,r"$\rm DMO$",fontsize=26,color='k')
fig.text(0.45,0.4,r"$dof = {0}$".format(dof),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(),100)
ax.plot(max_v,get_maxw(max_v,sigma_8)/N0,'#721c09',lw=1.6,label=r"$\rm 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(2,1,sharey=True,sharex=True,figsize=[10,10])
path = "/home/arturo/Documents/LAM/LAM2LUPM/HALOB/fdvs/fv_Eddington_halo_B_DM_baryons_Rmax=881.61kpc_3kpc_no_divergence.txt"
fdv_plot_chi2_max_edd(ax[0],path,2,4,index=0)
path = "/home/arturo/Documents/LAM/LAM2LUPM/HALOB/fdvs/fv_Eddington_halo_B_DMO_Rmax=785.94kpc_2kpc_no_divergence.txt"
fdv_plot_chi2_max_eddDMO(ax[1],path,1,3,index=0)
fig.tight_layout(h_pad=-0.4)
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)
plt.savefig("/home/arturo/Documents/LAM/LAM2LUPM/"+hydro.namenospace+"/VdfComparison/vdf_comparison_"+hydro.namenospace+"_3kpc.pdf")
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots(2,1,sharey=True,sharex=True,figsize=[10,10])
path = "/home/arturo/Documents/LAM/LAM2LUPM/HALOB/fdvs/fv_Eddington_halo_B_DM_baryons_Rmax=881.61kpc_8kpc_no_divergence.txt"
fdv_plot_chi2_max_edd(ax[0],path,7,9,index=1)
path = "/home/arturo/Documents/LAM/LAM2LUPM/HALOB/fdvs/fv_Eddington_halo_B_DMO_Rmax=785.94kpc_8kpc_no_divergence.txt"
fdv_plot_chi2_max_eddDMO(ax[1],path,7,9,index=1)
fig.tight_layout(h_pad=-0.4)
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)
plt.savefig("/home/arturo/Documents/LAM/LAM2LUPM/"+hydro.namenospace+"/VdfComparison/vdf_comparison_"+hydro.namenospace+"_8kpc.pdf")
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots(2,1,sharey=True,sharex=True,figsize=[10,10])
path = "/home/arturo/Documents/LAM/LAM2LUPM/HALOB/fdvs/fv_Eddington_halo_B_DM_baryons_Rmax=881.61kpc_20kpc_no_divergence.txt"
fdv_plot_chi2_max_edd(ax[0],path,19,21,index=2)
path = "/home/arturo/Documents/LAM/LAM2LUPM/HALOB/fdvs/fv_Eddington_halo_B_DMO_Rmax=785.94kpc_20kpc_no_divergence.txt"
fdv_plot_chi2_max_eddDMO(ax[1],path,19,21,index=2)
fig.tight_layout(h_pad=-0.4)
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)
plt.savefig("/home/arturo/Documents/LAM/LAM2LUPM/"+hydro.namenospace+"/VdfComparison/vdf_comparison_"+hydro.namenospace+"_20kpc.pdf")
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots(2,1,sharey=True,sharex=True,figsize=[10,10])
path = "/home/arturo/Documents/LAM/LAM2LUPM/HALOB/fdvs/fv_Eddington_halo_B_DM_baryons_Rmax=881.61kpc_50kpc_no_divergence.txt"
fdv_plot_chi2_max_edd(ax[0],path,49,51,index=3)
path = "/home/arturo/Documents/LAM/LAM2LUPM/HALOB/fdvs/fv_Eddington_halo_B_DMO_Rmax=785.94kpc_50kpc_no_divergence.txt"
fdv_plot_chi2_max_eddDMO(ax[1],path,49,51,index=3)
fig.tight_layout(h_pad=-0.4)
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)
plt.savefig("/home/arturo/Documents/LAM/LAM2LUPM/"+hydro.namenospace+"/VdfComparison/vdf_comparison_"+hydro.namenospace+"_50kpc.pdf")
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%% Cell type:code id: tags:
 
``` python
fig, ax = plt.subplots(2,1,sharey=True,sharex=True,figsize=[10,10])
path = "/home/arturo/Documents/LAM/LAM2LUPM/HALOB/fdvs/fv_Eddington_halo_B_DM_baryons_Rmax=881.61kpc_150kpc_no_divergence.txt"
fdv_plot_chi2_max_edd(ax[0],path,149,151,index=4)
path = "/home/arturo/Documents/LAM/LAM2LUPM/HALOB/fdvs/fv_Eddington_halo_B_DMO_Rmax=785.94kpc_150kpc_no_divergence.txt"
fdv_plot_chi2_max_eddDMO(ax[1],path,149,151,index=4)
fig.tight_layout(h_pad=-0.4)
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)
plt.savefig("/home/arturo/Documents/LAM/LAM2LUPM/"+hydro.namenospace+"/VdfComparison/vdf_comparison_"+hydro.namenospace+"_150kpc.pdf")
```
 
%%%% Output: display_data
 
 
%%%% Output: display_data
 
 
%% Cell type:code id: tags:
 
``` python
```
 
%% Cell type:code id: tags:
 
``` python
```
......
......@@ -101,7 +101,7 @@
},
{
"cell_type": "code",
"execution_count": 4,
"execution_count": 10,
"metadata": {
"collapsed": false
},
......@@ -110,12 +110,25 @@
"name": "stdout",
"output_type": "stream",
"text": [
"24303507 24303507\n"
"0.15167238\n",
"230812.75\n"
]
}
],
"source": [
"pos_dm = np.array(myhydro.dm.pos3d.reshape(len(myhydro.dm.pos3d)*3),dtype=np.float32)\n",
"print myhydro.gs.hsml.min()\n",
"print myhydro.dm.mass.min()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# pos_dm = np.array(myhydro.dm.pos3d.reshape(len(myhydro.dm.pos3d)*3),dtype=np.float32)\n",
"pos_gs = np.array(myhydro.gs.pos3d.reshape(len(myhydro.gs.pos3d)*3),dtype=np.float32)\n",
"pos_st = np.array(myhydro.st.pos3d.reshape(len(myhydro.st.pos3d)*3),dtype=np.float32)\n",
"pos = np.concatenate((pos_dm, pos_st, pos_gs))\n",
......@@ -173,7 +186,7 @@
},
{
"cell_type": "code",
"execution_count": 7,
"execution_count": 12,
"metadata": {
"collapsed": false
},
......@@ -182,12 +195,14 @@
"name": "stdout",
"output_type": "stream",
"text": [
"Halo C\n"
"176.3671875\n",
"275776.28\n"
]
}
],
"source": [
"print hydro.name"
"print myDMO.r200\n",
"print myDMO.dm.mass.min()"
]
},
{
......
%% 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.interpolate import interp1d
from matplotlib.patches import Circle
from scipy.special import gamma
import numpy as np
import emcee
from mpl_toolkits.axes_grid1 import make_axes_locatable
from numpy import exp, sqrt
from scipy.integrate import quad, dblquad, nquad
import matplotlib.patches as patches
from itertools import product
from scipy.integrate import quad
import scipy.optimize as optimize
import matplotlib.pyplot as plt
import matplotlib as mpl
from sklearn.neighbors import KDTree
import sys
import lmfit
from py_unsio import *
import pymc
import os
from pymodelfit import FunctionModel1DAuto
import wkbl
from wkbl.astro.halo_info import *
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm
import wkbl.astro.nbody_essentials as nbe
import cfalcon
CF =cfalcon.CFalcon()
import iminuit
from iminuit import Minuit, describe, Struct
import probfit
from matplotlib.colors import LogNorm
from matplotlib.ticker import FormatStrFormatter
import warnings
warnings.filterwarnings('ignore')
```
 
%% Cell type:markdown id: tags:
 
# Hydro
 
%% Cell type:code id: tags:
 
``` python
hydro = wkbl.astro.halo_info.HALOCHydro()
simname=hydro.name
myhydro = wkbl.Galaxy_Hound(hydro.path)
print myhydro.dm.pos3d[:,0].max()
myhydro.center_shift(hydro.c_dm_com)
myhydro.r_virial(600,n=1)
```
 
%%%% Output: stream
 
loading Dark matter..
loading Stars..
loading Gas..
19879.998
| r_200 = 182.23
| Diagonal matrix computed
| | 18, 0, 0|
| D =| 0, 13, 0|
| | 0, 0, 3|
 
%% Cell type:code id: tags:
 
``` python
pos_dm = np.array(myhydro.dm.pos3d.reshape(len(myhydro.dm.pos3d)*3),dtype=np.float32)
print myhydro.gs.hsml.min()
print myhydro.dm.mass.min()
```
%%%% Output: stream
0.15167238
230812.75
%% Cell type:code id: tags:
``` python
# pos_dm = np.array(myhydro.dm.pos3d.reshape(len(myhydro.dm.pos3d)*3),dtype=np.float32)
pos_gs = np.array(myhydro.gs.pos3d.reshape(len(myhydro.gs.pos3d)*3),dtype=np.float32)
pos_st = np.array(myhydro.st.pos3d.reshape(len(myhydro.st.pos3d)*3),dtype=np.float32)
pos = np.concatenate((pos_dm, pos_st, pos_gs))
phi_cord = np.concatenate((myhydro.dm.phi,myhydro.st.phi, myhydro.gs.phi))
 
mass = np.concatenate((myhydro.dm.mass,myhydro.st.mass,myhydro.gs.mass))
v = np.concatenate((myhydro.dm.v,myhydro.st.v,myhydro.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: stream
24303507 24303507
%% Cell type:markdown id: tags:
 
# DMO
 
%% Cell type:code id: tags:
 
``` python
dmo = wkbl.astro.halo_info.HALOCdmo()
myDMO = wkbl.Galaxy_Hound(dmo.path)
myDMO.center_shift(dmo.c_dm_com)
myDMO.r_virial(600,n=20)
```
 
%%%% Output: stream
 
loading Dark matter..
 
%% Cell type:code id: tags:
 
``` python
def eddingtong_from_file(path):
files = np.loadtxt(path)
return files[:,0], files[:,1]
 
```
 
%% Cell type:code id: tags:
 
``` python
print hydro.name
print myDMO.r200
print myDMO.dm.mass.min()
```
 
%%%% Output: stream
 
Halo C
176.3671875
275776.28
 
%% 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)
 
def eddingtong_from_file(path):
files = np.loadtxt(path)
return files[:,0], files[:,1]
 
def fdv_plot_chi2_max_edd(ax, path, limmin, limmax,index=0,spherical=False, save=False,outname="/home/arturo/Pictures/ploto.png",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 = (myhydro.dm.r>limmin)&(myhydro.dm.r<limmax)
if width == None:
width = limmin /6.
disc = (np.abs(myhydro.dm.pos3d[:,2])<width)
local_v_wc = myhydro.dm.v[shell_wc]
v_disc_in_the_shell = myhydro.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(myhydro.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 = N = test = means= binsize = np.array([])
for i in range(len(bins_wc)-1):
pop = local_v_wc[(local_v_wc>bins_wc[i])&(local_v_wc<=bins_wc[i+1])]
if len(pop)==0:
continue
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)
means, N = np.append(means,np.average(pop)), np.append(N,len(pop))
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])]
means, N = np.append(means,np.average(pop)), np.append(N,len(pop))
dof = len(test)
# normalizations factor for data histogram
ptot = np.sum(N)
Ntot = np.trapz(np.nan_to_num(N),x=np.nan_to_num(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(myhydro.p.G * enclosed_m / ((limmin + limmax) /2. )) * myhydro.p.kpctokm
# for gas
vc_gas = np.mean(myhydro.gs.vphi[(myhydro.gs.R>limmin)&(myhydro.gs.R<limmax)&(np.abs(myhydro.gs.pos3d[:,2])<2.)])
std_vc_gas = np.std(myhydro.gs.vphi[(myhydro.gs.R>limmin)&(myhydro.gs.R<limmax)&(np.abs(myhydro.gs.pos3d[:,2])<2.)])
vc_dm = np.nanmean(myhydro.dm.v[(myhydro.dm.R>limmin)&(myhydro.dm.R<limmax)&(np.abs(myhydro.dm.pos3d[:,2])<2.)])
#for stars
vc_stars = np.nanmean(myhydro.st.v[(myhydro.st.R>limmin)&(myhydro.st.R<limmax)&(np.abs(myhydro.st.pos3d[:,2])<2.)])
std_vc_stars = np.nanstd(myhydro.st.v[(myhydro.st.r>limmin)&(myhydro.st.r<limmax)])
leno = len(myhydro.st.v[(myhydro.st.r<limmin)&(myhydro.st.r<limmax)])
# this sigma is a parameter in the fdv not an error of any kind
sigma_8 = vc / np.sqrt(2.)
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
# eddington
chi2 = np.nansum((((N/Ntot) - fv_the)**2 / (sigma**2))[(fv_the>0)])
# maxwellian
chi2_max = np.nansum((((N/Ntot) - maxwellian)**2 / (sigma**2))[(fv_the>0)])
#for i in range(len(N)):
# print test[i],((((N/Ntot) - fv_the)**2 / (fv_the**2)))[i],((((N/Ntot) - maxwellian)**2 / (maxwellian**2)))[i]
# v escape from potential
if (spherical):
v_esc , sig_vesc = hydro.v_esc_sph[index], hydro.v_esc_sigma_sph[index]
else:
v_esc , sig_vesc = hydro.v_esc_iso[index], hydro.v_esc_sigma_iso[index]
# plotting
ax.set_xlim([0,1.2*v_esc])
ax.set_ylim([0,1.3*fv_the.max()])
x=np.linspace(v[0],v[-2],100)
ax.plot(x,f(x),'k-',lw=1.6,label=r"$\rm 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,fmt=".", label="data")
#ax.set_xlabel(r"$|\vec{v}|$ [km/s]",fontsize=16)
ax.set_ylabel(r"$f(v) $ ",fontsize=22)
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 = r"$\rm Maxwellian$"+"\n"
texto_max += r" $\chi^2_{red} = $"
texto_max += r"${0:.2f}$".format(chi2_max/dof)
texto_edd = r"$\rm Eddington$"+"\n"
texto_edd += r" $\chi^2_{red} = $"
texto_edd += r"${0:.2f}$".format(chi2/dof)
texto_vc = r"$v_c ^{stars} =$ "
texto_vc += r"${0:.2f}$".format(vc_stars,leno)
texto_vc += "\n"
texto_vc += "$dof = {0}$".format(dof)
 
 
 
fig.text(0.15,0.88,texto_max,fontsize=16,color='r')
fig.text(0.45,0.88,texto_edd,fontsize=16,color='k')
if (vc_stars>=0):
fig.text(0.15,0.80,texto_vc,fontsize=16,color='k')
fig.text(0.7,0.70,r"$\rm Hydro$",fontsize=26,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-',lw=1.6,label=r"$\rm 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)
 
 
def fdv_plot_chi2_max_eddDMO(ax, path, limmin, limmax,index=0, spherical=False, save=False,outname="/home/arturo/Pictures/ploto.png",width=None):