Commit da823473 authored by NUNEZ Arturo's avatar NUNEZ Arturo

Automatic commit vendredi 20 juillet 2018, 16:30:06 (UTC+0200)

parent c10161c5
{
"cells": [],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 1
}
......@@ -64,7 +64,7 @@
},
{
"cell_type": "code",
"execution_count": 3,
"execution_count": 4,
"metadata": {
"collapsed": false,
"hide_input": false,
......@@ -77,24 +77,28 @@
"text": [
"loading Dark matter..\n",
"loading Stars..\n",
"loading Gas..\n",
"36839.426\n",
"[14314.38720356 15227.74783398 15695.98061772]\n",
"| r_200 = 211.52\n",
"| Diagonal matrix computed \n",
"| | 19, 0, 0|\n",
"| D =| 0, 15, 0|\n",
"| | 0, 0, 4|\n"
"loading Gas..\n"
]
},
{
"ename": "LinAlgError",
"evalue": "Array must not contain infs or NaNs",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mLinAlgError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-4-806dbd19e21c>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0mmyhydro\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mwkbl\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mGalaxy_Hound\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mhydro\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpath\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0mmyhydro\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcenter_shift\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mhydro\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mc_dm_com\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0mmyhydro\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mr_virial\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m600\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m25\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;32m/home/arturo/Documents/git/WKBL/wkbl/astro/galaxy_peeker.pyc\u001b[0m in \u001b[0;36mr_virial\u001b[0;34m(self, r_max, r_min, rotate, n, bins, quiet)\u001b[0m\n\u001b[1;32m 82\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 83\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mrotate\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;32mand\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_sts\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 84\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrotate_galaxy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 85\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mredefine\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 86\u001b[0m \u001b[0mD\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmatrix_T\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmatrix_P\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtranspose\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmatrix_T\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/home/arturo/Documents/git/WKBL/wkbl/astro/galaxy_peeker.pyc\u001b[0m in \u001b[0;36mrotate_galaxy\u001b[0;34m(self, rmin, rmax)\u001b[0m\n\u001b[1;32m 120\u001b[0m \u001b[0msecond\u001b[0m \u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0maverage\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpos_ring\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0maverage\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpos_ring\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mj\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 121\u001b[0m \u001b[0mP\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mj\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfirst\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0msecond\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 122\u001b[0;31m \u001b[0meigen_values\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mevecs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meig\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mP\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 123\u001b[0m \u001b[0morder\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margsort\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mabs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0meigen_values\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 124\u001b[0m \u001b[0mT\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/usr/local/lib/python2.7/dist-packages/numpy/linalg/linalg.pyc\u001b[0m in \u001b[0;36meig\u001b[0;34m(a)\u001b[0m\n\u001b[1;32m 1141\u001b[0m \u001b[0m_assertRankAtLeast2\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1142\u001b[0m \u001b[0m_assertNdSquareness\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1143\u001b[0;31m \u001b[0m_assertFinite\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 1144\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mresult_t\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_commonType\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1145\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/usr/local/lib/python2.7/dist-packages/numpy/linalg/linalg.pyc\u001b[0m in \u001b[0;36m_assertFinite\u001b[0;34m(*arrays)\u001b[0m\n\u001b[1;32m 214\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0ma\u001b[0m \u001b[0;32min\u001b[0m \u001b[0marrays\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 215\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0misfinite\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mall\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 216\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mLinAlgError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Array must not contain infs or NaNs\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 217\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 218\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_isEmpty2d\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mLinAlgError\u001b[0m: Array must not contain infs or NaNs"
]
}
],
"source": [
"#halo = HALOBHydro(where=\"home\") \n",
"simname=\"Adicora\"\n",
"pathsim = \"/data/OWN/Adicora/SF0/Stable/output_00041\"\n",
"myhydro = wkbl.Galaxy_Hound(pathsim)\n",
"print myhydro.dm.pos3d[:,0].max()\n",
"zoom_reg = np.where(myhydro.dm.mass==myhydro.dm.mass.min())\n",
"hydro = wkbl.astro.halo_info.AdicoraHydro() \n",
"myhydro = wkbl.Galaxy_Hound(hydro.path)\n",
"nucenter = nbe.real_center(myhydro.dm.pos3d[zoom_reg], myhydro.dm.mass[zoom_reg])\n",
"print nucenter\n",
"myhydro.center_shift(nucenter)\n",
......@@ -109,6 +113,7 @@
},
"outputs": [],
"source": [
"simname = hydro.name\n",
"simname_nospace = list(simname)\n",
"for i in range(len(simname)):\n",
" if simname[i]==\" \":\n",
......@@ -2,21 +2,12 @@
"cells": [
{
"cell_type": "code",
"execution_count": 12,
"execution_count": 1,
"metadata": {
"collapsed": false,
"hide_input": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The autoreload extension is already loaded. To reload it, use:\n",
" %reload_ext autoreload\n"
]
}
],
"outputs": [],
"source": [
"%matplotlib notebook\n",
"%load_ext autoreload\n",
......@@ -73,7 +64,7 @@
},
{
"cell_type": "code",
"execution_count": 3,
"execution_count": null,
"metadata": {
"collapsed": false,
"hide_input": false,
......@@ -86,24 +77,17 @@
"text": [
"loading Dark matter..\n",
"loading Stars..\n",
"loading Gas..\n",
"density [9667.492 9865.035 9799.733]\n",
"| r_200 = 177.54\n",
"| Diagonal matrix computed \n",
"| | 19, 0, 0|\n",
"| D =| 0, 14, 0|\n",
"| | 0, 0, 2|\n"
"loading Gas..\n"
]
}
],
"source": [
"#halo = HALOBHydro(where=\"home\") \n",
"pathsim = \"/data/O\"\n",
"myhydro = wkbl.Galaxy_Hound(pathsim)\n",
"ok,rho,_= CF.getDensity(np.array(myhydro.st.pos3d.reshape(len(myhydro.st.pos3d)*3),dtype=np.float32), myhydro.st.mass)\n",
"centro_rho = myhydro.st.pos3d[np.where(rho == rho.max())][0]\n",
"print \"density\",centro_rho\n",
"myhydro.center_shift(centro_rho)\n",
"hydro = wkbl.astro.halo_info.HALOBHydro() \n",
"\n",
"\n",
"myhydro = wkbl.Galaxy_Hound(hydro.path)\n",
"\n",
"myhydro.center_shift(hydro.c_dm_com)\n",
"myhydro.r_virial(600,n=25)\n"
]
},
......@@ -2,7 +2,7 @@
"cells": [
{
"cell_type": "code",
"execution_count": null,
"execution_count": 1,
"metadata": {
"collapsed": false,
"hide_input": false
......@@ -64,7 +64,7 @@
},
{
"cell_type": "code",
"execution_count": 3,
"execution_count": null,
"metadata": {
"collapsed": false,
"hide_input": false,
......@@ -78,12 +78,7 @@
"loading Dark matter..\n",
"loading Stars..\n",
"loading Gas..\n",
"density [9868.701 9744.951 9766.655]\n",
"| r_200 = 182.23\n",
"| Diagonal matrix computed \n",
"| | 18, 0, 0|\n",
"| D =| 0, 13, 0|\n",
"| | 0, 0, 3|\n"
"density [9868.701 9744.951 9766.655]\n"
]
}
],
......@@ -100,20 +95,12 @@
},
{
"cell_type": "code",
"execution_count": 4,
"execution_count": null,
"metadata": {
"collapsed": false,
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"5.497e+11\n"
]
}
],
"outputs": [],
"source": [
"print \"{0:.3e}\".format(np.sum(myhydro.dm.mass[myhydro.dm.r<myhydro.r200]))"
]
......@@ -146,20 +133,12 @@
},
{
"cell_type": "code",
"execution_count": 7,
"execution_count": null,
"metadata": {
"collapsed": false,
"hide_input": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"168833547 168833547\n"
]
}
],
"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",
......@@ -185,7 +164,7 @@
},
{
"cell_type": "code",
"execution_count": 8,
"execution_count": null,
"metadata": {
"collapsed": false,
"hide_input": false
This source diff could not be displayed because it is too large. You can view the blob instead.
{
"cells": [],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 1
}
This diff is collapsed.
This diff is collapsed.
This source diff could not be displayed because it is too large. You can view the blob instead.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
-0.1696, 1.1696, 1.2471e+16, -2.1327e+16, 0.0720, 9.1355e+11, H
-0.1565, 1.1565, 1.2496e+16, -2.1610e+16, 0.0668, 9.1678e+11, H
-0.1565, 1.1565, 1.2496e+16, -2.1610e+16, 0.0668, 9.1678e+11, H
-0.1696, 1.1696, 1.2471e+16, -2.1327e+16, 0.0720, 9.1355e+11, H
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())
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')
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