Included number of radmc threads in mdl config file; fixed alignment between grid and image.

    with open(os.path.join('mdl', 'radmc3d_postprocessing.sh'), 'w') as outfile:

        outfile.write(f'cd {_radmc_input_path}\n')

        options_set = get_command_options(config_mdl)

        radmc_command = f'radmc3d image setthreads 4 {" ".join(options_set)}'

        radmc_command = f'radmc3d image {" ".join(options_set)}'

        outfile.write(radmc_command)

    engine = get_pg_engine(logger=logger)

        lines: Union[list, tuple],

        run_id: str,

        session: Session,

        is_isothermal: bool = False) -> Tuple[float, str, str, np.array, np.array]:

        is_isothermal: bool = False) -> Tuple[float, str, str, tuple, np.array, np.array]:

    # Select correct column for filtering

    if is_isothermal is True:

        dust_temperature_column = GridPars.dust_temperature

        dust_temperature_column = GridPars.dust_temperature_at_reference

    # Query database

    results = session.query(GridPars, RatioMaps, GridFiles) \

    results = session.query(GridPars, RatioMaps, GridFiles, ModelPars) \

        .join(GridFiles) \

        .join(StarsPars, isouter=True) \

        .join(ModelPars) \

        grid_values = fitsfile[0].data

    with fits.open(os.path.join('prs', 'fits', 'ratios', results[1].ratio_map_name)) as ratio_fitsfile:

        ratio_values = ratio_fitsfile[0].data

    return results[1].aggregated_ratio, results[1].ratio_map_name, results[2].fits_grid_name, grid_values, ratio_values

    assert (results[3].npix % results[0].grid_shape_1 == 0) and (results[3].npix % results[0].grid_shape_2 == 0)

    return results[1].aggregated_ratio, \

        results[1].ratio_map_name, \

        results[2].fits_grid_name, \

        (int(results[3].npix / results[0].grid_shape_1), int(results[3].npix / results[0].grid_shape_2)), \

        grid_values, \

        ratio_values

def main(run_id: str,

                                   'dust_temperature': dust_temperatures,

                                   'gas_density': central_densities

                               })

        density_correlation = xr.DataArray(np.empty(shape=[len(central_densities)]),

                                           dims=('gas_density'),

                                           coords={

                                               'gas_density': central_densities

                                           })

        # ratio_correlation =

        scatter_data = {}

        for lines in line_pairs:

            scatter_data[str(lines)] = {str(density): list() for density in central_densities}

            fig = plt.figure()

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

                aggregated_ratio = get_aggregated_ratio_from_db(dust_temperature=tdust,

                                                                gas_density=nh2,

                                                                session=session)

                results.loc[tdust, nh2] = aggregated_ratio

                logger.debug(f'The aggregated ratio for lines {lines}, using {nh2}, {tdust} is: {aggregated_ratio}')

                ratio, ratio_fitsname, grid_fitsname, zoom_ratios, density_grid, ratio_values = get_density_distribution(

                    dust_temperature=tdust,

                    gas_density=nh2,

                    lines=lines,

                    session=session,

                    run_id=run_id,

                    is_isothermal=is_isothermal

                )

                density_correlation.loc[nh2] = np.nanmean(np.where(density_grid == 0, np.nan, density_grid))

                avg_density_map = np.nansum(np.where(density_grid == 0, np.nan, density_grid ** 2), axis=2) / \

                                  np.nansum(np.where(density_grid == 0, np.nan, density_grid), axis=2)

                test = np.array([avg_density_map.flatten(), ratio_values[::zoom_ratios[0], ::zoom_ratios[1]].flatten()])

                scatter_data[str(lines)][str(nh2)].append(test)

                plt.scatter(test[0], test[1])

            plt.semilogx()

            plt.xlabel('avg(gas_density)')

            plt.ylabel('simulated ratio')

            plt.savefig(os.path.join('prs', 'output', f'ratio_vs_avg_density_los_{"-".join(lines)}.png'))

            plt.clf()

            results.plot(x='dust_temperature', y='gas_density', yscale='log')

            plt.savefig(os.path.join('prs', 'output', f'ratio_grid_lines_{"-".join(lines)}.png'))

            plt.clf()

if __name__ == '__main__':

    main(run_id='55d05c03-192a-47da-9dcf-c41df1882868', is_isothermal=False)

    # main(run_id='55d05c03-192a-47da-9dcf-c41df1882868', is_isothermal=False)

    main(run_id='ba7fd3ed-5947-4dc6-bcef-38f151f19b77', is_isothermal=False)

    if 'stars' in config:

        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')

            try:

                _threads = config_mdl['radmc_observation']['threads']

            except KeyError:

                _threads = 4

            os.system(f'radmc3d mctherm setthreads {str(_threads)}')

            populate_stars_table(config_stars=config['stars'],

                                 engine=engine,

                                 grid_zipfile=zip_filename,