From e967956c7037198a04e186ca11e2ce6cb13c29c6 Mon Sep 17 00:00:00 2001 From: kpolsen Date: Mon, 11 Sep 2017 12:29:39 -0700 Subject: [PATCH] edits to comments within the code --- sigame/GMC_module.py | 107 +++++++++++---------------------------- sigame/__init__.py | 7 +-- sigame/dif_module.py | 55 +++++++++----------- sigame/subgrid_module.py | 63 ++++++++++++++--------- 4 files changed, 97 insertions(+), 135 deletions(-) diff --git a/sigame/GMC_module.py b/sigame/GMC_module.py index f5af9fb..cffadeb 100644 --- a/sigame/GMC_module.py +++ b/sigame/GMC_module.py @@ -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 --------- @@ -41,19 +41,20 @@ 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') @@ -61,10 +62,10 @@ def create_GMCs(galname=galnames[0],zred=zreds[0],verbose=False): 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 @@ -79,11 +80,12 @@ 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] @@ -91,8 +93,8 @@ def create_GMCs(galname=galnames[0],zred=zreds[0],verbose=False): 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]) @@ -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'])) @@ -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']) @@ -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',\ @@ -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'])]) @@ -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') @@ -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') diff --git a/sigame/__init__.py b/sigame/__init__.py index 761a1cb..3d86fae 100644 --- a/sigame/__init__.py +++ b/sigame/__init__.py @@ -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' diff --git a/sigame/dif_module.py b/sigame/dif_module.py index 273387f..7dffb4c 100644 --- a/sigame/dif_module.py +++ b/sigame/dif_module.py @@ -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): ''' @@ -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' @@ -89,11 +97,10 @@ 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') @@ -101,10 +108,6 @@ def calc_line_emission(galname=galnames[0],zred=zreds[0],SFRsd=SFRsd_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 @@ -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') @@ -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) - diff --git a/sigame/subgrid_module.py b/sigame/subgrid_module.py index a2eff6a..e2c5fd2 100644 --- a/sigame/subgrid_module.py +++ b/sigame/subgrid_module.py @@ -12,28 +12,49 @@ import subprocess as sub import aux as aux import matplotlib.pyplot as plt -# From SIGAME submodules: -# import sigame.plot as plot params = np.load('sigame/temp/params.npy').item() for key,val in params.items(): exec(key + '=val') def subgrid(galname=galnames[0],zred=zreds[0]): - print('\n ** Use grid of starburst models to calculate local FUV field and pressure **') + ''' + Purpose + ------- + On subgrid scales calculate: + - local FUV field using grid of starburst99 models + - local pressure using velocity dispersion, surface densities of gas and stars + and hydrostatic equilibrium + + Arguments + --------- + galname: galaxy name - str + default = first galaxy name in galnames list from parameter file + + zred: redshift of galaxy - float/int + default = first redshift name in redshift list from parameter file + + ''' + + print('\n** Derive local FUV field and pressure **') + + # Declare these variables to be global for availability in other + # functions called by subprocess: global simgas,simgas1,simstar,L_FUV + + # Load simulation data for gas and stars simgas = pd.read_pickle(d_sim+'z'+'{:.2f}'.format(zred)+'_'+galname+'_sim0.gas') - simstar = pd.read_pickle(d_sim+'z'+'{:.2f}'.format(zred)+'_'+galname+'_sim0.star') - # simgas = simgas[0:100] - # simstar = simstar[0:100] - print('Number of gas elements: '+str(len(simgas))) + + if verbose: print('Number of gas elements: %s' % str(len(simgas))) print('\nFind local FUV field from stellar population synthesis! ') - grid = pd.read_pickle(d_t+'FUV/FUVgrid_'+z1+'_noneb') # read grid parameters for FUV grid of 1e4 stellar populations - FUV = pd.read_pickle(d_t+'FUV/FUV_'+z1+'_noneb') # [age,Z,L_FUV] dataframe + # Read grid parameters for FUV grid of 1e4 Msun stellar populations + grid = pd.read_pickle(d_t+'FUV/FUVgrid_'+z1+'_noneb') + # Read corresponding [age,Z,L_FUV] values + FUV = pd.read_pickle(d_t+'FUV/FUV_'+z1+'_noneb') l_FUV = FUV['L_FUV'].values l_FUV = l_FUV.reshape((len(grid['Ages']),len(grid['Zs']))) - print('Calculate FUV luminosity of each stellar particle') + print('First, calculate FUV luminosity of each stellar particle') part = 0.1 L_FUV = np.zeros(len(simstar)) for i in range(0,len(simstar)): @@ -43,31 +64,29 @@ def subgrid(galname=galnames[0],zred=zreds[0]): print(int(part*100),' % done!') part = part+0.1 simstar['L_FUV'] = L_FUV - print('Find FUV flux at gas particle positions') + print('Then, find FUV flux at gas particle positions') FUV = np.zeros(len(simgas)) - t1 = time.time() - print('Multiprocessing starting up!') + print('(Multiprocessing starting up!)') pool = mp.Pool(processes=3) # 8 processors on my Mac Pro, 16 on Betty results = [pool.apply_async(FUVfunc, args=(i,)) for i in range(0,len(simgas))]#len(simgas) FUV = [p.get() for p in results] - t2 = time.time() - dt2 = t2-t1 - # print('Time that it took: dt2') + print('(Multiprocessing done!)') simgas['FUV'] = FUV simgas['FUV'] = simgas['FUV']/(kpc2cm**2*FUV_ISM) - print('Assume that CR intensity follows FUV field for now') + print('Finally, scale local CR intensity to follow fluctuations in local FUV field') simgas['CR'] = simgas['FUV']*CR_ISM - print('\nFind local hydrostatic mid-plane pressure') - print('Do calculations of pressure using only star-forming gas!!') + print('\nFind local hydrostatic mid-plane pressure!') + # Extract star forming gas only: simgas1 = simgas.copy() simgas1 = simgas1[simgas1['SFR'] > 0].reset_index() global m_gas,m_star m_gas,m_star = simgas1['m'].values,simstar['m'].values - print('Multiprocessing starting up!') + print('(Multiprocessing starting up!)') pool = mp.Pool(processes=3) # 8 processors on my Mac Pro, 16 on Betty results = [pool.apply_async(Pfunc, args=(i,)) for i in range(0,len(simgas))]#len(simgas) res = [p.get() for p in results] + print('(Multiprocessing done!)') P_ext,surf_gas,surf_star,sigma_gas,sigma_star,vel_disp_gas = np.zeros(len(res)),np.zeros(len(res)),np.zeros(len(res)),np.zeros(len(res)),np.zeros(len(res)),np.zeros(len(res)) for i in range(0,len(res)): P_ext[i],surf_gas[i],surf_star[i],sigma_gas[i],sigma_star[i],vel_disp_gas[i] = res[i][0],res[i][1],res[i][2],res[i][3],res[i][4],res[i][5] @@ -75,13 +94,11 @@ def subgrid(galname=galnames[0],zred=zreds[0]): simgas['sigma_gas'] = sigma_gas simgas['sigma_star'] = sigma_star simgas['vel_disp_gas'] = vel_disp_gas - simgas['P_ext1'] = simgas['nH']*simgas['vel_disp_gas']**2.*1e6/kB # density [kg/cm^3] * disp [m/s]^2 / kB [J/K] -> K/cm^3 print('Min and max of log(P_ext): %s and %s' % (min(np.log10(simgas['P_ext'])),max(np.log10(simgas['P_ext'])))) - print('Min and max of alternative log(P_ext): %s and %s' % (min(np.log10(simgas['P_ext1'])),max(np.log10(simgas['P_ext1'])))) simgas['surf_gas'] = surf_gas simgas['surf_star'] = surf_star - # Save + # Save results simgas.to_pickle('sigame/temp/sim_FUV/z'+'{:.2f}'.format(zred)+'_'+galname+'_sim1.gas') simstar.to_pickle('sigame/temp/sim_FUV/z'+'{:.2f}'.format(zred)+'_'+galname+'_sim1.star')