Commit 360bf2ce authored by NUNEZ Arturo's avatar NUNEZ Arturo
Browse files

full fits for mochima 2

parent 206d93b5
This diff is collapsed.
......@@ -7,7 +7,9 @@
"collapsed": true
},
"outputs": [],
"source": []
"source": [
"import pymc, # load the model file\n"
]
}
],
"metadata": {
......
This diff is collapsed.
{
"cells": [],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 1
}
This diff is collapsed.
import numpy as np
from py_unsio import *
import pymc
import wkbl
import wkbl.astro.nbody_essentials as nbe
import cfalcon
CF =cfalcon.CFalcon()
import warnings
warnings.filterwarnings('ignore')
path = "/data/OWN/SF1test/AnH3/output_00041"
#path = "/media/arturo/ARTUROTECA/OUTPUTS/HaloB/output_00417"
myDMO = wkbl.Galaxy_Hound(path)
print "centering"
zoom_reg = np.where(myDMO.dm.mass == myDMO.dm.mass.min())
nucenter = nbe.real_center(myDMO.dm.pos3d[zoom_reg], myDMO.dm.mass[zoom_reg])
myDMO.center_shift(nucenter)
myDMO.r_virial(600)
print "done r200 = {0}".format(myDMO.r200)
ok,myDMO.dm.rho,_= CF.getDensity(np.array(myDMO.dm.pos3d.reshape(len(myDMO.dm.pos3d)*3),dtype=np.float32), myDMO.dm.mass)
Pcrit = myDMO.dm._p.rho_crit
Mdm = myDMO.dm.mass.min()
myradiuses = myDMO.dm.r[np.argsort(myDMO.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= myDMO.gs.hsml.min()
# R array logarithmic Bining
r_p = np.logspace(np.log10(hsml),np.log10(2.5*myDMO.r200),30)
n_dm,r = np.histogram(myDMO.dm.r,bins=r_p)
vol = np.array([])
r1,r2 =r[:-1],r[1:]
vol = 4.* np.pi * ((r2**3)-(r1**3)) / 3.
profileDMO = n_dm*myDMO.dm.mass[0]/vol
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 = np.array([])
std = np.array([])
n=np.array([])
for i in range(len(r_p)-1):
shell = np.where((myDMO.dm.r > r_p[i])&(myDMO.dm.r < r_p[i+1])&(myDMO.dm.r > hsml))
n = np.append(n,len(shell[0]))
mean = np.append(mean,np.mean(myDMO.dm.rho[shell]))
std = np.append(std,np.std(myDMO.dm.rho[shell]))
#################################################################################
def abg_profile(x,po,r_s,al,be,ga):
x = 10**x
power = (be - ga) / al
denominator = ((x/r_s)**ga) * ((1 + (x / r_s)**al)**power)
return (10**po) / denominator
#################################################################################
#priors
po = pymc.Uniform('po', 2.0, 15.0, value= 8.0)
r_s = pymc.Uniform('r_s', 1.0, 100.0, value= 10.0)
al = pymc.Uniform('al', 0.01, 5.0, value= 1.0)
be = pymc.Uniform('be', 0.01, 10.0, value= 3.0)
ga = pymc.Uniform('ga', 0.01, 5.0, value= 1.0)
#model
@pymc.deterministic(plot=False)
def abg_profile(x=np.log10(r),po=po,r_s=r_s,al=al,be=be,ga=ga):
power = (be - ga) / al
denominator = ((x/r_s)**ga) * ((1 + (x / r_s)**al)**power)
return np.log10((10**po) / denominator)
err = 1/np.sqrt(n_dm)
#likelihood
rho = pymc.Normal('rho', mu=abg_profile, tau=1, value=np.log10(profileDMO), observed=True)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment