Skip to content
GitLab
Commit dcffaa00 authored by NUNEZ Arturo's avatar NUNEZ Arturo
Browse files

Automatic commit vendredi 27 juillet 2018, 17:46:07 (UTC+0200)

parent 3e3de1e6
Hide whitespace changes
Inline Side-by-side
LAM2LUPM/S_T_beta_Mochima.ipynb
View file @ dcffaa00
%% Cell type:code id: tags:
``` python
%matplotlib notebook
%load_ext autoreload
%autoreload 2
```
%% Cell type:code id: tags:
``` python
from scipy.stats import rv_continuous
from scipy.interpolate import interp1d
from matplotlib.patches import Circle
from scipy.special import gamma
import numpy as np
import emcee
from mpl_toolkits.axes_grid1 import make_axes_locatable
from numpy import exp, sqrt
from scipy.integrate import quad, dblquad, nquad
import matplotlib.patches as patches
from itertools import product
from scipy.integrate import quad
import scipy.optimize as optimize
import matplotlib.pyplot as plt
import matplotlib as mpl
from sklearn.neighbors import KDTree
import sys
import lmfit
from py_unsio import *
import pymc
import os
from scipy.integrate import simps
from pymodelfit import FunctionModel1DAuto
import wkbl
from wkbl.astro.halo_info import *
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm
import wkbl.astro.nbody_essentials as nbe
import cfalcon
CF =cfalcon.CFalcon()
import iminuit
from iminuit import Minuit, describe, Struct
import probfit
from matplotlib.colors import LogNorm
from matplotlib.ticker import FormatStrFormatter
import warnings
warnings.filterwarnings('ignore')
```
%% Cell type:code id: tags:
``` python
hydro = wkbl.astro.halo_info.MochimaHydro()
```
%% Cell type:code id: tags:
``` python
#halo = HALOBHydro(where="home")
simname="Mochima"
pathsim = "/data/OWN/SF1test/SF0/mstar1_T3600/output_00041"
myhydro = wkbl.Galaxy_Hound(pathsim)
print myhydro.dm.pos3d[:,0].max()
zoom_reg = np.where(myhydro.dm.mass==myhydro.dm.mass.min())
nucenter = nbe.real_center(myhydro.dm.pos3d[zoom_reg], myhydro.dm.mass[zoom_reg])
print nucenter
myhydro.center_shift(nucenter)
myhydro.r_virial(600,n=3)
```
%% Output
loading Dark matter..
loading Stars..
loading Gas..
36844.594
[20418.88714131 17567.72033332 17124.40448217]
| r_200 = 212.70
| Diagonal matrix computed
| | 20, 0, 0|
| D =| 0, 14, 0|
| | 0, 0, 4|
%% Cell type:code id: tags:
``` python
simname_nospace = list(simname)
for i in range(len(simname)):
if simname[i]==" ":
simname_nospace[i]="_"
simname_nospace = "".join(simname_nospace)
```
%% Cell type:code id: tags:
``` python
pos_dm = np.array(myhydro.dm.pos3d.reshape(len(myhydro.dm.pos3d)*3),dtype=np.float32)
pos_gs = np.array(myhydro.gs.pos3d.reshape(len(myhydro.gs.pos3d)*3),dtype=np.float32)
pos_st = np.array(myhydro.st.pos3d.reshape(len(myhydro.st.pos3d)*3),dtype=np.float32)
pos = np.concatenate((pos_dm, pos_st, pos_gs))
phi_cord = np.concatenate((myhydro.dm.phi,myhydro.st.phi, myhydro.gs.phi))
mass = np.concatenate((myhydro.dm.mass,myhydro.st.mass,myhydro.gs.mass))
v = np.concatenate((myhydro.dm.v,myhydro.st.v,myhydro.gs.v))
print len(mass)*3, len(pos)
pos3d = pos.reshape(len(pos)/3,3)
r2 = pos3d[:,0]**2 + pos3d[:,1]**2 +pos3d[:,2]**2
```
%% Output
10194240 10194240
%% Cell type:markdown id: tags:
# Ellipticity T and S
%% Cell type:code id: tags:
``` python
get_mat = np.vectorize(nbe.m_matrix_for_r)
```
%% Cell type:code id: tags:
``` python
r = np.linspace(myhydro.gs.hsml.min(),2*myhydro.r200,300)
M_dm = np.array([nbe.m_matrix_for_r(myhydro,'halo',i) for i in r])
M_st = np.array([nbe.m_matrix_for_r(myhydro,'stars',i) for i in r])
a_dm, b_dm, c_dm = np.sqrt(M_dm[:,0,0]), np.sqrt(M_dm[:,1,1]), np.sqrt(M_dm[:,2,2])
a_st, b_st, c_st = np.sqrt(M_st[:,0,0]), np.sqrt(M_st[:,1,1]), np.sqrt(M_st[:,2,2])
S_dm = c_dm/a_dm
T_dm = ((a_dm**2) - (b_dm**2))/((a_dm**2) -(c_dm**2))
S_st = c_st/a_st
T_st = ((a_st**2) - (b_st**2))/((a_st**2) -(c_st**2))
```
%% Cell type:code id: tags:
``` python
outputing = open("../../datafiles/"+simname_nospace+"_S_and_T.txt","w")
outputing.write("# "+simname+" ellipticity parameters\n")
outputing.write("# r200 = {0:.2f} kpc\n".format(myhydro.r200))
outputing.write("# format:\n")
outputing.write("# r (kpc), S , T\n")
for i in range(len(T_dm)):
outputing.write("{0:.2f} {1:.6f} {2:.6f} \n".format(r[i],S_dm[i],T_dm[i]))
outputing.close()
```
%% Cell type:code id: tags:
``` python
fig, ax = plt.subplots(figsize=[10,4])
#ax.set_xlim(4,250)
font=15
ax.set_ylim(0,1.18)
ax.set_xlim(0,410)
ax.set_xlabel('r [kpc]',fontsize=font)
ax.plot([4,5],[1e3,2e3],color='gray',linestyle='--',lw=2,label="S")
ax.plot([4,5],[1e3,2e3],color='gray',linestyle='-',lw=2,label="T")
ax.plot(r,T_dm,'-',color='#155c9e',lw=2,label="Dark matter")
ax.plot(r,S_dm,'--',color='#155c9e',lw=2)
#ax.plot(r,T_st,'-',color='#a91e4f',lw=2,label='stars')
#ax.plot(r,np.sqrt(S_st),'--',color='#a91e4f',lw=2)
ax.axvline(x=myhydro.r200,color='k',linestyle='-.')
ax.text(myhydro.r200+2,0.2,r"R$_{200}$",fontsize=font)
ax.axvline(x=myhydro.r97,color='k',linestyle='-.')
ax.text(myhydro.r97+2,0.2,r"R$_{97}$",fontsize=font)
ax.axvline(x=8,color='k',linestyle='-.')
ax.text(10,0.2,r"R$_{\odot}$",fontsize=font)
ax.text(50,0.9,r"$\rm "+simname+"$",fontsize=1.5*font)
legend = ax.legend(loc='upper right', ncol=2, shadow=False, fontsize=font)
frame = legend.get_frame()
ax.tick_params(axis='both', which='major', labelsize=15, size=5,width=1.2)
ax.tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
fig.tight_layout()
```
%% Output
%% Cell type:markdown id: tags:
# Beta of r
%% Cell type:code id: tags:
``` python
point_num = 150
r_beta = np.logspace(-1,np.log10(4*myhydro.r200),point_num)
vphi = np.concatenate((myhydro.dm.vphi,myhydro.st.vphi,myhydro.gs.vphi))
vtheta = np.concatenate((myhydro.dm.vtheta,myhydro.st.vtheta,myhydro.gs.vtheta))
vr = np.concatenate((myhydro.dm.vr,myhydro.st.vr,myhydro.gs.vr))
```
%% Cell type:code id: tags:
``` python
def beta_param(i):
condition = np.where((r2>r_beta[i]**2)&(r2<=r_beta[i+1]**2))
v_r = vr[condition]
v_phi = vphi[condition]
v_theta = vtheta[condition]
#print (np.std(v_phi))**2 ,(np.std(v_theta))**2 , (np.std(v_r))**2
return 1 - ((np.std(v_phi))**2 +(np.std(v_theta))**2) / 2. / (np.std(v_r))**2
get_beta = np.vectorize(beta_param)
```
%% Cell type:code id: tags:
``` python
beta_r = get_beta(range(point_num-1))
```
%% Cell type:code id: tags:
``` python
outputing = open("../../datafiles/"+simname_nospace+"_Beta.txt","w")
outputing.write("# "+simname+" Beta parameters\n")
outputing.write("# r200 = {0:.2f} kpc\n".format(myhydro.r200))
outputing.write("# format:\n")
outputing.write("# r (kpc), Beta(r) \n")
for i in range(len(beta_r)):
outputing.write("{0:.2f} {1:.6f} \n".format(((r_beta[1:]+r_beta[:-1])/2.)[i],beta_r[i]))
outputing.close()
```
%% Cell type:code id: tags:
``` python
fig, ax = plt.subplots(figsize=[9,5])
r_2 = 6.95
ax.set_xscale('log')
#ax.set_title(r"$\beta(r)$ parameter", fontsize=1.5*font)
ax.set_xlabel(r"r [kpc]",fontsize=font)
ax.set_ylabel(r"$\beta(r)$",fontsize=1.5*font)
ax.set_ylim([-1.5,1])
ax.plot((r_beta[1:]+r_beta[:-1])/2.,beta_r,lw=1.5)
#ax.plot((r_beta[1:]+r_beta[:-1])/2/r_2,beta_r_2,lw=1.5)
#ax.plot((r_beta[1:]+r_beta[:-1])/2/r_2,beta_r_3,lw=1.5)
ax.axvline(x=myhydro.r200, color='k',lw=2,linestyle='--')
ax.axvline(x=myhydro.r97, color='gray',lw=2,linestyle='--')
ax.axvline(x=8/r_2, color='y',linestyle='--',lw=2,label=r'r$_{\odot}$')
fig.tight_layout()
ax.tick_params(axis='both', which='major', labelsize=15, size=5,width=1.2)
ax.tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
#plt.savefig("/home/arturo/Documents/git/LAMtoLUPM_latex/beta_r.png",dpi=300)
```
%% Output
%% Cell type:markdown id: tags:
# F(v) moments
%% Cell type:code id: tags:
``` python
def moments(rmin,rmax):
selection= np.where((myhydro.dm.r>rmin)&(myhydro.dm.r<=rmax))
fdv, vs = np.histogram(myhydro.dm.v[selection],bins=50,normed=1)
vcenter = (vs[:-1]+vs[1:])/2.
m1 = simps(vcenter*fdv,x=vcenter)
m2 = simps((vcenter**2)*fdv,x=vcenter)
m3 = simps((vcenter**3)*fdv,x=vcenter)
return m1,m2,m3
def moments2(rmin,rmax):
selection= np.where((myhydro.dm.r>rmin)&(myhydro.dm.r<=rmax))
fdv, vs = np.histogram(myhydro.dm.v[selection],bins=50,normed=1)
vcenter = (vs[:-1]+vs[1:])/2.
vw = vs[1:]-vs[:-1]
m1 = np.sum(vcenter*fdv*vw)
m2 = np.sum((vcenter**2)*fdv*vw)
m3 = np.sum((vcenter**3)*fdv*vw)
return m1,m2,m3
def moments3(rmin,rmax):
selection= np.where((myhydro.dm.r>rmin)&(myhydro.dm.r<=rmax))
npart = len(myhydro.dm.v[selection])
m_2 = np.sum(myhydro.dm.v[selection]**(-2))/npart
m_1 = np.sum(myhydro.dm.v[selection]**(-1))/npart
m1 = np.sum(myhydro.dm.v[selection])/npart
m2 = np.sum((myhydro.dm.v[selection])**2)/npart
m3 = np.sum((myhydro.dm.v[selection])**3)/npart
return m1,m2,m3
return m_2,m_1,m1, m2, np.sqrt(npart)
fdv_moments = np.vectorize(moments)
fdv_moments2 = np.vectorize(moments2)
fdv_moments3 = np.vectorize(moments3)
```
%% Cell type:code id: tags:
``` python
patho = "/home/arturo/Documents/LAM/LAM2LUPM/speed/"+hydro.namenospace
v_2_av = np.loadtxt(patho+"/Eddington/v_minus2_av_Eddington_Mochima2_DM_baryons_Rmax=2908.43kpc_no_divergence.txt")
v_1_av = np.loadtxt(patho+"/Eddington/v_minus1_av_Eddington_Mochima2_DM_baryons_Rmax=2908.43kpc_no_divergence.txt")
v_av = np.loadtxt(patho+"/Eddington/v_av_Eddington_Mochima2_DM_baryons_Rmax=2908.43kpc_no_divergence.txt")
v2_av = np.loadtxt(patho+"/Eddington/v_sq_av_Eddington_Mochima2_DM_baryons_Rmax=2908.43kpc_no_divergence.txt")
v_2_av_m = np.loadtxt("/home/arturo/Documents/LAM/LAM2LUPM/speed/Mochima2/Maxwellian/v_minus2_av_Mochima2_DM_baryons_Rmax=2908.43kpc_Maxwellian.txt")
v_1_av_m = np.loadtxt("/home/arturo/Documents/LAM/LAM2LUPM/speed/Mochima2/Maxwellian/v_minus1_av_Mochima2_DM_baryons_Rmax=2908.43kpc_Maxwellian.txt")
v_av_m = np.loadtxt("/home/arturo/Documents/LAM/LAM2LUPM/speed/Mochima2/Maxwellian/v_av_Mochima2_DM_baryons_Rmax=2908.43kpc_Maxwellian.txt")
v2_av_m = np.loadtxt("/home/arturo/Documents/LAM/LAM2LUPM/speed/Mochima2/Maxwellian/v_sq_av_Mochima2_DM_baryons_Rmax=2908.43kpc_Maxwellian.txt")
```
%% Cell type:code id: tags:
``` python
r_v = np.logspace(np.log10(4*myhydro.gs.hsml.min()),np.log10(3*myhydro.r200),150)
#m1,m2,m3 = fdv_moments(r_v[:-1],r_v[1:])
#mm1,mm2,mm3 = fdv_moments2(r_v[:-1],r_v[1:])
mmm1,mmm2,mmm3 = fdv_moments3(r_v[:-1],r_v[1:])
m_2, m_1,m1, m2, sigma = fdv_moments3(r_v[:-1],r_v[1:])
```
%% Cell type:code id: tags:
``` python
fig, ax = plt.subplots(3,1)
fig, ax = plt.subplots(4,1,figsize=[8,12])
ax[0].yaxis.get_major_formatter().set_powerlimits((0, 2))
ax[1].yaxis.get_major_formatter().set_powerlimits((0, 2))
ax[2].yaxis.get_major_formatter().set_powerlimits((0, 2))
ax[3].yaxis.get_major_formatter().set_powerlimits((0, 2))
ax[3].set_xlabel(r'$\rm [kpc]$',fontsize=18)
ax[0].set_ylabel(r'$\langle v^{-2} \rangle$$\rm \, [km/s]^{-2}$',fontsize=18)
ax[1].set_ylabel(r'$\langle v^{-1} \rangle$$\rm \,[km/s]^{-1}$',fontsize=18)
ax[2].set_ylabel(r'$\langle v \rangle$$\rm \,[km/s]$',fontsize=18)
ax[3].set_ylabel(r'$\langle v^{2} \rangle$$\rm \,[km/s]^2$',fontsize=18)
ax[2].set_xlabel('r [kpc]',fontsize=18)
ax[0].set_ylabel('[km/s]',fontsize=18)
ax[1].set_ylabel(r'[km/s]$^2$',fontsize=18)
ax[2].set_ylabel(r'[km/s]$^3$',fontsize=18)
ax[0].set_xscale('log')
ax[1].set_xscale('log')
ax[2].set_xscale('log')
ax[0].set_ylim([52,1.1*m1.max()])
ax[1].set_ylim([1e4,1.1*m2.max()])
ax[2].set_ylim([1e4,1.1*m3.max()])
ax[3].set_xscale('log')
ax[0].set_ylim([m_2.min(),1.1*m_2.max()])
ax[1].set_ylim([m_1.min(),1.1*m_1.max()])
ax[2].set_ylim([m1.min(),1.1*m1.max()])
ax[3].set_ylim([m2.min(),1.1*m2.max()])
ax[0].set_xlim([0.5,3*myhydro.r200])
ax[1].set_xlim([0.5,3*myhydro.r200])
ax[2].set_xlim([0.5,3*myhydro.r200])
ax[3].set_xlim([0.5,3*myhydro.r200])
ax[0].axvline(x=myhydro.r200,color='k',linestyle='--')
ax[1].axvline(x=myhydro.r200,color='k',linestyle='--')
ax[2].axvline(x=myhydro.r200,color='k',linestyle='--')
ax[3].axvline(x=myhydro.r200,color='k',linestyle='--')
r_v_c = (r_v[:-1]+r_v[1:])/2.
ax[0].plot(r_v_c,m1)
ax[0].plot(r_v_c,mm1)
ax[0].plot(r_v_c,mmm1)
ax[1].plot(r_v_c,m2)
ax[1].plot(r_v_c,mm2)
ax[1].plot(r_v_c,mmm2)
ax[2].plot(r_v_c,m3)
ax[2].plot(r_v_c,mm3)
ax[2].plot(r_v_c,mmm3)
ax[0].plot(r_v_c,m_2,lw=2)
ax[0].plot(v_2_av[:,0],v_2_av[:,1],"k",lw=1.5)
ax[0].plot(v_2_av_m[:,0],v_2_av_m[:,1],"r",lw=1.5)
fig.text(0.3,.9,r"$\rm"+ hydro.name +" $",fontsize=30)
fig.text(0.35,.87,r"$\rm Hydro $",fontsize=25)
ax[1].plot(r_v_c,m_1,lw=2)
ax[1].plot(v_1_av[:,0],v_1_av[:,1],"k",lw=1.5)
ax[1].plot(v_1_av_m[:,0],v_1_av_m[:,1],"r",lw=1.5)
ax[2].plot(r_v_c,m1,lw=2)
ax[2].plot(v_av[:,0],v_av[:,1],"k",lw=1.5)
ax[2].plot(v_av_m[:,0],v_av_m[:,1],"r",lw=1.5)
ax[3].plot(r_v_c,m2,lw=2,label=r"$\rm Data$")
ax[3].plot(v2_av[:,0],v2_av[:,1],"k",lw=1.5,label=r"$\rm Eddington$")
ax[3].plot(v2_av_m[:,0],v2_av_m[:,1],"r",lw=1.5,label=r"$\rm Maxwellian$")
legend = ax[3].legend(loc='upper right', ncol=2, shadow=False, fontsize=font)
frame = legend.get_frame()
fig.tight_layout(h_pad=-2.7)
ax[0].text(22,250,r'$\int_0 ^{v_{max}} v f(v) dv$',fontsize=17)
ax[1].text(22,70000,r'$\int_0 ^{v_{max}} v^2 f(v) dv$',fontsize=17)
ax[2].text(22,3e7,r'$\int_0 ^{v_{max}} v^3 f(v) dv$',fontsize=17)
ax[2].text(170,1e7,r'R$_{200}$',fontsize=17)
ax[0].tick_params(axis='y', which='major', labelsize=15, size=5,width=1.2)
ax[0].tick_params(axis='x', which='major', labelsize=0, size=5,width=1.2)
ax[0].tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
ax[1].tick_params(axis='y', which='major', labelsize=15, size=5,width=1.2)
ax[1].tick_params(axis='x', which='major', labelsize=0, size=5,width=1.2)
ax[1].tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
ax[2].tick_params(axis='both', which='major', labelsize=15, size=5,width=1.2)
ax[2].tick_params(axis='y', which='major', labelsize=15, size=5,width=1.2)
ax[2].tick_params(axis='x', which='major', labelsize=0, size=5,width=1.2)
ax[2].tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
ax[3].tick_params(axis='both', which='major', labelsize=15, size=5,width=1.2)
ax[3].tick_params(axis='both', which='minor', labelsize=15, size=3,width=1.2)
plt.savefig("/home/arturo/Documents/git/LAMtoLUPM_latex/moments.png",dpi=300)
```
%% Output
%% Cell type:code id: tags:
``` python
outputing = open("../../datafiles/"+simname_nospace+"_Moments.txt","w")
outputing.write("# "+simname+" Statistical moments of dm velocity\n")
outputing.write("# r200 = {0:.2f} kpc\n".format(myhydro.r200))
outputing.write("# format:\n")
outputing.write("# r (kpc), 1stmoment, 2ndmoment, 3rdmoment\n")
for i in range(len(r_v_c)):
outputing.write("{0:.2f} {1:.6e} {2:.6e} {3:.6e}\n".format(r_v_c[i],m1[i],m2[i],m3[i]))
outputing.close()
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
......
Simulations/LUMINOSITyinge.ipynb
View file @ dcffaa00
This source diff could not be displayed because it is too large. You can view the blob instead.
Simulations/SFR_SIMUSSS.ipynb
View file @ dcffaa00