I wrote a program intended to calculate the plasma beta of a hydrostatic solar atmosphere:
B_0 = 0.02 # base magnetic field [T]
p_0 = 0.015 # base pressure [J/m^-3]
h_D = 7.5E7 # dipole depth [m]
R_☉ = 6.96E8 # solar radius [m]
M_☉ = 1.9891E30 # solar mass [kg]
G = 6.673E-11 # gravitational constant [m^3 kg^-1 s^-2]
μ = 0.61 # mean molecular weight of solar corona
m_H = 1.6726E-27 # mass of H particle [kg]
μ_0 = (4*pi)*10^(-7) # permeability [H/m]
k = 1.3806E-23 # boltzmann constant [J/K]
T = 1E6 # temperature [K]
g = (G*M_☉)/(R_☉^2)
lambda_p = (k*T)/(μ*m_H*g) # coronal scale height for 1MK [m]
p(h) = p_0*exp(-h/lambda_p)
B(h) = B_0*(1 + (h/h_D))^(-3)
I am not very familiar with using iterative/numerical methods (let alone implement them with Julia) and was wondering if there was a way to determine the value of h when beta = 1?
You can use the classical Newton method.
I was wondering specifically if there was a Julia-esque way, or a way that uses a certain Julia package? Otherwise, I might as well learn Newton method.
The simplest answer is probably to formulate this as a root finding problem (ie find the root of f(h) = p(h)/(B(h)^2 / (2*μ_0) )-beta and use Roots.jl
1 Like
So, something like this? I am not familiar with find_zero, specifically the second argument.
f = h -> p(h)/(B(h)^2 / (2*μ_0) ) - beta(h)
Roots.find_zero(f,2*h_D)
The second argument is an initial guess. You can pass anything you want, but if you happen to have an idea of what the root should be, it can significantly improve performance.
Understood, thanks! Anyway, is this an appropriate approach for determining h(beta = 1)? And do I make the result print more decimals?
it should print all the decimals it can compute. You should be able to get more by doing your calculation with BigFloats (but it will be ~50x slower).
Being an ex-MATLAB et al., I really miss the simplicity of doing such things symbolically. Even if a symbolic solution was not found, the algorithm will automatically fallback to a numerical solution without exposing the user to advanced details of the solution algorithm used.
>> h = solve(2*mu0*p/B^2 - beta)
Warning: Unable to solve symbolically. Returning a numeric solution using vpasolve.
h =
-234392042.24591558379314753229238
If you want a symbolic solution, you can try Symbolics.jl. I’m not sure how it works though.