diff --git a/src/construct_parker.py b/src/construct_parker.py index 649a39a..81e9284 100644 --- a/src/construct_parker.py +++ b/src/construct_parker.py @@ -189,7 +189,7 @@ def save_plain_parker_profile(planet, Mdot, T, spectrum, h_fraction=0.9, def save_temp_parker_profile(planet, Mdot, T, spectrum, zdict, pdir, - mu_bar=None, mu_struc=None, no_tidal=False, altmax=20): + mu_bar, mu_struc=None, no_tidal=False, altmax=20): """ Uses the p-winds code (dos Santos et al. 2022) Runs p_winds and saves a 'pprof' txt file with the r, rho, v, mu structure. @@ -218,9 +218,8 @@ def save_temp_parker_profile(planet, Mdot, T, spectrum, zdict, pdir, where the isothermal parker wind density and velocity profiles are saved. Different folders may exist there for a given planet, to separate for example profiles with different assumptions such as stellar SED/semi-major axis/composition. - mu_bar : float, optional + mu_bar : float Weighted mean of the mean particle mass. Based on Eq. A.3 of Lampon et al. (2020). - If None, p-winds will calculate mu(r) and the associated mu_bar. By default None. mu_struc : numpy.ndarray, optional Mean particle mass profile, must be provided if mu_bar is None. Typically, this is a mu(r)-profile as given by Cloudy. By default None. @@ -235,10 +234,6 @@ def save_temp_parker_profile(planet, Mdot, T, spectrum, zdict, pdir, ------- save_name : str Full path + filename of the saved Parker wind profile file. - mu_bar : float - Weighted mean of the mean particle mass. Based on Eq. A.3 of Lampon et al. (2020). - If the input mu_bar was None, this will return the value as calculated by p-winds. - If the input mu_bar was not None, this will return that same value. launch_velocity : float Velocity at the planet radius in units of the sonic speed. If it is larger than 1, the wind is "launched" already supersonic, and hence the assumption of a transonic @@ -254,33 +249,8 @@ def save_temp_parker_profile(planet, Mdot, T, spectrum, zdict, pdir, m_dot = 10 ** Mdot # Total atmospheric escape rate in g / s r = np.logspace(0, np.log10(altmax), 1000) # Radial distance profile in unit of planetary radii - - if mu_bar is None: #if not given by a Cloudy run, let p-winds calculate it (used the first iteration) - #pretend that the metals don't exist and just calculate the h_fraction with only H and He abundances - abundances = tools.get_abundances(zdict) #solar abundances - h_fraction = abundances['H'] / (abundances['H'] + abundances['He']) #approximate it by this for now, later Cloudy will give mu - - # A few assumptions about the planet's atmosphere - he_fraction = 1 - h_fraction # He number fraction - he_h_fraction = he_fraction / h_fraction - mean_f_ion = 0.0 # Mean ionization fraction (will be self-consistently calculated later) - mu_0 = (1 + 4 * he_h_fraction) / (1 + he_h_fraction + mean_f_ion) - # mu_0 is the constant mean molecular weight (assumed for now, will be updated later) - - initial_f_ion = 0. - - f_r, mu_bar = pw_hydrogen.ion_fraction(r, R_pl, T, h_fraction, - m_dot, M_pl, mu_0, - spectrum_at_planet=spectrum, exact_phi=True, - initial_f_ion=initial_f_ion, relax_solution=True, - return_mu=True, atol=1e-8, rtol=1e-5, - convergence=0.0001, max_n_relax=30) #I personally think we can use more than 0.01 convergence - - mu_array = ((1-h_fraction)*4.0 + h_fraction)/(h_fraction*(1+f_r)+(1-h_fraction)) #this assumes no Helium ionization - - else: #used later iterations - assert np.abs(mu_struc[0,0] - 1.) < 0.03 and np.abs(mu_struc[-1,0] - altmax) < 0.0001, "Looks like Cloudy didn't simulate to 1Rp: "+str(mu_struc[0,0]) #ensure safe extrapolation - mu_array = interp1d(mu_struc[:,0], mu_struc[:,1], fill_value='extrapolate')(r) + assert np.abs(mu_struc[0,0] - 1.) < 0.03 and np.abs(mu_struc[-1,0] - altmax) < 0.0001, "Looks like Cloudy didn't simulate to 1Rp: "+str(mu_struc[0,0]) #ensure safe extrapolation + mu_array = interp1d(mu_struc[:,0], mu_struc[:,1], fill_value='extrapolate')(r) vs = pw_parker.sound_speed(T, mu_bar) # Speed of sound (km/s, assumed to be constant) if no_tidal: @@ -305,7 +275,7 @@ def save_temp_parker_profile(planet, Mdot, T, spectrum, zdict, pdir, launch_velocity = v_array[0] #velocity at Rp in units of sonic speed - return save_name, mu_bar, launch_velocity + return save_name, launch_velocity def run_parker_with_cloudy(filename, T, planet, zdict): @@ -457,9 +427,10 @@ def save_cloudy_parker_profile(planet, Mdot, T, spectrum, zdict, pdir, tools.verbose_print("Making initial parker profile while assuming a completely neutral mu_bar...", verbose=verbose) neutral_mu_bar = calc_neutral_mu(zdict) neutral_mu_struc = np.array([[1., neutral_mu_bar], [altmax, neutral_mu_bar]]) #set up an array with constant mu(r) at the neutral value - filename, previous_mu_bar, launch_velocity = save_temp_parker_profile(planet, Mdot, T, spectrum, zdict, pdir, - mu_bar=neutral_mu_bar, mu_struc=neutral_mu_struc, no_tidal=no_tidal, altmax=altmax) - tools.verbose_print(f"Saved temp parker profile with neutral mu_bar: {previous_mu_bar}" , verbose=verbose) + filename, launch_velocity = save_temp_parker_profile(planet, Mdot, T, spectrum, zdict, pdir, + neutral_mu_bar, mu_struc=neutral_mu_struc, no_tidal=no_tidal, altmax=altmax) + tools.verbose_print(f"Saved temp parker profile with neutral mu_bar: {neutral_mu_bar}" , verbose=verbose) + previous_mu_bar = neutral_mu_bar for itno in range(maxit): tools.verbose_print(f"Iteration number: {itno+1}", verbose=verbose) @@ -474,8 +445,8 @@ def save_cloudy_parker_profile(planet, Mdot, T, spectrum, zdict, pdir, mu_bar = calc_mu_bar(sim) tools.verbose_print(f"Making new parker profile with p-winds based on Cloudy's reported mu_bar: {mu_bar}", verbose=verbose) mu_struc = np.column_stack((sim.ovr.alt.values[::-1]/planet.R, sim.ovr.mu[::-1].values)) #pass Cloudy's mu structure to save in the pprof - filename, mu_bar, launch_velocity = save_temp_parker_profile(planet, Mdot, T, spectrum, zdict, pdir, - mu_bar=mu_bar, mu_struc=mu_struc, no_tidal=no_tidal, altmax=altmax) + filename, launch_velocity = save_temp_parker_profile(planet, Mdot, T, spectrum, zdict, pdir, + mu_bar, mu_struc=mu_struc, no_tidal=no_tidal, altmax=altmax) tools.verbose_print("Saved temp parker profile.", verbose=verbose) if np.abs(mu_bar - previous_mu_bar)/previous_mu_bar < convergence: