Revised figures, added figures requested by referee, and revised infall rate

import numpy as np

import matplotlib.pyplot as plt

from astropy import units as u

from astropy.constants import m_p

# Selection of the source to compute

source = 'G332' # Options: 'G28' or 'G332'

source = 'G332'

if source == 'G28':

    # Computation carried out for G28

    # Data for AGAL028.273-00.167

    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

    # Explicitly defining the DIAMETERS as measured on the sky

    diam_major_proj = 1.0 * u.pc

    diam_minor_proj = 0.5 * u.pc

    depth_los = 0.2 * u.pc # Average of 0.15-0.2 pc to match paper's H_LOS

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

else:

    # Tentative estimate for G332

    # Data for AGAL332.969-00.029

    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

    diam_major_proj = 0.7 * u.pc

    diam_minor_proj = 0.2 * u.pc

    depth_los = 0.2 * u.pc

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

# Inclination angles

# gamma = major axis inclination, theta = minor axis inclination

los_angle_filament = 30 * u.deg

los_angle_sheet = 45 * u.deg

mean_molecular_weight = 2.33 # Kauffmann 2008

# Physical constants and assumptions

mean_molecular_weight = 2.33 # Consistent with high-mass clump studies

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

# \tau_{acc} = \frac{W_p}{\cos \theta} \left( \frac{\cos \theta}{V_{grad}} \frac{\cos \theta}{W_p} \right)

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

# 1. Calculate Deprojected Volume (V)

# We must use the semi-axes (radii) for the area of the ellipse!

# V = pi * (L_p / 2) * (W_p / 2) * H_LOS

radius_major = diam_major_proj / 2

radius_minor = diam_minor_proj / 2

volume = np.pi * radius_major * radius_minor * depth_los

# 2. Calculate Total Mass (M = rho * V)

mass = (volume * mean_molecular_weight * m_p * number_density).to(u.Msun)

# 3. Calculate Corrected Accretion Timescale

# tau = tan(theta) / grad

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

# 4. Calculate Final Infall Rate (M_dot = M / tau)

infall_rate = (mass / accretion_timescale).to(u.Msun / u.yr)

print(f"--- Results for {source} ---")

print(f"Total Mass in Sheet: {mass:.2f}")

print(f"Accretion Timescale: {accretion_timescale:.2e}")

print(f"Infall Rate: {infall_rate / 1e-4:.2f} x 1e-4 Msun/yr")

print(f"Infall Rate: {infall_rate * 1e6:.2f} Msun/Myr")

# --- Plotting the Geometrical Factor ---

# Define range for angles in degrees

gamma_deg = np.linspace(10, 80, 200)

theta_deg = np.linspace(10, 80, 200)

G, T = np.meshgrid(gamma_deg, theta_deg)

# Calculate geometrical factor (Now strictly cot(theta))

factor = 1.0 / np.tan(np.radians(T))

# Set up the plot

plt.figure(figsize=(8, 6))

# Adjusted levels for cot(theta) range in 10-80 degrees

levels = np.linspace(0, 4, 21)

cp = plt.contourf(G, T, factor, levels=levels, cmap='viridis', extend='max')

cbar = plt.colorbar(cp)

cbar.set_label(r'Geometrical Factor $f(\theta) = \cot\theta$', fontsize=12)

# Highlight the 30-60 degree "realistic" region

plt.plot([30, 60, 60, 30, 30], [30, 30, 60, 60, 30], color='red', lw=2, linestyle='--', label=r'Realistic Range ($30^\circ-60^\circ$)')

# Add a marker at the nominal angles (30, 45)

nominal_factor = 1.0 / np.tan(np.radians(45))

plt.plot(30, 45, 'w*', markersize=12, markeredgecolor='black', label=f'Nominal ($30^\circ, 45^\circ$): $f \\approx {nominal_factor:.2f}$')

plt.title('Combined Geometrical Factor Affecting $\\dot{M}$', fontsize=14)

plt.xlabel(r'Major Axis Inclination $\gamma$ (degrees)', fontsize=12)

plt.ylabel(r'Minor Axis Inclination $\theta$ (degrees)', fontsize=12)

plt.legend(loc='upper right')

plt.grid(alpha=0.3)

plt.savefig('geom_factor.png', dpi=300, bbox_inches='tight')

 No newline at end of file

import pickle

import os

import numpy as np

import matplotlib.pyplot as plt

from sklearn.mixture import BayesianGaussianMixture

from scipy.optimize import curve_fit

from mpl_toolkits.axes_grid1 import make_axes_locatable

from assets.constants import config_vis3d as config, source_properties, datadir

fontsize_small = 14

fontsize_medium = 16

fontsize_large = 18

def plane_3d(xy, ra_grad, dec_grad, offset):

    """3D plane function matching your fitting routine."""

    return ra_grad * xy[0] - dec_grad * xy[1] + offset

def plot_2d_clustered_kinematics(source, species):

    print(f"Generating 2D clustered map for {source} - {species}...")

    # 1. Load the data from the pickles generated by your 3D script

    try:

        with open(os.path.join('prs', 'pickles', f'flattened_velocity_fields_{source}_{species}.pickle'), 'rb') as f:

            flattened_velocity_fields = pickle.load(f)

        with open(os.path.join('prs', 'pickles', f'mask_velocity_fields_{source}_{species}.pickle'), 'rb') as f:

            mask = pickle.load(f)

    except FileNotFoundError:

        print(f"Pickle files not found for {source} {species}. Run the 3D script first.")

        return

    # Extract valid data points

    valid_data = flattened_velocity_fields[species][mask[species]]

    velocities = valid_data[:, 0]

    x_pix = valid_data[:, 1]

    y_pix = valid_data[:, 2]

    x_coords = valid_data[:, 3]

    y_coords = valid_data[:, 4]

    # 2. Re-run the Bayesian Gaussian Mixture clustering

    velocity_array = valid_data[:, 0:3]  # [vel, x_pix, y_pix]

    model = BayesianGaussianMixture(n_components=config[source]['n_clusters'],

                                    random_state=config['random_seed'])

    model.fit(velocity_array)

    y_hat = model.predict(velocity_array)

    # Find the surface cluster ID based on the config surface point

    surface_point = config[source]['surface_point']

    sp_mask = (velocity_array[:, 1] == surface_point[0]) & \

              (velocity_array[:, 2] == surface_point[1]) & \

              (np.isclose(velocity_array[:, 0], surface_point[2]))

    if not np.any(sp_mask):

        print(f"Warning: Surface point {surface_point} not found in valid data. Skipping.")

        return

    surface_cluster = model.predict(velocity_array[sp_mask])[0]

    # 3. Create masks for Main (planar) vs Secondary (LOS) components

    min_vel = config[source]['min_velocity']

    if isinstance(min_vel, str):

        if min_vel == '-inf':

            min_vel = -np.inf

        else:

            min_vel = float(min_vel)

    main_mask = (y_hat == surface_cluster) & (velocities > min_vel)

    secondary_mask = ~main_mask

    # 4. Fit the plane to the main component to extract the gradient vector

    pixel_size_pc = source_properties[source]['distance'] * source_properties[source]['pixel_size'] / (

                360 * 3600 / (2 * np.pi))

    x_fit = (x_pix[main_mask] - 255) * pixel_size_pc

    y_fit = (y_pix[main_mask] - 255) * pixel_size_pc

    v_fit = velocities[main_mask]

    guiding_plane, _ = curve_fit(plane_3d, xdata=[x_fit, y_fit], ydata=v_fit)

    ra_grad, dec_grad, offset = guiding_plane

    # 5. Create the 2D Plot

    fig, ax = plt.subplots(figsize=(8, 6))

    # Plot secondary LOS points in transparent grey (squares, pushed to background)

    ax.scatter(x_coords[secondary_mask], y_coords[secondary_mask],

               c='grey', alpha=0.2, s=20, marker='s', zorder=1,

               label='Secondary LOS Components', edgecolors='none')

    ax.tick_params(labelsize=fontsize_small)

    # Plot the main component

    # Center the colormap exactly on the median velocity of the main cluster

    v_mid = np.median(v_fit)

    v_range = max(np.abs(np.percentile(v_fit, 95) - v_mid),

                  np.abs(np.percentile(v_fit, 5) - v_mid))

    # Use square markers and remove black edges to tile it cleanly like an image

    scatter = ax.scatter(x_coords[main_mask], y_coords[main_mask],

                         c=v_fit, cmap='viridis', s=20, marker='s', zorder=2,

                         vmin=v_mid - v_range, vmax=v_mid + v_range,

                         edgecolors='none',

                         label='Main Planar Component')

    divider = make_axes_locatable(ax)

    cax = divider.append_axes("right", size="5%", pad=0.1)

    cbar = plt.colorbar(scatter, cax=cax)

    cbar.set_label(r'$V_{LSR}$ (km s$^{-1}$)', fontsize=fontsize_medium)

    cbar.ax.tick_params(labelsize=fontsize_small)

    # 6. Add the Gradient Direction Arrow

    # Plane eq: V = ra_grad * x - dec_grad * y + offset

    # Vector of steepest velocity increase: (ra_grad, -dec_grad)

    arrow_length = np.max(np.abs(x_coords)) * 0.5  # Scale arrow to the plot size

    grad_norm = np.sqrt(ra_grad ** 2 + (-dec_grad) ** 2)

    dx = (ra_grad / grad_norm) * arrow_length

    dy = (-dec_grad / grad_norm) * arrow_length

    # Start the arrow near the center of the main cluster

    center_x, center_y = np.median(x_fit), np.median(y_fit)

    # We subtract dx/2 to center the arrow visually on the core

    ax.arrow(center_x - dx / 2, center_y - dy / 2, dx, dy,

             color='black', width=0.005, head_width=0.03,

             length_includes_head=True, zorder=5)

    # Add a dummy line for the legend to explain the arrow

    ax.plot([], [], color='black', marker=r'$\rightarrow$', markersize=15,

            linestyle='None', label='Gradient Direction')

    # 7. Formatting

    species_map = {

        'CH3OH': r'CH$_3$OH',

        'N2H+': r'N$_2$H$^+$'

    }

    ax.set_xlabel('dRA [pc]', fontsize=fontsize_medium)

    ax.set_ylabel('dDec [pc]', fontsize=fontsize_medium)

    ax.set_title(f'{source} ({species_map[species]})', fontsize=fontsize_large)

    ax.set_aspect('equal')

    # FORCE identical plot box sizes by setting a fixed physical window for all sources

    plot_radius_pc = 0.8  # Adjust this value so it covers your largest clump

    ax.set_xlim(plot_radius_pc, -plot_radius_pc)  # Standard astronomical convention (RA increases left)

    ax.set_ylim(-plot_radius_pc, plot_radius_pc)

    ax.legend(loc='upper right', fontsize=fontsize_small, framealpha=0.9)

    ax.grid(alpha=0.3)

    plt.tight_layout()

    output_filename = f'projected_2d_kinematics_{source}_{species}.png'

    plt.savefig(output_filename, dpi=300, bbox_inches='tight')

    print(f"Saved plot to {output_filename}\n")

    # plt.show() # Uncomment if you want to view them interactively

if __name__ == '__main__':

    # You can change these to loop over your specific sources/species

    sources = ['G28', 'G332.96', 'G14']

    species_list = ['CH3OH', 'N2H+']

    for src in sources:

        for spc in species_list:

            plot_2d_clustered_kinematics(src, spc)

 No newline at end of file

import matplotlib.pyplot as plt

import numpy as np

fontsize_small = 14

fontsize_medium = 16

fontsize_large = 18

def create_clean_2d_schematic():

    # Set math font to match LaTeX style

    plt.rcParams['mathtext.fontset'] = 'cm'

    fig, axs = plt.subplots(2, 2, figsize=(14, 12))

    axA, axB, axC, axD = axs[0, 0], axs[0, 1], axs[1, 0], axs[1, 1]

    # --- Parameters ---

    L, W = 1.0, 0.5

    gamma = np.radians(30)

    theta = np.radians(45)

    phi = np.radians(45)

    t_int = 0.15  # Intrinsic thickness (H_int)

    # ==========================================

    # Panel A: Physical Geometry (Face-on)

    # ==========================================

    axA.set_title("A) Physical Face-on View", fontsize=fontsize_large)

    t_pts = np.linspace(0, 2 * np.pi, 100)

    axA.plot((L / 2) * np.cos(t_pts), (W / 2) * np.sin(t_pts), color='darkgreen', lw=2)

    axA.fill((L / 2) * np.cos(t_pts), (W / 2) * np.sin(t_pts), color='darkgreen', alpha=0.1)

    axA.annotate('', xy=(L / 2, 0), xytext=(-L / 2, 0), arrowprops=dict(arrowstyle='<->', color='black'))

    axA.text(-0.24, 0.04, r'Physical Major Axis ($L$)', ha='center', fontsize=fontsize_small)

    axA.annotate('', xy=(0, W / 2), xytext=(0, -W / 2), arrowprops=dict(arrowstyle='<->', color='black'))

    axA.text(0.04, 0, r'Physical Minor Axis ($W$)', rotation=90, va='center', fontsize=fontsize_small)

    axA.set_aspect('equal');

    axA.set_xlim(-0.6, 0.6);

    axA.set_ylim(-0.6, 0.6);

    axA.axis('off')

    # Helper function to draw 2D cross-section (slab)

    def draw_2d_slab(ax, length, angle):

        # Create unrotated rectangle centered at origin

        xs = np.array([-length / 2, length / 2, length / 2, -length / 2, -length / 2])

        ys = np.array([t_int / 2, t_int / 2, -t_int / 2, -t_int / 2, t_int / 2])

        # Rotate by angle

        x_rot = xs * np.cos(angle) - ys * np.sin(angle)

        y_rot = xs * np.sin(angle) + ys * np.cos(angle)

        ax.fill(x_rot, y_rot, color='darkgreen', alpha=0.2)

        # Top and bottom parallel lines

        ax.plot(x_rot[:2], y_rot[:2], color='darkgreen', lw=2.5)  # Top

        ax.plot(x_rot[2:4], y_rot[2:4], color='darkgreen', lw=2.5)  # Bottom

        return x_rot, y_rot

    # ==========================================

    # Panel B: Side View along Major Axis (\gamma)

    # ==========================================

    axB.set_title(r"B) Side View: Major Axis Inclination ($\gamma$)", fontsize=fontsize_large)

    axB.axhline(0, color='blue', lw=2)

    axB.text(0.4, 0.02, "Sky Plane", color='blue', fontsize=fontsize_small)

    x_rotB, y_rotB = draw_2d_slab(axB, L, gamma)

    Lp = L * np.cos(gamma)

    # Projections

    axB.plot([Lp / 2, Lp / 2], [0, (L / 2) * np.sin(gamma)], color='gray', linestyle='--')

    axB.plot([-Lp / 2, -Lp / 2], [0, -(L / 2) * np.sin(gamma)], color='gray', linestyle='--')

    axB.annotate('', xy=(Lp / 2, -0.05), xytext=(-Lp / 2, -0.05), arrowprops=dict(arrowstyle='<->', color='red'))

    axB.text(0, -0.12, r'Projected Length ($L_p = L \cos \gamma$)', color='red', ha='center', fontsize=fontsize_small)

    # Gamma angle

    arc_gamma = np.linspace(0, gamma, 20)

    axB.plot(0.2 * np.cos(arc_gamma), 0.2 * np.sin(arc_gamma), color='purple')

    axB.text(0.22, 0.05, r'$\gamma$', color='purple', fontsize=fontsize_medium)

    # LOS Arrow

    axB.annotate('', xy=(-0.4, 0.4), xytext=(-0.4, 0.1), arrowprops=dict(arrowstyle='<-', color='black', lw=2))

    axB.text(-0.43, 0.4, "Line of Sight", rotation=90, va='top', fontsize=fontsize_small)

    axB.set_aspect('equal');

    axB.set_xlim(-0.6, 0.6);

    axB.set_ylim(-0.3, 0.5);

    axB.axis('off')

    # ==========================================

    # Panel C: Side View along Minor Axis (\theta) & Thickness

    # ==========================================

    axC.set_title(r"C) Side View: Minor Axis ($\theta$) & Thickness", fontsize=fontsize_large)

    axC.axhline(0, color='blue', lw=2)

    axC.text(0.4, 0.02, "Sky Plane", color='blue', fontsize=fontsize_small)

    x_rotC, y_rotC = draw_2d_slab(axC, W, theta)

    Wp = W * np.cos(theta)

    # Projections

    axC.plot([Wp / 2, Wp / 2], [0, (W / 2) * np.sin(theta)], color='gray', linestyle='--')

    axC.plot([-Wp / 2, -Wp / 2], [0, -(W / 2) * np.sin(theta)], color='gray', linestyle='--')

    axC.annotate('', xy=(Wp / 2, -0.1), xytext=(-Wp / 2, -0.1), arrowprops=dict(arrowstyle='<->', color='red'))

    axC.text(0, -0.16, r'Projected Width ($W_p = W \cos \theta$)', color='red', ha='center', fontsize=fontsize_small)

    # Theta angle

    arc_theta = np.linspace(0, theta, 20)

    axC.plot(0.15 * np.cos(arc_theta), 0.15 * np.sin(arc_theta), color='purple')

    axC.text(0.14, 0.07, r'$\theta$', color='purple', fontsize=fontsize_medium)

    # --- Thickness Annotations ---

    # H_int (Perpendicular)

    axC.annotate('', xy=(-0.08 * np.sin(theta), 0.08 * np.cos(theta)),

                 xytext=((t_int - 0.07) * np.sin(theta), -(t_int - 0.07) * np.cos(theta)),

                 arrowprops=dict(arrowstyle='<->', color='purple', lw=2))

    axC.text(-0.12, 0.08, r'$H_{\mathrm{int}}$', color='purple', fontsize=fontsize_medium)

    # H_LOS (Vertical)

    H_LOS = t_int / np.cos(theta)

    axC.plot([0, 0], [-H_LOS / 2, H_LOS / 2], color='black', linestyle=':', lw=2)

    axC.annotate('', xy=(0, H_LOS / 2), xytext=(0, -H_LOS / 2),

                 arrowprops=dict(arrowstyle='<->', color='black', lw=1.5))

    axC.text(-0.10, -0.03, r'$H_{\mathrm{LOS}}$', va='center', fontsize=fontsize_medium)

    # --- Velocity Vectors ---

    v_len = 0.25

    x_edge = (W / 2) * np.cos(theta)

    y_edge = (W / 2) * np.sin(theta)

    # v_inf (parallel to slab)

    axC.annotate('', xy=(x_edge - v_len * np.cos(theta), y_edge - v_len * np.sin(theta)), xytext=(x_edge, y_edge),

                 arrowprops=dict(arrowstyle='->', color='darkgreen', lw=2))

    axC.text(x_edge - 0.17, y_edge - 0.1, r"$v_{\mathrm{inf}}$", color='darkgreen', fontsize=fontsize_medium)

    # v_rad (vertical projection of v_inf)

    axC.annotate('', xy=(x_edge, y_edge - v_len * np.sin(theta)), xytext=(x_edge, y_edge),

                 arrowprops=dict(arrowstyle='->', color='orange', lw=2))

    axC.text(x_edge + 0.02, y_edge - 0.15, r"$v_{\mathrm{rad}}$", color='orange', fontsize=fontsize_medium)

    axC.set_aspect('equal');

    axC.set_xlim(-0.6, 0.6);

    axC.set_ylim(-0.3, 0.5);

    axC.axis('off')

    # ==========================================

    # Panel D: Sky View (Position Angle \phi)

    # ==========================================

    axD.set_title(r"D) Sky View: Position Angle ($\phi$)", fontsize=fontsize_large)

    xD = (Lp / 2) * np.cos(t_pts)

    yD = (Wp / 2) * np.sin(t_pts)

    x_rot = xD * np.cos(phi) - yD * np.sin(phi)

    y_rot = xD * np.sin(phi) + yD * np.cos(phi)

    axD.plot(x_rot, y_rot, color='red', lw=2)

    axD.fill(x_rot, y_rot, color='red', alpha=0.1)

    axD.axhline(0, color='black', alpha=0.3, linestyle=':')

    axD.axvline(0, color='black', alpha=0.3, linestyle=':')

    axD.text(0.5, 0.02, "x (RA)", alpha=0.7, fontsize=fontsize_small)

    axD.text(0.02, 0.5, "y (Dec)", alpha=0.7, fontsize=fontsize_small)

    maj_x = np.array([-Lp / 2, Lp / 2]) * np.cos(phi)

    maj_y = np.array([-Lp / 2, Lp / 2]) * np.sin(phi)

    axD.plot(maj_x, maj_y, color='black', linestyle='--')

    arc_phi = np.linspace(0, phi, 20)

    axD.plot(0.2 * np.cos(arc_phi), 0.2 * np.sin(arc_phi), color='purple')

    axD.text(0.22, 0.1, r'$\phi$', color='purple', fontsize=fontsize_medium)