cumulative_distributions.py 文件源码

def calc_cumulative_dist(momlim=0.3,region_list=['L1688','NGC1333','B18','OrionA'],
                         file_extension='DR1_rebase3'):
    projDir = '/media/DATAPART/projects/GAS/testing/'
    min_bin = 20.5 # log 10 N(H2) - set by B18
    max_bin = 26. # log 10 N(H2) - set by Orion A
    bin_size = 0.05
    nbins = np.int((max_bin - min_bin)/bin_size)
    data_bin_vals = [min_bin + x * bin_size for x in range(nbins+1)]
    # Loop over regions
    for region_i in range(len(region_list)):
        region = region_list[region_i]
        nh3ImFits  = projDir + '{0}/{0}_NH3_11_{1}_mom0_QA_trim.fits'.format(region,file_extension)
        herColFits = 'nh2_regridded/{0}_NH2_regrid.fits'.format(region)
        nh3_fits = fits.open(nh3ImFits)
        nh2_regrid_fits = fits.open(herColFits)
        h2_data = np.log10(nh2_regrid_fits[0].data)
        nh3_data = nh3_fits[0].data
        h2_mean_array = np.zeros(nbins)
        h2_std_array = np.zeros(nbins)
        nh3_frac_array = np.zeros(nbins)
        for bin_i in range(nbins-1):
            bin_h2_indices = np.where(np.logical_and(h2_data >= data_bin_vals[bin_i],
                                                     h2_data < data_bin_vals[bin_i+1]))
            bin_h2_data = h2_data[bin_h2_indices]
            bin_nh3_data = nh3_data[bin_h2_indices]
            if np.count_nonzero(bin_nh3_data) != 0:
                frac_above_mom = np.count_nonzero(bin_nh3_data > momlim)/(1.0*np.count_nonzero(bin_nh3_data))
            else:
                frac_above_mom = 0
            bin_mean_h2 = np.nanmean(bin_h2_data)
            bin_std_h2 = np.nanstd(bin_h2_data)
            h2_mean_array[bin_i] = bin_mean_h2
            nh3_frac_array[bin_i] = frac_above_mom
            h2_std_array[bin_i] = bin_std_h2
            # Write out bins
            np.savetxt('cumulative/{0}_cumulative.txt'.format(region),
                       np.transpose([h2_mean_array,nh3_frac_array,h2_std_array]))