Commit 36e79b86 authored by Andrea Giannetti's avatar Andrea Giannetti
Browse files

Added aperture spectra fits

parent b49f74ac

etl/assets/commons/__init__.py

+33 −3
Changes for etl/assets/commons/__init__.py: 33 added lines, 3 removed lines.
Original line number Diff line number Diff line
@@ -66,8 +66,33 @@ def load_config_file(config_file_path: str, 
    return config





def convert_ratio_to_density(ratio):

    return 10 ** ((np.arctanh((ratio - 0.523) / 1.009) + 6.086) ** 0.971)

def convert_ratio_to_density(ratio_value, ratio_name=None):

    """

    Convert the line ratio of methanol to a guess for the number density of H2



    :param ratio_value: The line ratio

    :param ratio_name: The name of the line ratio

    :return: The guessed density, in cm^-3



    This function is an empirical fit to the relationship between line ratio and number density of H2.

    It is based on the relationship between the line ratio and the number density, derived from the results of a grid of models.

    """

    _ratio_name = validate_parameter(ratio_name, default='87-86')

    fit_pars = {

        '87-86': {

            'a': 0.523,

            'b': 6.086,

            'scale': 1.009,

            'nu': 0.971

        },

        '88-86': {

            'a': 0.520,

            'b': 6.728,

            'scale': 1.065,

            'nu': 0.976

        }

    }

    return 10 ** ((np.arctanh((ratio_value - fit_pars[_ratio_name]['a']) / fit_pars[_ratio_name]['scale']) + fit_pars[_ratio_name]['b']) ** fit_pars[_ratio_name]['nu'])





def gradient_plane_3d(xy, ra_grad, dec_grad, offset):
@@ -100,8 +125,13 @@ def get_cube(filename, datadir, frequency_limits=None, smooth_data=False): 
def fill_guesses(line_peak, doppler_adjusted_frequencies, shift_hz, linewidth_guess_hz, rms,

                 peak_scaling_factors=(1, 1, 0.5, 0.3)):

    guesses = []

    if isinstance(line_peak, list) or isinstance(line_peak, tuple):

        assert len(line_peak) == len(doppler_adjusted_frequencies)

        _line_peaks = line_peak

    else:

        _line_peaks = [peak_scaling_factors[idx] * line_peak for idx in range(len(peak_scaling_factors))]

    for idx, freq in enumerate(doppler_adjusted_frequencies):

        guesses += [np.max((rms, peak_scaling_factors[idx] * line_peak)),

        guesses += [np.max((rms, _line_peaks[idx])),

                    freq.to(u.Hz).value + shift_hz,

                    linewidth_guess_hz]

    return guesses

etl/assets/data/source_properties.yml

+1 −1
Changes for etl/assets/data/source_properties.yml: 1 added line, 1 removed line.
Original line number Diff line number Diff line
@@ -25,7 +25,7 @@ G28: { 
    'data': 'G28_CH3OH_3mm_feather.fits',

    'linewidth_guess_km_s': 1.5,

    'rms': 0.003,

    'frequency_limits_ghz': [ 96.6995, 96.740 ] },

    'frequency_limits_ghz': [ 96.6995, 96.760 ] },



  N2H+:{

    'data': 'G28_N2H+_3mm_feather.fits',

etl/config/ch3oh_regions.yml

+41 −1
Changes for etl/config/ch3oh_regions.yml: 41 added lines, 1 removed line.
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],

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

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

  'envelope_03': [349, 393, 25, 1],

  'envelope_04': [200, 300, 25, 1],

  '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],

}

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],

    '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],

}

G332.96: {

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

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

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

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

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

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

    'envelope_04': [210, 204, 30, 1],

    'envelope_05_global': [285, 215, 55, 3],

}
 
 No newline at end of file

etl/config/config_visfit.yml

+3 −3
Changes for etl/config/config_visfit.yml: 3 added lines, 3 removed lines.
Original line number Diff line number Diff line
@@ -5,7 +5,7 @@ G332.96: { 
  min_velocity: -68,

  surface_point: [291, 216, -66.43255],

  plot_velocity_limits: [-70, -62],

  plot_species: ['CH3OH'],

  plot_species: ['N2H+','CH3OH'],

  continuum_limits: [0.1, 1]

}

G28: {
@@ -13,7 +13,7 @@ G28: { 
  min_velocity: -inf,

  surface_point: [272, 278, 79.2268463],

  plot_velocity_limits: [74, 81],

  plot_species: ['CH3OH'],

  plot_species: ['N2H+', 'CH3OH'],

  continuum_limits: [0.1, 3]

}

G14: {
@@ -21,6 +21,6 @@ G14: { 
  min_velocity: -inf,

  surface_point: [197, 291, 40.2304702],

  plot_velocity_limits: [37.5, 42],

  plot_species: ['N2H+'],

  plot_species: ['N2H+', 'CH3OH'],

  continuum_limits: [0.1, 10]

}

etl/mdl/mdl_extract_ch3oh_spectra.py

+85 −29
Changes for etl/mdl/mdl_extract_ch3oh_spectra.py: 85 added lines, 29 removed lines.
Original line number Diff line number Diff line
@@ -5,8 +5,14 @@ import matplotlib 
import matplotlib.pyplot as plt

import numpy as np

from astropy import units as u

from assets.commons import get_cube, convert_dvelocity_to_hz, fill_guesses

from astropy.io import fits

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 typing import Union





def fill_tied_pars(tied_pars):
@@ -15,8 +21,31 @@ def fill_tied_pars(tied_pars): 
    return _tied_pars_list





source = 'G28'

def create_limits(par_number, channelwidth_hz, rms, max_linewidth: Union[float, None] = None):

    limits = []

    _max_linewidth = validate_parameter(max_linewidth, 1e7)

    for idx in range(par_number):

        if idx % 3 == 0:

            limits.append((rms, 1))

        elif idx % 3 == 2:

            limits.append((channelwidth_hz / 2, _max_linewidth))

        else:

            limits.append((0, 0))

    limited = [(True, True) if (idx % 3 == 0) or (idx % 3 == 2) else (False, False) for idx in range(par_number)]

    return limits, limited





sources = ['G14', ]

species = 'CH3OH'

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)



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

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



for source in sources:

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

    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
@@ -32,11 +61,34 @@ delta_freq = rest_frequencies[1:] - rest_frequencies[0] 
    pars_per_component = 3 * len(rest_frequencies)



    tied = []

guesses = []

peak_guess = 0.03

manual_components = 1

for components in range(manual_components):

    idx_shift = components * pars_per_component

    line_guesses = []

    matplotlib.use('TkAgg')



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

    cube = get_cube(filename=filename,

                    datadir=datadir,

                    frequency_limits=frequency_limits,

                    smooth_data=config_mdl['smooth_data'])

    channelwidth_hz = cube.header['CDELT3']

    max_linewidth = source_properties[source][species]['max_linewidth_km_s'] if \

        'max_linewidth_km_s' in source_properties[source][species].keys() else None

    if max_linewidth is not None:

        max_linewidth_hz = convert_dvelocity_to_hz(dvelocity=max_linewidth,

                                                   rest_frequencies=rest_frequencies).value / 2.355

    else:

        max_linewidth_hz = None



    for region in regions[source]:

        print(region)

        _region_coordinates = regions[source][region][0:3]

        components = regions[source][region][3]

        region_guesses = regions[source][region]

        apspec = cube.get_apspec(_region_coordinates)

        line_guesses = []

        tied = []



        for component in range(components):

            idx_shift = component * pars_per_component



            if len(rest_frequencies) > 1:

                tied_pars = {
@@ -48,32 +100,36 @@ for components in range(manual_components): 
                }

                tied += fill_tied_pars(tied_pars=tied_pars)



        guesses += fill_guesses(line_peak=peak_guess,

                line_guesses += fill_guesses(line_peak=guesses[source][region][0][component],

                                             doppler_adjusted_frequencies=doppler_adjusted_frequencies,

                                shift_hz=linewidth_guess_hz / 2 * idx_shift,

                                linewidth_guess_hz=linewidth_guess_hz,

                                             shift_hz=guesses[source][region][1][component],

                                             linewidth_guess_hz=guesses[source][region][2][component],

                                             rms=rms)



matplotlib.use('TkAgg')



source = 'G28'

species = 'CH3OH'

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

cube = get_cube(filename=filename,

                datadir=datadir,

                frequency_limits=frequency_limits,

                smooth_data=config_mdl['smooth_data'])

                limits, limited = create_limits(par_number=len(line_guesses),

                                                channelwidth_hz=np.abs(channelwidth_hz),

                                                rms=rms / np.sqrt(np.pi * _region_coordinates[2] ** 2),

                                                max_linewidth=max_linewidth_hz)



regions = {

    'centre': (256, 256, 20)

}

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

for region in regions:

    apspec = cube.get_apspec(regions[region])

    apspec.plotter()

        apspec.specfit(fittype='gaussian',

                   guesses=guesses,

                       guesses=line_guesses,

                       tied=tied,

                   clear=True,)

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

                       limits=limits,

                       limited=limited,

                       clear=True)

        apspec.plotter()

        apspec.specfit.plot_fit(show_components=True)

        plt.show()



        # 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

        else:

            ratio2 = np.nan



        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