source terms module for gravitational acceleration due to two pointmasses More...
Data Types | |
| type | gravity_binary |
Functions/Subroutines | |
| subroutine | initgravity_binary (this, Mesh, Physics, config, IO) |
| Constructor of gravity binary module. More... | |
| subroutine | infogravity (this) |
| write binary parameters to screen More... | |
| subroutine | setoutput (this, Mesh, Physics, config, IO) |
| subroutine | updategravity_single (this, Mesh, Physics, Fluxes, pvar, time, dt) |
| calculates the resultant gravitational accerleration of a binary system More... | |
| subroutine | updatepositions (this, time) |
| compute new positions for both components of the binary system More... | |
| pure subroutine | calcdiskheight_single (this, Mesh, Physics, pvar, bccsound, h_ext, height) |
| compute the pressure scale height of a disk for the binary potential More... | |
| real function | getmass_primary (this, time) |
| get mass of the primary component, eventually scaled by switch on factor More... | |
| real function | getmass_secondary (this, time) |
| get mass of the primary component, eventually scaled by switch on factor More... | |
| subroutine | finalize (this) |
| Closes the binary source term. More... | |
| pure subroutine | func (x, fx, plist) |
| find root of the Kepler equation More... | |
Detailed Description
source terms module for gravitational acceleration due to two pointmasses
Function/Subroutine Documentation
◆ calcdiskheight_single()
| pure subroutine gravity_binary_mod::calcdiskheight_single | ( | class(gravity_binary), intent(inout) | this, |
| class(mesh_base), intent(in) | Mesh, | ||
| class(physics_base), intent(in) | Physics, | ||
| class(marray_compound), intent(inout) | pvar, | ||
| type(marray_base), intent(inout) | bccsound, | ||
| type(marray_base), intent(inout) | h_ext, | ||
| type(marray_base), intent(inout) | height | ||
| ) |
compute the pressure scale height of a disk for the binary potential
Since the potential of a binary system is given by the superposition of the potential of two individual stars we get
\[ h = c_s \left(\sum_{i=1}^2\frac{\partial^2\Phi_i}{\partial z^2}\right)^{-1/2} = c_s \left(\Omega_1^2 + \Omega_2^2\right)^{-1/2} \]
where \( \Phi_i = -GM_i / r_i \) is the pointmass potential in the equatorial plane with respect to the \( i \) th component of the binary system and \( \Omega_i \) is the associated local Keplerian angular velocity.
Definition at line 511 of file gravity_binary.f90.
◆ finalize()
| subroutine gravity_binary_mod::finalize | ( | class(gravity_binary), intent(inout) | this | ) |
Closes the binary source term.
Definition at line 560 of file gravity_binary.f90.
◆ func()
|
private |
find root of the Kepler equation
Definition at line 576 of file gravity_binary.f90.
◆ getmass_primary()
|
private |
get mass of the primary component, eventually scaled by switch on factor
Definition at line 528 of file gravity_binary.f90.
◆ getmass_secondary()
|
private |
get mass of the primary component, eventually scaled by switch on factor
Definition at line 544 of file gravity_binary.f90.
◆ infogravity()
| subroutine gravity_binary_mod::infogravity | ( | class(gravity_binary), intent(in) | this | ) |
write binary parameters to screen
Definition at line 232 of file gravity_binary.f90.
◆ initgravity_binary()
| subroutine gravity_binary_mod::initgravity_binary | ( | class(gravity_binary), intent(inout) | this, |
| class(mesh_base), intent(in) | Mesh, | ||
| class(physics_base), intent(in) | Physics, | ||
| type(dict_typ), pointer | config, | ||
| type(dict_typ), pointer | IO | ||
| ) |
Constructor of gravity binary module.
This subroutine reads the necessary config data for the binary. It initializes the gravity type and various mesh data arrays. Some of those are marked for output.
- Parameters
-
[in,out] this [in,out] this all gravity data [in] mesh [in] Mesh mesh type [in] physics [in] Physics physics type config [in,out] config sub-dictionary with binary configuration data io [in,out] IO output sub-dictionary
Definition at line 114 of file gravity_binary.f90.
◆ setoutput()
|
private |
Definition at line 254 of file gravity_binary.f90.
◆ updategravity_single()
| subroutine gravity_binary_mod::updategravity_single | ( | class(gravity_binary), intent(inout) | this, |
| class(mesh_base), intent(in) | Mesh, | ||
| class(physics_base), intent(in) | Physics, | ||
| class(fluxes_base), intent(in) | Fluxes, | ||
| class(marray_compound), intent(inout) | pvar, | ||
| real, intent(in) | time, | ||
| real, intent(in) | dt | ||
| ) |
calculates the resultant gravitational accerleration of a binary system
The acceleration is computed with respect to the curvilinear local frame of reference. Each component of the binary system contributes according to
\[ \vec{g} = -\frac{d\Phi}{d r} = -r \Omega^2 \hat{e}_r = -\Omega^2 \vec{r} \]
to the total acceleration. Thereby \( \Omega^2 = GM/r^3 \) is the square of the local Keplerian velocity.
Definition at line 284 of file gravity_binary.f90.
◆ updatepositions()
|
private |
compute new positions for both components of the binary system
This subroutine solves the Kepler equation to get the positions at a given time. A detailed explanation and a derivation of the relevant equations can be found in [taff1985] .
First we map the dimensionless orbital time into the interval \([0,2\pi]\) to calculate the mean anomaly \( \tau \). Then we have to solve the Kepler equation:
\[ \tau= E-\epsilon\sin{E} \]
to calculate the eccentric anomaly \( E \). \( \epsilon \) is the eccentricity of the orbit which is given as an input parameter of the binary. We can solve the equation using the function roots::getroot_regulafalsi . To do this, we have to know an Intervall that includes the correct solution for \( E \). This should be as small as possible to speed up the calculation. Looking at the Kepler equation we can see, that if \( \tau\) is in the interval \([0,\pi]\) \(E\) also has a value in \([0,\pi]\). The same is true for the other possible intervall of values for \( \tau \) and \( E \) \(]\pi,2\pi[\). We can further constrain the Intervall for the correct solution for \( E \) by looking at these separate cases, keeping in mind that the eccentricity \( \epsilon \) can only have values of:
\[ 0 \leq \epsilon \leq 1 \]
Case 1 : \( E,\tau \in [0,\pi] \)
\begin{eqnarray*} E &=& \tau + \epsilon\sin{E} \leq \tau + \epsilon \\ E &=& \tau + \epsilon\sin{E} \geq \tau \\ \end{eqnarray*}
Case 2 : \( E,\tau \in ]\pi,2\pi[ \)
\begin{eqnarray*} E &=& \tau + \epsilon\sin{E} < \tau \\ E &=& \tau + \epsilon\sin{E} > \tau - \epsilon \\ \end{eqnarray*}
If we know the mean anomaly \( E \) we can calculate the distance between the stars
\[ r = a (1-\epsilon \cos{E}). \]
\( a \) is the semimajor axis of the binary. We can also calculate the angle \( \phi \) between the semimajoraxis of the ellipse and the line connecting the primary component to the center of mass
\[ \phi = 2\arctan\left(\sqrt{\frac{1+\epsilon}{1-\epsilon}} \tan{\frac{E}{2}} \right) \]
This equation can be written as
\[ \phi = 2\arctan\left(\pm\sqrt{\frac{1+\epsilon}{1-\epsilon} \frac{1-\cos{E}}{1+\cos{E}}} \right) \]
to avoid one evaluation of \( \tan \). To choose the coorect sign in this expression we have to look at the sign of \( \tan{(E/2)} \) which is positive for \( E\in[0,\pi] \) and negative for \( E\in[\pi,2\pi]\). When we know the angle \( \phi \) and the distance between the stars \( r \) we can calculate the positions of the stars \( \vec{r}_1,\vec{r}_2 \) using the equations:
\begin{eqnarray*} \vec{r}_1 &=& \frac{m_2}{m_1+m_2} \vec{r}\\ \vec{r}_2 &=& \vec{r}-\vec{r}_1 \\ \end{eqnarray*}
with the masses of the binary components \( m_1, m_2 \) and \( r_1=|\vec{r}_1|\).
Definition at line 437 of file gravity_binary.f90.
