import numpy as np import matplotlib.pyplot as plt import sys import csv import warnings warnings.simplefilter('ignore') """ Usage: [1] from gasmass import * [2] gasmass? (for help) [3] gasmass(f13=2,f18=1,d=100,Jup=3,fcont=100) """ def gasmass(f13=None, f18=None, e13=None, e18=None, d=None, Jup=None, fcont=None, pdf=False): """ Plot C18O vs 13CO flux, compare with the WB grid, show histogram of gas mass values, and print out the mean mass optionally calculate dust mass from continuum and derive gas/dust ratio f13 = 13CO flux in Jy km/s f18 = C18O flux in Jy km/s d = distance to source in pc Jup = upper J level of transition [1,2,3] e13 = 13CO error in Jy km/s (optional) e18 = C18O error in Jy km/s (optional) fcont = continuum flux in mJy [optional] pdf = True for hardcopy (optional) """ bad = 0 if f13 == None: print "*** Error: Must enter 13CO flux, f13" bad = 1 if f18 == None: print "*** Error: Must enter C18O flux, f18" bad = 1 if d == None: print "*** Error: Must enter distance, d" bad = 1 if Jup == None: print "*** Error: Must enter upper J level, Jup" bad = 1 if bad == 1: sys.exit() if e13 == None or e13 == 0: print "Warning: No 13CO error given (e13), setting to 20% of 13CO flux" e13 = 0.2*f13 if e18 == None or e18 == 0: print "Warning: No C18O error given (e18), setting to 20% of C18O flux" e18 = 0.2*f18 if Jup < 2 or Jup > 3: sys.exit("*** Error: Only implemented for Jup = 1, 2, 3") if Jup == 1: xline = 'f_1-0_13co' xlab = "F($^{13}$CO 1-0) [Jy km/s]" yline1 = 'f_1-0_c18o' yline2 = 'f_1-0_c18o_low' ylab = "F(C$^{18}$O 1-0) [Jy km/s]" nu = 110.0 if Jup == 2: xline = 'f_2-1_13co' xlab = "F($^{13}$CO 2-1) [Jy km/s]" yline1 = 'f_2-1_c18o' yline2 = 'f_2-1_c18o_low' ylab = "F(C$^{18}$O 2-1) [Jy km/s]" nu = 220.0 if Jup == 3: xline = 'f_3-2_13co' xlab = "F($^{13}$CO 3-2) [Jy km/s]" yline1 = 'f_3-2_c18o' yline2 = 'f_3-2_c18o_low' ylab = "F(C$^{18}$O 3-2) [Jy km/s]" nu = 330.0 # read in WB grid and scale fluxes # from model calculation at 140 pc to user input header,data = readcsv('gasgrid.csv') scale = (140.0/d)**2 xplot = scale * data[header.index(xline),:] yplot1 = scale * data[header.index(yline1),:] yplot2 = scale * data[header.index(yline2),:] if f13 > 0: xlim = (0.1*f13, 10*f13) else: xlim = (-0.3*f13, -30*f13) if f18 > 0: ylim = (0.1*f18, 10*f18) else: ylim = (-0.3*f18, -30*f18) mgas = data[1,:] mgas_values = np.sort(list(set(mgas))) mgas_list = ['1e-1','3e-2','1e-2', '3e-3', '1e-3', '3e-4', '1e-4'] mgas_col = ['orange','red','green','mediumseagreen','c','b','k'] # plot results xsize = 16.0 ysize = 10.0 # definitions for the axes left = 0.1 right = 0.95 bottom = 0.15 height = 0.75 dpx = 0.08 width1 = height * ysize/xsize width2 = right - (left + width1 + dpx) pos1 = [left, bottom, width1, height] pos2 = [left+width1+dpx, bottom, width2, height] # define axes plt.figure(figsize=(xsize,ysize)) axScatter = plt.axes(pos1) axHist = plt.axes(pos2) # logarithmic axes axScatter.set_xlim(xlim) axScatter.set_xscale('log') axScatter.set_ylim(ylim) axScatter.set_yscale('log') axScatter.set_xlabel(xlab, fontsize=18, weight='ultralight') axScatter.set_ylabel(ylab, fontsize=18, weight='ultralight') axScatter.tick_params(length=8, width=1.5, which='major') axScatter.tick_params(length=4, width=1.5, which='minor') [i.set_linewidth(1.5) for i in axScatter.spines.itervalues()] for tick in (axScatter.get_xticklabels() + axScatter.get_yticklabels()): tick.set_fontsize(14) tick.set_weight('bold') i = 0 for m in mgas_list: indx = (mgas == float(m)) x = xplot[indx] y1 = yplot1[indx] y2 = yplot2[indx] c = mgas_col[i] axScatter.plot(x, y1, 'o', alpha=0.5, mec='none', color=c, ms=4, rasterized=True) axScatter.plot(x, y2, 'o', alpha=0.5, fillstyle='none', color=c, ms=4, rasterized=True) axScatter.scatter(x[0]+9999, y1[0]+9999, alpha=1.0, edgecolor='none', color=c, label='$'+m.replace('e-','\\times 10^{-')+'} M_\odot$') i += 1 if (f13 > 0 and f18 > 0): axScatter.plot([f13,f13], [f18-e18,f18+e18], 'k-_', lw=3, mew=2) axScatter.plot([f13-e13,f13+e13], [f18,f18], 'k-|', lw=3, mew=2) axScatter.plot(f13, f18, 'ko', ms=10) ifit13 = (abs(f13-xplot) < e13) ifit18 = (abs(f18-yplot1) < e18) + (abs(f18-yplot2) < e18) if (f13 > 0 and f18 < 0): f18 = -3*f18 axScatter.plot([f13-e13,f13+e13], [f18,f18], 'k-|', lw=2, mew=2) axScatter.plot([f13,f13], [f18,0.7*f18], 'k-', lw=2) axScatter.plot([f13,f13], [0.7*f18,0.7*f18], 'kv', lw=2) axScatter.plot(f13, f18, 'ko', ms=10) ifit13 = (abs(f13-xplot) < e13) ifit18 = (yplot1 < f18) + (yplot2 < f18) if (f13 < 0 and f18 < 0): f13 = -3*f13 f18 = -3*f18 axScatter.plot([f13,f13], [f18,0.7*f18], 'k-', lw=3) axScatter.plot([f13,f13], [0.7*f18,0.7*f18], 'kv', lw=3) axScatter.plot([f13,0.7*f13], [f18,f18], 'k-', lw=3) axScatter.plot([0.7*f13,0.7*f13], [f18,f18], 'k<', lw=3) axScatter.plot(f13, f18, 'ko', ms=10) ifit13 = (xplot < f13) ifit18 = (yplot1 < f18) + (yplot2 < f18) axScatter.legend(loc='upper left', scatterpoints=1, markerscale=2.0, borderaxespad=1.2, title='$M_{gas}$', fontsize='medium') # plot histogram of mass values ifit_both = ifit13 & ifit18 nfit13 = np.count_nonzero(ifit13) nfit_both = np.count_nonzero(ifit_both) if nfit13 == 0: sys.exit("*** no solutions found in grid!") mgas_fit13 = data[header.index("M_gas"),ifit13] mgas_fit_both = data[header.index("M_gas"),ifit_both] log_mgas_bins = np.log10([1e-4,3e-4,1e-3,3e-3,1e-2,3e-2,1e-1,3e-1])-0.2386 xlim = [-4.5,-0.5] xticks = [-4.0,-3.5,-3.0,-2.5,-2.0,-1.5,-1.0] xlabel = '$log_{10}(M_{gas})$' hist13, bins = np.histogram(np.log10(mgas_fit13), bins=log_mgas_bins) hist_both, bins = np.histogram(np.log10(mgas_fit_both), bins=log_mgas_bins) norm13 = float(np.sum(hist13)) hist13 = hist13/norm13 hist_both = hist_both/norm13 center = (bins[:-1]+bins[1:])/2 log_mgas_median13 = center[hist13.argmax()] log_mgas_mean13 = np.sum(center*hist13) log_mgas_median_both = center[hist_both.argmax()] log_mgas_mean_both = np.sum(center*hist_both)/np.sum(hist_both) #print "Median gas mass = {0:.2e} Msun".format(10**log_mgas_median_both) #print "Mean gas mass = {0:.2e} Msun".format(10**log_mgas_mean_both) Mgas = 10**log_mgas_mean_both axHist.annotate("$Mgas = {0:.1e} M_\odot$".format(Mgas), xy=(0.07,0.94), xycoords='axes fraction', fontsize=12) col1 = 'c' alp1 = 0.7 col2 = 'b' alp2 = 0.7 width = 0.95*(bins[1]-bins[0]) axHist.bar(center, hist13, align = 'center', width = width, color=col1, alpha=alp1, linewidth=0) axHist.bar(center, hist_both, align = 'center', width = width, color=col2, alpha=alp2, linewidth=0) axHist.set_xlabel(xlabel, fontsize=20, weight='semibold', labelpad=10) axHist.set_ylabel('$Normalized$ $fraction$', fontsize=20, weight='normal', labelpad=10) axHist.set_xlim(xlim) axHist.set_xticks(xticks) axHist.tick_params(length=8, width=1.5) [i.set_linewidth(1.5) for i in axHist.spines.itervalues()] for item in (axHist.get_xticklabels() + axHist.get_yticklabels()): item.set_fontsize(14) item.set_weight('medium') # plot legend with the bold (not transparent) colors axHist.bar(center, 0*hist13, color=col1, linewidth=0, label='13CO') axHist.bar(center, 0*hist_both, color=col2, linewidth=0, label='13CO/C18O') axHist.legend(loc=1, borderaxespad=1.5) # if a continuum flux is given, calculate the dust mass # kappa is the Beckwith dust opacity (not including gas!) # assume dust temperature of 20 K for this simple calculation if fcont != None and fcont > 0: kappa_nu = 10.0 * (nu / (1000)) # cm2/g Tdust = 20.0 # K B_nu = 1.4745e-20 * nu**3 / (np.exp(4.80e-2*nu/Tdust) - 1) # cgs Mdust = 4.798e-23 * fcont * d**2 / (kappa_nu * B_nu) # Msun #print "Dust mass = {0:.1e Msun".format(Mdust) #print "Gas/dust ratio = {0:.1f}".format(Mgas/Mdust) axHist.annotate("$Mdust = {0:.1e} M_\odot$".format(Mdust), xy=(0.07,0.91), xycoords='axes fraction', fontsize=12) axHist.annotate("$Gas/dust$ $ratio = {0}$".format(int(np.round(Mgas/Mdust))), xy=(0.07,0.88), xycoords='axes fraction', fontsize=12) if pdf: plt.savefig('gasmass.pdf') plt.show() def readcsv(filename): # get number of rows d = open(filename, "rU") r = csv.reader(d) nrow = -1 for row in r: nrow = nrow + 1 d.close() # now read the data into a header and data array d = open(filename, "rU") r = csv.reader(d) rownum = 0 for row in r: if rownum == 0: header = row ncol = len(header) data = np.zeros((ncol,nrow)) else: colnum = 0 for col in row: data[colnum,rownum-1] = float(col) colnum += 1 rownum += 1 d.close() data = data[:,0:rownum-1] return header,data