I think I may have spotted a small mistake in the average_molecular_weight() function in the parker.py module. The function takes a "v_profile" argument, and multiplies that with 1000 to presumably get the velocity structure in units of m / s. However, when the function is called from within ion_fraction() in the hydrogen.py module (both the first time and in the relax_solution loop), the velocity is passed in units of the sound speed. In the average_molecular_weight() function, the sound speed is thus effectively always assumed to be 1 km/s. This difference in velocity leads to a slightly different mu_bar and hence a slightly different final parker wind structure. Similarly, the average_molecular_weight() function also takes an "r_profile" argument, which is multiplied by a Jupiter radius, but the function is called with an r-grid in units of the planet radius, so the planet is effectively always assumed to be the size of Jupiter, which has a slight impact on mu_bar as well.
I must say that I have not carefully looked through all of the p_winds source code, so excuse me if I am misinterpreting your code and if this is actually not a bug! Also, maybe the average_molecular_weight() function is also called from other places than just the ion_fraction() function, and if that's the case, I'm not sure whether the units are correct there.
Hi Leonardo,
Thanks for your great work on p_winds!
I think I may have spotted a small mistake in the average_molecular_weight() function in the parker.py module. The function takes a "v_profile" argument, and multiplies that with 1000 to presumably get the velocity structure in units of m / s. However, when the function is called from within ion_fraction() in the hydrogen.py module (both the first time and in the relax_solution loop), the velocity is passed in units of the sound speed. In the average_molecular_weight() function, the sound speed is thus effectively always assumed to be 1 km/s. This difference in velocity leads to a slightly different mu_bar and hence a slightly different final parker wind structure. Similarly, the average_molecular_weight() function also takes an "r_profile" argument, which is multiplied by a Jupiter radius, but the function is called with an r-grid in units of the planet radius, so the planet is effectively always assumed to be the size of Jupiter, which has a slight impact on mu_bar as well.
I must say that I have not carefully looked through all of the p_winds source code, so excuse me if I am misinterpreting your code and if this is actually not a bug! Also, maybe the average_molecular_weight() function is also called from other places than just the ion_fraction() function, and if that's the case, I'm not sure whether the units are correct there.
Thanks again!