First commit with RADMC3d input generator.

#!/bin/bash

# Set PYCHARM_HOME, INSTALL_DEPENDENCIES environment variable

PYCHARM_HOME="${PYCHARM_HOME:=~/PycharmProjects}"

PYCHARM_HOME=${PYCHARM_HOME:=$HOME/PycharmProjects}

INSTALL_DEPENDENCIES="${INSTALL_DEPENDENCIES:=false}"

echo $PYCHARM_HOME

cd $PYCHARM_HOME

#!/bin/bash

# Set PYCHARM_HOME, INSTALL_DEPENDENCIES environment variable

PYCHARM_HOME=${PYCHARM_HOME:=$HOME/PycharmProjects}

INSTALL_DEPENDENCIES="${INSTALL_DEPENDENCIES:=false}"

echo $PYCHARM_HOME

cd $PYCHARM_HOME

export RADMC_HOME=$PYCHARM_HOME/radmc3d-2.0

if [[ ! -d "$RADMC_HOME" ]]

then

    echo RADMC not found, cloning into "$(pwd)"

    git clone https://github.com/dullemond/radmc3d-2.0.git

fi

if $INSTALL_DEPENDENCIES

then

    sudo apt-get install gfortran

fi

cd $RADMC_HOME/src

make

make install

# Add this to thee .bashrc configuration file (or equivalent)

export PATH=$HOME/bin:$PATH

export PYTHONPATH=$HOME/bin/python:$PYTHONPATH

cd $RADMC_HOME/python/radmc3dPy

python setup.py install

import yaml

import numpy as np

import os

import urllib.request

from astropy import units as u

from assets.constants import (radmc_grid_map,

                              radmc_coord_system_map,

                              leiden_url_mapping)

def load_config_file(config_file_path: str) -> dict:

    with open(config_file_path) as config_file:

        config = yaml.load(config_file, Loader=yaml.FullLoader)

    return config

def compute_power_law_radial_profile(

        central_value: float,

        power_law_index: float,

        distance_matrix: np.array,

        distance_reference: float = 1.0,

        fill_reference_pixel: bool = True

) -> np.array:

    _distance_matrix = np.where(distance_matrix == 0, 1,

                                distance_matrix) if fill_reference_pixel is True else distance_matrix

    return central_value * (_distance_matrix / distance_reference) ** power_law_index

def get_grid_properties(grid_config: dict,

                        keyword: str) -> list:

    grid_property_list = [grid_config['dim1'][keyword]]

    for axis_idx in range(2):

        try:

            grid_property_list.append(grid_config[f'dim{axis_idx + 2}'][keyword])

        except KeyError:

            grid_property_list.append(grid_config['dim1'][keyword])

    return grid_property_list

def compute_cartesian_coordinate_grid(indices: np.array,

                                      physical_px_size: np.array) -> np.array:

    _physical_px_size_cm = [px_size.to(u.cm).value for px_size in physical_px_size]

    return indices * _physical_px_size_cm

def extract_grid_metadata(config):

    grid_config = config['grid']

    grid_metadata = grid_config.copy()

    grid_metadata['grid_type'] = radmc_grid_map[grid_config['grid_type']]

    grid_metadata['coordinate_system'] = radmc_coord_system_map[grid_config['coordinate_system']]

    grid_metadata['grid_shape'] = get_grid_properties(grid_config=grid_config,

                                                      keyword='steps')

    empty_grid = np.ones(grid_metadata['grid_shape'])

    grid_metadata['refpix'] = get_grid_properties(grid_config=grid_config,

                                                  keyword='refpix')

    grid_metadata['grid_size'] = get_grid_properties(grid_config=grid_config,

                                                     keyword='size')

    grid_metadata['grid_size_units'] = get_grid_properties(grid_config=grid_config,

                                                           keyword='unit')

    grid_metadata['physical_px_size'] = get_physical_px_size(grid_metadata)

    indices = np.indices(empty_grid.shape).T - grid_metadata['refpix']

    grid_metadata['coordinate_grid'] = compute_cartesian_coordinate_grid(indices=indices,

                                                                         physical_px_size=grid_metadata[

                                                                             'physical_px_size'])

    grid_metadata['distance_matrix'] = get_distance_matrix(grid_metadata, indices)

    grid_metadata['grid_edges'] = get_grid_edges(grid_metadata)

    return grid_metadata

def get_grid_edges(grid_metadata):

    grid_edges = []

    for axis_idx in range(len(grid_metadata['grid_shape'])):

        # valid for regular grid

        grid_edges.append(

            (np.arange(grid_metadata['grid_shape'][axis_idx] + 1) - grid_metadata['refpix'][axis_idx] - 0.5) *

            grid_metadata['physical_px_size'][axis_idx])

    # Transposed for consistent flattening

    return np.array(grid_edges).T

def get_distance_matrix(grid_metadata, indices):

    distance_matrix = np.zeros(grid_metadata['grid_shape'])

    for axis_idx in range(len(grid_metadata['grid_shape'])):

        distance_matrix += (indices[:, :, :, axis_idx] * grid_metadata['physical_px_size'][axis_idx].value) ** 2

    return np.sqrt(distance_matrix)

def get_physical_px_size(grid_metadata):

    physical_px_size = []

    for axis_idx in range(len(grid_metadata['grid_shape'])):

        physical_px_size.append(

            (grid_metadata['grid_size'][axis_idx] * u.Unit(grid_metadata['grid_size_units'][axis_idx])).to(u.cm) /

            grid_metadata['grid_shape'][axis_idx])

    return physical_px_size

def get_moldata(species_names: list):

    for species in species_names:

        data = urllib.request.urlopen(leiden_url_mapping[species]).read().decode()

        with open(os.path.join('mdl', 'radmc_input_files', f'molecule_{species}.inp'), 'w') as outfile:

            outfile.writelines(data)

from .physical_constants import mean_molecular_mass

polaris_grid_map = {

    'spherical': 30

}

polaris_data_map = {

    "gas_number_density": 0,

    "dust_temperature": 2,

    "gas_temperature": 3,

    "molecular_abundance_mass_ratio": 17

}

radmc_grid_map = {

    'regular': 0,

    'spherical': 100

}

radmc_coord_system_map = {

    'cartesian': 0,

    'polar': 100

}

radmc_input_headers = {

    'amr_grid.inp': ['iformat', 'grid_type', 'coordinate_system', 'grid_info', 'active_axes', 'ncells_per_axis'],

    'dust_density.inp': ['iformat', 'ncells', 'dust_species'],

    'dust_temperature.dat': ['iformat', 'ncells', 'dust_species'],

    'microturbulence.inp': ['iformat', 'ncells'],

    'gas_velocity.inp': ['iformat', 'ncells'],

    'gas_temperature.inp': ['iformat', 'ncells'],

    'wavelength_micron.inp': ['continuum_lambdas'],

    'numberdens_mol.inp': ['iformat', 'ncells'],

}

leiden_url_mapping = {

    'co': 'https://home.strw.leidenuniv.nl/~moldata/datafiles/co.dat',

    'e-ch3oh': 'https://home.strw.leidenuniv.nl/~moldata/datafiles/e-ch3oh.dat',

    'a-ch3oh': 'https://home.strw.leidenuniv.nl/~moldata/datafiles/ch3oh_a.dat'

}

mean_molecular_mass = 2.34