<commit_before>#!/usr/bin/env python

import os
import numpy as np
import scipy.stats as stats
import scipy.constants as const
from astropy.io import fits

from .sqbase import datadir

# shorthands
exp,sqrt = np.exp,np.sqrt
c_kms = const.c/1e3
sqrt_pi = sqrt(np.pi)
sigma_c = 6.33e-18 # cm^-2
fourpi = 4*np.pi

def _getlinelistdata():
	# Line list obtained from Prochaska's XIDL code
	# https://svn.ucolick.org/xidl/trunk/Spec/Lines/all_lin.fits
	linelist = fits.getdata(datadir+'all_lin.fits')
	Hlines = np.array([i for i in range(linelist.size) 
	                       if 'HI' in linelist.ION[i]])
	transitionParams = {}
	for n,idx in enumerate(Hlines[::-1],start=2):
		transitionParams[n] = (linelist.WREST[idx],
		                       linelist.F[idx],
		                       linelist.GAMMA[idx])
	return transitionParams

transitionParams = _getlinelistdata()

# default is to go up to 32->1
default_lymanseries_range = (2,33)

Fan99_model = {
  'forest':{'zrange':(0.0,6.0),
            'logNHrange':(13.0,17.3),
            'N0':50.3,
            'gamma':2.3,
            'beta':1.41,
            'b':30.0},
     'LLS':{'zrange':(0.0,6.0),
            'logNHrange':(17.3,20.5),
            'N0':0.27,
            'gamma':1.55,
            'beta':1.25,
            'b':70.0},
     'DLA':{'zrange':(0.0,6.0),
            'logNHrange':(20.5,22.0),
            'N0':0.04,
            'gamma':1.3,
            'beta':1.48,
            'b':70.0},
}

WP11_model = {
 'forest0':{'zrange':(0.0,1.5),
            'logNHrange':(12.0,19.0),
            'gamma':0.2,
            'beta':1.55,
            'B':0.0170,
            'N0':340.,
            'brange':(10.,100.),
            'bsig':24.0},
 'forest1':{'zrange':(1.5,4.6),
            'logNHrange':(12.0,14.5),
            'gamma':2.04,
            'beta':1.50,
            'B':0.0062,
            'N0':102.0,
            'brange':(10.,100.),
            'bsig':24.0},
 'forest2':{'zrange':(1.5,4.6),
            'logNHrange':(14.5,17.5),
            'gamma':2.04,
            'beta':1.80,
            'B':0.0062,
            'N0':4.05,
            'brange':(10.,100.),
            'bsig':24.0},
 'forest3':{'zrange':(1.5,4.6),
            'logNHrange':(17.5,19.0),
            'gamma':2.04,
            'beta':0.90,
            'B':0.0062,
            'N0':0.051,
            'brange':(10.,100.),
            'bsig':24.0},
    'SLLS':{'zrange':(0.0,4.6),
            'logNHrange':(19.0,20.3),
            'N0':0.0660,
            'gamma':1.70,
            'beta':1.40,
            'brange':(10.,100.),
            'bsig':24.0},
     'DLA':{'zrange':(0.0,4.6),
            'logNHrange':(20.3,22.0),
            'N0':0.0440,
            'gamma':1.27,
            'beta':2.00,
            'brange':(10.,100.),
            'bsig':24.0},
}

forestModels = {'Fan1999':Fan99_model,
                'Worseck&Prochaska2011':WP11_model}

def generate_los(model,zmin,zmax):
	'''Given a model for the distribution of absorption systems, generate
	   a random line-of-sight populated with absorbers.
	   returns (z,logNHI,b) for each absorption system.
	'''
	abs_dtype = [('z',np.float32),('logNHI',np.float32),('b',np.float32)]
	absorbers = []
	for component,p in model.items():
		if zmin > p['zrange'][1] or zmax < p['zrange'][0]:
			# outside the redshift range of this forest component
			continue
		# parameters for the forest component (LLS, etc.) absorber distribution
		NHImin,NHImax = p['logNHrange']
		NHImin,NHImax = 10**NHImin,10**NHImax
		z1 = max(zmin,p['zrange'][0])
		z2 = min(zmax,p['zrange'][1])
		beta = p['beta'] 
		mbeta1 = -beta+1
		gamma1 = p['gamma'] + 1
		# expectation for the number of absorbers at this redshift
		#  (inverting n(z) = N0*(1+z)^gamma)
		N = (p['N0']/gamma1) * ( (1+z2)**gamma1 - (1+z1)**gamma1 )
		# sample from a Poisson distribution for <N>
		n = stats.poisson.rvs(N,size=1)[0]
		# invert the dN/dz CDF to get the sample redshifts
		x = np.random.random_sample(n)
		z = (1+z1)*((((1+z2)/(1+z1))**gamma1 - 1)*x + 1)**(1/gamma1) - 1
		# invert the NHI CDF to get the sample column densities
		x = np.random.random_sample(n)<commit_msg>add forest model used in McGreer+2013 for z<~6<commit_after>#!/usr/bin/env python

import os
import numpy as np
import scipy.stats as stats
import scipy.constants as const
from astropy.io import fits

from .sqbase import datadir

# shorthands
exp,sqrt = np.exp,np.sqrt
c_kms = const.c/1e3
sqrt_pi = sqrt(np.pi)
sigma_c = 6.33e-18 # cm^-2
fourpi = 4*np.pi

def _getlinelistdata():
	# Line list obtained from Prochaska's XIDL code
	# https://svn.ucolick.org/xidl/trunk/Spec/Lines/all_lin.fits
	linelist = fits.getdata(datadir+'all_lin.fits')
	Hlines = np.array([i for i in range(linelist.size) 
	                       if 'HI' in linelist.ION[i]])
	transitionParams = {}
	for n,idx in enumerate(Hlines[::-1],start=2):
		transitionParams[n] = (linelist.WREST[idx],
		                       linelist.F[idx],
		                       linelist.GAMMA[idx])
	return transitionParams

transitionParams = _getlinelistdata()

# default is to go up to 32->1
default_lymanseries_range = (2,33)

Fan99_model = {
  'forest':{'zrange':(0.0,6.0),
            'logNHrange':(13.0,17.3),
            'N0':50.3,
            'gamma':2.3,
            'beta':1.41,
            'b':30.0},
     'LLS':{'zrange':(0.0,6.0),
            'logNHrange':(17.3,20.5),
            'N0':0.27,
            'gamma':1.55,
            'beta':1.25,
            'b':70.0},
     'DLA':{'zrange':(0.0,6.0),
            'logNHrange':(20.5,22.0),
            'N0':0.04,
            'gamma':1.3,
            'beta':1.48,
            'b':70.0},
}

WP11_model = {
 'forest0':{'zrange':(0.0,1.5),
            'logNHrange':(12.0,19.0),
            'gamma':0.2,
            'beta':1.55,
            'B':0.0170,
            'N0':340.,
            'brange':(10.,100.),
            'bsig':24.0},
 'forest1':{'zrange':(1.5,4.6),
            'logNHrange':(12.0,14.5),
            'gamma':2.04,
            'beta':1.50,
            'B':0.0062,
            'N0':102.0,
            'brange':(10.,100.),
            'bsig':24.0},
 'forest2':{'zrange':(1.5,4.6),
            'logNHrange':(14.5,17.5),
            'gamma':2.04,
            'beta':1.80,
            'B':0.0062,
            'N0':4.05,
            'brange':(10.,100.),
            'bsig':24.0},
 'forest3':{'zrange':(1.5,4.6),
            'logNHrange':(17.5,19.0),
            'gamma':2.04,
            'beta':0.90,
            'B':0.0062,
            'N0':0.051,
            'brange':(10.,100.),
            'bsig':24.0},
    'SLLS':{'zrange':(0.0,4.6),
            'logNHrange':(19.0,20.3),
            'N0':0.0660,
            'gamma':1.70,
            'beta':1.40,
            'brange':(10.,100.),
            'bsig':24.0},
     'DLA':{'zrange':(0.0,4.6),
            'logNHrange':(20.3,22.0),
            'N0':0.0440,
            'gamma':1.27,
            'beta':2.00,
            'brange':(10.,100.),
            'bsig':24.0},
}

McG13hiz_model = {
 'forest1':{'zrange':(1.5,6.1),
            'logNHrange':(12.0,14.5),
            'gamma':3.5,
            'beta':1.50,
            'N0':8.5 * 1.1,
            'brange':(10.,100.),
            'bsig':24.0},
 'forest2':{'zrange':(1.5,6.1),
            'logNHrange':(14.5,17.2),
            'gamma':3.5,
            'beta':1.70,
            'N0':0.33 * 1.1,
            'brange':(10.,100.),
            'bsig':24.0},
     'LLS':{'zrange':(1.5,6.1),
            'logNHrange':(17.2,20.3),
            'gamma':2.0,
            'beta':1.3,
            'N0':0.13 * 1.1,
            'brange':(10.,100.),
            'bsig':24.0},
  'subDLA':{'zrange':(0.0,6.1),
            'logNHrange':(20.3,21.0),
            'N0':0.13 / 7.5 * 1.1,
            'gamma':1.70,
            'beta':1.28,
            'brange':(10.,100.),
            'bsig':24.0},
     'DLA':{'zrange':(0.0,6.1),
            'logNHrange':(21.0,22.0),
            'N0':0.13 / 33 * 1.1,
            'gamma':2.0,
            'beta':1.40,
            'brange':(10.,100.),
            'bsig':24.0},
}

forestModels = {'Fan1999':Fan99_model,
                'Worseck&Prochaska2011':WP11_model,
                'McGreer+2013':McG13hiz_model}

def generate_los(model,zmin,zmax):
	'''Given a model for the distribution of absorption systems, generate
	   a random line-of-sight populated with absorbers.
	   returns (z,logNHI,b) for each absorption system.
	'''
	abs_dtype = [('z',np.float32),('logNHI',np.float32),('b',np.float32)]
	absorbers = []
	for component,p in model.items():
		if zmin > p['zrange'][1] or zmax < p['zrange'][0]:
			# outside the redshift range of this forest component
			continue
		# parameters for the forest component (LLS, etc.) absorber distribution
		NHImin,NHImax = p['logNHrange']
		NHImin,NHImax = 10**NHImin,10**NHImax
		z1 = max(zmin,p['zrange'][0])
		z2 = min(zmax,p['zrange'][1])
		beta = p['beta'] 
		mbeta1 = -beta+1
		gamma1 = p['gamma'] + 1
		# expectation for the number of absorbers at this redshift
		#  (inverting n(z) = N0*(1+z)^gamma)
		N = (p['N0']/gamma1) * ( (1+z2)**gamma1 - (1+z1)**gamma1 )
		# sample from a Poisson distribution for <N>
		n = stats.poisson.rvs(N,size=1)[0]
		# invert the dN/dz CDF to get the sample redshifts
		x = np.random.random_sample(n)
		z = (1+z1)*((((1+z2)/(1+z1))**gamma1 - 1)*x + 1)**(1/gamma1) - 1
		# invert the NHI CDF to get the sample column densities
		x = np.random.random_sample(n)
{'fraction_non_alphanumeric': 0.20008514261387825, 'fraction_numerical': 0.0899318859088974, 'mean_word_length': 1.9634184799747714, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 9, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7204662', 'n_tokens_mistral': 4169, 'n_tokens_neox': 3726, 'n_words': 594}