Added example figure of the profiles.

import matplotlib.pyplot as plt

import numpy as np

import os

# plt.rcParams.update({'font.size': 12})

p = -1.5

q = -0.5

n_0 = 1e5

T_0 = 25

radius = np.linspace(0, 0.9, 100)

density = n_0 * (radius / 0.5) ** p

density[0] = 1e8

temperature = T_0 * (radius / 0.5) ** q

temperature[0] = 2000

plt.clf()

fig, ax = plt.subplots(layout='constrained')

ln1 = ax.plot(radius, density, label=r'n(H$_2$)')

ax.semilogy()

ax.set_ylabel(r'Number density [cm$^{-3}$]')

ax.set_xlabel(r'Radius [pc]')

ax2 = ax.twinx()

ln2 = ax2.plot(radius, temperature, color='red', label='T')

lns = ln1 + ln2

labs = [l.get_label() for l in lns]

ax.legend(lns, labs, loc=0)

ax2.semilogy()

ax2.set_ylabel('Temperature [K]')

plt.savefig(os.path.join('prs', 'output', 'density_temperature_profiles.png'))