q_parametre_simplifie.py 2.62 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
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]
myDMO = wkbl.Galaxy_Hound(path)
print "centering recursively on Zoom Region COM"
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])
print nucenter
myDMO.center_shift(nucenter)
myDMO.r_virial(600,n=3)
myDMO.r200
myDMO.redefine(25)
print "recentering"
centro = nbe.real_center(myDMO.dm.pos3d[myDMO.dm.r<50],myDMO.dm.mass[myDMO.dm.r<50])
print centro
myDMO.center_shift(centro)
myDMO.r_virial(600,n=1)
myDMO.r200
myDMO.redefine(1)

print "r200 = {0}".format(myDMO.r200)
print "M_dm_200 = {0}".format(np.sum(myDMO.dm.mass[myDMO.dm.r<myDMO.r200]))

myDMO.redefine(1)
myGkm = 6.673e-11*(1e-3**3)*myDMO.p.msuntokg#km^ 3 Msun^-1 s^-2


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

r_sorted = np.argsort(myDMO.dm.r)
M_i = np.cumsum(myDMO.dm.mass[r_sorted]) - myDMO.dm.mass[r_sorted]
m_i = myDMO.dm.mass[r_sorted]
r_i = myDMO.dm.r[r_sorted]*(1e-2*myDMO.p.pctocm)# in km
U =  np.sum(-myGkm*M_i*m_i/r_i)
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((K/U)+1)     
print "| q parameter a la Shapiro2004                |"
print "|  (q = 2Ekin / |Epot|)                       |"
print "| q = {0:.5f}                                 |".format(K/np.abs(U)) 
print "|_____________________________________________|"
pos = np.array(myDMO.dm.pos3d.reshape(len(myDMO.dm.pos3d)*3),dtype=np.float32)
mass = np.array(myDMO.dm.mass,dtype=np.float32)
ok, acc, myDMO.dm.Phy = CF.getGravity(pos,mass,0.09,G=myDMO.p.G)
grav_c = myDMO.dm.pos3d[myDMO.dm.Phy==myDMO.dm.Phy.min()][0] 
grav_mod = np.sqrt(grav_c[0]**2+grav_c[1]**2+grav_c[2]**2)
line =  "{0:.4f}, {1:.4f}, {2:.4e}, {3:.4e}, {4:.4f}, {5:.4e}, ".format((K/U)+1,(K/np.abs(U)),K/2.,U,grav_mod,myDMO.dm.total_m)
line += simname
os.system('echo '+line+'>> Qdata.txt')