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, public | kepler (time, period, excent, semax, r, phi) |
| solve Kepler-Equation and return new distance r and position angle phi 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 440 of file gravity_binary.f90.
◆ finalize()
| subroutine gravity_binary_mod::finalize | ( | type(gravity_binary), intent(inout) | this | ) |
Closes the binary source term.
Definition at line 489 of file gravity_binary.f90.
◆ func()
| pure subroutine gravity_binary_mod::func | ( | real, intent(in) | x, |
| real, intent(out) | fx, | ||
| real, dimension(:), intent(in), optional | plist | ||
| ) |
find root of the Kepler equation
Definition at line 606 of file gravity_binary.f90.
◆ getmass_primary()
|
private |
get mass of the primary component, eventually scaled by switch on factor
Definition at line 457 of file gravity_binary.f90.
◆ getmass_secondary()
|
private |
get mass of the primary component, eventually scaled by switch on factor
Definition at line 473 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 236 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 115 of file gravity_binary.f90.
◆ kepler()
| pure subroutine, public gravity_binary_mod::kepler | ( | real, intent(in) | time, |
| real, intent(in) | period, | ||
| real, intent(in) | excent, | ||
| real, intent(in) | semax, | ||
| real, intent(out) | r, | ||
| real, intent(out) | phi | ||
| ) |
solve Kepler-Equation and return new distance r and position angle phi
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 [52] .
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]\).
- Parameters
-
[in] time orbital time [in] period orbital period [in] excent excentricity [in] semax semi-major axis [out] r mean anomaly [out] phi excentricity
Definition at line 556 of file gravity_binary.f90.
◆ setoutput()
|
private |
Definition at line 263 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 293 of file gravity_binary.f90.
◆ updatepositions()
|
private |
compute new positions for both components of the binary system
When we know the angle \( \phi \) and the distance between the stars \( r \) we can calculate the positions of the two 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 400 of file gravity_binary.f90.
