Skip to content

Commit

Permalink
Merge pull request #7 from kpolsen/KPO-branch
Browse files Browse the repository at this point in the history 
edits to comments within the code
  • Loading branch information
Karen Pardos Olsen authored Sep 11, 2017
2 parents 2fb7ea5 + e967956 commit 172ea02
Show file tree
Hide file tree
Showing 4 changed files with 97 additions and 135 deletions.
107 changes: 30 additions & 77 deletions sigame/GMC_module.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,8 +26,8 @@
def create_GMCs(galname=galnames[0],zred=zreds[0],verbose=False):
'''
Purpose
---------
Splits the molecular part of each SPH particle into GMCs (called by main.run)
-------
Splits the molecular part of each SPH particle into GMCs (called by main.run)
Arguments
---------
Expand All @@ -41,30 +41,31 @@ def create_GMCs(galname=galnames[0],zred=zreds[0],verbose=False):
default = False
'''
print('\n ** Split gas elements into GMCs **')

print('\n ** Create GMCs **')

plt.close('all') # close all windows
# Read in raw simulation data files (gas and stars) for galname simulation

# Load simulation data for gas and stars
simgas = pd.read_pickle('sigame/temp/sim_FUV/z'+'{:.2f}'.format(zred)+'_'+galname+'_sim1.gas')
simstar = pd.read_pickle('sigame/temp/sim_FUV/z'+'{:.2f}'.format(zred)+'_'+galname+'_sim1.star')
if verbose:
print(galname)
print('# gas element particles: ',len(simgas))
print('# star element particles: ',len(simstar))
# Neutral gas mass:
# Mneu = simgas['m'].values*simgas['f_H2'].values
print('Using f_H2 from simulation')
Mneu = simgas['m'].values*simgas['f_H21'].values

if verbose: print('Number of gas elements: %s' % str(len(simgas)))
if verbose: print('Number of stars: %s' % str(len(simstar)))

print('Using f_H2 from simulation to estimate neutral gas mass per gas element')
Mneu = simgas['m'].values*simgas['f_H2'].values
print('Total dense gas mass: '+str(np.sum(Mneu))+' Msun')
print('Dense gas mass fraction out of total ISM mass: '+str(np.sum(Mneu)/np.sum(simgas['m'])*100.)+' %')
print('Min and max particle dense gas mass: '+str(np.min(Mneu))+' '+str(np.max(Mneu))+' Msun')
print('Throwing away this much gas mass: '+str(np.sum(simgas.loc[Mneu < 1e4]['m']))+' Msun')
print('In percent: '+str(np.sum(simgas.loc[Mneu < 1e4]['m'])/np.sum(simgas['m']))+' %')

# Create mass spectrum (probability function = dN/dM normalized)
b = 1.8
if ext_DENSE == '_b3p0': b = 3.0 # powerlaw slope [Blitz+07]
if ext_DENSE == '_b1p5': b = 1.5 # powerlaw slope [Blitz+07]
print('b used is %s' % b)
b = 1.8 # MW powerlaw slope [Blitz+07]
if ext_DENSE == '_b3p0': b = 3.0 # outer MW and M33 powerlaw slope [Rosolowsky+05, Blitz+07]
if ext_DENSE == '_b1p5': b = 1.5 # inner MW powerlaw slope [Rosolowsky+05, Blitz+07]
print('Powerlaw slope for mass spectrum (beta) used is %s' % b)
Mmin = 1.e4 # min mass of GMC
Mmax = 1.e6 # max mass of GMC
tol = Mmin # precision in reaching total mass
Expand All @@ -79,20 +80,21 @@ def create_GMCs(galname=galnames[0],zred=zreds[0],verbose=False):
print('# Gas particles with Mmol > 2*1e4 Msun: %s out of %s' % (np.size(Mneu[Mneu > 2*Mmin]),np.size(Mneu)))
simgas.index = range(0,len(simgas))
j = 0
print('Starting up multiprocessing')
pool = mp.Pool(processes=5) # 8 processors on my Mac Pro, 16 on Betty
print('(Multiprocessing starting up!)')
pool = mp.Pool(processes=3) # 8 processors on my Mac Pro, 16 on Betty
np.random.seed(len(GMCgas)) # so we don't end up with the same random numbers for every galaxy
results = [pool.apply_async(GMC_generator, args=(i,Mneu,h,Mmin,Mmax,b,tol,nn,)) for i in range(0,len(simgas))]#
GMCs = [p.get() for p in results]
print('(Multiprocessing done!)')
GMCgas = simgas.iloc[0:GMCs[0][0]]
Mgmc = GMCs[0][1]
newx = simgas.loc[0]['x']+GMCs[0][2]
newy = simgas.loc[0]['y']+GMCs[0][3]
newz = simgas.loc[0]['z']+GMCs[0][4]
i1 = 0
part = 0.1
print('Make new (GMC) dataframe')
for i in range(1,len(simgas)):#
print('Append results to new dataframe')
for i in range(1,len(simgas)):
for ii in range (0,GMCs[i][0]):
GMCgas = pd.DataFrame.append(GMCgas,simgas.loc[i],ignore_index=True)
Mgmc = np.append(Mgmc,GMCs[i][1])
Expand Down Expand Up @@ -185,13 +187,18 @@ def calc_line_emission(galname=galnames[0],zred=zreds[0],SFRsd=SFRsd_MW):
default = SFR surface density of the MW
'''
print('\n ** Calculate total line emission etc. from GMCs **')

print('\n ** Calculate total line emission from GMCs **')

plt.close('all') # close all windows

ext_DENSE0 = ext_DENSE

# Load dataframe with GMCs
GMCgas = pd.read_pickle('sigame/temp/GMC/z'+'{:.2f}'.format(zred)+'_'+galname+'_GMC.gas')
if ext_DENSE0 == '_b1p5': GMCgas = pd.read_pickle('sigame/temp/GMC/z'+'{:.2f}'.format(zred)+'_'+galname+'_GMC_b1p5.gas')
if ext_DENSE0 == '_b3p0': GMCgas = pd.read_pickle('sigame/temp/GMC/z'+'{:.2f}'.format(zred)+'_'+galname+'_GMC_b3p0.gas')

nGMC = len(GMCgas)
print('Number of GMCs: %s' % nGMC)
print('Mass of GMCs: %s Msun' % np.sum(GMCgas['m']))
Expand All @@ -211,12 +218,6 @@ def calc_line_emission(galname=galnames[0],zred=zreds[0],SFRsd=SFRsd_MW):
cloudy_grid_param = pd.read_pickle('cloudy_models/GMC/grids/GMCgrid'+ext_DENSE0+'_'+z1+'.param')
cloudy_grid = pd.read_pickle('cloudy_models/GMC/grids/GMCgrid'+ext_DENSE0+'_'+z1+'_CII.models')

print('Z parameters: ')
print(cloudy_grid_param['Zs'])
print('Mass parameters: ')
print(cloudy_grid_param['Mgmcs'])
print('Min and max of actual clouds: %s and %s' % (min(GMCgas['Z']),max(GMCgas['Z'])))

# make sure we don't go outside of grid:
GMCgas1 = np.log10(GMCgas) # Going to edit a bit in the DataFrame
GMCgas1.ix[GMCgas1['m'] < min(cloudy_grid_param['Mgmcs']),'m'] = min(cloudy_grid_param['Mgmcs'])
Expand All @@ -229,7 +230,7 @@ def calc_line_emission(galname=galnames[0],zred=zreds[0],SFRsd=SFRsd_MW):
GMCgas1.ix[GMCgas1['P_ext'] > max(cloudy_grid_param['P_exts']),'P_ext'] = max(cloudy_grid_param['P_exts'])
# Values used for interpolation in cloudy models:
GMCs = np.column_stack((GMCgas1['m'].values,GMCgas1['FUV'].values,GMCgas1['Z'].values,GMCgas1['P_ext'].values))

# List of target items that we are interpolating for:
target_list = ['Mgmc_fit','FUV_fit','Z_fit','P_ext_fit',\
'f_H2','f_HI','m_dust','m_H','m_H2','m_HI','m_HII',\
Expand All @@ -242,10 +243,6 @@ def calc_line_emission(galname=galnames[0],zred=zreds[0],SFRsd=SFRsd_MW):
cloudy_grid['Z_fit'] = cloudy_grid['Z']
cloudy_grid['P_ext_fit'] = cloudy_grid['P_ext']

print('Check:')
print(np.log10(GMCgas['FUV'][0:10].values))


for target in target_list:
# Make grid of corresponding target grid values:
target_grid = np.zeros([len(cloudy_grid_param['Mgmcs']),len(cloudy_grid_param['FUVs']),len(cloudy_grid_param['Zs']),len(cloudy_grid_param['P_exts'])])
Expand Down Expand Up @@ -286,43 +283,6 @@ def calc_line_emission(galname=galnames[0],zred=zreds[0],SFRsd=SFRsd_MW):
i_P_ext = np.argmin(abs(cloudy_grid_param['P_exts']-lP_ext))
GMCgas['closest_model_i'][i] = int(model_numbers[i_Mgmc,i_FUV,i_Z,i_P_ext])

print('Check sum of L_[CII] using closest model number!')
L_tot = 0.
for i in range(0,nGMC):
L_tot += cloudy_grid['L_CII'][GMCgas['closest_model_i'][i]]
print(L_tot)
print('To sum of interpolation: %s ' % np.sum(GMCgas['L_CII']))

# # Histogram of all model L_[CII] values
# fig = plt.figure(0)
# ax1 = fig.add_subplot(1,1,1)
# n, bins, patches = plt.hist(np.log10(cloudy_grid['L_CII'][cloudy_grid['L_CII']>0]), \
# 50, normed=1, facecolor='green', alpha=0.75)
# plt.show(block=False)

# Quick check!
# print('\n')
# print('Model numbers with very high FUV:')
# very_high = cloudy_grid['index'][cloudy_grid['FUV']>4].values
# print(very_high)
# for model_number in very_high:
# model_number = int(model_number)
# print('Model %s has L_CII: %s Lsun ' % (model_number,cloudy_grid['L_CII'][model_number]))
# print('With log Mgmc: %s, FUV: %s, Z: %s, P_ext: %s' % (cloudy_grid['Mgmc'][model_number],cloudy_grid['FUV'][model_number],cloudy_grid['Z'][model_number],cloudy_grid['P_ext'][model_number]))
# dex = 0.1
# GMCgas1['L_CII'] = GMCgas['L_CII']
# # pdb.set_trace()
# GMCs = aux.within_dex(GMCgas1,cloudy_grid['Mgmc'][model_number],'m',dex)
# GMCs = aux.within_dex(GMCs,cloudy_grid['FUV'][model_number],'FUV',dex)
# GMCs = aux.within_dex(GMCs,cloudy_grid['Z'][model_number],'Z',dex)
# GMCs = aux.within_dex(GMCs,cloudy_grid['P_ext'][model_number],'P_ext',0.5)
# print('GMCs with parameters within 1 %s of these: %s' % (dex,len(GMCs)))
# print('Their [CII] luminosities: ' )
# print(GMCs['L_CII'].values)
# print('\n')

# pdb.set_trace()

print('Total GMC mass available: '+str(sum(GMCgas['m'])/1e8)+' x 10^8 Msun')
print('Mass of selected models: '+str(sum(GMCgas['m_H'])/1e8)+' x 10^8 Msun')
print('Total mass of CII: '+str(sum(GMCgas['m_CII'])/1e8)+' x 10^8 Msun')
Expand All @@ -341,16 +301,9 @@ def calc_line_emission(galname=galnames[0],zred=zreds[0],SFRsd=SFRsd_MW):
for line in lines:
print('Total L_'+line+': %s x 10^8 L_sun' % (sum(GMCgas['L_'+line])/1e8))

print('Which models are "most popular"? Printing the first 20 (model number: # occurences)')
GMCs = np.column_stack((GMCgas1['m'].values,GMCgas1['FUV'].values,GMCgas1['Z'].values,GMCgas1['P_ext'].values))
model_numbers = np.arange(0,len(cloudy_grid))
close_model_numbers = griddata((cloudy_grid['Mgmc'],cloudy_grid['FUV'],cloudy_grid['Z'],cloudy_grid['P_ext']),model_numbers,GMCs,method='nearest')
counter = collections.Counter(close_model_numbers)
print(counter.most_common(20))

print('Add mass-weighted density and G0 to GMCgas dataframe for diagnostic line ratio plot')

sigma_vs = 1.2*(GMCgas['P_ext'].values/1e4)**(1.0/4.0)*(GMCgas['Rgmc'].values)**(1.0/2.0)
sigma_vs = 1.2*(GMCgas['P_ext'].values/1e4)**(1.0/4.0)*(GMCgas['Rgmc'].values)**(1.0/2.0)
GMCgas['vel_disp_gas'] = sigma_vs
GMCgas.to_pickle(GMC_path+'z'+'{:.2f}'.format(zred)+'_'+galname+'_GMC'+ext_DENSE+'_em.gas')

Expand Down
7 changes: 2 additions & 5 deletions sigame/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,12 +28,9 @@
# '-- A code to simulate the far-IR emission lines of the ISM --'+'\n'+\
# '----------- in galaxies from hydrodynamical codes ------------'+'\n'+\
# '--- for the interpretation and prediction of observations. ---'+'\n'+\
# '-- Contact: Karen Pardos Olsen, kpolsen (at) asu.edu (2016) --'+'\n'+\
# '-- Contact: Karen Pardos Olsen, kpolsen (at) asu.edu (2017) --'+'\n'+\
# ''+'\n')
# style from http://www.kammerl.de/ascii/AsciiSignature.php
# (alternatives: epic, roman, blocks, varsity, pepples, soft, standard, starwars)

################ Some instructions ######################
# the font style above is from http://www.kammerl.de/ascii/AsciiSignature.php

# Remember to select and edit a parameter file:
params_file = 'parameters_z6.txt'
Expand Down
55 changes: 25 additions & 30 deletions sigame/dif_module.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,22 +36,29 @@ def create_dif(galname=galnames[0],zred=zreds[0]):
'''

print('\n ** Identify diffuse gas containing hot, ionized medium (DIG) and warm, neutral medium (DNG)**')
plt.close('all') # close all windows
simgas = pd.read_pickle('sigame/temp/sim_FUV/z'+'{:.2f}'.format(zred)+'_'+galname+'_sim1.gas')
difgas = simgas.copy()
# difgas['m'] = difgas['m'].values*(1.-difgas['f_HI'].values*difgas['f_H2'].values)
difgas['m'] = difgas['m'].values*(1.-difgas['f_H2'].values)
print('Total gas mass: '+str(sum(simgas['m'])))
print('Diffuse gas mass: '+str(sum(difgas['m'])))

# Load simulation data for gas
simgas = pd.read_pickle('sigame/temp/sim_FUV/z'+'{:.2f}'.format(zred)+'_'+galname+'_sim1.gas')

# Start new dataframe with only the diffuse gas
difgas = simgas.copy()
difgas['m'] = difgas['m'].values*(1.-difgas['f_H2'].values)
print('Total gas mass in galaxy: '+str(sum(simgas['m'])))
print('Diffuse gas mass in galaxy: '+str(sum(difgas['m'])))
print('in percent: '+str(sum(difgas['m'])/sum(simgas['m'])*100.)+'%')
difgas['R'] = difgas['h']
difgas['nH'] = 0.75*np.array(difgas['m'],dtype=np.float64)/(4/3.*np.pi*np.array(difgas['R'],dtype=np.float64)**3.)*Msun/mH/kpc2cm**3 # Hydrogen atoms per cm^3
difgas['R'][difgas['m'] == 0] = 0
difgas['nH'][difgas['m'] == 0] = 0
dir_gas = 'sigame/temp/dif/'

# Calculate radius of diffuse gas clouds
difgas['R'] = difgas['h']
difgas['R'][difgas['m'] == 0] = 0

# Calculate density of diffuse gas clouds
difgas['nH'] = 0.75*np.array(difgas['m'],dtype=np.float64)/(4/3.*np.pi*np.array(difgas['R'],dtype=np.float64)**3.)*Msun/mH/kpc2cm**3 # Hydrogen atoms per cm^3
difgas['nH'][difgas['m'] == 0] = 0

# Save results
dir_gas = 'sigame/temp/dif/'
difgas.to_pickle(dir_gas+'z'+'{:.2f}'.format(zred)+'_'+galname+'_dif.gas')
# pdb.set_trace()

def calc_line_emission(galname=galnames[0],zred=zreds[0],SFRsd=SFRsd_MW):
'''
Expand All @@ -67,18 +74,19 @@ def calc_line_emission(galname=galnames[0],zred=zreds[0],SFRsd=SFRsd_MW):
zred: redshift of galaxy - float/int
default = first redshift name in redshift list from parameter file
SFRsd: SFR surface density of this galaxy - float
default = SFR surface density of the MW
'''

printe('\n ** Calculate total line emission from diffuse gas **')

plt.close('all') # close all windows

ext_DIFFUSE0 = ext_DIFFUSE

# Load galaxy
# Load dataframe with diffuse gas
difgas = pd.read_pickle('sigame/temp/dif/z'+'{:.2f}'.format(zred)+'_'+galname+'_dif.gas')
# TEST!!!
if ext_DIFFUSE0 == '_Z0p05': difgas['Z'],ext_DIFFUSE0 = difgas['Z'].values*0.+0.05,'_Ztest'
if ext_DIFFUSE0 == '_Z1': difgas['Z'],ext_DIFFUSE0 = difgas['Z'].values*0.+1.,'_Ztest'
if ext_DIFFUSE0 == '_Zx3': difgas['Z'],ext_DIFFUSE0 = difgas['Z'].values*3.,'_highZ'
Expand All @@ -89,22 +97,17 @@ def calc_line_emission(galname=galnames[0],zred=zreds[0],SFRsd=SFRsd_MW):
difgas = pd.DataFrame.reset_index(difgas)
ndif = len(difgas)

# Load the grid that best represents this Sigma_SFR
# Load the grid that best represents the SFR surface density of this galaxy
UV = [5,35]
if ext_DIFFUSE == '_FUV': UV,ext_DIFFUSE0 = [5,15,25,35,45,120],'_ism'
UV1 = str(int(UV[np.argmin(np.abs(np.array(UV)-SFRsd/SFRsd_MW))]))
# UV1 = '120'
cloudy_grid_param = pd.read_pickle('cloudy_models/dif/grids/difgrid_'+UV1+'UV'+ext_DIFFUSE0+'_'+z1+'.param')
cloudy_grid = pd.read_pickle('cloudy_models/dif/grids/difgrid_'+UV1+'UV'+ext_DIFFUSE0+'_'+z1+'_CII.models')

print('SFR surface density is '+str(SFRsd/SFRsd_MW)+' x that of MW ')
print('Using grid at '+UV1+' x ISM FUV field')
difgas['FUV'] = difgas['FUV']*0.+UV[np.argmin(np.abs(np.array(UV)-SFRsd/SFRsd_MW))]

print('Z parameters: ')
print(cloudy_grid_param['Zs'])
print('Min and max of actual clouds: %s and %s' % (min(difgas['Z']),max(difgas['Z'])))

# make sure we don't go outside of grid:
difgas1 = difgas.copy()
difgas1 = np.log10(difgas1) # Going to edit a bit in the DataFrame
Expand Down Expand Up @@ -232,13 +235,6 @@ def calc_line_emission(galname=galnames[0],zred=zreds[0],SFRsd=SFRsd_MW):
print('Z: %s > 0' % (np.mean(np.log10(difgas1['Z'][difgas1['Z']>0]))))
print('Tk: %s' % (np.mean(np.log10(difgas1['Tk']))))

# print('Which models are "most popular"? Printing the first 20 (model number: # occurences)')
# difs = np.column_stack((difgas1['nH'].values,difgas1['R'].values,difgas1['Z'].values,difgas1['Tk'].values))
# model_numbers = np.arange(0,len(cloudy_grid))
# close_model_numbers = griddata((cloudy_grid['nH'],cloudy_grid['R'],cloudy_grid['Z'],cloudy_grid['Tk']),model_numbers,difs,method='nearest')
# counter = collections.Counter(close_model_numbers)
# print(counter.most_common(20))

print('Total [CII] luminosity from diffuse gas (DIG/DNG): '+str(sum(difgas1['L_CII_DNG']+difgas1['L_CII_DIG'])/1e8)+' x 10^8 L_sun')

difgas1.to_pickle(dif_path+'z'+'{:.2f}'.format(zred)+'_'+galname+'_dif'+ext_DIFFUSE0+'_em.gas')
Expand All @@ -261,4 +257,3 @@ def calc_line_emission(galname=galnames[0],zred=zreds[0],SFRsd=SFRsd_MW):
data.to_pickle('sigame/temp/line_profiles/'+galname+'.diffuse')

return(dif_results)

Loading

0 comments on commit 172ea02

Please sign in to comment.