Added solid body rotation. (b0d62ef3) · Commits · Andrea Giannetti / swiss_army_knife_stable

etl/assets/commons/init.py

+3 −4

Original line number	Diff line number	Diff line
		@@ -112,13 +112,12 @@ def extract_grid_metadata(config: dict) -> dict:

		grid_metadata['physical_px_size'] = get_physical_px_size(grid_metadata)

		indices = get_centered_indices(grid_metadata)

		grid_metadata['coordinate_grid'] = compute_cartesian_coordinate_grid(indices=indices,
		grid_metadata['centered_indices'] = get_centered_indices(grid_metadata)
		grid_metadata['coordinate_grid'] = compute_cartesian_coordinate_grid(indices=grid_metadata['centered_indices'],
		physical_px_size=grid_metadata[
		'physical_px_size'])

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

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

etl/main.py

+1 −2

Original line number	Diff line number	Diff line
		@@ -10,7 +10,7 @@ from prs.prs_compute_integrated_fluxes_and_ratios import main as prs_main
		dust_temperatures = np.arange(10, 25, 5)
		central_densities = np.array([1e4, 1e5, 1e6])
		lines = [2, 3]
		# loop override definition

		for (tdust, nH2) in product(dust_temperatures, central_densities):
		overrides = {
		'grid': {
		@@ -32,4 +32,3 @@ 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')
		# call ETL routines
		No newline at end of file

etl/stg/config/config.yml

+8 −3

Original line number	Diff line number	Diff line
		@@ -10,14 +10,19 @@ grid:
		microturbulence: 5
		microturbulence_unit: 'km/s'
		dim1: {"size":0.1, "size_units": "pc", "shape": 5, "refpix": 2}
		velocity_field: 'solid'
		velocity_gradient: 5
		velocity_gradient_unit: "km/s/pc"

		lines:
		species_to_include: ['co']
		molecular_abundances: {
		"co": 1e-7
		"co": 1e-7,
		"o-h2": 0.75,
		"p-h2": 0.25,
		}
		lines_mode: 1
		collision_partners: []
		lines_mode: 2
		collision_partners: ['o-h2', 'p-h2']

		radmc:
		nphotons: 1000000,

etl/stg/stg_radmc_input_generator.py

+19 −3

Original line number	Diff line number	Diff line
		@@ -61,10 +61,24 @@ def write_radmc_lines_input(line_config: dict,
		outfile.write(f'{coll_partner}\n')


		def get_solid_body_rotation_y(grid_metadata):
		distance_xz = grid_metadata['distance_matrix'][:, grid_metadata['grid_refpix'][1], :]
		gradient = float(grid_metadata['velocity_gradient']) * u.Unit(grid_metadata['velocity_gradient_unit'])
		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 = 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):
		profiles_mapping = {
		'gas_number_density': {
		'central_value': float(grid_metadata['central_density']),
		'central_value': (float(grid_metadata['central_density']) * u.Unit(grid_metadata['density_unit']))
		.to(u.cm ** -3).value,
		'power_law_index': 0,
		},
		'dust_temperature': {
		@@ -99,6 +113,8 @@ def get_profiles(grid_metadata):
		distance_matrix=grid_metadata['distance_matrix']
		)

		if grid_metadata['velocity_field'] == 'solid':
		profiles['velocity_x'], profiles['velocity_z'] = get_solid_body_rotation_y(grid_metadata=grid_metadata)
		profiles['velocity_field'] = np.array([profiles['velocity_x'], profiles['velocity_y'], profiles['velocity_z']])
		return profiles

		@@ -140,7 +156,7 @@ def write_grid_input_files(profiles, grid_metadata):
		def write_molecular_number_density_profiles(profiles: dict,
		line_config: dict,
		grid_metadata: dict):
		for species in line_config['species_to_include']:
		for species in line_config['species_to_include'] + line_config['collision_partners']:
		write_radmc_input(filename=f'numberdens_{species}.inp',
		quantity=profiles['gas_number_density'] * float(line_config['molecular_abundances'][species]),
		path=os.path.join('mdl', 'radmc_files'),
		@@ -165,7 +181,7 @@ def main(override_config: Union[dict, None] = None):
		profiles = get_profiles(grid_metadata=grid_metadata)
		write_grid_input_files(grid_metadata=grid_metadata,
		profiles=profiles)
		get_moldata(species_names=config['lines']['species_to_include'] + config['lines']['collision_partners'])
		get_moldata(species_names=config['lines']['species_to_include'])
		write_radmc_lines_input(line_config=config['lines'],
		path=os.path.join('mdl', 'radmc_files'))
		write_molecular_number_density_profiles(profiles=profiles,

etl/tests/test_commons.py

+2 −4

Original line number	Diff line number	Diff line
		@@ -108,8 +108,6 @@ class Test(TestCase):
		'central_density': '1e6',
		'density_unit': 'cm^-3',
		'dim1': {"size": grid_size, "size_units": "pc", "shape": npix, "refpix": 2},
		'source_distance': '1000',
		'source_distance_unit': 'pc'
		}
		expected_result = [(grid_size * u.pc).to(u.cm) / npix] * 3
		create_test_config(config_dict=config_dict,
		@@ -127,7 +125,7 @@ class Test(TestCase):
		physical_px_size=physical_px_size)
		self.assertTrue(np.array_equal(expected_results, computed_grid))

		def test_compute_cartesian_coordinate_grid(self):
		def test_compute_cartesian_coordinate_grid_2d(self):
		indices = np.indices([3, 4])
		physical_px_size = [1 * u.cm, 1 * u.cm]
		expected_results = indices