Commit 1f6864bc authored by Andrea Giannetti's avatar Andrea Giannetti
Browse files

Added possibility to define a hot core abundance profile.

parent 23c25514

etl/stg/stg_radmc_input_generator.py

+31 −2
Original line number Diff line number Diff line
@@ -269,14 +269,43 @@ def write_molecular_number_density_profiles(profiles: dict, 
                                            grid_metadata: dict,

                                            path: Union[str, None] = None):

    _path = validate_parameter(path, default=os.path.join('mdl', 'radmc_files'))

    hot_core_specs = read_abundance_variation_schema(line_config=line_config)

    for species in line_config['species_to_include'] + line_config['collision_partners']:

        write_radmc_input(filename=f'numberdens_{species}.inp',

                          quantity=profiles['gas_number_density'] * float(

                              line_config['molecular_abundances'][species]),

                          quantity=compute_molecular_number_density_hot_core(

                              gas_number_density_profile=profiles['gas_number_density'],

                              abundance=float(line_config['molecular_abundances'][species]),

                              temperature_profile=profiles['dust_temperature'],

                              threshold=float(hot_core_specs[species]['threshold']),

                              abundance_jump=float(hot_core_specs[species]['abundance_jump'])),

                          path=_path,

                          grid_metadata=grid_metadata)





def read_abundance_variation_schema(line_config):

    species_list = list(line_config['molecular_abundances'].keys())

    results_dict = {}

    for species in species_list:

        if species in line_config['hot_core_specs']:

            results_dict[species] = line_config['hot_core_specs'][species]

        else:

            results_dict[species] = {

                'threshold': 20,

                'abundance_jump': 1

            }

    return results_dict





def compute_molecular_number_density_hot_core(gas_number_density_profile: np.array,

                                              abundance: float,

                                              temperature_profile: np.array,

                                              threshold: float,

                                              abundance_jump: Union[float, int]) -> np.array:

    return np.where(temperature_profile < threshold,

                    gas_number_density_profile * abundance,

                    gas_number_density_profile * abundance * abundance_jump)





def save_fits_grid_profile(quantity: np.array,

                           filename: str,

                           path: str = None):

etl/tests/test_stg_radmc_input_generator.py

+80 −37
Original line number Diff line number Diff line 
import os

import numpy as np

from stg.stg_radmc_input_generator import (get_solid_body_rotation_y,

                                           get_grid_name)

                                           get_grid_name,

                                           read_abundance_variation_schema,

                                           compute_molecular_number_density_hot_core)

from astropy import units as u

from unittest import TestCase

from assets.commons import (load_config_file,
@@ -20,6 +22,47 @@ class Test(TestCase): 
    def setUp(self):

        self.config_filename = 'config.yml'



    def test_read_abundance_variation_schema(self):

        line_config = {

            'species_to_include': ['e-ch3oh'],

            'molecular_abundances': {

                "e-ch3oh": 1e-9,

                "p-h2": 1,

            },

            'hot_core_specs': {

                "e-ch3oh": {

                    'threshold': 90,

                    'abundance_jump': 100

                }

            },

            'lines_mode': 'lvg',

            'collision_partners': ['p-h2']

        }

        result_dict = read_abundance_variation_schema(line_config=line_config)

        expected_dict = {

            'e-ch3oh': {

                'threshold': 90,

                'abundance_jump': 100

            },

            'p-h2': {

                'threshold': 20,

                'abundance_jump': 1

            },

        }

        self.assertDictEqual(result_dict, expected_dict)



    def test_compute_molecular_number_density_hot_core(self):

        gas_number_density_profile = np.array([1, 1, 2, 3])

        temperature_profile = np.array([10, 100, 200, 20])

        abundance_array = compute_molecular_number_density_hot_core(

            gas_number_density_profile=gas_number_density_profile,

            abundance=1e-9,

            temperature_profile=temperature_profile,

            abundance_jump=100,

            threshold=90)

        expected_results = np.array([1e-9, 1e-7, 2e-7, 3e-9])

        self.assertTrue(np.allclose(abundance_array, expected_results, 1e-7))



    def test_get_grid_filename_composite(self):

        grid_name = get_grid_name(method='composite_grid',

                                  zip_filename='abc.def.zip',
@@ -28,10 +71,10 @@ class Test(TestCase): 


    def test_get_grid_filename_missing_info(self):

        with self.assertRaises(AssertionError):

            grid_name = get_grid_name(method='composite_grid',

            _ = get_grid_name(method='composite_grid',

                              quantity_name='h2_density')

        with self.assertRaises(AssertionError):

            grid_name = get_grid_name(method='composite_grid',

            _ = get_grid_name(method='composite_grid',

                              zip_filename='abc.def.zip')



    def test_get_grid_filename_undefined_method(self):