Commit adf6edeb authored by Andrea Giannetti's avatar Andrea Giannetti
Browse files

Completed density estimates in aperture spectra

parent 9653e2dd

etl/assets/commons/__init__.py

+2 −2
Changes for etl/assets/commons/__init__.py: 2 added lines, 2 removed lines.
Original line number Diff line number Diff line
@@ -82,13 +82,13 @@ def convert_ratio_to_density(ratio_value, ratio_name=None): 
        '87-86': {

            'a': 0.523,

            'b': 6.086,

            'scale': 1.009,

            'scale': 1.009 / 2,

            'nu': 0.971

        },

        '88-86': {

            'a': 0.520,

            'b': 6.728,

            'scale': 1.065,

            'scale': 1.065 / 2,

            'nu': 0.976

        }

    }

etl/config/ch3oh_regions.yml

+25 −24
Changes for etl/config/ch3oh_regions.yml: 25 added lines, 24 removed lines.
Original line number Diff line number Diff line 
G28: {

  'cont_core_01': [229, 256, 10, 1],

  'cont_core_02': [239, 278, 10, 1],

  'cont_core_03': [254, 255, 10, 1],

  'cont_core_04': [278, 254, 10, 3],

  'cont_core_05': [248, 265, 10, 1],

  'cont_core_06': [159, 186, 10, 2],

  'cont_core_07': [261, 176, 10, 2],

  'cont_core_08': [291, 190, 10, 4],

  'core_cont_01': [229, 256, 10, 1],

  'core_cont_02': [239, 278, 10, 1],

  'core_cont_03': [254, 255, 10, 1],

  'core_cont_04': [278, 254, 10, 3],

  'core_cont_05': [248, 265, 10, 1],

  'core_cont_06': [159, 186, 10, 2],

  'core_cont_07': [261, 176, 10, 2],

  'core_cont_08': [291, 190, 10, 4],

  'envelope_01': [204, 241, 15, 1],

  'envelope_02': [307, 341, 25, 1],

  'envelope_03': [349, 393, 25, 1],
@@ -14,27 +14,28 @@ G28: { 
  'envelope_05': [317, 259, 25, 1],

  'envelope_06': [200, 194, 25, 1],

  'envelope_07': [167, 336, 25, 3],

  'core_envelope_01': [150, 200, 10, 1],

  'envelope_core_01': [150, 200, 10, 1],

}

G14: {

    'cont_core_01': [220, 260, 8, 1],

    'cont_core_02': [301, 251, 8, 2],

    'cont_core_03': [252, 237, 8, 3],

    'cont_core_04': [292, 271, 8, 1],

    'cont_core_05': [198, 289, 8, 3],

    'cont_core_06': [284, 245, 8, 2],

    'cont_core_07': [305, 268, 8, 2],

    'core_cont_01': [220, 260, 8, 1],

    'core_cont_02': [301, 251, 8, 2],

    'core_cont_03': [252, 237, 8, 3],

    'core_cont_04': [292, 271, 8, 1],

    'core_cont_05': [198, 289, 8, 3],

    'core_cont_06': [284, 245, 8, 2],

    'core_cont_07': [305, 268, 8, 2],

    'envelope_01': [251, 292, 25, 2],

    'envelope_02': [181, 249, 25, 2],

    'envelope_03': [290, 335, 25, 1],

#    'envelope_04': [345, 309, 25, 1],

    'ch3oh_core_01': [251, 213, 25, 3],

    'ch3oh_core_02': [301, 224, 8, 1],

    'envelope_02': [181, 249, 25, 1],

    'envelope_03': [288, 338, 25, 1],

    'envelope_04': [328, 318, 20, 1],

    'envelope_05': [185, 197, 25, 1],

    'core_ch3oh_01': [251, 208, 20, 3],

    'core_ch3oh_02': [301, 224, 8, 1],

}

G332.96: {

    'cont_core_01': [137, 272, 15, 1],

    'cont_core_02': [249, 303, 25, 1],

    'cont_core_03': [303, 200, 25, 1],

    'core_cont_01': [137, 272, 15, 1],

    'core_cont_02': [249, 303, 25, 1],

    'core_cont_03': [303, 200, 25, 1],

    'envelope_01': [276, 244, 25, 2],

    'envelope_02': [328, 242, 25, 3],

    'envelope_03': [334, 157, 30, 1],

etl/config/guesses.yml

+34 −29
Changes for etl/config/guesses.yml: 34 added lines, 29 removed lines.
Original line number Diff line number Diff line 
G28: {

  cont_core_01: [

  core_cont_01: [

    [ [ 0.07, 0.1, 0.02, 0.01 ], ],

    [ 0, ],

    [ 5.0e+4, ]

  ],

  cont_core_02: [

  core_cont_02: [

    [ [ 0.12, 0.15, 0.05, 0.02 ] ],

    [ 0, ],

    [ 2.5e+5 ]

  ],

  cont_core_03: [

  core_cont_03: [

    [ [ 0.05, 0.07, 0.03, 0.01 ] ],

    [ 0, ],

    [ 2.5e+5, ]

  ],

  cont_core_04: [

  core_cont_04: [

    [ [ 0.025, 0.03, 0.001, 0.0005 ], [ 0.025, 0.03, 0.001, 0.0003 ], [ 0.02, 0.02, 0.001, 0.0003 ] ],

    [ 0, -10.0e+5, 10.0e+5 ],

    [ 1.5e+5, 1.5e+5, 1.5e+5 ]

  ],

  cont_core_05: [

  core_cont_05: [

    [ [ 0.05, 0.07, 0.03, 0.01 ] ],

    [ 0, ],

    [ 2.5e+5, ]

  ],

  cont_core_06: [

  core_cont_06: [

    [ [ 0.05, 0.07, 0.05, 0.01 ], [ 0.02, 0.05, 0.001, 0.0003 ] ],

    [ 0, -5.0e+5 ],

    [ 2.0e+5, 2.0e+5 ]

  ],

  cont_core_07: [

  core_cont_07: [

    [ [ 0.015, 0.02, 0.003, 0.001 ], [ 0.012, 0.017, 0.001, 0.0003 ] ],

    [ -5.0e+5, 5.0e+5 ],

    [ 2.0e+5, 2.0e+5 ]

  ],

  cont_core_08: [

  core_cont_08: [

    [ [ 0.015, 0.02, 0.005, 0.001 ], [ 0.02, 0.03, 0.01, 0.003 ], [ 0.012, 0.013, 0.002, 0.001 ], [ 0.001, 0.012, 0.001, 0.0005 ] ],

    [ -4.0e+5, 0.0e+5, 9.5e+5, 35.0e+5 ],

    [ 1.9e+5, 1.9e+5, 1.8e+5, 3.0e+5 ]
@@ -74,7 +74,7 @@ G28: { 
    [ 6.0e+5, -6.0e+5, -2.0e+5 ],

    [ 2.0e+5, 2.0e+5, 2.0e+5 ]

  ],

  core_envelope_01: [

  envelope_core_01: [

    [ [ 0.025, 0.035, 0.003, 0.001 ], ],

    [ 0, ],

    [ 2.5e+5, ]
@@ -82,37 +82,37 @@ G28: { 


}

G14: {

  cont_core_01: [

  core_cont_01: [

    [ [ 0.07, 0.1, 0.02, 0.01 ], ],

    [ 0, ],

    [ 5.0e+4, ]

  ],

  cont_core_02: [

  core_cont_02: [

    [ [ 0.12, 0.15, 0.05, 0.02 ], [ 0.015, 0.02, 0.01, 0.003 ] ],

    [ 0, 10.0e+5 ],

    [ 2.5e+5, 10.0e+5 ]

  ],

  cont_core_03: [

  core_cont_03: [

    [ [ 0.05, 0.07, 0.03, 0.01 ], [ 0.01, 0.02, 0.01, 0.003 ], [ 0.01, 0.02, 0.01, 0.003 ] ],

    [ 0, 10.0e+5, -30.0e+5 ],

    [ 2.5e+5, 2.5e+5, 10.0e+5 ]

  ],

  cont_core_04: [

  core_cont_04: [

    [ [ 0.06, 0.08, 0.03, 0.01 ], ],

    [ 0, ],

    [ 2.5e+5, ]

  ],

  cont_core_05: [

  core_cont_05: [

    [ [ 0.05, 0.1, 0.03, 0.01 ], [ 0.02, 0.03, 0.01, 0.003 ], [ 0.03, 0.04, 0.02, 0.003 ] ],

    [ 0, -10.0e+5, 5.0e+5 ],

    [ 2.5e+5, 10.0e+5, 2.5e+5 ]

  ],

  cont_core_06: [

  core_cont_06: [

    [ [ 0.05, 0.07, 0.05, 0.01 ], [ 0.02, 0.05, 0.001, 0.0003 ] ],

    [ -1.0e+5, -30.0e+5 ],

    [ 2.0e+5, 20.0e+5 ]

  ],

  cont_core_07: [

  core_cont_07: [

    [ [ 0.05, 0.07, 0.03, 0.01 ], [ 0.01, 0.02, 0.01, 0.003 ] ],

    [ 0, 10.0e+5 ],

    [ 2.5e+5, 10.0e+5 ]
@@ -123,43 +123,48 @@ G14: { 
    [ 2.5e+5, 10.0e+5 ]

  ],

  envelope_02: [

    [ [ 0.05, 0.07, 0.01, 0.003 ], [ 0.01, 0.02, 0.01, 0.003 ] ],

    [ 0, 12.0e+5 ],

    [ 2.5e+5, 4.0e+5 ]

    [ [ 0.05, 0.07, 0.02, 0.003 ] ],

    [ 6.0e+5, ],

    [ 2.0e+5, ]

  ],

  envelope_03: [

    [ [ 0.025, 0.035, 0.003, 0.001 ], ],

    [ 0, ],

    [ 2.5e+5, ]

  ],

  #  envelope_04: [

  #    [ [ 0.06, 0.08, 0.03, 0.01 ], ],

  #    [ 0, ],

  #    [ 2.5e+5, ]

  #  ],

  ch3oh_core_01: [

  envelope_04: [

    [ [ 0.02, 0.03, 0.003, 0.001 ], ],

    [ 0, ],

    [ 2.5e+5, ]

  ],

  envelope_05: [

    [ [ 0.02, 0.03, 0.003, 0.001 ], ],

    [ 6.0e+5, ],

    [ 2.5e+5, ]

  ],

  core_ch3oh_01: [

    [ [ 0.05, 0.07, 0.03, 0.01 ], [ 0.01, 0.02, 0.01, 0.003 ], [ 0.001, 0.002, 0.001, 0.0003 ] ],

    [ 0, 10.0e+5, -15.0e+6 ],

    [ 2.5e+5, 2.0e+5, 2.0e+6 ]

  ],

  ch3oh_core_02: [

  core_ch3oh_02: [

    [ [ 0.06, 0.08, 0.03, 0.01 ], ],

    [ 0, ],

    [ 2.5e+5, ]

  ],

}

G332.96: {

  cont_core_01: [

  core_cont_01: [

    [ [ 0.015, 0.02, 0.003, 0.001 ], ],

    [ 0, ],

    [ 2.5e+5, ]

  ],

  cont_core_02: [

  core_cont_02: [

    [ [ 0.006, 0.008, 0.001, 0.001 ], ],

    [ -1.0e+5, ],

    [ 1.8e+5, ]

  ],

  cont_core_03: [

  core_cont_03: [

    [ [ 0.012, 0.015, 0.001, 0.0005 ], ],

    [ 0, ],

    [ 2.5e+5, ]

etl/mdl/mdl_extract_ch3oh_spectra.py

+145 −15
Changes for etl/mdl/mdl_extract_ch3oh_spectra.py: 145 added lines, 15 removed lines.
Original line number Diff line number Diff line 
import pyspeckit

import yaml

import os

import matplotlib

import matplotlib.pyplot as plt

import numpy as np

from sklearn.neighbors import KernelDensity

from astropy import units as u

from astropy.io import fits

from astropy.wcs import WCS

from matplotlib.patches import Circle

from assets.commons import (get_cube,

                            convert_dvelocity_to_hz,

                            fill_guesses,

                            validate_parameter,

                            convert_ratio_to_density)

from assets.constants import source_properties, config_mdl, datadir

from matplotlib.collections import PatchCollection

from assets.constants import source_properties, config_mdl, datadir, continuum_fits

from typing import Union
@@ -35,8 +38,16 @@ def create_limits(par_number, channelwidth_hz, rms, max_linewidth: Union[float, 
    return limits, limited





sources = ['G14', ]

def filter_dict_by_key(d, key_root):

    return {k: v for k, v in d.items() if k.startswith(key_root)}





sources = ['G14', 'G28', 'G332.96']

sources = ['G28', ]

species = 'CH3OH'

samples = 100

matplotlib.use('TkAgg')



plotdir = os.path.join('prs', 'figures', 'aperture_spectra_fits')

with open(os.path.join('config', 'ch3oh_regions.yml'), 'r') as infile:

    regions = yaml.load(infile, Loader=yaml.FullLoader)
@@ -44,8 +55,20 @@ with open(os.path.join('config', 'ch3oh_regions.yml'), 'r') as infile: 
with open(os.path.join('config', 'guesses.yml'), 'r') as infile:

    guesses = yaml.load(infile, Loader=yaml.FullLoader)





def compute_density_pdf(density, density_axis=None):

    _density_axis = validate_parameter(density_axis, np.linspace(2.5, 7.5, 500))

    kde = KernelDensity(bandwidth=0.1, kernel='epanechnikov')

    _density = density[~np.isnan(density)]

    kde.fit(np.log10(_density).reshape(-1, 1))

    pdf = np.exp(kde.score_samples(list(_density_axis.reshape(-1, 1))))

    norm = np.trapz(pdf, _density_axis)

    return _density_axis, pdf / norm





for source in sources:

    moment_map = fits.open(f'moment_zero_{source}_{species}.fits')

    moment_map = fits.open(f'moment_zero_{source}_{species}.fits')[0]

    wcs = WCS(moment_map.header)

    vlsr = source_properties[source]['vlsr'] * u.km / u.s

    rms = source_properties[source][species]['rms']

    linewidth_guess = source_properties[source][species]['linewidth_guess_km_s'] * u.km / u.s
@@ -62,7 +85,6 @@ for source in sources: 


    tied = []

    line_guesses = []

    matplotlib.use('TkAgg')



    filename = source_properties[source][species]['data']

    cube = get_cube(filename=filename,
@@ -78,6 +100,11 @@ for source in sources: 
    else:

        max_linewidth_hz = None



    hdu_cont = fits.open(os.path.join(datadir, continuum_fits[source]))[0]

    continuum_data = np.where(hdu_cont.data >= 2.5e-4, hdu_cont.data, np.nan)

    wcs_cont = WCS(hdu_cont.header)



    densities = {}

    for region in regions[source]:

        print(region)

        _region_coordinates = regions[source][region][0:3]
@@ -87,6 +114,7 @@ for source in sources: 
        line_guesses = []

        tied = []



        densities[region] = []

        for component in range(components):

            idx_shift = component * pars_per_component
@@ -117,19 +145,121 @@ for source in sources: 
                       limits=limits,

                       limited=limited,

                       clear=True)

        apspec.plotter()

        apspec.specfit.plot_fit(show_components=True)

        plt.show()



        fig = plt.figure(figsize=(15,5))

        ax_s = fig.add_subplot(122)

        apspec.plotter(title=f'{source} {region}', figure=fig)

        # plt.annotate('a', xy=(0.8, 0.75), xycoords='figure fraction')

        apspec.specfit.plot_fit(show_components=True, annotate=False)

        # apspec.plotter.savefig(os.path.join(plotdir, f'aperture_spectra_{source}_{species}_{region}.svg'))



        spec_rms = np.nanstd(apspec.specfit.residuals)

        ratio = apspec.specfit.parinfo['AMPLITUDE2'].value / apspec.specfit.parinfo['AMPLITUDE0'].value



        if apspec.specfit.parinfo['AMPLITUDE3'].value > 3 * spec_rms:

            ratio2 = apspec.specfit.parinfo['AMPLITUDE3'].value / apspec.specfit.parinfo['AMPLITUDE0'].value

        for vel_component in range(components):

            if ((source == 'G332.96') or ((source == 'G14') and (region == 'envelope_01'))) and (vel_component > 0):

                pass

            else:

            ratio2 = np.nan

                idx_shift = vel_component * 4

                line_1 = np.random.default_rng().normal(scale=spec_rms, size=samples) + apspec.specfit.parinfo[

                    f'AMPLITUDE{idx_shift}'].value

                line_2 = np.random.default_rng().normal(scale=spec_rms, size=samples) + apspec.specfit.parinfo[

                    f'AMPLITUDE{idx_shift + 2}'].value

                line_3 = np.random.default_rng().normal(scale=spec_rms, size=samples) + apspec.specfit.parinfo[

                    f'AMPLITUDE{idx_shift + 3}'].value

                ratio_87 = np.where(line_2 <= 0, 0.3 * spec_rms, line_2) / np.where(line_1 <= 0, spec_rms, line_1)

                ratio_88 = np.where(line_3 <= 0, 0.3 * spec_rms, line_3) / np.where(line_1 <= 0, spec_rms, line_1)



            # if apspec.specfit.parinfo[f'AMPLITUDE{idx_shift+3}'].value > 1.5 * spec_rms:

            #     ratio2 = apspec.specfit.parinfo[f'AMPLITUDE{idx_shift+3}'].value / apspec.specfit.parinfo[f'AMPLITUDE{idx_shift}'].value

            # else:

            #     ratio2 = np.nan





            density_pdfs = {}

            for ratio_string, ratio in zip(['87-86', '88-86'], [ratio_87, ratio_88]):

                density = convert_ratio_to_density(ratio, ratio_name=ratio_string)

                if (np.isnan(density).sum() == len(density)) or (np.nanmean(ratio) > 0.95):

                    pdf = np.ones(500)

                    density_axis = np.linspace(2.5, 7.5, 500)

                else:

                    density_axis, pdf = compute_density_pdf(density)

                density_pdfs[ratio_string] = {'density': density_axis, 'pdf': pdf}

            final_density_pdf = density_pdfs['87-86']['pdf'] * density_pdfs['88-86']['pdf']

            norm = np.trapz(final_density_pdf, density_axis)

            final_density_pdf /= norm

            densities[region].append(density_axis[np.argmax(final_density_pdf)])



            # plt.clf()

            # plt.plot(density_axis, final_density_pdf)

            # plt.plot(density_axis, density_pdfs['87-86']['pdf'])

            # plt.plot(density_axis, density_pdfs['88-86']['pdf'])

            # plt.show()



            # fig = plt.figure()

            # density_axis, pdf_87 = compute_density_pdf(density)

            # _, pdf_88 = compute_density_pdf(density)

            # kde2 = KernelDensity(bandwidth=0.1, kernel='epanechnikov')

            # kde2.fit(np.log10(density2).reshape(-1, 1))

            # z2 = np.exp(kde2.score_samples(list(density_axis.reshape(-1, 1))))

            # plt.plot(density_axis, z * z2)

            # sns.kdeplot(np.log10(density), figure=fig)

            # sns.kdeplot(np.log10(density2), figure=fig)

            # plt.show()

            # print(np.nanmean(ratio), np.nanmean(ratio2), np.nanmean(density), np.nanstd(density),

            #       np.nanmean(density2), np.nanstd(density2), spec_rms)



        # plt.show()



        # plt.clf()

        # hdu_cont = fits.open(os.path.join(datadir, continuum_fits[source]))[0]

        # continuum_data = np.where(hdu_cont.data >= 2.5e-4, hdu_cont.data, np.nan)

        # wcs = WCS(hdu_cont.header)

        # fig = plt.figure()



        ax = fig.add_subplot(121, projection=wcs)

        plt.imshow(np.log10(moment_map.data), origin='lower')

        plt.colorbar(label=r'log(K km s$^{-1}$)')

        patch = Circle(_region_coordinates[0:2], _region_coordinates[2], edgecolor='black', facecolor='none')

        ax.set_autoscale_on(False)

        ax.add_patch(patch)

        plt.contour(continuum_data,

                    levels=np.logspace(np.log10(2.5e-4), np.log10(0.9 * np.nanmax(continuum_data)), num=5),

                    transform=ax.get_transform(wcs_cont),

                    colors='red',

                    origin='lower')

        plt.ylabel('Dec[J2000]')

        plt.xlabel('RA[J2000]')

        plt.tight_layout()

        plt.savefig(os.path.join(plotdir, f'aperture_spectra_extraction_{source}_{region}_{species}.svg'))



    region_types = ['envelope', 'core']

    fig = plt.figure(figsize=(7, 5))



        density = convert_ratio_to_density(ratio)

        print(ratio, ratio2, density, convert_ratio_to_density(ratio2, ratio_name='88-86'))
 
 No newline at end of file
 
    for rtype in region_types:

        plt.clf()

        ax = plt.subplot(projection=wcs)

        plt.imshow(np.log10(moment_map.data), origin='lower')

        plt.colorbar(label=r'log(K km s$^{-1}$)')

        ax.set_autoscale_on(False)

        regions_to_plot = filter_dict_by_key(regions[source], rtype)

        patches = []

        density_values = []

        for reg in regions_to_plot:

            patch = Circle(regions_to_plot[reg][0:2], regions_to_plot[reg][2], edgecolor='black', facecolor='none',

                           zorder=10)

            patches.append(patch)

            density_values.append(np.max(densities[reg]))

        pc = PatchCollection(patches, cmap=matplotlib.cm.hot, zorder=10)

        pc.set_array(density_values)

        pc.set_facecolor('none')

        pc.set_clim(4.8, 5.8)

        ax.add_collection(pc)

        plt.contour(continuum_data,

                    levels=np.logspace(np.log10(2.5e-4), np.log10(0.9 * np.nanmax(continuum_data)), num=5),

                    transform=ax.get_transform(wcs_cont),

                    colors='grey',

                    origin='lower')

        plt.ylabel('Dec[J2000]')

        plt.xlabel('RA[J2000]')

        plt.colorbar(pc, label=r'log($n($H$_2)$ [cm$^{-3}$])')

        plt.tight_layout()

        plt.savefig(os.path.join(plotdir, f'aperture_spectra_extraction_{source}_{rtype}_{species}.svg'))

etl/prs/prs_mom0_dust.py

+1 −6
Changes for etl/prs/prs_mom0_dust.py: 1 added line, 6 removed lines.
Original line number Diff line number Diff line
@@ -7,6 +7,7 @@ import yaml 
import pyspeckit

import aplpy

from scipy.optimize import curve_fit

from assets.constants import continuum_fits

from astropy.io import fits

from astropy import units as u

from astropy.constants import c as speed_of_light
@@ -23,12 +24,6 @@ source_names = ['G332.96', 'G28', 'G14', 'G12'] 
source_names = ['G332.96', 'G28', 'G14', ]

# source = 'G332.96'



continuum_fits = {

    'G332.96': 'tG332.96.cont_1.33mm.12m7m.multiscale.image.fits',

    'G28': 'tG028.27.cont_1.33mm.12m7m.multiscale.image.fits',

    'G14': 'tG014.49.cont_1.33mm.12m7m.multiscale.image.fits',

}



marker_opacity = 0.05

marker_size = 1.5

marker_color = 'DarkSlateGrey'