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

Added script to visualize results.

parent 47c18865
Loading
Loading
Loading
Loading
+14 −0
Original line number Diff line number Diff line
@@ -55,6 +55,20 @@ def load_config_file(config_file_path: str,
    return config


def parse_grid_overrides(par_name: str,
                         config: dict) -> np.array:
    config_overrides = config['overrides']
    if config_overrides[f'{par_name}_grid_type'] == 'linear':
        grid_values = np.arange(float(config_overrides[f'{par_name}_limits'][0]),
                                float(config_overrides[f'{par_name}_limits'][1]),
                                float(config_overrides[f'{par_name}_step']))
    else:
        grid_values = 10 ** np.arange(np.log10(float(config_overrides[f'{par_name}_limits'][0])),
                                      np.log10(float(config_overrides[f'{par_name}_limits'][1])),
                                      np.log10(float(config_overrides[f'{par_name}_step'])))
    return grid_values


def compute_power_law_radial_profile(
        central_value: float,
        power_law_index: float,
+1 −1
Original line number Diff line number Diff line
@@ -6,4 +6,4 @@ overrides:
    gas_density_limits: [1e4, 1e8]
    gas_density_step: 10
    gas_density_unit: cm^-3
    lines_to_process: [87, 86]
    lines_to_process: ['87', '86']
+33 −41
Original line number Diff line number Diff line
import os
import numpy as np
from itertools import product
from assets.commons import load_config_file
from assets.commons import (load_config_file,
                            parse_grid_overrides)
from stg.stg_radmc_input_generator import main as stg_main
from mdl.mdl_prepare_radmc_command import main as create_radmc_script
from mdl.mdl_save_results_as_fits import main as save_results
from prs.prs_compute_integrated_fluxes_and_ratios import main as prs_main

config = load_config_file(os.path.join('config', 'config.yml'))


def parse_grid_overrides(par_name: str,
                         config: dict) -> np.array:
    config_overrides = config['overrides']
    if config_overrides[f'{par_name}_grid_type'] == 'linear':
        grid_values = np.arange(float(config_overrides[f'{par_name}_limits'][0]),
                                float(config_overrides[f'{par_name}_limits'][1]),
                                float(config_overrides[f'{par_name}_step']))
    else:
        grid_values = 10 ** np.arange(np.log10(float(config_overrides[f'{par_name}_limits'][0])),
                                      np.log10(float(config_overrides[f'{par_name}_limits'][1])),
                                      np.log10(float(config_overrides[f'{par_name}_step'])))
    return grid_values

def build_model_grid():
    config = load_config_file(os.path.join('config', 'config.yml'))

    # grid definition
    dust_temperatures = parse_grid_overrides(par_name='dust_temperature',
                                             config=config)
    central_densities = parse_grid_overrides(par_name='gas_density',
                                             config=config)
lines = config['lines_to_process']
    lines = config['overrides']['lines_to_process']

    for (tdust, nH2) in product(dust_temperatures, central_densities):
        overrides = {
@@ -52,3 +40,7 @@ for (tdust, nH2) in product(dust_temperatures, central_densities):
        prs_main(cube_fits1=f'td_{tdust}_nh2_{nH2}_line_{lines[0]}.fits',
                 cube_fits2=f'td_{tdust}_nh2_{nH2}_line_{lines[1]}.fits',
                 ratio_fits_name=f'ratio_td_{tdust}_nh2_{nH2}_lines_{lines[0]}-{lines[1]}.fits')


if __name__ == '__main__':
    build_model_grid()
 No newline at end of file
+16 −0
Original line number Diff line number Diff line
import pyspeckit
from astropy.io import fits


fitsfile = '/home/andrea/PycharmProjects/swiss_army_knife/etl/prs/fits/cubes/td_15_nh2_100000.0_line_86.fits'
cube = fits.open(fitsfile)

data = cube[0].data
header = cube[0].header

header['CUNIT3'] = 'Hz'

fits.writeto('fixed_cube.fits', data=data, header=header, overwrite=True)

psk_cube = pyspeckit.Cube('fixed_cube.fits')
psk_cube.mapplot()
+51 −0
Original line number Diff line number Diff line
import numpy as np
import os
import xarray as xr
from typing import Union
from itertools import product
from astropy.io import fits
from assets.commons import (load_config_file,
                            parse_grid_overrides,
                            setup_logger)


def get_fitsfile_name(dust_temperature: float,
                      gas_density: float,
                      lines: Union[list, tuple]) -> str:
    return f'ratio_td_{str(int(dust_temperature))}_nh2_{str(round(gas_density, 1))}_lines_{"-".join(lines)}.fits'


def aggregate_image(data:np.array):
    return np.nanmean(data)


logger = setup_logger(name='PRS - INSPECT')
config = load_config_file(os.path.join('config', 'config.yml'))

# grid definition
dust_temperatures = parse_grid_overrides(par_name='dust_temperature',
                                         config=config)
central_densities = parse_grid_overrides(par_name='gas_density',
                                         config=config)
lines = config['overrides']['lines_to_process']

results = xr.DataArray(np.empty(shape=[len(dust_temperatures), len(central_densities)]),
                       dims=('dust_temperature', 'gas_density'),
                       coords={
                           'dust_temperature': dust_temperatures,
                           'gas_density': central_densities
                       })
for (tdust, nH2) in product(dust_temperatures, central_densities):
    overrides = {
        'grid': {
            'dust_temperature': tdust,
            'central_density': nH2,
        }
    }
    filename = get_fitsfile_name(dust_temperature=tdust, gas_density=nH2, lines=lines)
    logger.debug(get_fitsfile_name(dust_temperature=tdust, gas_density=nH2, lines=lines))
    hdu = fits.open(os.path.join('prs', 'fits', 'ratios', filename))
    aggregated_ratio = aggregate_image(hdu[0].data)
    results.loc[tdust, nH2] = aggregated_ratio
    logger.debug(aggregated_ratio)
results.plot(x='dust_temperature', y='gas_density', yscale='log')
Loading