q_parametre_hydro.py 3.55 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93
import numpy as np
import sys
from py_unsio import *
import os
import wkbl
import wkbl.astro.nbody_essentials as nbe
import cfalcon
CF =cfalcon.CFalcon()
import warnings
warnings.filterwarnings('ignore')


path = sys.argv[1]
simname = sys.argv[2]
"""
myhalo = wkbl.Galaxy_Hound(path,"halo,gas,stars")
print "centering recursively on Zoom Region COM"
ok,rho,_= CF.getDensity(np.array(myhalo.st.pos3d.reshape(len(myhalo.st.pos3d)*3),dtype=np.float32), myhalo.st.mass)
centro_rho = myhalo.st.pos3d[np.where(rho == rho.max())][0]
print "density",centro_rho
myhalo.center_shift(centro_rho)
myhalo.r_virial(600,n=2)
"""
myhalo= wkbl.Galaxy_Hound(path)
print "loaded"

zoomreg = np.where(myhalo.dm.mass==myhalo.dm.mass.min())
centro = nbe.real_center(myhalo.dm.pos3d[zoomreg],myhalo.dm.mass[zoomreg])
print centro
myhalo.center_shift(centro)
myhalo.r_virial(600,n=2)
myhalo.r200
myhalo.redefine(25)
print "recentering"
centro = nbe.real_center(myhalo.dm.pos3d[myhalo.dm.r<50],myhalo.dm.mass[myhalo.dm.r<50])
print centro
myhalo.center_shift(centro)
myhalo.r_virial(600,n=2)
myGkm = 6.673e-11*(1e-3**3)*myhalo.p.msuntokg#km^ 3 Msun^-1 s^-2

#"""

print "r200 = {0}".format(myhalo.r200)
print "M_dm_200 = {0:.5e}".format(np.sum(myhalo.dm.mass[myhalo.dm.r<myhalo.r200]))
print "M_st_200 = {0:.5e}".format(np.sum(myhalo.st.mass[myhalo.st.r<myhalo.r200]))
print "M_dm_0.1 = {0:.5e}".format(np.sum(myhalo.st.mass[myhalo.st.r<myhalo.r200*0.1]))
print "M_st_Fire = {0:.5e}".format(myhalo.st.fire_m)
pos_dm = np.array(myhalo.dm.pos3d.reshape(len(myhalo.dm.pos3d)*3),dtype=np.float32)
pos_gs = np.array(myhalo.gs.pos3d.reshape(len(myhalo.gs.pos3d)*3),dtype=np.float32)
pos_st = np.array(myhalo.st.pos3d.reshape(len(myhalo.st.pos3d)*3),dtype=np.float32)
pos = np.concatenate((pos_dm, pos_st))#, pos_gs))
mass = np.concatenate((myhalo.dm.mass,myhalo.st.mass))#,myhalo.gs.mass))
v = np.concatenate((myhalo.dm.v,myhalo.st.v))#,myhalo.gs.v))
print "getting potential..."
pos3d = pos.reshape(len(pos)/3,3)
r = np.sqrt(pos3d[:,0]**2 + pos3d[:,1]**2 +pos3d[:,2]**2)
inside_halo = np.where(r<(myhalo.r200))
r_sorted = np.argsort(r)
M_i = np.cumsum(mass[r_sorted]) - mass[r_sorted]
m_i = mass[r_sorted]
r_i = r[r_sorted]*(1e-2*myhalo.p.pctocm)# in km
U =  np.sum(-myGkm*M_i*m_i/r_i)
K_thermal_gas = np.sum(myhalo.gs.mass[myhalo.gs.r<myhalo.r200]*myhalo.gs.temp[myhalo.gs.r<myhalo.r200])*(3./2.)
K_dm_st = np.sum(mass[inside_halo]*(v[inside_halo])**2)/2.
K = K_dm_st #+ K_thermal_gas
print K_dm_st, K_thermal_gas




kmtokpc = 1 / 3.08567758128e+16 
q = (2. * K / U ) + 1. 

print "\n\n\n"
print "______________________________________________"
print "|                                             |"
#K = np.sum(mass[inside_halo]*(v[inside_halo]*kmtokpc)**2)
print "|     2E_kin = {0:1.3e} m_sun kms^1 s^-2     |".format(K)
print "|                                             |"

print "| First simplification from Shapiro2004       |"
print "|                                             |"
print "|    E_pot = {0:1.3e} m_sun km^2 s^-2       |".format(U)
print "|                                             |"
print "| q parameter a la Zjupa2016                  |"
print "|  (q = 1+ (2Ekin/Epot2))                     |"
print "|                                             |"
print "| q = {0:.5f}                                |".format(q) 
print "| q parameter a la Shapiro2004                |"
print "|  (q = 2Ekin / |Epot|)                       |"
print "| q = {0:.5f}                                 |".format(K/np.abs(U)) 
print "|_____________________________________________|"
#potmin = np.where(Phy[r2<myhalo.r200**2]==Phy[r2<myhalo.r200**2].min())