Source code for lenstools.simulations.amiga

from __future__ import division
import re

from .nbody import NbodySnapshot

import numpy as np

from astropy.cosmology import w0waCDM
import astropy.constants as cnst
import astropy.units as u

###############################################
#########AmigaHalos class######################
###############################################

[docs]class AmigaHalos(NbodySnapshot): """ A class that handles the Amiga Halo Finder (AHF) halo output files. Inherits from the abstract NbodySnapshot :warning: Not tested yet! """ ############################################################################################### #########################Abstract methods implementation####################################### ############################################################################################### @classmethod def buildFilename(cls,root,pool,**kwargs): #Make sure this is an AHF halo catalog assert root.endswith(".AHF_halos") #Substitute xxxx with the task number if pool is not None: return re.sub(r"\.[0-9]{4}\.",".{0:04d}.".format(pool.rank),root) else: return root @classmethod def int2root(cls,name,n): return re.sub(r"\.[0-9]+\.",".{0:04d}.".format(n),name) def getHeader(self,h=0.72,w0=-1.,wa=0.): #Create header dictionary header = dict() header["h"] = h header["w0"] = w0 header["wa"] = wa #Read the redshift from the filename header["redshift"] = float(re.search(r"\.z([0-9\.]+)\.",self.fp.name).groups()[0]) header["scale_factor"] = 1./(1.+header["redshift"]) #Look for the log file that contains all the information we need in the header task_number = int(re.search(r"\.([0-9]{4})\.",self.fp.name).groups()[0]) logfile = re.sub(r"\.[0-9]{4}\..*",".{0:02d}.log".format(task_number),self.fp.name) #Read the log file and fill the header with open(logfile,"r") as logfp: ahflog = logfp.read() #Cosmological parameters and box size Om0,Ode0,box_size = re.search(r"simu\.omega0\s*:\s*([0-9\.]+)\nsimu\.lambda0\s*:\s*([0-9\.]+)\nsimu\.boxsize\s*:\s*([0-9\.]+)",ahflog).groups() header["Om0"] = float(Om0) header["Ode0"] = float(Ode0) header["box_size"] = float(box_size)*1.0e3 #These are not important header["masses"] = np.zeros(6) header["num_particles_file"] = 1. header["num_particles_total"] = 1. header["num_files"] = 1 #Finally compute the comoving distance header["comoving_distance"] = w0waCDM(H0=h*100,Om0=header["Om0"],Ode0=header["Ode0"],w0=header["w0"],wa=header["wa"]).comoving_distance(header["redshift"]).to(u.kpc).value / h #Return to user return header def getPositions(self,first=None,last=None,save=True): #Matter density today rhoM = self.cosmology.critical_density0 * self.cosmology.Om0 if first is None: first = 0 if last is None: m,x,y,z,rv,c = np.loadtxt(self.fp,usecols=(3,5,6,7,11,42))[first:].T else: m,x,y,z,rv,c = np.loadtxt(self.fp,usecols=(3,5,6,7,11,42))[first:last].T positions = np.array((x,y,z)).astype(np.float32).T * self.kpc_over_h self.virial_radius = rv * self.kpc_over_h self.concentration = c self.weights = ((1./(4*np.pi)) * (c**3/(np.log(1.+c)-c/(1.+c))) * m*(u.Msun/self.header["h"]) / (rhoM*(self.virial_radius**3))).decompose().value if save: self.positions = positions return self.positions return positions ############################################ ###########These are not necessary########## ############################################ def getVelocities(self,first=None,last=None,save=True): raise NotImplementedError def getID(self,first=None,last=None,save=True): raise NotImplementedError def write(self,filename,files=1): raise NotImplementedError