diff --git a/analyse_spectra.py b/analyse_spectra.py index 06ccffe..5345bd4 100755 --- a/analyse_spectra.py +++ b/analyse_spectra.py @@ -137,6 +137,7 @@ def read_spectra(): em_std = np.zeros(len(results_array)) sigma_tau = np.zeros(len(results_array)) i = 0 + has_emission = False for row in results_array: opacity = row['opacity'] velocities[i] = row['velocity'] / 1000.0 @@ -146,6 +147,7 @@ def read_spectra(): sigma_tau[i] = row['sigma_tau'] if 'em_mean' in results_array.dtype.names: em_temps[i] = row['em_mean'] + has_emission = True if 'em_std' in results_array.dtype.names: em_std[i] = row['em_std'] i += 1 @@ -180,6 +182,7 @@ def read_spectra(): spectrum.em_temps = em_temps spectrum.em_std = em_std spectrum.sigma_tau = sigma_tau + spectrum.has_emission = has_emission spectra.append(spectrum) return spectra @@ -429,7 +432,8 @@ def plot_lv_image(x, y, c, filename): v_max * 1000) # Add an outline of the emission from GASS III - plt.contour(np.log10(np.flipud(gass_subset)), 1, cmap='Pastel2') + contour_set = plt.contour(np.log10(np.flipud(gass_subset)), 1, cmap='Pastel2') + #plt.clabel(contour_set) # Set the axis ticks and scales x_step = int(20 * (gass_subset.shape[1] / (l_max - l_min))) @@ -504,6 +508,15 @@ def is_resolved(day, field_name, island_id, source_id, beam_area, islands): return isle['area'] > isle['beam_area'] +def set_field_metadata(field, ucd, unit, description): + if ucd: + field.ucd = ucd + if unit: + field.unit = unit + if description: + field.description = description + + def output_spectra_catalogue(spectra, isle_day_map): """ Output the list of spectrum stats to a VOTable file magmo-spectra.vot @@ -531,10 +544,13 @@ def output_spectra_catalogue(spectra, isle_day_map): continuum_temp = np.zeros(rows) max_s_max_n = np.zeros(rows) max_em_std = np.zeros(rows) + max_emission = np.zeros(rows) + min_emission = np.zeros(rows) rating = np.empty(rows, dtype=object) used = np.empty(rows, dtype=bool) resolved = np.empty(rows, dtype=bool) duplicate = np.empty(rows, dtype=bool) + has_emission = np.empty(rows, dtype=bool) filenames = np.empty(rows, dtype=object) local_paths = np.empty(rows, dtype=object) local_emission_paths = np.empty(rows, dtype=object) @@ -557,6 +573,9 @@ def output_spectra_catalogue(spectra, isle_day_map): min_velocity[i] = np.min(spectrum.velocity) max_velocity[i] = np.max(spectrum.velocity) max_em_std[i] = np.max(spectrum.em_std) + if spectrum.has_emission: + max_emission[i] = np.max(spectrum.em_temps) + min_emission[i] = np.min(spectrum.em_temps) opacity_range[i] = spectrum.opacity_range max_s_max_n[i] = spectrum.max_s_max_n @@ -568,9 +587,9 @@ def output_spectra_catalogue(spectra, isle_day_map): isle_day_map[spectrum.day] if spectrum.day in isle_day_map else None) duplicate[i] = spectrum.duplicate + has_emission[i] = spectrum.has_emission used[i] = not spectrum.low_sn - prefix = 'day' + spectrum.day + '/' + spectrum.field_name + \ - "_src" + spectrum.src_id + prefix = 'day' + spectrum.day + '/' + spectrum.field_name + "_src" + spectrum.src_id filenames[i] = prefix + "_plot.png" em_filename = prefix + "_emission.png" local_paths[i] = base_path + '/' + filenames[i] @@ -579,17 +598,23 @@ def output_spectra_catalogue(spectra, isle_day_map): spectra_table = Table( [ids, days, fields, sources, longitudes, latitudes, eq_ras, eq_decs, max_flux, min_opacity, - max_opacity, rms_opacity, min_velocity, max_velocity, used, continuum_temp, - opacity_range, max_s_max_n, continuum_sd, max_em_std, rating, resolved, duplicate, + max_opacity, rms_opacity, min_emission, max_emission, min_velocity, max_velocity, used, continuum_temp, + opacity_range, max_s_max_n, continuum_sd, max_em_std, rating, resolved, duplicate, has_emission, filenames, local_paths, local_emission_paths], names=['Name', 'Day', 'Field', 'Source', 'Longitude', 'Latitude', 'RA', 'Dec', 'Max_Flux', - 'Min_Opacity', 'Max_Opacity', 'RMS_Opacity', 'Min_Velocity', + 'Min_Opacity', 'Max_Opacity', 'RMS_Opacity', 'Min_Emission', 'Max_Emission', 'Min_Velocity', 'Max_Velocity', 'Used', 'Continuum_Temp', 'Opacity_Range', 'Max_S_Max_N', - 'Continuum_SD', 'max_em_std', 'Rating', 'Resolved', 'Duplicate', + 'Continuum_SD', 'max_em_std', 'Rating', 'Resolved', 'Duplicate', 'Has_Emission', 'Filename', 'Local_Path', 'Local_Emission_Path'], meta={'ID': 'magmo_spectra', 'name': 'MAGMO Spectra ' + str(datetime.date.today())}) votable = from_table(spectra_table) + table = votable.get_first_table() + set_field_metadata(table.get_field_by_id('Min_Emission'), 'stat.min', 'K', + 'Minimum average emission') + set_field_metadata(table.get_field_by_id('Max_Emission'), 'stat.max', 'K', + 'Maximum average emission') + filename = "magmo-spectra.vot" writeto(votable, filename) print("Wrote out", i, "spectra to", filename) @@ -890,7 +915,8 @@ def filter_duplicate_sources(spectra, field_map): # match2.rating, # match1.continuum_sd, # match2.continuum_sd)) - if match1.rating > match2.rating or match1.continuum_sd > match2.continuum_sd: + if match1.rating > match2.rating or (not match1.has_emission and match2.has_emission) or \ + match1.continuum_sd > match2.continuum_sd: match1.duplicate = True print("Flagging {} rating {} as duplicate. Dist {:.2f}".format(match1.name, match1.rating, dist_2d[i].to(u.arcsec))) diff --git a/decompose.py b/decompose.py index c07365a..c0b48e2 100755 --- a/decompose.py +++ b/decompose.py @@ -196,7 +196,7 @@ def decompose(spectra, out_filename, alpha1, alpha2, snr_thresh, data_filename): (len(spectra), time.strftime('%Y-%m-%d %H:%M:%S', time.localtime(end)), (end - end_read))) -def output_component_catalogue(spectra, data, data_decomposed): +def output_component_catalogue(spectra, data, data_decomposed, folder): names = [] comp_names = [] days = [] @@ -212,7 +212,6 @@ def output_component_catalogue(spectra, data, data_decomposed): fwhms_fit_errs = [] means_fit_errs = [] - num_no_comps = {} for i in range(len(data_decomposed['fwhms_fit'])): @@ -252,6 +251,7 @@ def output_component_catalogue(spectra, data, data_decomposed): else: rating = spectrum['Rating'] num_no_comps[rating] = num_no_comps.get(rating, 0) + 1 + print ("Unable to find components for ") temp_table = Table( [comp_names, names, days, field_names, sources, longitudes, latitudes, amps, fwhms, means, best_fit_rchi2s, @@ -263,6 +263,8 @@ def output_component_catalogue(spectra, data, data_decomposed): votable = from_table(temp_table) filename = "magmo-components.vot" writeto(votable, filename) + filename = folder + "/magmo-components.vot" + writeto(votable, filename) total_nc = 0 for rating, count in num_no_comps.items(): @@ -351,19 +353,38 @@ def plot_single_spectrum(ax, velo, opacity, fit_amps, fit_fwhms, fit_means, name return residual +def find_bounds(velo, fit_fwhms, fit_means): + buff = 25 # km/s + + if len(fit_fwhms) == 0: + return 0, len(velo)-1 + + means = np.array(fit_means) + fwhms = np.array(fit_fwhms) + lowest_vals = means - fwhms + highest_vals = means + fwhms + low_velo = np.min(lowest_vals) - buff + high_velo = np.max(highest_vals) + buff + low_idx = (np.abs(velo - low_velo)).argmin() + high_idx = (np.abs(velo - high_velo)).argmin() + return low_idx, high_idx + + def plot_spectrum(velo, opacity, fit_amps, fit_fwhms, fit_means, name, filename, formats): fig = plt.figure(figsize=(8, 6)) gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1]) ax = fig.add_subplot(gs[0]) y = convert_from_ratio(opacity) - residual = plot_single_spectrum(ax, velo, y, fit_amps, fit_fwhms, fit_means, name) + min_bound, max_bound = find_bounds(velo, fit_fwhms, fit_means) + residual = plot_single_spectrum(ax, velo[min_bound:max_bound], y[min_bound:max_bound], + fit_amps, fit_fwhms, fit_means, name) residual_rms = np.sqrt(np.mean(np.square(residual))) ax.set_ylabel('$e^{-\\tau}$') # Residual plot ax = fig.add_subplot(gs[1]) - ax.plot(velo, residual, 'or', markerfacecolor='None', markersize=2, markeredgecolor='blue') + ax.plot(velo[min_bound:max_bound], residual, 'or', markerfacecolor='None', markersize=2, markeredgecolor='blue') ax.grid() plt.xlabel('LSR Velocity (km/s) n=%d rms=%.4f' % (len(fit_amps), residual_rms)) @@ -483,7 +504,7 @@ def output_decomposition(spectra, out_filename, folder, data_filename, alpha1, a plt.savefig(folder + "/magmo-decomp.pdf") plt.close() - output_component_catalogue(spectra, data, data_decomposed) + output_component_catalogue(spectra, data, data_decomposed, folder) output_decomposition_catalogue(folder, spectra, data, data_decomposed, alpha1, alpha2) plot_spectra(spectra, data, data_decomposed, alpha1, alpha2, folder=folder) diff --git a/examine_gas.py b/examine_gas.py index e8b5a9a..2a03485 100755 --- a/examine_gas.py +++ b/examine_gas.py @@ -12,6 +12,7 @@ import argparse import csv import datetime +import glob import math import os import re @@ -20,16 +21,20 @@ from astropy.coordinates import SkyCoord, matching from astropy.io.votable import parse, from_table, writeto -from astropy.table import Table, Column, hstack +from astropy.table import Table, Column, hstack, vstack from astropy.io import ascii import astropy.units as u import matplotlib import matplotlib.pyplot as plt +from matplotlib.ticker import LogFormatter import numpy as np +import numpy.ma as ma from scipy import interpolate import scipy.stats as stats import seaborn as sns +import magmo + class Gas(object): def __init__(self, day, field, src): @@ -73,6 +78,7 @@ def get_spectra_key(day, field, source): def read_spectra(filename): spectra_array = read_votable_results(filename) + spectra_array = spectra_array[spectra_array['Duplicate'] == False] spectra_map = {} for row in spectra_array: key = get_spectra_key(row['Day'], row['Field'], row['Source']) @@ -114,45 +120,64 @@ def get_emission_filename(day, field, source): def get_temp(emission, comp_vel): velocities = emission['velocity'] - temp = 0 for i in range(0, len(velocities)): if velocities[i] >= comp_vel: temp = emission['em_mean'][i] - return temp, velocities[i] - return 0, 0 + temp_err = emission['em_std'][i] + return temp, temp_err, velocities[i] + return 0, 0, 0 def analyse_components(components, spectra_map, mmb_map): all_gas = [] for component in components: + + comp_vel = component['Mean'] + comp_width = component['FWHM'] + comp_amp = component['Amplitude'] + optical_depth = 1 - comp_amp # From 1-e^-tau to e^-tau + + gas = Gas(component['Day'], component['Field'], component['Source']) + gas.comp_vel = comp_vel + gas.comp_width = math.fabs(comp_width) + gas.optical_depth = optical_depth + gas.longitude = component['Longitude'] + gas.latitude = component['Latitude'] + gas.tau = -1 * np.log(np.maximum(optical_depth, 1e-16)) + gas.t_kmax = 21.866 * comp_width ** 2 + print("Gas T_kmax=", gas.t_kmax, "comp_width=", comp_width) + gas.t_off = None + gas.t_s = None + gas.em_vel = None + gas.name = component['Comp_Name'] + gas.spectra_name = component['Spectra_Name'] + + spectrum = spectra_map[get_spectra_key(component['Day'], component['Field'], component['Source'])] + gas.rating = spectrum['Rating'] + gas.continuum_sd = spectrum['Continuum_SD'] + + loc = SkyCoord(gas.longitude, gas.latitude, frame='galactic', unit="deg") + gas.loc = loc + gas.ra = loc.icrs.ra.degree + gas.dec = loc.icrs.dec.degree + + maser = mmb_map.get(component['Field']) + if maser is None: + print("unable to find maser for " + component['Field']) + else: + gas.maser_vel_low = maser['VL'] + gas.maser_vel_high = maser['VH'] + gas.maser_loc = SkyCoord(maser['RAJ2000'], maser['DEJ2000'], frame='fk5', unit="deg") + # Load emission data + t_off = 0 emission_filename = get_emission_filename(component['Day'], component['Field'], component['Source']) if os.path.exists(emission_filename): emission = read_emission(emission_filename) + t_off, t_off_err, em_vel = get_temp(emission, comp_vel * 1000) - comp_vel = component['Mean'] - comp_width = component['FWHM'] - comp_amp = component['Amplitude'] - t_off, em_vel = get_temp(emission, comp_vel * 1000) - optical_depth = 1 - comp_amp # From 1-e^-tau to e^-tau - - gas = Gas(component['Day'], component['Field'], component['Source']) - gas.comp_vel = comp_vel - gas.comp_width = math.fabs(comp_width) - gas.optical_depth = optical_depth gas.em_vel = em_vel - gas.longitude = component['Longitude'] - gas.latitude = component['Latitude'] - gas.tau = -1 * np.log(np.maximum(optical_depth, 1e-16)) - gas.t_kmax = 21.866 * comp_width**2 - print ("Gas T_kmax=", gas.t_kmax, "comp_width=", comp_width) - gas.t_off = None - gas.t_s = None - gas.name = component['Comp_Name'] - gas.spectra_name = component['Spectra_Name'] - - spectrum = spectra_map[get_spectra_key(component['Day'], component['Field'], component['Source'])] - gas.rating = spectrum['Rating'] + # Validate the component velocity if not spectrum['Min_Velocity'] <= comp_vel <= spectrum['Max_Velocity']: @@ -160,32 +185,23 @@ def analyse_components(components, spectra_map, mmb_map): spectrum['Min_Velocity'], spectrum['Max_Velocity'], comp_vel)) continue - loc = SkyCoord(gas.longitude, gas.latitude, frame='galactic', unit="deg") - gas.loc = loc - gas.ra = loc.icrs.ra.degree - gas.dec = loc.icrs.dec.degree - - maser = mmb_map.get(component['Field']) - if maser is None: - print("unable to find maser for " + component['Field']) - else: - gas.maser_vel_low = maser['VL'] - gas.maser_vel_high = maser['VH'] - gas.maser_loc = SkyCoord(maser['RAJ2000'], maser['DEJ2000'], frame='fk5', unit="deg") - - all_gas.append(gas) - - if t_off > 0 and optical_depth < 0.9: - # Calculate spin temperature and column density - #denominator = math.min(1-comp_amp) - t_s = t_off / (1-np.exp(-gas.tau)) # (tb / (1-e^-tau) - - # Record - gas.t_off = t_off - gas.t_s = t_s - # component['t_s '] = t_s - print("Component %s at velocity %.4f has t_s %.3f (%.3f/%.3f)" % ( - gas.name, comp_vel, t_s, t_off, comp_amp)) + all_gas.append(gas) + + if t_off > 0 and optical_depth < 0.9: + # Calculate spin temperature and column density + #denominator = math.min(1-comp_amp) + t_s = t_off / (1-np.exp(-gas.tau)) # (tb / (1-e^-tau) + delta_t_s = t_s * ((t_off_err ** 2 / t_off ** 2) + + (gas.continuum_sd ** 2 / (1 - np.exp(-gas.tau)) ** 2))**0.5 + + # Record + gas.t_off = t_off + gas.delta_t_off = t_off_err + gas.t_s = t_s + gas.delta_t_s = delta_t_s + # component['t_s '] = t_s + print("Component %s at velocity %.4f has t_s %.3f (%.3f/%.3f)" % ( + gas.name, comp_vel, t_s, t_off, comp_amp)) return all_gas @@ -217,13 +233,16 @@ def output_gas_catalogue(all_gas): ras = [] decs = [] velocities = np.zeros(num_gas) - em_velocities = np.zeros(num_gas) + em_velocities = np.ma.array(np.zeros(num_gas)) optical_depths = np.zeros(num_gas) comp_widths = np.zeros(num_gas) temps_off = np.ma.array(np.zeros(num_gas)) + delta_temps_off = np.ma.array(np.zeros(num_gas)) temps_spin = np.ma.array(np.zeros(num_gas)) + delta_temps_spin = np.ma.array(np.zeros(num_gas)) temps_kmax = np.zeros(num_gas) tau = np.zeros(num_gas) + continuum_sd = np.zeros(num_gas) maser_region = np.empty(num_gas, dtype=bool) ratings = np.empty(num_gas, dtype=object) filenames = np.empty(num_gas, dtype=object) @@ -244,21 +263,26 @@ def output_gas_catalogue(all_gas): ras.append(gas.ra) decs.append(gas.dec) velocities[i] = gas.comp_vel - em_velocities[i] = gas.em_vel / 1000 + em_velocities[i] = gas.em_vel / 1000 if gas.em_vel else np.ma.masked optical_depths[i] = gas.optical_depth comp_widths[i] = gas.comp_width temps_kmax[i] = gas.t_kmax if gas.t_off is None: temps_off[i] = np.ma.masked + delta_temps_off[i] = np.ma.masked else: temps_off[i] = gas.t_off + delta_temps_off[i] = gas.delta_t_off if gas.t_s is None: temps_spin[i] = np.ma.masked + delta_temps_spin[i] = np.ma.masked else: temps_spin[i] = gas.t_s + delta_temps_spin[i] = gas.delta_t_s tau[i] = gas.tau maser_region[i] = is_gas_near_maser(gas) ratings[i] = gas.rating + continuum_sd[i] = gas.continuum_sd # Need to read in spectra to get rating and include it in the catalogue and # link to the fit preview: e.g. plots/A/012.909-0.260_19_src4-1_fit prefix = 'day' + str(gas.day) + '/' + gas.field + \ @@ -272,25 +296,28 @@ def output_gas_catalogue(all_gas): # bulk calc fields vel_diff = np.abs(velocities - em_velocities) - equiv_width = np.abs((optical_depths) * comp_widths) - #tau = -np.log(optical_depths) - column_density = tau * np.abs(comp_widths) * temps_spin * 1.83E18 + equiv_width = np.abs((1-optical_depths) * comp_widths) fwhm = np.abs(comp_widths) + column_density = tau * fwhm * temps_spin * 1.823E18 * 1.064 sigma = fwhm / (2 * math.sqrt(2 * math.log(2))) + mach_num = np.sqrt(4.2*((temps_kmax/temps_spin)-1)) + n_wnm = ((column_density * temps_spin) / 50 ) - column_density temp_table = Table( [names, days, field_names, sources, velocities, em_velocities, optical_depths, temps_off, temps_spin, temps_kmax, longitudes, latitudes, ras, decs, fwhm, sigma, vel_diff, equiv_width, tau, maser_region, column_density, - ratings, filenames, local_paths, local_emission_paths, local_spectra_paths], + mach_num, n_wnm, ratings, delta_temps_spin, delta_temps_off, continuum_sd, + filenames, local_paths, local_emission_paths, local_spectra_paths], names=['Comp_Name', 'Day', 'Field', 'Source', 'Velocity', 'em_velocity', 'optical_depth', 'temp_off', - 'temp_spin', 'temp_kmax', 'longitude', 'latitude', 'ra', 'dec', 'fwhm', 'sigma', 'vel_diff', 'equiv_width', 'tau', - 'near_maser', 'column_density', 'Rating', 'Filename', 'Local_Path', 'Local_Emission_Path', - 'Local_Spectrum_Path'], + 'temp_spin', 'temp_kmax', 'longitude', 'latitude', 'ra', 'dec', 'fwhm', 'sigma', 'vel_diff', + 'equiv_width', 'tau', 'near_maser', 'column_density', 'mach', 'n_wnm', 'Rating', + 'delta_temp_spin', 'delta_temp_off', 'delta_optical_depth', + 'Filename', 'Local_Path', 'Local_Emission_Path', 'Local_Spectrum_Path'], meta={'ID': 'magmo_gas', 'name': 'MAGMO Gas ' + str(datetime.date.today())}) votable = from_table(temp_table) table = votable.get_first_table() - set_field_metadata(table.get_field_by_id('longitude'), 'pos.galactic.long', 'deg', + set_field_metadata(table.get_field_by_id('longitude'), 'pos.galactic.lon', 'deg', 'Galactic longitude of the background source') set_field_metadata(table.get_field_by_id('latitude'), 'pos.galactic.lat', 'deg', 'Galactic latitude of the background source') @@ -306,58 +333,202 @@ def output_gas_catalogue(all_gas): 'The excitation or spin temperature of the gas') set_field_metadata(table.get_field_by_id('optical_depth'), '', '', 'The peak optical depth of the component ($e^(-\\tau)$)') + set_field_metadata(table.get_field_by_id('column_density'), '', 'cm-2', + 'The density of Cold HI gas in the Gaussian component') + set_field_metadata(table.get_field_by_id('n_wnm'), '', 'cm-2', + 'The density of Warm HI gas in the Gaussian component') + set_field_metadata(table.get_field_by_id('mach'), '', '', + 'The turbulent mach number of the gas in the Gaussian component') filename = "magmo-gas.vot" writeto(votable, filename) return table -def plot_equiv_width_lv(gas_table): - values = gas_table.array - cm = plt.cm.get_cmap('RdYlBu_r') - sc = plt.scatter(values['longitude'], values['Velocity'], c=values['equiv_width'], s=35, cmap=cm) - cb = plt.colorbar(sc, norm=matplotlib.colors.LogNorm()) +def plot_equiv_width_lv(values, ax, key, label, max_clip=None): + scale_vals = values[key] + if max_clip: + scale_vals = np.clip(scale_vals, 0, max_clip) + base = np.array([values['longitude'], values['Velocity'], scale_vals], + dtype=[('l', float), ('b', float), ('scale', float)]) + + sample = base # np.sort(base, axis=0, order='ew') + #vmax = np.max(sample[2,:]) + #cm = plt.cm.get_cmap('RdYlBu_r') + cm = plt.cm.get_cmap('viridis') + x = sample[0,:] + y = sample[1,:] + z = sample[2,:] + idx = z.argsort() + x, y, z = x[idx], y[idx], z[idx] + + # equiv width plot needs log, min 3 + sc = ax.scatter(x, y, c=z, s=25, cmap=cm, norm=matplotlib.colors.LogNorm(vmin=3)) + #sc = ax.scatter(sample[0,:], sample[1,:], c=sample[2,:], s=25, cmap=cm, norm=matplotlib.colors.PowerNorm(0.5, vmin=3)) + #cb = plt.colorbar(sc, norm=matplotlib.colors.LogNorm()) + formatter = LogFormatter(10, labelOnlyBase=False) + ticks = [0, 5, 10, 20, 30, 40, 60, 80] + cb = plt.colorbar(sc, ticks=ticks, format=formatter) + + #cb.set_ticklabels([0, 10, 20, 30, 40 , 50 ,60 , 70, 80]) + #cb.ax.get_yaxis().set_ticks([]) + #for j in range(0, 80, 10): + # cb.ax.text(.5, (2 * j + 1) / 8.0, str(j), ha='center', va='center') - ax = plt.gca() ax.set_xlim(values['longitude'].max() + 5, values['longitude'].min() - 5) - plt.title("Equivalent Width of Fitted Gas Components") - plt.xlabel('Galactic longitude (deg)') - plt.ylabel('LSR Velocity (km/s)') - cb.set_label('Equivalent Width (km/s)') + ax.set_xlabel('Galactic longitude (deg)') + ax.set_ylabel('LSR Velocity (km/s)') + cb.set_label(label) - filename = 'magmo-equiv-width-lv.pdf' - # plt.show() - plt.savefig(filename) - plt.close() return None +def plot_spin_temp_lv(values, ax, key, label, max_clip=None): + scale_vals = values[key] + dataset = values[~ma.getmask(scale_vals)] + if max_clip: + scale_vals = np.clip(dataset[key], 0, max_clip) + base = np.array([dataset['longitude'], dataset['Velocity'], scale_vals], + dtype=[('l', float), ('b', float), ('scale', float)]) + + sample = base + #vmax = np.max(sample[2,:]) + #cm = plt.cm.get_cmap('RdYlBu_r') + cm = plt.cm.get_cmap('viridis') + + x = sample[0,:] + y = sample[1,:] + z = sample[2,:] + idx = z.argsort() + x, y, z = x[idx], y[idx], z[idx] + + # spin temp plot needs log, min 3 + sc = ax.scatter(x, y, c=z, s=25, cmap=cm, norm=matplotlib.colors.PowerNorm(0.5, vmin=3)) + formatter = LogFormatter(10, labelOnlyBase=False) + ticks = [0, 5, 10, 20, 40, 70, 100, 150, 200, 300] + cb = plt.colorbar(sc, ticks=ticks, format=formatter) + + ax.set_xlim(values['longitude'].max() + 5, values['longitude'].min() - 5) + + ax.set_xlabel('Galactic longitude (deg)') + ax.set_ylabel('LSR Velocity (km/s)') + cb.set_label(label) + + return None + + +def plot_spectra(spectra): + print("## Plotting spectra stats ##") + + fig = plt.figure(figsize=(4, 3)) + gs = matplotlib.gridspec.GridSpec(2, 2) + labels = ('A', 'B', 'C', 'D') + ratings = 'ABCD' + ind = np.arange(len(ratings)) + + # Baseline noise box chart + ax1 = fig.add_subplot(gs[0, 0]) + noise = [] + for rating in ratings: + noise.append(spectra['Continuum_SD'][spectra['Rating'] == rating].compressed()) + ax1.boxplot(noise, labels=labels, sym='k.') + ax1.axhline(y=1 / 3, linestyle='--', color='darkgray') + ax1.semilogy() + ax1.set_ylabel('$\sigma_{continuum}$') + ax1.set_title('Baseline Noise', fontsize=10) + + ax2 = fig.add_subplot(gs[0, 1]) + maxsn = [] + for rating in ratings: + maxsn.append(spectra['Max_S_Max_N'][spectra['Rating'] == rating].compressed()) + ax2.boxplot(maxsn, labels=labels, sym='k.') + ax2.axhline(y=3, linestyle='--', color='darkgray') + ax2.semilogy() + ax2.set_title('Max Signal:Max Noise', fontsize=10) + + # Quality histogram + ax3 = fig.add_subplot(gs[1, :]) + counts = [] + for rating in ratings: + counts.append(len(spectra['Rating'][spectra['Rating'] == rating])) + ax3.bar(ind, counts, width=1, align='center') + ax3.set_xticks(ind) + ax3.set_xticklabels(labels) + ax3.set_xlim(-0.5, len(ind) - 0.5) + ax3.set_xlabel('Quality Rating') + ax3.set_ylabel('Number') + + gs.update(wspace=0.5, hspace=0.5) + filename = 'magmo-spectra-stats.pdf' + plt.savefig(filename, bbox_inches="tight") + plt.close() + + def plot_maser_comparison(gas_table): print("## Comparing gas near and away from masers ##") #df = pd.DataFrame(np.ma.filled(gas_array)) gas_array = gas_table.array - df = gas_table.to_table().to_pandas() - grid = sns.FacetGrid(df, col="near_maser", margin_titles=True, sharey=False) - bins = np.logspace(0.01, 3, 21) - grid.map(plt.hist, "temp_spin", color="steelblue", lw=0, bins=bins) - plt.savefig('near_maser_temp.pdf') - plt.close() + #df = gas_table.to_table().to_pandas() + #grid = sns.FacetGrid(df, col="near_maser", margin_titles=True, sharey=False) + #bins = np.logspace(0.01, 3, 21) + #grid.map(plt.hist, "temp_spin", color="steelblue", lw=0, bins=bins) + #plt.savefig('near_maser_temp.pdf') + #plt.close() ts = gas_array['temp_spin'] - indexes = [~np.isnan(ts)] + indexes = [(0 < ts)] ts_gas = gas_array[indexes] - - #print(ts[0:10], ts_gas[0:10]) near_maser = ts_gas[ts_gas['near_maser']] away_maser = ts_gas[ts_gas['near_maser'] == False] - print(near_maser[0:10], away_maser[0:10]) - print("Median near:{} v away:{} , sd n:{} v a:{}".format(np.median(near_maser['temp_spin']), + #print(near_maser[0:10], away_maser[0:10]) + + fig = plt.figure(figsize=(7.5, 3)) + gs = matplotlib.gridspec.GridSpec(1, 2) + + # Spin Temperature + ax1 = fig.add_subplot(gs[0, 0]) + away_sample = away_maser['temp_spin'] + bins = np.linspace(0, 450, 19) + hist, edges = build_hist_fraction(away_sample, bins, 450) + ax1.step(edges, hist) + + near_sample = near_maser['temp_spin'] + hist, edges = build_hist_fraction(near_sample, bins, 450) + ax1.step(edges, hist, color='black', ls='--') + + ax1.set_xlabel('Spin Temperature (K)') + ax1.set_ylabel('Fraction of components') + + # Optical Depth + ax2 = fig.add_subplot(gs[0, 1]) + away_sample = away_maser['tau'] + bins = np.linspace(0, 4, 21) + hist, edges = build_hist_fraction(away_sample, bins, 4) + ax2.step(edges, hist) + + near_sample = near_maser['tau'] + hist, edges = build_hist_fraction(near_sample, bins, 4) + ax2.step(edges, hist, color='black', ls='--') + + ax2.set_xlabel('Optical depth $(\\tau)$') + ax2.set_ylabel('Fraction of components') + + gs.update(wspace=0.5, hspace=0.5) + plt.savefig('near_maser_temp.pdf', bbox_inches="tight") + plt.close() + + + print("Median near:{} v away:{} , sd n:{} v a:{}, count n:{} v a:{}".format(np.median(near_maser['temp_spin']), np.median(away_maser['temp_spin']), np.std(near_maser['temp_spin']), - np.std(away_maser['temp_spin']))) + np.std(away_maser['temp_spin']), + len(near_maser['temp_spin']), + len(away_maser['temp_spin']))) statistic, p_value = stats.ks_2samp(np.ma.filled(near_maser['temp_spin']), np.ma.filled(away_maser['temp_spin'])) - print ('Population similarity p_value={}'.format(p_value)) + print ('Population T_S similarity K-s={} p_value={}'.format(statistic, p_value)) + statistic, p_value = stats.ks_2samp(np.ma.filled(near_maser['tau']), np.ma.filled(away_maser['tau'])) + print ('Population tau similarity K-s={} p_value={}'.format(statistic, p_value)) def plot_single_hist(gas_array, field, min, max, label=None, sample=None): @@ -381,23 +552,194 @@ def plot_single_hist(gas_array, field, min, max, label=None, sample=None): outliers, max)) +def plot_observed(gas_array): + fig = plt.figure(figsize=(7.5, 3)) + gs = matplotlib.gridspec.GridSpec(1, 2) + + # Optical depth + ax1 = fig.add_subplot(gs[0, 0]) + sample = gas_array['tau'] + bins = np.linspace(0,4,21) + hist, edges = np.histogram(sample, bins=bins, range=(0, 4)) + outliers = len(sample[sample > 4]) + hist[-1] += outliers + hist = np.append(hist, hist[-1]) + #print ("tau edges",edges) + ax1.step(edges, hist, where='post') # , width=edges[1]) + ax1.set_xlabel('Optical depth $(\\tau)$') + ax1.set_ylabel('Number of components') + + # FWHM + ax2 = fig.add_subplot(gs[0, 1]) + sample = gas_array['fwhm'] + bins = np.linspace(0,60,25) + hist, edges = np.histogram(sample, bins=bins, range=(0, 60)) + outliers = len(sample[sample > 60]) + hist[-1] += outliers + hist = np.append(hist[0], hist) + ax2.step(edges, hist) # , width=edges[1]) + ax2.set_xlabel('FWHM (km/s)') + ax2.set_ylabel('Number of components') + + gs.update(wspace=0.5, hspace=0.5) + filename = 'magmo-hist-observed.pdf' + plt.savefig(filename, bbox_inches="tight") + plt.close() + + +def plot_equiv_width(gas_array): + fig = plt.figure(figsize=(7.5, 3)) + gs = matplotlib.gridspec.GridSpec(1, 2) + + # Histogram - all + ax1 = fig.add_subplot(gs[0, 0]) + sample = gas_array['equiv_width'] + bins = np.linspace(0,40,41) + hist, edges = np.histogram(sample, bins=bins) + outliers = len(sample[sample > 40]) + hist[-1] += outliers + hist = np.append(hist[0], hist) + ax1.step(edges, hist) # , width=edges[1]) + ax1.set_xlabel('Equivalent Width (km/s)') + ax1.set_ylabel('Number of components') + + # Histogram - FWHM < 50 + ax2 = fig.add_subplot(gs[0, 1]) + sample = sample[gas_array['fwhm'] < 50] + hist, edges = np.histogram(sample, bins=bins) + outliers = len(sample[sample > 40]) + hist[-1] += outliers + hist = np.append(hist[0], hist) + ax2.step(edges, hist) # , width=edges[1]) + ax2.set_xlabel('Equivalent Width (km/s)') + ax2.set_ylabel('Number of components') + + gs.update(wspace=0.5, hspace=0.5) + filename = 'magmo-equiv-width.pdf' + plt.savefig(filename, bbox_inches="tight") + plt.close() + + +def plot_derived(gas_array): + fig = plt.figure(figsize=(7.5, 3)) + gs = matplotlib.gridspec.GridSpec(1, 2) + + # Spin Temperature + ax1 = fig.add_subplot(gs[0, 0]) + sample = np.ma.array(gas_array['temp_spin']).compressed() + bins = np.linspace(0,450,19) + hist, edges = np.histogram(sample, bins=bins) + outliers = len(sample[sample > 450]) + hist[-1] += outliers + hist = np.append(hist[0], hist) + ax1.step(edges, hist) # , width=edges[1]) + ax1.set_xlabel('Spin Temperature (K)') + ax1.set_ylabel('Number of components') + + # Column Density + ax2 = fig.add_subplot(gs[0, 1]) + sample = np.log10(np.array(gas_array['column_density'])) + bins = np.linspace(19, 24, 21) + hist, edges = np.histogram(sample, bins=bins) + outliers = len(sample[sample > 24]) + hist[-1] += outliers + hist = np.append(hist[0], hist) + + ## Call out the saturated values + sat_index = gas_array['tau'] > 3 + sat_sample = np.log10(np.array(gas_array['column_density'][sat_index])) + sat_hist, sat_edges = np.histogram(sat_sample, bins=bins) + outliers = len(sat_sample[sat_sample > 24]) + sat_hist[-1] += outliers + sat_hist = np.append(sat_hist[0], sat_hist) + + ax2.step(edges, hist) #, width=edges[1]-edges[0]) + ax2.step(edges, sat_hist, color='r') # , width=edges[1]-edges[0] + label = 'Column Density $\\log_{10}(N_{H}$) (cm$^{-2}$)' + ax2.set_xlabel(label) + ax2.set_ylabel('Number of components') + + gs.update(wspace=0.5, hspace=0.5) + filename = 'magmo-hist-derived.pdf' + plt.savefig(filename, bbox_inches="tight") + plt.close() + + +def plot_derived2(gas_array): + fig = plt.figure(figsize=(7.5, 3)) + gs = matplotlib.gridspec.GridSpec(1, 2) + + # Mach number + ax1 = fig.add_subplot(gs[0, 0]) + sample = np.ma.array(gas_array['mach']).compressed() + bins = np.linspace(0,50,21) + hist, edges = np.histogram(sample, bins=bins) + outliers = len(sample[sample > 50]) + hist[-1] += outliers + hist = np.append(hist[0], hist) + ax1.step(edges, hist) # , width=edges[1]) + ax1.set_xlabel('Turbulent Mach number') + ax1.set_ylabel('Number of components') + + # Cold fraction + ax2 = fig.add_subplot(gs[0, 1]) + sample = 48 / np.ma.array(gas_array['temp_spin']).compressed() + bins = np.linspace(0, 1, 21) + hist, edges = np.histogram(sample, bins=bins) + outliers = len(sample[sample > 1]) + hist[-1] += outliers + hist = np.append(hist[0], hist) + ax2.step(edges, hist) # , width=edges[1]-edges[0]) + label = 'Fraction of cold gas' + ax2.set_xlabel(label) + ax2.set_ylabel('Number of components') + + gs.update(wspace=0.5, hspace=0.5) + filename = 'magmo-hist-derived2.pdf' + plt.savefig(filename, bbox_inches="tight") + plt.close() + + def plot_histograms(gas_table): gas_array = gas_table.array + + plot_observed(gas_array) + plot_equiv_width(gas_array) + plot_derived(gas_array) + plot_derived2(gas_array) + + fig = plt.figure(figsize=(7.5, 3)) + ax = plt.gca() + plot_equiv_width_lv(gas_array, ax, 'equiv_width', 'Equivalent Width (km/s)') + plt.savefig('magmo-equiv-width-lv.pdf', bbox_inches="tight") + plt.close() + + fig = plt.figure(figsize=(7.5, 3)) + ax = plt.gca() + plot_spin_temp_lv(gas_array, ax, 'temp_spin', 'Spin Temperature (K)', max_clip=300) + plt.savefig('magmo-spin-temp-lv.pdf', bbox_inches="tight") + plt.close() + + plot_single_hist(gas_array, 'optical_depth', 0, 1.2, label='$e^{-\\tau}$') plot_single_hist(gas_array, 'tau', 0, 8, label='$\\tau$') plot_single_hist(gas_array, 'fwhm', 0, 60, label='FWHM (km/s)') - plot_single_hist(gas_array, 'equiv_width', 0, 60, label='Equivalent Width (km/s)') - plot_single_hist(gas_array, 'temp_spin', 0, 450, label='Spin Temperature (K)') + #plot_single_hist(gas_array, 'mach', 0, 50, label='Mach Number (M$_t$)') + plot_single_hist(gas_array, 'temp_spin', 0, 450, label='Spin Temperature (K)', sample = np.ma.array(gas_array['temp_spin']).compressed()) # Column density is special - sample = np.array(gas_array['column_density']) / 1e20 - bins = 10 ** np.linspace(np.log10(np.min(sample)), np.log10(np.max(sample)), 50) + sample = np.log10(np.array(gas_array['column_density'])) # / 20.0 + bins = np.linspace(19, 24, 21) + #bins = 10 ** np.linspace(np.log10(np.min(sample)), np.log10(np.max(sample)), 50) hist, edges = np.histogram(sample, bins=bins) - plt.bar(edges[:-1], hist, width=edges[1]) - label = 'Column Density $N_{H20}$ (1E20 g/cm^3)' + outliers = len(sample[sample > 24]) + hist[-1] += outliers + + plt.bar(edges[:-1], hist, width=edges[1]-edges[0]) + label = 'Column Density $\\log_{10}(N_{H}$) (g/cm^3)' plt.xlabel(label) plt.ylabel('Number of components') filename = 'magmo-hist-column_density.pdf' - plt.savefig(filename) + plt.savefig(filename, bbox_inches="tight") plt.close() #sample_no_outliers = sample[sample <= max] print("{} had {} values, mean {:0.3f}, median {:0.3f}, sd {:0.3f} outliers".format('column_density', len(sample), @@ -423,6 +765,52 @@ def plot_histograms(gas_table): # np.std(tau), outliers)) +def plot_uncertainty(gas_table): + fig = plt.figure(figsize=(7.5, 3)) + gs = matplotlib.gridspec.GridSpec(1, 2) + + # Uncertainty spread + gas_array = gas_table.array + x = np.ma.array(gas_array['temp_spin']).compressed() + y = np.ma.array(gas_array['delta_temp_spin']).compressed() + + # Calculate the point density + xy = np.vstack([x, y]) + z = stats.gaussian_kde(xy)(xy) + idx = z.argsort() + x, y, z = x[idx], y[idx], z[idx] + + ax1 = fig.add_subplot(gs[0, 0]) + ax1.scatter(x, y, c=z, s=20, edgecolor='', cmap='viridis') + ax1.set_ylim((0, 400)) + ax1.set_xlim(0) + ax1.set_xlabel("Spin Temperature (K)") + ax1.set_ylabel("$\Delta$ Spin Temperature (K)") + + # Uncertainty by fraction + x = np.ma.array(gas_array['temp_spin']).compressed() + y = np.ma.array(gas_array['delta_temp_spin']).compressed() / x + print("Median T_S uncertainty fracton is", np.median(y)) + # Calculate the point density + xy = np.vstack([x, y]) + z = stats.gaussian_kde(xy)(xy) + idx = z.argsort() + x, y, z = x[idx], y[idx], z[idx] + + #fig, ax = plt.subplots() + ax2 = fig.add_subplot(gs[0, 1]) + ax2.scatter(x, y, c=z, s=20, edgecolor='', cmap='viridis') + ax2.set_ylim((0,2)) + ax2.set_xlim(0) + # plt.yscale('log') + ax2.set_xlabel("Spin Temperature (K)") + ax2.set_ylabel("$\Delta T_S / T_S$") + #plt.savefig('magmo-spin-temp-uncertainty-fraction.pdf', bbox_inches="tight") + #plt.close() + plt.savefig('magmo-spin-temp-uncertainty.pdf', bbox_inches="tight") + plt.close() + + def find_best_matches(magmo_coords, other_coords, max_dist, magmo_table): """ Match the best MAGMO spectra within a defined distance to each source in a list of locations. @@ -585,6 +973,8 @@ def compare_brown_2014(magmo_coords, magmo_table): """ print("## Comparing with Brown et al 2014 ##") + magmo.ensure_dir_exists('figures') + brown_table = ascii.read('../Brown-2014-HII-Regions.dat', format='cds') brown_coords = SkyCoord(brown_table['GLON'], brown_table['GLAT'], frame='galactic', unit="deg") @@ -596,9 +986,9 @@ def compare_brown_2014(magmo_coords, magmo_table): combined.add_column(dist_col) combined.sort('GLON') noisy_count = np.zeros(len(combined)) - dist_col = Column(name='Noisy_Count', data=noisy_count, + noisy_col = Column(name='Noisy_Count', data=noisy_count, description='Count of optical depth measurements more than 4 sigma difference between MAGMO and SGPS') - combined.add_column(dist_col) + combined.add_column(noisy_col) # Produce plots and a latex template i = 1 @@ -643,7 +1033,7 @@ def compare_brown_2014(magmo_coords, magmo_table): writer.writerow( [combined[i]['SIMBAD'], combined[i]['Name'], "{:0.0f}".format(combined[i]['Separation']*3600), "{:0.0f}".format(combined[i]['Noisy_Count']), combined[i]['Qual'], combined[i]['Rating'], - "{:0.3f}".format(math.log(1 - combined[i]['Continuum_SD']) * -1), + "{:0.2f}".format(math.log(1 - combined[i]['Continuum_SD']) * -1), "{:0.4f}".format(combined[i]['GLON']), "{:0.4f}".format(combined[i]['GLAT'])]) return len(matches) @@ -675,6 +1065,139 @@ def compare_dickey_2003(magmo_coords, magmo_table): combined.write('dm-combined.vot', format='votable', overwrite=True) +def build_hist_fraction(sample, bins, clip_val): + hist, edges = np.histogram(sample, bins=bins) + outliers = len(sample[sample > clip_val]) + hist[-1] += outliers + hist = hist / np.sum(hist) + hist = np.append(hist[0], hist) + return hist, edges + + +def compare_heiles_2003(magmo_gas): + print("## Comparing with Heiles & Troland 2003 ##") + + # Read in ht03 data + rdr = ascii.get_reader(Reader=ascii.Csv) + ht03_table = rdr.read('../Millennium_data.csv') + # filter for just the CNM data |b| < 10 + cnm = np.array(ht03_table['CNM']) + sample = ht03_table[cnm >= '0'] + abs_lat = np.absolute(np.array(sample['GLAT'])) + ht03_low_cnm = sample[abs_lat <= 10] + spin_temp = np.array(ht03_low_cnm['Ts']) + print("Sample had {} values, mean {:0.3f}, median {:0.3f}, sd {:0.3f}".format(len(spin_temp), + np.mean(spin_temp), + np.median(spin_temp), + np.std(spin_temp))) + + votable = from_table(ht03_low_cnm) + writeto(votable, 'millenium_spin.vot') + + # comparative histogram of the two CNM spin temp sets + fig = plt.figure(figsize=(7.5, 3)) + gs = matplotlib.gridspec.GridSpec(1, 2) + gas_array = magmo_gas + + # Spin Temperature + ax1 = fig.add_subplot(gs[0, 0]) + sample = np.ma.array(gas_array['temp_spin']).compressed() + bins = np.linspace(0,450,19) + hist, edges = build_hist_fraction(sample, bins, 450) + ax1.step(edges, hist) + + ht03_sample = np.array(ht03_low_cnm['Ts']) + hist, edges = build_hist_fraction(ht03_sample, bins, 450) + ax1.step(edges, hist, color='black', ls='--') + + ax1.set_xlabel('Spin Temperature (K)') + ax1.set_ylabel('Fraction of components') + + statistic, p_value = stats.ks_2samp(np.ma.filled(sample), np.ma.filled(ht03_sample)) + print ('Spin temp population similarity p_value={}'.format(p_value)) + + + # Column Density + ax2 = fig.add_subplot(gs[0, 1]) + sample = np.ma.array(gas_array['column_density']).compressed() + sample = np.log10(sample) + bins = np.linspace(19, 24, 21) + hist, edges = build_hist_fraction(sample, bins, 24) + + sample = np.array(ht03_low_cnm['NHI']) * 1E20 + sample = np.log10(sample[sample > 0]) + ht03_hist, edges = build_hist_fraction(sample, bins, 24) + + ax2.step(edges, hist) #, width=edges[1]-edges[0]) + ax2.step(edges, ht03_hist, color='black', ls=':') # , width=edges[1]-edges[0] + label = 'Column Density $\\log_{10}(N_{H}$) (cm$^{-2}$)' + ax2.set_xlabel(label) + ax2.set_ylabel('Fraction of components') + + statistic, p_value = stats.ks_2samp(np.ma.filled(gas_array['column_density']), np.ma.filled(np.array(ht03_low_cnm['NHI']) * 1E20)) + print ('Column density population similarity p_value={}'.format(p_value)) + + gs.update(wspace=0.5, hspace=0.5) + filename = 'magmo-heiles_2003_comp.pdf' + plt.savefig(filename, bbox_inches="tight") + plt.close() + + return + + +def compare_strasser_2007(magmo_gas): + print("## Comparing with Strasser et al 2007 ##") + + # Read in s07 data + s07_files = glob.glob('../strasser2007/*.dat') + s07_data = None + for file_name in s07_files: + data = ascii.read(file_name) + spin_temp = np.array(data['Ts']) + subset = data[spin_temp > 0] + if s07_data is None: + s07_data = subset + else: + s07_data = vstack([s07_data, subset]) + + s07_spin_temp = s07_data['Ts'] + s07_spin_temp = s07_spin_temp[s07_spin_temp<1E4] + print("Sample had {} values, mean {:0.3f}, median {:0.3f}, sd {:0.3f}".format(len(s07_spin_temp), + np.mean(s07_spin_temp), + np.median(s07_spin_temp), + np.std(s07_spin_temp))) + votable = from_table(s07_data) + writeto(votable, 'strasser_2007_spin.vot') + + fig = plt.figure(figsize=(7.5, 3)) + gs = matplotlib.gridspec.GridSpec(1, 2) + gas_array = magmo_gas + + # Spin Temperature + ax1 = fig.add_subplot(gs[0, 0]) + sample = np.ma.array(gas_array['temp_spin']).compressed() + bins = np.linspace(0,450,19) + hist, edges = build_hist_fraction(sample, bins, 450) + ax1.step(edges, hist) + + s07_sample = np.array(s07_data['Ts']) + s07_sample = s07_sample[s07_sample<1E4] + hist, edges = build_hist_fraction(s07_sample, bins, 450) + ax1.step(edges, hist, color='black', ls='--') + + ax1.set_xlabel('Spin Temperature (K)') + ax1.set_ylabel('Fraction of components') + + statistic, p_value = stats.ks_2samp(np.ma.filled(sample), np.ma.filled(s07_sample)) + print ('Spin temp population similarity p_value={}'.format(p_value)) + + gs.update(wspace=0.5, hspace=0.5) + filename = 'magmo-strasser_2007_comp.pdf' + plt.savefig(filename, bbox_inches="tight") + plt.close() + return + + def report_close_neighbours(magmo_coords, magmo_table): idx, sep2d, dist3d = matching.match_coordinates_sky(magmo_coords, magmo_coords, 2) idx_orig = np.arange(len(magmo_coords)) @@ -713,19 +1236,24 @@ def main(): # Output a catalogue gas_table = output_gas_catalogue(all_gas) - plot_equiv_width_lv(gas_table) # Plots of MAGMO specific plot_maser_comparison(gas_table) plot_histograms(gas_table) + plot_spectra(spectra_array) + plot_uncertainty(gas_table) # Comparisons num_matches_brown = 0 if not args.no_compare: magmo_coords, magmo_table = get_magmo_table() + votable = parse("magmo-gas.vot", pedantic=False) + magmo_gas = votable.get_first_table().to_table() num_matches_brown = compare_brown_2014(magmo_coords, magmo_table) compare_dickey_2003(magmo_coords, magmo_table) report_close_neighbours(magmo_coords, magmo_table) + compare_heiles_2003(magmo_gas) + compare_strasser_2007(magmo_gas) # Report end = time.time()