Commit 9f4dfd6d authored by Andrea Giannetti's avatar Andrea Giannetti
Browse files

Added central stellar heating source.

parent 069b51f7
Loading
Loading
Loading
Loading
+1 −0
Original line number Diff line number Diff line
@@ -220,6 +220,7 @@ def extract_grid_metadata(config: dict) -> dict:
    grid_metadata['distance_matrix'] = get_distance_matrix(grid_metadata, grid_metadata['centered_indices'])

    grid_metadata['grid_edges'] = get_grid_edges(grid_metadata)
    grid_metadata['continuum_lambdas'] = 100
    return grid_metadata


+1 −0
Original line number Diff line number Diff line
@@ -32,6 +32,7 @@ radmc_input_headers = {
    'wavelength_micron.inp': ['continuum_lambdas'],
    'numberdens_mol.inp': ['iformat', 'ncells'],
    'escprob_lengthscale.inp': ['iformat', 'ncells'],
    'stars.inp': ['iformat', 'nstars', 'continuum_lambdas', 'stars_properties', 'lambdas']
}


+12 −1
Original line number Diff line number Diff line
@@ -3,6 +3,7 @@ import sys
from itertools import product
from typing import Union
import sqlalchemy
import shutil
from astropy import units as u
from datetime import datetime
from radmc3dPy import image
@@ -190,7 +191,8 @@ def populate_line_table(config_lines: dict,
def main(grid_tarfile: str,
         run_id: str,
         override_config: Union[dict, None] = None,
         radmc_input_path: Union[str, None] = None) -> str:
         radmc_input_path: Union[str, None] = None,
         compute_dust_temperature: bool = True) -> str:
    # This is necessary, because the lines_mode is needed both in the lines.inp and radmc3d.inp files
    # The reason for splitting the main input file from the rest is that some parameters can be changed
    # independently of the grid for the modeling. The mdl hash should depend on all the mdl parameters, not a subset
@@ -215,8 +217,17 @@ def main(grid_tarfile: str,
    engine = get_pg_engine(logger=logger)

    # Execute radmc
    logger.debug(f'Executing command: {radmc_command}')
    execution_dir = os.getcwd()
    os.chdir(_radmc_input_path)
    if 'stars' in config_stg:
        if compute_dust_temperature is True:
            logger.info('Computing dust temperature distribution using the stars in the configuration file')
            os.system('radmc3d mctherm setthreads 16')
        else:
            logger.info('Using cached dust temperature distribution')
            shutil.copy(os.path.join(execution_dir, 'model', 'data', 'dust_temperature.dat'),
                        os.path.join('.', 'dust_temperature.dat'))
    os.system(radmc_command)
    os.chdir(execution_dir)
    logger.debug(f'Checking presence of file: {os.path.join(_radmc_input_path, "image.out")}')
+10 −0
Original line number Diff line number Diff line
@@ -20,6 +20,16 @@ grid:
    velocity_gradient: 2
    velocity_gradient_unit: "km/(s pc)"

stars:
    nstars: 1
    rstars: [6.96e10, ]
    mstars: [1.99e36, ]
    star_positions: [[0, 0, 0]]
    star_fluxes: [[-32000,], ]
    nlambdas: 300
    spacing: 'log'
    lambdas_micron_limits: [0.001, 1300]

lines:
    species_to_include: ['e-ch3oh']
    molecular_abundances: {
+51 −8
Original line number Diff line number Diff line
@@ -34,7 +34,8 @@ def write_radmc_input(filename,
                      quantity,
                      grid_metadata,
                      path: Union[None, str] = None,
                      override_defaults: Union[None, dict] = None):
                      override_defaults: Union[None, dict] = None,
                      flatten_style: Union[None, str] = None):
    rt_metadata_default = {
        'iformat': 1,
        'grid_type': grid_metadata['grid_type'],
@@ -45,10 +46,11 @@ def write_radmc_input(filename,
        'ncells_per_axis': ' '.join([str(axis_size) for axis_size in grid_metadata['grid_shape']]),
        'dust_species': 1,  ##
        'active_axes': ' '.join(['1'] * len(grid_metadata['grid_shape'])),
        'continuum_lambdas': 100,  ##
        'continuum_lambdas': grid_metadata['continuum_lambdas'],  ##
    }
    _path = validate_parameter(path, default='.')
    _filename = 'numberdens_mol.inp' if filename.startswith('numberdens_') else filename
    _flatten_style = validate_parameter(flatten_style, default='F')
    with open(os.path.join(_path, filename), "w") as outfile:
        for key in radmc_input_headers[_filename]:
            try:
@@ -56,7 +58,7 @@ def write_radmc_input(filename,
            except (KeyError, TypeError):
                header_line = rt_metadata_default[key]
            outfile.write(f'{header_line}\n')
        outfile.write('\n'.join(quantity.flatten('F').astype(str)))
        outfile.write('\n'.join(quantity.flatten(_flatten_style).astype(str)))


def write_radmc_lines_input(line_config: dict,
@@ -217,12 +219,9 @@ def get_reference_value(grid_metadata: dict,

def write_grid_input_files(profiles: dict,
                           grid_metadata: dict,
                           wavelengths_micron: np.array,
                           path: Union[str, None] = None):
    _path = validate_parameter(path, default=os.path.join('mdl', 'radmc_files'))
    n_wavelengths = 100
    start_wavelength = 5
    end_wavelength = 1300
    wavelengths_micron = np.logspace(np.log10(start_wavelength), np.log10(end_wavelength), n_wavelengths)

    quantity_mapping = {
        'amr_grid.inp': {
@@ -402,6 +401,27 @@ def populate_grid_files(quantity: str,
    )


def write_stellar_input_file(stars_metadata: dict,
                             grid_metadata: dict,
                             path: str,
                             wavelengths_micron: np.array):
    star_properties = [' '.join([str(rstar), str(mstar), str(pos[0]), str(pos[1]), str(pos[2])]) for rstar, mstar, pos in
                       zip(stars_metadata['rstars'], stars_metadata['mstars'], stars_metadata['star_positions'])]
    override_defaults = {
        'iformat': 2,
        'nstars': stars_metadata['nstars'],
        'continuum_lambdas': stars_metadata['nlambdas'],
        'stars_properties': '\n'.join(star_properties),
        'lambdas': ' \n'.join(wavelengths_micron.astype(str)),
    }
    write_radmc_input(filename='stars.inp',
                      path=path,
                      quantity=np.array(stars_metadata['star_fluxes']),
                      grid_metadata=grid_metadata,
                      override_defaults=override_defaults,
                      flatten_style='C')


def main(run_id: str,
         override_config: Union[dict, None] = None,
         path_radmc_files: Union[str, None] = None) -> str:
@@ -415,10 +435,33 @@ def main(run_id: str,
    copy_additional_files(dst_path=input_files_dir)

    grid_metadata = extract_grid_metadata(config=config)
    if 'stars' in config:
        stars_metadata = config['stars']
        if stars_metadata['spacing'] == 'log':
            wavelengths_micron = np.logspace(np.log10(stars_metadata['lambdas_micron_limits'][0]),
                                             np.log10(stars_metadata['lambdas_micron_limits'][1]),
                                             stars_metadata['nlambdas'])
        elif stars_metadata['spacing'] == 'linear':
            wavelengths_micron = np.linspace(stars_metadata['lambdas_micron_limits'][0],
                                             stars_metadata['lambdas_micron_limits'][1],
                                             stars_metadata['nlambdas'])
        else:
            raise (NotImplemented('Spacing not defined. Choose between {linear, log}'))
        write_stellar_input_file(stars_metadata=stars_metadata,
                                 grid_metadata=grid_metadata,
                                 path=input_files_dir,
                                 wavelengths_micron=wavelengths_micron)
    else:
        wavelengths_micron = np.logspace(np.log10(5),
                                         np.log10(1300),
                                         100)
    grid_metadata['continuum_lambdas'] = len(wavelengths_micron)

    profiles = get_profiles(grid_metadata=grid_metadata)
    write_grid_input_files(grid_metadata=grid_metadata,
                           profiles=profiles,
                           path=input_files_dir)
                           path=input_files_dir,
                           wavelengths_micron=wavelengths_micron)
    get_moldata(species_names=config['lines']['species_to_include'],
                logger=logger,
                path=input_files_dir)