# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# ## Introduction to the Interstellar Medium
# ### <NAME>

# ### Figures 9.1, 9.2, and 9.3: Bonnor-Ebert profiles and mass

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate
import scipy.interpolate as interpolate
# %matplotlib inline

def lane_emden_integrate(x):
    # solve Lane-Emden equation
    nsteps = x.size
    y = np.zeros(nsteps)
    yp = np.zeros(nsteps)
    yp2 = np.zeros(nsteps)

    # initial condition on d2y/dx2
    # (verified that solutions are insensitive to this beyond x = 2)
    yp2[0] = 1/3

    # integrate outwards step by step
    # (logarithmic steps)
    for i in np.arange(1,nsteps):
        dx = x[i] - x[i-1]
        y[i] = y[i-1] + yp[i-1]*dx + yp2[i-1]*dx**2/2
        yp[i] = yp[i-1] + yp2[i-1]*dx
        yp2[i] = np.exp(-y[i]) - 2*yp[i]/x[i]

    return(y,yp)


def plot_profiles():
    # plot Bonnor-Ebert density profile
    nsteps = 1000
    xmax = 1e4
    x = np.logspace(-2, np.log10(xmax), nsteps)
    y,yp = lane_emden_integrate(x)

    # scale for various physical parameters
    r0 = 1.243e3        # radial scale factor in pc

    fig = plt.figure(figsize=(6,4))
    ax = fig.add_subplot(111)

    ax.set_xlim(0.002,1.0)
    ax.set_ylim(1e8,1e13)
    ax.set_xscale('log')
    ax.set_yscale('log')
    ax.set_xlabel(r'${\rm Radius}\ {\rm (pc)}$', fontsize=14)
    ax.set_ylabel(r'${\rm H_2\ density}\ {\rm (m^{-3})}$', fontsize=14)

    T = 10               # isothermal temperature (K)
    n_ext = 8e9/T        # lower density limit from pressure equilibrium

    n0 = np.array([1,0.2,5,25,125])*14.2*n_ext
    ls = ['-','-','--','--','--']
    lw = [2,2,2,2,2]
    alpha = [1,0.3,0.3,0.3,0.3]

    for i in range(len(n0)):
        r = x * r0 * np.sqrt(T/n0[i])
        n = n0[i] / np.exp(y)
        if i == 0:
            ax.plot(r, n, linestyle=ls[i], color='k', lw=lw[i], alpha=alpha[i], label='Critical')
        else:
            ax.plot(r, n, linestyle=ls[i], color='k', lw=lw[i], alpha=alpha[i])

    # singular isothermal sphere
    r = np.logspace(-3,1,2)
    ax.plot(r,3.09e6*T/r**2, 'k--', lw=2, label='Singular')

    ax.plot([0.2,10], [n_ext,n_ext], 'k:', label='Ambient density')
    ax.legend()

    ax.text(0.0027, 2.7e9, 'Stable', fontsize=10)
    ax.text(0.0027, 6.9e10, 'Unstable', fontsize=10)

    x_labels = ['0.01','0.1','1']
    x_loc = np.array([float(x) for x in x_labels])
    ax.set_xticks(x_loc)
    ax.set_xticklabels(x_labels)

    fig.tight_layout() 
    plt.savefig('bonnor_ebert_profiles.pdf')


def plot_mass():
    # plot mass for given P_ext
    nsteps = 10000
    xmax = 1e4
    x = np.logspace(-4, np.log10(xmax), nsteps)
    y,yp = lane_emden_integrate(x)

    T = 10              # isothermal temperature (K)
    r0 = 1.243e3        # radial scale factor in pc
    n_ext = 8e9/T       # exterior density in m-3

    n0 = np.logspace(np.log10(1.1*n_ext),12,300)
    ndens = n0.size
    r_ext = np.zeros(ndens)
    m_ext = np.zeros(ndens)
    m_tot = np.zeros(ndens)
    for i in range(ndens):
        y_ext = np.log(n0[i]/n_ext)
        j = np.where(np.abs(y/y_ext - 1) < 0.1)[0]
        ycubic = interpolate.UnivariateSpline(x[j],y[j]-y_ext)
        x_ext = ycubic.roots()[0]
        k = np.where(x < x_ext)[0]
        m_ext[i] = 1.19e3 * integrate.simps(x[k]**2 / np.exp(y[k]), x[k]) * np.sqrt(T**3/n0[i])

    # max pressure contrast
    Pratio = n0/n_ext
    imax = m_ext.argmax()
    m_ext_max = m_ext[imax]
    Pratio_max = Pratio[imax]

    fig = plt.figure(figsize=(6,4))
    ax1 = fig.add_subplot(111)

    ax1.set_xlim(1,3e2)
    ax1.set_xscale('log')
    #ax1.set_yscale('log')
    ax1.set_xlabel(r'$\rho_{\rm cen}/\rho_{\rm amb}$', fontsize=14)

    ax1.set_ylim(0,6.5)
    #ax1.set_yscale('log')
    ax1.set_ylabel(r'${\rm Mass}\ (M_\odot)$', fontsize=14)
    #mplot = ax1.plot(Pratio, m_ext, 'k-', lw=3, label='Mass')

    ax1.plot(Pratio[0:imax-1], m_ext[0:imax-1], 'k-', lw=2, alpha=0.3, zorder=99)
    ax1.plot(Pratio[imax+1:], m_ext[imax+1:], 'k--', lw=2, alpha=0.3, zorder=99)
    ax1.plot(Pratio_max, m_ext_max, 'ko', markersize=4, zorder=999)
    ax1.text(2.05, 3.2, 'Stable', fontsize=12, rotation=58, backgroundcolor='white', zorder=2)
    ax1.text(50, 4.6, 'Unstable', fontsize=12, rotation=-21, zorder=2)
    ax1.text(9.5, m_ext_max+0.15, r'$M_{\rm BE}$', fontsize=12)

    # SIS
    m_SIS = 1.06 * np.sqrt(1e10/n_ext) * (T/10)**1.5
    ax1.plot([1,300], [m_SIS,m_SIS], 'k:', zorder=1)
    ax1.text(150, m_SIS-0.33, r'$M_{\rm SIS}$', fontsize=12)
    print('      M_SIS = {0:5.2f} Msun'.format(m_SIS))

    print('      M_max = {0:5.2f} Msun'.format(m_ext_max))
    print('M_max/M_SIS = {0:4.2f}'.format(m_ext_max/m_SIS))
    print('  P_0/P_ext = {0:5.2f}'.format(Pratio_max))
    ax1.plot([Pratio_max,Pratio_max], [0,10], 'k:')

    #x_labels = ['1','10','100']
    x_labels = ['1','3','10','30','100','300']
    x_loc = np.array([float(x) for x in x_labels])
    ax1.set_xticks(x_loc)
    ax1.set_xticklabels(x_labels)

    fig.tight_layout() 
    plt.savefig('bonnor_ebert_mass.pdf')


def plot_b68():
    fig = plt.figure(figsize=(6,4))
    ax = fig.add_subplot(111)

    # observed profile
    # data from Alves et al. Nature 2001
    # Figure 2 digitized using https://apps.automeris.io/wpd/
    r, Av = np.genfromtxt('Alves_Av.txt', unpack=True, delimiter=',')
    ax.plot(r, Av, 'ko', markersize=3, label='Observations')

    nsteps = 10000
    xmax = 10
    x = np.logspace(-2, np.log10(xmax), nsteps)
    y, yp = lane_emden_integrate(x)

    # set outer boundary
    # (note that I find a value a bit lower than in Alves et al...)
    xmax = 4.5
    y[x > xmax] = 10

    b = np.logspace(-2, np.log10(xmax)+0.5, 1000)
    Av = np.zeros(b.size)
    yinterp = interpolate.interp1d(x, y, kind='cubic',
        bounds_error=False, fill_value='extrapolate')
    for i in range(b.size):
        b1 = b[i]
        xpath = np.sqrt(x**2 + b1**2)
        Av[i] = integrate.simps(np.exp(-yinterp(xpath)), xpath)

    # manually scale axes to match Av
    # this has physical significance but that's the point of the paper
    # (this illustrative plot is only to show that an observed core does indeed look like a BE sphere)
    Ascale = 35/Av.max()
    Av *= Ascale
    b *= 26
    ax.plot(b, Av, 'k-', lw=2, alpha=0.5, label='Bonnor Ebert profile')

    ax.set_xlim(8,150)
    ax.set_ylim(0.3,45)
    ax.set_xscale('log')
    ax.set_yscale('log')
    ax.set_xlabel("Projected radius ('')", fontsize=14)
    ax.set_ylabel(r"${\rm A_V\ (mag)}$", fontsize=14)
    ax.legend(loc=3, bbox_to_anchor=(0.04, 0.05))
    ax.text(0.24, 0.24, 'B68 visual extinction', fontsize=12, ha='center', transform = ax.transAxes)

    x_labels = ['10','30','100']
    x_loc = np.array([float(x) for x in x_labels])
    ax.set_xticks(x_loc)
    ax.set_xticklabels(x_labels)
    y_labels = ['1','3','10','30']
    y_loc = np.array([float(y) for y in y_labels])
    ax.set_yticks(y_loc)
    ax.set_yticklabels(y_labels)

    fig.tight_layout() 
    plt.savefig('b68_profile.pdf')


# Figure 9.1
plot_profiles()

# Figure 9.2
plot_mass()

# Figure 9.3
plot_b68()



StarcoderdataJupyterScriptsDedupFiltered