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

Amended input for LVG; started refining spherical object creation.

parent b0d62ef3

etl/assets/commons/__init__.py

+6 −2
Original line number Diff line number Diff line
@@ -38,7 +38,7 @@ def compute_power_law_radial_profile( 
        power_law_index: float,

        distance_matrix: np.array,

        value_at_reference: Union[float, None] = None,

        distance_reference: float = 1.0,

        distance_reference: Union[float, u.Quantity] = 1.0,

        fill_reference_pixel: bool = True) -> np.array:

    """

    Compute a power law distribution over the specified grid
@@ -54,7 +54,11 @@ def compute_power_law_radial_profile( 
    _value_at_reference = validate_parameter(value_at_reference, default=central_value)

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

                                distance_matrix) if fill_reference_pixel is True else distance_matrix

    profile = _value_at_reference * (_distance_matrix / distance_reference) ** power_law_index

    try:

        _distance_reference = distance_reference.to(u.cm).value

    except AttributeError:

        _distance_reference = distance_reference

    profile = _value_at_reference * (_distance_matrix / _distance_reference) ** power_law_index

    # If the routine fills the 0-distance point (the centre), it fixes the profile making the central value the maximum

    if fill_reference_pixel is True:

        profile = np.where(profile > central_value, central_value, profile)

etl/main.py

+1 −1
Original line number Diff line number Diff line
@@ -9,7 +9,7 @@ from prs.prs_compute_integrated_fluxes_and_ratios import main as prs_main 
# grid definition

dust_temperatures = np.arange(10, 25, 5)

central_densities = np.array([1e4, 1e5, 1e6])

lines = [2, 3]

lines = [86, 87]



for (tdust, nH2) in product(dust_temperatures, central_densities):

    overrides = {

etl/mdl/config/config.yml

+1 −1
Original line number Diff line number Diff line
@@ -4,6 +4,6 @@ radmc: 
    iline: 2

    central_frequency: 230.538

    frequency_units: GHz

    width_kms: 30

    width_kms: 15

    nchannels: 100

    npix: 200

etl/stg/config/config.yml

+8 −6
Original line number Diff line number Diff line
@@ -4,6 +4,9 @@ grid: 
    central_density: 1e6

    density_unit: "cm^-3"

    density_powerlaw_idx: 0

#    density_value_at_reference: 1e6

#    distance_reference: 0.5

#    distance_reference_unit: 'pc'

    dust_temperature: 15

    dust_temperature_unit: 'K'

    dust_temperature_powerlaw_idx: 0
@@ -11,18 +14,17 @@ grid: 
    microturbulence_unit: 'km/s'

    dim1: {"size":0.1, "size_units": "pc", "shape": 5, "refpix": 2}

    velocity_field: 'solid'

    velocity_gradient: 5

    velocity_gradient: 0.1

    velocity_gradient_unit: "km/s/pc"



lines:

    species_to_include: ['co']

    species_to_include: ['e-ch3oh']

    molecular_abundances: {

                              "co": 1e-7,

                              "o-h2": 0.75,

                              "e-ch3oh": 1e-7,

                              "p-h2": 0.25,

    }

    lines_mode: 2

    collision_partners: ['o-h2', 'p-h2']

    lines_mode: 3

    collision_partners: ['p-h2']



radmc:

    nphotons: 1000000,

etl/stg/stg_radmc_input_generator.py

+34 −5
Original line number Diff line number Diff line
@@ -67,40 +67,67 @@ def get_solid_body_rotation_y(grid_metadata): 
    radial_velocity = (distance_xz * u.cm * gradient).to(u.cm / u.second)

    angle_in_radians = np.arctan(

        abs(grid_metadata['centered_indices'][2, ...] / grid_metadata['centered_indices'][0, ...]))

    velocity_x = radial_velocity[:, np.newaxis, :] * -np.sin(angle_in_radians) * np.sign(grid_metadata['centered_indices'][2, ...])

    velocity_z = radial_velocity[:, np.newaxis, :] * np.cos(angle_in_radians) * np.sign(grid_metadata['centered_indices'][0, ...])

    velocity_x = radial_velocity[:, np.newaxis, :] * -np.sin(angle_in_radians) * np.sign(

        grid_metadata['centered_indices'][2, ...])

    velocity_z = radial_velocity[:, np.newaxis, :] * np.cos(angle_in_radians) * np.sign(

        grid_metadata['centered_indices'][0, ...])

    velocity_x = np.where(np.isnan(velocity_x), 0, velocity_x)

    velocity_z = np.where(np.isnan(velocity_z), 0, velocity_z)

    return velocity_x, velocity_z





def get_profiles(grid_metadata):

    try:

        density_ref = {

            'value_at_reference': (

                (float(grid_metadata['density_value_at_reference']) * u.Unit(grid_metadata['density_unit']))

                .to(u.cm ** -3)).value,

            'distance_reference': (float(grid_metadata['distance_reference']) * u.Unit(

                grid_metadata['distance_reference_unit']).to(u.cm)).value

        }

    except KeyError:

        density_ref = {

            'value_at_reference': None,

            'distance_reference': 1.0

        }

    profiles_mapping = {

        'gas_number_density': {

            'central_value': (float(grid_metadata['central_density']) * u.Unit(grid_metadata['density_unit']))

            .to(u.cm ** -3).value,

            'power_law_index': 0,

            'power_law_index': float(grid_metadata['density_powerlaw_idx']),

            'value_at_reference': density_ref['value_at_reference'],

            'distance_reference': density_ref['distance_reference']

    },

        'dust_temperature': {

            'central_value': float(grid_metadata['dust_temperature']),

            'power_law_index': 0,

            'value_at_reference': None,

            'distance_reference': 1.0

        },

        'microturbulence': {

            'central_value': (float(grid_metadata['microturbulence']) * u.Unit(grid_metadata['microturbulence_unit']))

            .to(u.cm / u.second).value,

            'power_law_index': 0,

            'value_at_reference': None,

            'distance_reference': 1.0

        },

        'velocity_x': {

            'central_value': 0,

            'power_law_index': 0,

            'value_at_reference': None,

            'distance_reference': 1.0

        },

        'velocity_y': {

            'central_value': 0,

            'power_law_index': 0,

            'value_at_reference': None,

            'distance_reference': 1.0

        },

        'velocity_z': {

            'central_value': 0,

            'power_law_index': 0,

            'value_at_reference': None,

            'distance_reference': 1.0

        },

    }
@@ -110,7 +137,9 @@ def get_profiles(grid_metadata): 
        profiles[profile] = compute_power_law_radial_profile(

            central_value=profiles_mapping[profile]['central_value'],

            power_law_index=profiles_mapping[profile]['power_law_index'],

            distance_matrix=grid_metadata['distance_matrix']

            distance_matrix=grid_metadata['distance_matrix'],

            value_at_reference=profiles_mapping[profile]['value_at_reference'],

            distance_reference=profiles_mapping[profile]['distance_reference'],

        )



    if grid_metadata['velocity_field'] == 'solid':