import numpy as np
from astropy.constants import M_sun
from astropy import units as u
from scipy.integrate import solve_ivp
# class for the Weaver et al. (1977) stellar wind solution
# initiated with the
# terminal velocity v_inf
# mass loss rate M_dot
# background density rho_0
# which can be used to calculate the density and velocity at time t.
[docs]
class Weaver:
def __init__(self, v_inf, M_dot, rho_0, p_0, num_xi = 100, gamma = 5/3):
# input wind parameters
## terminal velocity of the wind
self.v_inf = v_inf
## mass loss rate of the wind
self.M_dot = M_dot
## background density of the ISM
self.rho_0 = rho_0
## background pressure of the ISM
self.p_0 = p_0
# derived wind parameters
## mechanical luminosity of the wind
self.L_w = 0.5 * M_dot * v_inf ** 2
# constants
self.xi_crit = 0.86
self.alpha = 0.88
self.gamma = gamma
# number of points to calculate the shell profiles
self.num_xi = num_xi
# calculate the shell profiles
self.calculate_shell_profiles()
[docs]
def calculate_shell_profiles(self):
# integrate equations 6 to 7 in weaver II, with
# 3 * (U - xi) * U' - 2U + 3 P' / G = 0
# (U - xi) * G' / G + U' + 2U / xi = 0
# 3 * (U - xi) * P' - 3 * gamma * P * (U - xi) * G' / G - 4P = 0
# where ' denotes \partial_\xi
# from xi = 1 to xi = xi_crit with G(1) = 4, U(1) = 3/4, P(1) = 3/4
# r = R_2 * xi; rho = rho_0 * G(xi)
# v = V_2 * U(xi); p = rho_0 * V_2^2 * P(xi)
gamma = self.gamma
# the equations are given by
def shell_equation(xi, y):
G, U, P = y
G_prime = 2*(-2*xi**2*G*U + 5*xi*G*U**2 + 2*xi*P - 3*G*U**3)*G/(3*xi*(gamma*xi*P - gamma*P*U - xi**3*G + 3*xi**2*G*U - 3*xi*G*U**2 + G*U**3))
U_prime = 2*(-3*gamma*P*U + xi**2*G*U - xi*G*U**2 + 2*xi*P)/(3*xi*(gamma*P - xi**2*G + 2*xi*G*U - G*U**2))
P_prime = 2*(-2*gamma*xi*U + 3*gamma*U**2 + 2*xi**2 - 2*xi*U)*G*P/(3*xi*(gamma*P - xi**2*G + 2*xi*G*U - G*U**2))
return [G_prime, U_prime, P_prime]
# initial conditions
G1 = 4
U1 = 3/4
P1 = 3/4
sol = solve_ivp(shell_equation, [1, self.xi_crit], [G1, U1, P1], t_eval=np.linspace(1, self.xi_crit, self.num_xi))
self.xi = np.flip(sol.t)
self.U = np.flip(sol.y[1])
self.G = np.flip(sol.y[0])
self.P = np.flip(sol.y[2])
# returns the inner shock radius R_1
[docs]
def get_inner_shock_radius(self, t):
return 0.9 * self.alpha ** 1.5 * (1 / self.rho_0 * self.M_dot) ** (3/10) * self.v_inf ** (1/10) * t ** (2/5)
# returns the outer shock radius R_2
[docs]
def get_outer_shock_radius(self, t):
return self.alpha * (self.L_w * t ** 3 / self.rho_0) ** (1/5)
# returns the critical radius R_c
[docs]
def get_critical_radius(self, t):
return self.xi_crit * self.get_outer_shock_radius(t)
[docs]
def get_radial_range_wind_interior(self, delta_R, t):
R_1 = self.get_inner_shock_radius(t)
R_c = self.get_critical_radius(t)
return np.arange(R_1.to(u.parsec).value, R_c.to(u.parsec).value, delta_R.to(u.parsec).value) * u.parsec
[docs]
def get_radial_range_free_wind(self, delta_R, t):
R_1 = self.get_inner_shock_radius(t)
return np.arange(delta_R.to(u.parsec).value, R_1.to(u.parsec).value, delta_R.to(u.parsec).value) * u.parsec
[docs]
def get_radial_range_undisturbed_ism(self, delta_R, R_max, t):
R_2 = self.get_outer_shock_radius(t)
return np.arange(R_2.to(u.parsec).value, R_max.to(u.parsec).value, delta_R.to(u.parsec).value) * u.parsec
# get inner pressure
[docs]
def get_pressure_profile(self, delta_R, R_max, t):
# free wind profile (?????)
# Rs_free_wind = self.get_radial_range_free_wind(delta_R, t)
# pressures_free_wind = 1/20 * self.M_dot * self.v_inf / (4 * np.pi * Rs_free_wind ** 2)
# wind interior
Rs_wind_interior = np.array([self.get_inner_shock_radius(t).to(u.parsec).value, self.get_critical_radius(t).to(u.parsec).value]) * u.parsec
pressure_wind_interior = 5 / (22 * np.pi * (self.xi_crit * self.alpha) ** 3) * (self.L_w ** 2 * self.rho_0 ** 3) ** (1/5) * t ** (-4/5) * np.array([1, 1])
# shell
# as of the RH-jump conditions at a contact discontinuity,
# the pressure is continuous across the contact discontinuity
Rs_shell = self.xi * self.get_outer_shock_radius(t)
V2 = 15 / 25 * self.get_critical_radius(t) / t
pressure_shell = self.rho_0 * V2 ** 2 * self.P
pressure_shell = pressure_shell / pressure_shell[0] * pressure_wind_interior[1] # why is this necessary? sth wrong with the solution??
# undisturbed ISM
Rs_undisturbed_ism = self.get_radial_range_undisturbed_ism(delta_R, R_max, t)
pressure_undisturbed_ism = self.p_0 * np.ones(len(Rs_undisturbed_ism))
return np.concatenate((Rs_wind_interior, Rs_shell, Rs_undisturbed_ism)), np.concatenate((pressure_wind_interior, pressure_shell, pressure_undisturbed_ism))
# get velocity profile
[docs]
def get_velocity_profile(self, delta_R, R_max, t):
# free wind profile
Rs_free_wind = self.get_radial_range_free_wind(delta_R, t)
velocities_free_wind = self.v_inf * np.ones(len(Rs_free_wind))
# wind interior
Rs_wind_interior = self.get_radial_range_wind_interior(delta_R, t)
R_c = self.get_critical_radius(t)
velocities_wind_interior = 11 / 25 * R_c ** 3 / (Rs_wind_interior ** 2 * t) + 4 / 25 * Rs_wind_interior / t
# shell
Rs_shell = self.xi * self.get_outer_shock_radius(t)
# by the RH-jump conditions at a contact discontinuity,
# the velocity at the beginning of the shell is the
# velocity at the end of the wind interior, so
V2 = 15 / 25 * R_c / t
velocities_shell = V2 * self.U / self.U[0]
# undisturbed ISM
# here the velocity is 0
Rs_undisturbed_ism = self.get_radial_range_undisturbed_ism(delta_R, R_max, t)
velocities_unisturbed_ism = np.zeros(len(Rs_undisturbed_ism))
return np.concatenate((Rs_free_wind, Rs_wind_interior, Rs_shell, Rs_undisturbed_ism)), np.concatenate((velocities_free_wind, velocities_wind_interior, velocities_shell, velocities_unisturbed_ism))
# get density profile
[docs]
def get_density_profile(self, delta_R, R_max, t):
# free wind profile
Rs_free_wind = self.get_radial_range_free_wind(delta_R, t)
densities_free_wind = self.M_dot / (4 * np.pi * Rs_free_wind ** 2 * self.v_inf)
# profile in the wind interior
Rs_wind_interior = self.get_radial_range_wind_interior(delta_R, t)
R_c = self.get_critical_radius(t)
densities_wind_interior = 0.628 * (self.M_dot ** 2 * self.rho_0 ** 3 * self.v_inf ** (-6)) ** (1 / 5) * t ** (-4 / 5) * (1 - (Rs_wind_interior / R_c) ** 3) ** (-8 / 33)
# profile in the shell
Rs_shell = self.xi * self.get_outer_shock_radius(t)
densities_shell = self.rho_0 * self.G
# profile in the undisturbed ISM
Rs_undisturbed_ism = self.get_radial_range_undisturbed_ism(delta_R, R_max, t)
densities_unisturbed_ism = self.rho_0 * np.ones(len(Rs_undisturbed_ism))
return np.concatenate((Rs_free_wind, Rs_wind_interior, Rs_shell, Rs_undisturbed_ism)), np.concatenate((densities_free_wind, densities_wind_interior, densities_shell, densities_unisturbed_ism))
