Commit 756a790c authored by NUNEZ Arturo's avatar NUNEZ Arturo
Browse files

Automatic commit mercredi 6 décembre 2017, 16:30:01 (UTC+0100)

parent 90734f5a
{
"cells": [
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%matplotlib notebook\n",
"from scipy.stats import rv_continuous\n",
"from scipy.special import gamma\n",
"import numpy as np\n",
"import emcee\n",
"from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
"from numpy import exp, sqrt\n",
"from scipy.integrate import quad, dblquad\n",
"import matplotlib.patches as patches\n",
"from scipy.integrate import quad\n",
"import scipy.optimize as optimize\n",
"import matplotlib.pyplot as plt\n",
"import matplotlib as mpl\n",
"from sklearn.neighbors import KDTree\n",
"import sys\n",
"import lmfit\n",
"from py_unsio import *\n",
"import pymc\n",
"import os\n",
"from pymodelfit import FunctionModel1DAuto\n",
"import wkbl\n",
"from mpl_toolkits.mplot3d import axes3d\n",
"from matplotlib import cm\n",
"import wkbl.astro.nbody_essentials as nbe\n",
"import cfalcon\n",
"CF =cfalcon.CFalcon()\n",
"import iminuit\n",
"from iminuit import Minuit, describe, Struct\n",
"import probfit\n",
"import warnings\n",
"warnings.filterwarnings('ignore')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [
{
"ename": "TypeError",
"evalue": "__init__() got multiple values for keyword argument 'getcen'",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-4-2b1f21952a39>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0;31m#path = \"/data/POL/HALOA/output_01274\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0;31m#path = \"/data/MANU/anunez/output_00690\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0mmyhalo\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[0mpath\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\"halo,gas,stars\"\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mgetcen\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;31mTypeError\u001b[0m: __init__() got multiple values for keyword argument 'getcen'"
]
}
],
"source": [
"path = \"/data/POL/HALOB/output_00417\"\n",
"#path = \"/data/POL/HALOA/output_01274\"\n",
"#path = \"/data/MANU/anunez/output_00690\"\n",
"myhalo = wkbl.Galaxy_Hound(path,\"halo,gas,stars\",getcen=False)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"ok,rho,_= CF.getDensity(np.array(myhalo.st.pos3d.reshape(len(myhalo.st.pos3d)*3),dtype=np.float32), myhalo.st.mass)\n",
"centro_rho = myhalo.st.pos3d[np.where(rho == rho.max())][0]\n",
"print \"density\",centro_rho\n",
"myhalo.center_shift(centro_rho)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"myhalo.r_virial(600,n=4)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"pos_dm = np.array(myhalo.dm.pos3d.reshape(len(myhalo.dm.pos3d)*3),dtype=np.float32)\n",
"pos_gs = np.array(myhalo.gs.pos3d.reshape(len(myhalo.gs.pos3d)*3),dtype=np.float32)\n",
"pos_st = np.array(myhalo.st.pos3d.reshape(len(myhalo.st.pos3d)*3),dtype=np.float32)\n",
"pos = np.concatenate((pos_dm, pos_st, pos_gs))\n",
"mass = np.concatenate((myhalo.dm.mass,myhalo.st.mass,myhalo.gs.mass))\n",
"v = np.concatenate((myhalo.dm.v,myhalo.st.v,myhalo.gs.v))\n",
"print len(mass)*3, len(pos)\n",
"pos3d = pos.reshape(len(pos)/3,3)\n",
"r2 = pos3d[:,0]**2 + pos3d[:,1]**2 +pos3d[:,2]**2"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"myhalo.gs.hsml.min()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"ok, acc, Phy = CF.getGravity(pos,mass,0.15,G=myhalo.p.G)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"#r_array = np.logspace(np.log10(min_R),2*np.log10(4*myhalo.r200),100)"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"# Potential and Rmax "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"\n",
"bin_num = 400\n",
"\n",
"pot_sph, bins_pot = np.histogram(r2,bins=bin_num,\n",
" weights=Phy)\n",
"n, _ = np.histogram(r2,bins=bin_num)\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig, ax = plt.subplots()\n",
"ax.set_title(\"spherical shells\",fontsize=30)\n",
"#ax.set_ylim([pot_sph.min(),np.abs(pot_sph.min())/6.])\n",
"ax.set_ylabel(r'$\\Phi$ [kpc $^2$ M$_{\\odot}$ s$^{-2}$ ] ', fontsize=15)\n",
"#ax.set_xscale('log')\n",
"ax.set_xlabel(r'r [kpc] ', fontsize=15)\n",
"ax.plot(np.sqrt(bins_pot[1:]),pot_sph/n,'bo-')\n",
"\n",
"ax.add_patch(\n",
" patches.Rectangle((450, -1.4e-29),\n",
" 100,\n",
" 0.9e-29,\n",
" fill=False # remove background\n",
" )\n",
")\n",
"\n",
"\n",
"left, bottom, width, height = [0.35, 0.23, 0.4, 0.4]\n",
"ax2 = fig.add_axes([left, bottom, width, height])\n",
"ax2.set_xlim([450,550])\n",
"ax2.set_ylim([-1.2e-29,-0.6e-29])\n",
"\n",
"ax2.set_ylabel(r'$\\Phi$ [kpc $^2$ M$_{\\odot}$ s$^{-2}$ ] ', fontsize=15)\n",
"ax2.set_xlabel(r'r [kpc] ', fontsize=15)\n",
"ax2.plot(np.sqrt(bins_pot[1:]),pot_sph/n,'bo-')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"pot_sph, bins_pot = np.histogram(r2,bins=np.logspace(np.log10(2),2*np.log10(4*myhalo.r200),1000),\n",
" weights=Phy)\n",
"n, _ = np.histogram(r2,bins=np.logspace(np.log10(2),2*np.log10(4*myhalo.r200),1000))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"hide_input": true
},
"outputs": [],
"source": [
"fig, ax = plt.subplots()\n",
"ax.set_title(\"spherical shells\",fontsize=30)\n",
"#ax.set_ylim([pot_sph.min(),np.abs(pot_sph.min())/6.])\n",
"ax.set_ylabel(r'$\\Phi$ [kpc $^2$ M$_{\\odot}$ s$^{-2}$ ] ', fontsize=15)\n",
"#ax.set_xscale('log')\n",
"ax.set_xlabel(r'r [kpc] ', fontsize=15)\n",
"ax.plot(np.sqrt(bins_pot[1:]),pot_sph/n,'bo-')\n",
"\n",
"ax.add_patch(\n",
" patches.Rectangle((450, -1.4e-29),\n",
" 100,\n",
" 0.9e-29,\n",
" fill=False # remove background\n",
" )\n",
")\n",
"\n",
"\n",
"left, bottom, width, height = [0.35, 0.23, 0.4, 0.4]\n",
"ax2 = fig.add_axes([left, bottom, width, height])\n",
"ax2.set_xlim([450,550])\n",
"ax2.set_ylim([-1.2e-29,-0.6e-29])\n",
"\n",
"ax2.set_ylabel(r'$\\Phi$ [kpc $^2$ M$_{\\odot}$ s$^{-2}$ ] ', fontsize=15)\n",
"ax2.set_xlabel(r'r [kpc] ', fontsize=15)\n",
"ax2.plot(np.sqrt(bins_pot[1:]),pot_sph/n,'bo-')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# l.o.s story "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"ok,rho_dm,_= CF.getDensity(np.array(myhalo.dm.pos3d.reshape(len(myhalo.dm.pos3d)*3),dtype=np.float32), myhalo.dm.mass)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# center of big neighbour\n",
"neighbourg = rho_dm[np.where((myhalo.dm.r<520)&(myhalo.dm.r>480))].max()\n",
"myneighbourg = myhalo.dm.pos3d[np.where(rho_dm==neighbourg)][0]\n",
"# its radius squared\n",
"r_nei2 = myneighbourg[0]**2 + myneighbourg[1]**2 + myneighbourg[2]**2\n",
"# distance to every point proyected to the line conecting center to neigblourg \n",
"adyacent = (pos3d[:,0]*myneighbourg[0] + pos3d[:,1]*myneighbourg[1] + pos3d[:,2]*myneighbourg[2]) / np.sqrt(r_nei2)\n",
"# distance to each point\n",
"hipoteneuse = np.sqrt(pos3d[:,0]**2 + pos3d[:,1]**2 + pos3d[:,2]**2)\n",
"# angle of cone\n",
"alpha = np.radians(5)\n",
"cos_alpha = np.cos(alpha)\n",
"# cosine of all particles respective to their angle to the l.o.s\n",
"cos_all = adyacent / hipoteneuse \n",
"# final selection of cone\n",
"my_cone = pos3d[np.where((cos_all)>cos_alpha)]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"hide_input": true
},
"outputs": [],
"source": [
"fig, ax = plt.subplots(figsize=[6,6])\n",
"ax.set_xlim([-600,0])\n",
"ax.set_ylim([-320,300])\n",
"\n",
"ax.scatter(my_cone[:,0], my_cone[:,2],c='r',lw=0,s=0.1,alpha=0.5)\n",
"ax.scatter(myneighbourg[0],myneighbourg[2],s=20)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Virial radius and mass of neighbour \n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# radii from center of neighbour\n",
"r_neig = np.sqrt((pos3d[:,0] - myneighbourg[0])**2 +(pos3d[:,1] - myneighbourg[1])**2 +(pos3d[:,2] - myneighbourg[2])**2 )\n",
"# calculating the R200 of neighbout\n",
"mhist, rhist = np.histogram(r_neig,range=(0.0,np.sqrt(r_nei2)),bins=512, weights=mass )\n",
"vol_bin = (4./3.)*np.pi*(rhist[:-1]**3)\n",
"r_bin = rhist[:-1]+ 0.5*(rhist[2]-rhist[1])\n",
"rho_s = np.cumsum(mhist) / vol_bin\n",
"r200_neigh = r_bin[np.argmin(np.abs(rho_s - (200 * myhalo.p.rho_crit)))]\n",
"# mass of neighbour\n",
"mass_neigh = np.sum(mass[np.where(r2<r200_neigh**2)])\n",
"print \"the biggest neighbourg has a mass of {0:.4e} Msun and a r200 of {1:.4f} kpc\".format(mass_neigh,r200_neigh)\n",
"print \"and it is located at {0:.4f} kpc from de center of our halo \".format(np.sqrt(r_nei2))\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"hide_input": true,
"scrolled": true
},
"outputs": [],
"source": [
"fig, ax = plt.subplots(figsize=[6,6])\n",
"ax.set_xlim([-600,0])\n",
"ax.set_ylim([-320,300])\n",
"\n",
"ax.scatter(my_cone[:,0], my_cone[:,2],c='r',lw=0,s=0.1,alpha=0.5)\n",
"ax.scatter(myneighbourg[0],myneighbourg[2],s=20)\n",
"circle = plt.Circle(((0),(0)),myhalo.r200, color='k',lw=1, fill=False)\n",
"ax.add_patch( circle )\n",
"circle = plt.Circle(((myneighbourg[0]),(myneighbourg[2])),r200_neigh, color='k',lw=1, fill=False)\n",
"ax.add_patch( circle )"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"los_condition = np.where((cos_all)>cos_alpha)\n",
"cone_pos = my_cone\n",
"cone_mass = mass[los_condition]\n",
"cone_Phy = Phy[los_condition]\n",
"cone_r2 = my_cone[:,0]**2 + my_cone[:,1]**2 + my_cone[:,2]**2\n",
"# now the histogram\n",
"bin_num = 1000\n",
"pot_los, bins_pot_los = np.histogram(cone_r2,bins=bin_num,\n",
" weights=cone_Phy)\n",
"n_los, _ = np.histogram(cone_r2,bins=bin_num)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"hide_input": true
},
"outputs": [],
"source": [
"fig, ax = plt.subplots()\n",
"ax.set_title(\"l.o.s\",fontsize=30)\n",
"#ax.set_ylim([pot_sph.min(),np.abs(pot_sph.min())/6.])\n",
"ax.set_ylabel(r'$\\Phi$ [kpc $^2$ M$_{\\odot}$ s$^{-2}$ ] ', fontsize=15)\n",
"#ax.set_xscale('log')\n",
"ax.set_xlabel(r'r [kpc] ', fontsize=15)\n",
"ax.plot(np.sqrt(bins_pot_los[1:]),pot_los/n_los,'bo-')\n",
"\n",
"ax.add_patch(\n",
" patches.Rectangle((450, -1.4e-29),\n",
" 100,\n",
" 0.9e-29,\n",
" fill=False # remove background\n",
" )\n",
")\n",
"\n",
"\n",
"left, bottom, width, height = [0.35, 0.23, 0.4, 0.4]\n",
"ax2 = fig.add_axes([left, bottom, width, height])\n",
"ax2.set_xlim([480,545])\n",
"ax2.set_ylim([-1e-29,-0.8e-29])\n",
"\n",
"ax2.set_ylabel(r'$\\Phi$ [kpc $^2$ M$_{\\odot}$ s$^{-2}$ ] ', fontsize=15)\n",
"ax2.set_xlabel(r'r [kpc] ', fontsize=15)\n",
"ax2.plot(np.sqrt(bins_pot_los[1:]),pot_los/n_los,'bo-')\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# virial ratio q\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"inside_halo = np.where(r2<(myhalo.r200)**2)\n",
"kmtokpc = 1 / 3.08567758128e+16 \n",
"q = (np.sum(mass[inside_halo]*(v[inside_halo]*kmtokpc)**2) / np.sum(mass[inside_halo]*Phy[inside_halo])) + 1 "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"print \"q = \",q,r\", numerator = \",np.sum(mass[inside_halo]*(v[inside_halo]*kmtokpc)**2),\", denominator = \",np.sum(mass[inside_halo]*Phy[inside_halo])\n",
"print \"ratio = \",(np.sum(mass[inside_halo]*(v[inside_halo]*kmtokpc)**2) / np.sum(mass[inside_halo]*Phy[inside_halo]))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"dm = len(myhalo.dm.mass)+len(myhalo.st.mass)\n",
"inside_halo = np.where(r2[:dm]<(myhalo.r200)**2)\n",
"kmtokpc = 1 / 3.08567758128e+16 \n",
"q = (np.sum(mass[inside_halo]*(v[inside_halo]*kmtokpc)**2) / np.sum(mass[inside_halo]*Phy[inside_halo])) + 1 "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"print \"q = \",q,r\", numerator = \",np.sum(mass[inside_halo]*(v[inside_halo]*kmtokpc)**2),\", denominator = \",np.sum(mass[inside_halo]*Phy[inside_halo])\n",
"print \"ratio = \",(np.sum(mass[inside_halo]*(v[inside_halo]*kmtokpc)**2) / np.sum(mass[inside_halo]*Phy[inside_halo]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# $\\beta$ parameter\n",
"\n",
"$$\\beta = 1 - \\frac{ \\sigma_{\\phi}^2 + \\sigma_{\\theta}^2}{2 \\sigma_{r}^2}$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"point_num = 300\n",
"r_beta = np.logspace(-1,np.log10(4*myhalo.r200),point_num)\n",
"vphi = np.concatenate((myhalo.dm.vphi,myhalo.st.vphi,myhalo.gs.vphi))\n",
"vtheta = np.concatenate((myhalo.dm.vtheta,myhalo.st.vtheta,myhalo.gs.vtheta))\n",
"vr = np.concatenate((myhalo.dm.vr,myhalo.st.vr,myhalo.gs.vr))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def beta_param(i):\n",
" condition = np.where((r2>r_beta[i]**2)&(r2<=r_beta[i+1]**2))\n",
" v_r = vr[condition]\n",
" v_phi = vphi[condition]\n",
" v_theta = vtheta[condition]\n",
" #print (np.std(v_phi))**2 ,(np.std(v_theta))**2 , (np.std(v_r))**2 \n",
" return 1 - ((np.std(v_phi))**2 +(np.std(v_theta))**2) / 2. / (np.std(v_r))**2 \n",
"get_beta = np.vectorize(beta_param)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"beta_r = get_beta(range(point_num-1))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"print len(r_beta[:-1]), len(beta_r)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig, ax = plt.subplots()\n",
"ax.set_xscale('log')\n",
"ax.set_xlabel(\"r [kpc]\",fontsize=15)\n",
"ax.set_ylabel(r\"$\\beta(r)$\",fontsize=18)\n",
"ax.plot((r_beta[1:]+r_beta[:-1])/2/10**0.8,beta_r)\n",
"#ax.axvline(x=myhalo.r200, color='k',linestyle='--',label=r'r$_{200}$')\n",
"#ax.axvline(x=myhalo.r97, color='gray',linestyle='--',label=r'r$_{97}$')\n",
"#ax.axvline(x=np.sqrt(r_nei2), color='g',linestyle='--',label=r'first neighbourg')\n",
"#legend = ax.legend(loc='lower right', ncol=1, shadow=False, fontsize=14)\n",
"#frame = legend.get_frame()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.12"
},
"latex_envs": {
"bibliofile": "biblio.bib",
"cite_by": "apalike",
"current_citInitial": 1,
"eqLabelWithNumbers": true,
"eqNumInitial": 0
}
},
"nbformat": 4,
"nbformat_minor": 1
}
This diff is collapsed.
This diff is collapsed.
{
"cells": [],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 1
}