Commit 045b6a99 authored by Andrea Giannetti's avatar Andrea Giannetti
Browse files

Added scalebar to main figure and script to compute the gravitational velocity

parent 6ccbd4f5

etl/mdl/mdl_compute_gravitational_velo.py

0 → 100644
+39 −0
Changes for etl/mdl/mdl_compute_gravitational_velo.py: 39 added lines, 0 removed lines.
Original line number Diff line number Diff line 
import numpy as np



def gravitational_escape_velocity(mass_solar_masses, radius_parsecs):

    """

    Computes the gravitational escape velocity of a clump.



    Parameters:

    - mass_solar_masses (float): The mass of the clump in solar masses.

    - radius_parsecs (float): The radius of the clump in parsecs.



    Returns:

    - float: The escape velocity in kilometers per second (km/s).

    """

    # Define physical constants

    G = 6.67430e-11  # Gravitational constant in m^3 kg^-1 s^-2

    M_sun = 1.989e30  # Solar mass in kilograms (kg)

    pc = 3.086e16  # Parsec in meters (m)



    # Convert units to SI (meters and kilograms)

    mass_kg = mass_solar_masses * M_sun

    radius_m = radius_parsecs * pc



    # Calculate escape velocity using the formula: v_e = sqrt(2GM/R)

    v_escape_m_s = np.sqrt((G * mass_kg) / radius_m)



    # Convert velocity from m/s to km/s for more common use

    v_escape_km_s = v_escape_m_s / 1000



    return v_escape_km_s



# Example usage:

# Let's say we have a clump with a mass of 1000 solar masses and a radius of 1 parsec.

masses = [1330, 4170, 2050]

radius = 1



for mass in masses:

    velocity = gravitational_escape_velocity(mass, radius)



    print(f"The gravitational escape velocity for a clump with a mass of {mass} solar masses and a radius of {radius} parsecs is: {velocity:.2f} km/s")
 
 No newline at end of file

etl/prs/prs_mom0_dust.py

+17 −3
Changes for etl/prs/prs_mom0_dust.py: 17 added lines, 3 removed lines.
Original line number Diff line number Diff line
@@ -10,6 +10,7 @@ from astropy import units as u 
from astropy.constants import c as speed_of_light

from astropy.wcs import WCS

from astropy.coordinates import SkyCoord

from astropy.visualization.wcsaxes import add_scalebar

import pyspeckit

matplotlib.use('TkAgg')
@@ -30,7 +31,12 @@ levels = { 
    'G28': [1, 2, 4, 8, 12],

    'G14': [1, 8, 20, 50, 100],

}



# Distances in parsecs for each source

distances_pc = {

    'G14': 3100 * u.pc,

    'G28': 4500 * u.pc,

    'G332.96': 4000 * u.pc

}



# Utility functions

def load_yaml(path):
@@ -106,11 +112,19 @@ def plot_swapped(moment0, moment_hdr, continuum, cont_wcs, source, species, core 


    if species == 'CH3OH':

        plt.figure()

        ax = plt.subplot(projection=wcs2d)

        im = ax.imshow(np.log10(continuum), origin='lower', transform=ax.get_transform(cont_wcs))

        ax = plt.subplot(projection=cont_wcs)

        im = ax.imshow(np.log10(continuum), origin='lower')

        plt.colorbar(im, label=r'log(Cont. Intensity [Jy beam$^{-1}$])', shrink=0.9)

        ax.scatter(coord.ra.deg, coord.dec.deg, transform=ax.get_transform('world'), s=50, marker='*',

                   edgecolor='gainsboro', facecolor='none', zorder=2, linewidths=0.7)



        # Add 0.1 pc ruler

        dist = distances_pc[source]

        ang_size_01pc = (0.5 * u.pc / dist * u.rad).to(u.arcsec, equivalencies=u.parallax())

        # The add_scalebar function from astropy.visualization.wcsaxes is a more robust way to do this

        add_scalebar(ax, length=ang_size_01pc, label='0.5 pc',

                     corner='bottom right', frame=False, pad=0.5, color='black')



        ax.set_xlabel('RA')

        ax.set_ylabel('Dec')

        plt.tight_layout()