Minor changes to scripts; added core-to-core velocity variations and Jeans mass

G332.96: {

  'vlsr': -67.6,

  tdust: 12.2,

  mass: 1330,

  pixel_size: 0.2,

  distance: 4000,

  rotation_angle: 45,

  reff: 0.42,

  CH3OH: {

    'data': 'G332.96_CH3OH_3mm_feather.fits',

    'rms': 0.0034,

}

G28: {

  'vlsr': 79.6,

  tdust: 11.8,

#  mass: 4170,

  mass: 2250,

  pixel_size: 0.2,

  distance: 4500,

  rotation_angle: -41.3785,

#  reff: 0.6,

  reff: 0.4,

  CH3OH: {

    'data': 'G28_CH3OH_3mm_feather.fits',

    'linewidth_guess_km_s': 1.5,

}

G14: {

  'vlsr': 40.5,

  tdust: 16.6,

  mass: 2100,

  pixel_size: 0.2,

  distance: 3100,

  rotation_angle: -45,

  reff: 0.26,

  CH3OH: {

    'data': 'G14_CH3OH_3mm_feather.fits',

    'rms': 0.0032,

    'frequency_limits_ghz': [ 93.1627, 93.184 ],

    'max_linewidth_km_s': 3 }

}

G12: {

  'vlsr': 34.4,

  CH3OH: {

    'data': 'G12_CH3OH_3mm_ACA_MA.fits',

    'linewidth_guess_km_s': 3,

    'rms': 0.002,

    'frequency_limits_ghz': [ 96.692, 96.786 ] }

}

 No newline at end of file

import yaml

import os

import numpy as np

from assets.constants import source_properties

def isothermal_sound_speed(T, mu=2.37):

    """

    Calculate the isothermal sound speed in km/s.

    Parameters:

    T : float

        Temperature in Kelvin.

    mu : float

        Mean molecular weight, defaults to 2.33 for molecular gas.

    Returns:

    c_s : float

        Sound speed in km/s.

    """

    # Physical constants

    k_B = 1.380649e-23  # J/K

    m_p = 1.6726219e-27  # kg

    # Calculate sound speed in m/s

    c_s_m = np.sqrt(k_B * T / (mu * m_p))

    # Convert m/s to km/s

    c_s_km_s = c_s_m / 1000

    return c_s_km_s

def density_from_mass_radius(M_solar, radius_cm, mu=2.37):

    """

    Calculate particle number density in particles per cubic centimeter.

    Parameters:

    M_solar : float

        Total mass in solar masses.

    radius_cm : float

        Radius in centimeters.

    mu : float

        Mean molecular weight (default is 2.33 for molecular hydrogen gas).

    Returns:

    n_H2 : float

        Number density in particles per cubic centimeter.

    """

    # Constants

    M_sun = 1.9885e33  # grams

    m_p = 1.6726219e-24  # grams (proton mass)

    # Convert mass to grams

    M = M_solar * M_sun

    # Calculate volume in cm^3

    volume = (4 / 3) * np.pi * radius_cm ** 3

    # Calculate number density

    n_H2 = (M / volume) / (mu * m_p)

    return n_H2

pc_in_cm = 3.0857e18

# Define the parameters for each clump

clumps = {

    "G332.96": {

        "sound_speed": 0.22,  # km/s

        "density": 4.8e4       # cm^-3

    },

    "G28": {

        "sound_speed": 0.22,

        "density": 5.3e4

    },

    "G14": {

        "sound_speed": 0.258,

        "density": 10.8e4

    }

}

def calculate_jeans_mass(sound_speed, density):

    """

    Calculate Jeans mass using your formula:

    M_J = 1.578 * (sound_speed/0.188)^3 * (density/10^5)^(-1/2)

    """

    mass = 1.578 * (sound_speed / 0.188) ** 3 * (density / 1e5) ** (-0.5)

    # mass = 0.6285 * (temperature / 10) ** (3/2) * (density / 1e5) ** (-0.5) # From Palau et al 2015

    return mass

# Calculate and output results

for name, params in source_properties.items():

    with open(os.path.join('mdl', 'inferred_densities', f'clump_densities_{name}.yml'), 'r') as f:

        data = yaml.safe_load(f)

    # s_speed = isothermal_sound_speed(T=source_properties[name]['tdust'])

    s_speed = clumps[name]['sound_speed']

    dens = 10 ** data['clump']['best_fit'][0]

    dens = clumps[name]['density']

    dens = density_from_mass_radius(source_properties[name]['mass'], source_properties[name]['reff'] * pc_in_cm)

    M_j = calculate_jeans_mass(s_speed, dens)

    # M_j = calculate_jeans_mass(source_properties[name]['tdust'], dens)

    total_mass = source_properties[name]['mass']

    if total_mass:

        number_jeans_masses = total_mass / M_j

        print(f"Clump {name}:")

        print(f"  Jeans mass (M_J) = {M_j:.2f} M_sun")

        print(f"  Td = {source_properties[name]['tdust']:.2f} K")

        print(f"  Sound speed = {s_speed:.2f} km/s")

        print(f"  Density = {dens:.2f} cm^(-3)")

        print(f"  Number of Jeans masses = {number_jeans_masses:.0f}\n")

    else:

        print(f"Clump {name}:")

        print(f"  Jeans mass (M_J) = {M_j:.2f} M_sun")

        print(f"  Td = {source_properties[name]['tdust']:.2f} K")

        print(f"  Sound speed = {s_speed:.2f} km/s")

        print(f"  Density = {dens:.2f} cm^(-3)")

        print("  Total mass not specified, cannot compute number of Jeans masses.\n")

 No newline at end of file

import numpy as np

from assets.constants import source_properties

def gravitational_escape_velocity(mass_solar_masses, radius_parsecs):

    """

    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)

    v_escape_m_s = np.sqrt(2 * (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

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

for source in source_properties:

    velocity = gravitational_escape_velocity(source_properties[source]['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

    print(f"The gravitational escape velocity for a clump with a mass of {source_properties[source]['mass']} solar masses and a radius of {radius} parsecs is: {velocity:.2f} km/s")

 No newline at end of file

from astropy.constants import m_p

source = 'G332'

if source == 'G28':

    # Computation carried out for G28

    velocity_gradient = 2.96 * u.km / u.s / u.pc

    # These are radii

    proj_length = 0.5 * u.pc

    proj_width = 0.25 * u.pc

    depth = 0.1 * u.pc

    # Clump density

    number_density = 5.3e4 * u.cm**-3

else:

    # Tentative estimate for G332

    velocity_gradient = 3.6 * u.km / u.s / u.pc

    proj_length = 0.35 * u.pc

    proj_width = 0.1 * u.pc

    depth = 0.1 * u.pc

    number_density = 4.8e4 * u.cm**-3

los_angle_filament = 30 * u.deg

los_angle_sheet = 45 * u.deg

mean_molecular_weight = 2.33 # Kauffmann 2008

number_density = 5e4 * u.cm**-3

volume = np.pi * proj_length / np.cos(los_angle_filament) * proj_width / np.cos(los_angle_sheet) * depth

mass = volume * mean_molecular_weight * m_p * number_density

accretion_timescale = (np.cos(los_angle_sheet) / velocity_gradient).to(u.yr)

infall_rate = (mass / accretion_timescale).to(u.Msun / u.yr) # Same accretion on both sides

# Along the filament it would be around one dex lower, as the velocity gradient is 0.2-0.3 km/s/pc

print(f'{infall_rate / 1e-4} x 1e-4 Msun/yr or {infall_rate * 1e6} Msun/Myr')

 No newline at end of file

import os

import yaml

import pandas as pd

import glob

from assets.constants import source_properties

def analyze_core_velocities(directory_path):

    """

    Reads all YAML files from a specified directory, extracts core velocity

    data, and computes the standard deviation of velocities grouped by source.

    Args:

        directory_path (str): The path to the directory containing the YAML files.

    Returns:

        pandas.Series: A pandas Series containing the standard deviation of

                       the 'vlsr' for each source file. Returns None if the

                       directory is not found or is empty.

    """

    if not os.path.isdir(directory_path):

        print(f"Error: Directory not found at '{directory_path}'")

        return None

    all_records = []

    sources = []

    # Iterate over each file in the given directory

    for filename in glob.glob(os.path.join(directory_path, 'densities_*.yml')):

        # Use the filename (without extension) as the source identifier

        source_name = os.path.splitext(os.path.basename(filename))[0].replace('densities_', '')

        sources.append(source_name)

        try:

            with open(filename, 'r') as f:

                data = yaml.safe_load(f)

            if not data:

                print(f"Warning: File '{filename}' is empty. Skipping.")

                continue

            # Extract records from the YAML data

            for component, values in data.items():

                # Ensure 'vlsr' key exists and is not empty

                if 'vlsr' in values and values['vlsr']:

                    for velocity in values['vlsr']:

                        record = {

                            'source': source_name,

                            'component': component,

                            'vlsr': velocity  # vlsr is a list with one value

                        }

                    if component.startswith('core'):

                        all_records.append(record)

                else:

                    print(f"Warning: 'vlsr' key missing or empty for component '{component}' in '{filename}'.")

        except yaml.YAMLError as e:

            print(f"Error parsing YAML file '{filename}': {e}")

        except Exception as e:

            print(f"An unexpected error occurred with file '{filename}': {e}")

    if not all_records:

        print("No valid data could be processed from the YAML files.")

        return None

    # Convert the list of records into a pandas DataFrame

    df = pd.DataFrame(all_records)

    # --- Per-source Velocity Filtering ---

    # Create a list to hold the filtered data for each source

    filtered_dataframes = []

    # Iterate through each unique source found in the data

    for source_name in df['source'].unique():

        if source_name in source_properties:

            # Isolate the data for the current source

            source_df = df[df['source'] == source_name]

            # Get the central velocity for this source from our constants

            center_velocity = source_properties[source_name]['vlsr']

            # Define the velocity range for filtering (+/- 5 km/s)

            min_vlsr = center_velocity - 5

            max_vlsr = center_velocity + 5

            # Apply the filter and keep rows within the velocity range

            filtered_source_df = source_df[source_df['vlsr'].between(min_vlsr, max_vlsr)]

            filtered_dataframes.append(filtered_source_df)

        else:

            print(f"Warning: No velocity properties found for source '{source_name}'. It will be excluded from filtered analysis.")

    if not filtered_dataframes:

        print("No data remains after applying velocity filters.")

        return None

    # Combine the filtered dataframes back into a single dataframe

    df_filtered = pd.concat(filtered_dataframes)

    # Group by the 'source' and compute the standard deviation of 'vlsr' on the filtered data

    velocity_std_dev = df_filtered.groupby('source')['vlsr'].agg(['count', 'std'])

    return velocity_std_dev

if __name__ == '__main__':

    # Run the analysis function on the directory

    std_dev_results = analyze_core_velocities(os.path.join('mdl', 'inferred_densities'))

    # --- Output ---

    if std_dev_results is not None:

        print("Standard Deviation of Core Velocities (vlsr) by Source:")

        print(std_dev_results)