☰Double Star Calculator - Documentation

Content

Abstract

Core Astrometric and Kinematic Functions

Angular Separation
Position Angle
Spatial Separation
Proper Motion and Velocity Analysis
Proper Motion Check
Combined Physical Binding Evidence

Auxiliary Astrometric Functions

Distance from Parallax
Projected Separation
Absolute Magnitude
Proper Motion Angular Difference
Total Proper Motion Difference
Velocity Differences
Spectral Type Estimation
Mass–Luminosity Relation

Orbit Determination

Optimization Summary

Abstract

This document provides a detailed methodological description of the computational procedures implemented in the Double Star Calculator, a tool designed for the analysis of visual and astrometric double stars based primarily on data from the Gaia mission.

The calculator derives a comprehensive set of astrometric and kinematic quantities for individual double star systems, including angular and spatial separations, position angles, proper motion vectors, distance estimates, tangential and spatial velocities, and derived physical parameters such as absolute magnitudes, luminosities, and mass estimates.

All computations are based on established geometrical relations, classical error propagation using partial derivatives, and standard astronomical conventions. The mathematical formulations underlying each calculation are explicitly documented to ensure transparency, reproducibility, and independent verification of the results.

This documentation is intended to serve as a technical reference for users of the Double Star Calculator, as well as a methodological supplement to scientific publications that make use of its results.

Core Astrometric and Kinematic Functions

This section describes the core astrometric and kinematic computations performed for individual double star systems. These calculations form the basis of the physical interpretation and are directly used in scientific analyses.

Angular Separation

This chapter describes the mathematical formulation used to compute the angular separation (ρ) between two stellar components based on Gaia astrometric data.

Input Parameters

Right ascension: α_A, α_B (degrees)
Declination: δ_A, δ_B (degrees)
Astrometric uncertainties: σ_{α_A}, σ_{δ_A}, σ_{α_B}, σ_{δ_B} (milliarcseconds)

Output Parameters

The function returns the angular separation ρ (radians) and its associated uncertainty σρ (radians).

ρ: angular separation between the two stars (radians)
σρ: uncertainty of the angular separation (radians)

Angular Separation

The angular separation 𝜌 between the two components was computed using spherical trigonometry based on their Gaia DR3 equatorial coordinates. The separation is given by:

The angular separation ρ is calculated using the spherical law of cosines:

ρ = \arccos (\sin δ_{A} \cdot \sin δ_{B} + \cos δ_{A} \cdot \cos δ_{B} \cdot \cos (α_{A} - α_{B}))

Error Propagation

The uncertainty of the angular separation σ_ρ is computed using first-order Gaussian error propagation:
The uncertainty of the separation was derived using first-order Gaussian error propagation, assuming uncorrelated uncertainties in right ascension and declination:

σ_{ρ} = \sqrt{\sum_{i}^{4} {(\frac{\partial ρ}{\partial x_{i}} \cdot σ_{x_{i}})}^{2}}

Where x_i ∈ {α_A, δ_A, α_B, δ_B}

All astrometric uncertainties were converted from milliarcseconds to radians prior to the calculation.

Assumptions

Astrometric uncertainties are treated as independent
Covariances between Gaia parameters are not considered
The calculation is performed in the ICRS reference frame

Position Angle

This function computes the position angle θ of a double star measured from North towards East (0°–360°), based on Gaia DR3 coordinates of the two components. The associated uncertainty is calculated using first-order Gaussian error propagation.

Input Parameters

Right ascension: α_A, α_B (degrees)
Declination: δ_A, δ_B (degrees)
Astrometric uncertainties: σ_αA, σ_δA, σ_αB, σ_δB (milliarcseconds)

Position Angle

The position angle θ is calculated using the arctangent of the differential coordinates, where atan2(y, x) denotes the quadrant-correct two-argument arctangent, ensuring the correct quadrant of the position angle.

θ = atan2 (\cos δ_{B} \cdot \sin (α_{B} - α_{A}), \cos δ_{A} \cdot \sin δ_{B} - \sin δ_{A} \cdot \cos δ_{B} \cdot \cos (α_{B} - α_{A}))

The resulting angle is converted to degrees and shifted to the range 0°–360°.

Error Propagation

The uncertainty σ_θ is calculated using first-order Gaussian error propagation, based on the partial derivatives of atan2:

σ_{θ} = \sqrt{\sum_{i}^{4} {(\frac{\partial θ}{\partial x_{i}} \cdot σ_{x_{i}})}^{2}}

Where x_i ∈ {α_A, δ_A, α_B, δ_B}

Output Parameters

θ: position angle from North → East → South → West (degrees, 0°–360°)
σθ: uncertainty of the position angle (degrees)

Assumptions

Astrometric uncertainties are treated as independent.
Covariances between Gaia parameters are not considered.
The calculation is performed in the ICRS reference frame.

Spatial Separation

This function calculates the spatial separation between two stars in parsecs, using their distances and angular separation. It also computes the uncertainty of the spatial separation via error propagation.

Spatial separation

s = \sqrt{{d_{A}}^{2} + {d_{B}}^{2} - 2 \cdot d_{A} \cdot d_{B} \cdot \cos (θ)}

where $d_{A}$ and $d_{B}$ are the distances to the two stars and $θ$ is their angular separation in radians.

Error propagation

The uncertainty of the spatial separation is calculated using partial derivatives:

σ_{s} = \sqrt{{(\frac{∂s}{{∂d}_{A}} \cdot σ d_{A})}^{2} + {(\frac{∂s}{{∂d}_{B}} \cdot σ d_{B})}^{2} + {(\frac{∂s}{∂θ} \cdot σ θ_{rad})}^{2}}

Output Parameters

s – spatial separation in parsecs
σ_s – uncertainty of spatial separation

Proper Motion and Velocity Analysis

This function performs a kinematic analysis of a single star based on Gaia astrometric data. It derives the direction and magnitude of the proper motion vector, the tangential velocity, and, if available, the total space velocity including the radial component. All quantities are accompanied by uncertainty estimates.

Proper Motion Vector Angle

The position angle of the proper motion vector is calculated relative to the north direction, measured counterclockwise from north through east (0°–360°). The quadrant-correct two-argument arctangent is used.

θ = atan2 (μ α, μ δ)

Here, $μ α$ and $μ δ$ denote the proper motion components in right ascension and declination, respectively. The function atan2(y, x) ensures correct quadrant assignment.

Uncertainty of the Proper Motion Angle

The uncertainty of the proper motion angle is estimated using a first-order approximation based on the relative uncertainties of the proper motion components. This approach provides numerically stable results for typical Gaia data but does not represent a full analytical error propagation of the atan2 function.

σ_{θ} \approx θ \cdot \sqrt{{(\frac{σ μ_{α}}{μ_{δ}})}^{2} + {(\frac{σ μ_{δ}}{μ_{α}})}^{2}}

Total Proper Motion

The total proper motion is calculated as the magnitude of the proper motion vector:

μ = \sqrt{{μ_{α}}^{2} + {μ_{δ}}^{2}}

Uncertainty of Total Proper Motion

σ_{μ} = \sqrt{{(\frac{μ_{α}}{μ} \cdot σ_{μ_{α}})}^{2} + {(\frac{μ_{δ}}{μ} \cdot σ_{μ_{δ}})}^{2}}

Tangential Velocity

The tangential velocity is derived from the total proper motion and the parallax:

V_{t} = 4.74057 \cdot \frac{μ}{π}

The numerical factor 4.74057 converts proper motion in milliarcseconds per year and distance in parsecs into velocity in km s⁻¹.

Uncertainty of Tangential Velocity

σ_{V_{t}} = \sqrt{{(\frac{4.74057}{π} \cdot σ μ)}^{2} + {(\frac{4.74057}{π^{2}} \cdot μ \cdot σ π)}^{2}}

Space Velocity

If a radial velocity measurement is available, the total space velocity is computed as:

V = \sqrt{{V_{t}}^{2} + {V_{r}}^{2}}

Uncertainty of Space Velocity

σ_{V} = \sqrt{{(\frac{V_{t}}{V} \cdot σ V_{t})}^{2} + {(\frac{V_{r}}{V} \cdot σ V_{r})}^{2}}

Output Parameters

Proper motion angle (degrees)
Uncertainty of proper motion angle (degrees)
Total proper motion (mas/yr)
Uncertainty of total proper motion
Tangential velocity (km/s)
Uncertainty of tangential velocity
Total space velocity (km/s)
Uncertainty of total space velocity

Proper Motion Check

The probability of a pair being a CPM (Common Proper Motion) pair is estimated based on the following measured quantities and their respective weights:

Proper Motion Angular Difference (weight 0.3)
Total Proper Motion Difference (weight 0.3)
Tangential Velocity Difference (weight 0.1)
Radial Velocity Difference (weight 0.1)
Space Velocity Difference (weight 0.2)

The Double Star Calculator computes a weighted sum of the normalized absolute differences listed above. This sum is then mapped through an exponential decay function to obtain the final CPM probability. The result is scaled to the range 0–100 % and rounded.

P_{CPM} = 100 \cdot \exp (- k \cdot S^{1.5})

with k:

k = \frac{\ln (2)}{50}

Here, 𝑆 denotes the weighted sum of the normalized absolute differences of the proper-motion and velocity parameters. The exponential decay ensures a rapid decrease of the CPM probability for increasing kinematic discrepancies, while small differences result in high CPM probabilities.

The resulting value represents an estimated likelihood that the two stars form a CPM pair. It should be emphasized that this approach is not intended as a rigorous statistical model, but rather as a practical indicator to assess whether a given pair exhibits common proper motion characteristics.

CPM Evidence factor

The calculated CPM Evidence factor (range 0.0 … 1.0) is derived from the estimated CPM probability using a power-law transformation. This transformation is intended to reduce dominance and to reflect the fact that common proper motion is a necessary, but not sufficient, condition for a physical (gravitationally bound) system.

The CPM Evidence factor is computed according to:

CPM evidence = {(\frac{CPM probability}{100.0})}^{2.0}

The exponent of 2.0 was chosen conservatively in order to prevent the CPM Evidence factor from becoming dominant in the overall physical-binding assessment. With this formulation, very similar proper motion values result in a high CPM Evidence factor, while intermediate cases are attenuated and clearly inconsistent proper motions yield very low evidence values.

Further testing may reveal that the current exponent is overly conservative. In such cases, a reduced exponent (e.g. 1.5) could be considered to moderately increase sensitivity.

Interpretation of the CPM Evidence factor:

    cpm_evidence  Meaning
    ------------------------------------------
    ≥ 0.8         strong CPM support
    0.4 – 0.8     weak to moderate CPM support
    < 0.4         limited CPM support
    < 0.2         practically no CPM support

Combined Physical Binding Evidence

To assess the overall likelihood that a double star system is physically bound, the individual evidential indicators are combined into a single Combined Physical Binding Evidence factor.

Currently, two independent evidence factors are considered:

Separation-Based Evidence factor (sep_evidence), reflecting the geometric plausibility of a bound system based on spatial separation and its uncertainty.
CPM-Based Evidence factor (cpm_evidence), reflecting the kinematic similarity of the two components based on proper motion and velocity differences.

Both factors are normalized to the range 0.0 … 1.0 and represent evidential support, not probabilities. The combined evidence is as follows computed:

\begin{array}{l} w(sep) = 4 \cdot k \cdot {| sep - 0.5 |}^{1.7} + (1 - k) \\ comb = w(sep) \cdot sep + (1 - w(sep)) \cdot cpm \end{array}

Factor Combination Logic: Adaptive Dominance

The resulting Combined Physical Binding Evidence factor (comb) is derived from a weighted combination of Separation Evidence and CPM Evidence factors. Unlike a simple linear average, this model employs a non-linear weighting function to account for the increasing dominance of the Separation factor at its extreme values.

Key Characteristics:

Extreme SEP Dominance: When SEP is near 0.0 or 1.0, it becomes the primary driver of the result. In these regions, the influence of CPM factor is significantly suppressed.
Balanced Mid-Range: When SEP is in the central range (around 0.5), both factors contribute with approximately equal relevance to the final COMB value.
Non-Linear Transition: By using a power-based weight distribution, the transition from balanced weighting to SEP dominance is accelerated, creating a "sharp" response curve that better aligns with empirical measurements.

Combined Evidence Factor — Illustration 3: Combined Physical Binding Evidence factor

The following table provides a preliminary interpretation guideline. The numerical boundaries are heuristic and may be refined based on empirical validation using systems with well-determined orbits and known optical pairs.

    Combined Evidence   Interpretation
    ---------------------------------------------------------------
    ≥ 0.85              Very strong evidence for a physical pair
    0.65 – 0.85         Strong evidence; likely physical association
    0.45 – 0.65         Moderate evidence; candidate physical pair
    0.25 – 0.45         Weak evidence; ambiguous, requires scrutiny
    < 0.25              Very limited evidence; likely optical pair

Notes:

The combined evidence is intended as a ranking and filtering metric, not as a strict classifier.
Systems with moderate combined evidence may still represent physical binaries, especially in cases of incomplete or noisy data.
Final assessment should consider additional indicators (e.g. orbital motion, astrophysical consistency, long-term monitoring).

Auxiliary Astrometric Functions

This chapter describes auxiliary functions used in the analysis of stellar kinematics and physical association. These functions support the main astrometric routines and provide derived quantities and comparison metrics.

Distance from Parallax

This function converts the stellar parallax into distance, expressed in parsecs and light-years, including the propagation of the parallax uncertainty.

The calculation follows the standard astronomical relation between parallax and distance and assumes Gaussian error propagation.

Distance in parsecs

d_{pc} = \frac{1000}{ϖ}

where $ϖ$ is the parallax in milliarcseconds (mas).

Uncertainty of the distance

σ d_{pc} = \frac{1000}{ϖ^{2}} \cdot σ ϖ

The uncertainty is derived using standard Gaussian error propagation.

Conversion to light-years

d_{ly} = d_{pc} \cdot C_{pc→ly}

σ d_{ly} = σ d_{pc} \cdot C_{pc→ly}

with the conversion constant $C_{pc→ly} = 3.26156$ .

Output Parameters

d_pc – distance in parsecs
σ_d,pc – distance uncertainty in parsecs
d_ly – distance in light-years
σ_d,ly – distance uncertainty in light-years

Projected Separation

Calculates a conservative lower bound of the projected physical separation between two stars by evaluating the projected separation at the distance of the nearer component.

s = d_{near} \cdot \tan (θ)

where d_near is the smaller of the two distances and θ is the angular separation in radians.

Returns: projected separation in parsec

Absolute Magnitude

This function calculates the absolute magnitude of a star from its apparent Gaia G-band mean magnitude and its distance in parsecs. The uncertainty of the absolute magnitude is derived using the exact propagation of uncertainty, correctly accounting for the logarithmic dependence on distance.

Absolute magnitude

M = m - 5 \cdot \log_{10} (d_{pc}) + 5

where $m$ is the apparent Gaia G-band mean magnitude, and $d_{pc}$ is the distance in parsecs.

Uncertainty of absolute magnitude

σ_{M} = \frac{5}{\ln 10} \cdot \frac{σ d_{pc}}{d_{pc}}

This expression uses exact error propagation, taking into account the logarithmic dependence on distance.

Output Parameters

M – absolute magnitude
σ_M – uncertainty of absolute magnitude (exact propagation)

Proper Motion Angular Difference

Computes the absolute angular difference between two proper motion position angles, normalized to the range 0°–180°.

Δθ = | θ_{B} - θ_{A} |

Uncertainty:

σ_{Δθ} = \sqrt{{σ_{A}}^{2} + {σ_{B}}^{2}}

Total Proper Motion Difference

Computes the absolute difference in total proper motion between two stars.

Δμ = | μ_{A} - μ_{B} |

Uncertainty:

σ_{Δμ} = \sqrt{{σ_{A}}^{2} + {σ_{B}}^{2}}

Velocity Differences

The following auxiliary functions compute absolute differences in tangential, radial, and total space velocity using identical mathematical formulations.

Generic Equation:

Δv = | v_{A} - v_{B} |

Uncertainty:

σ_{Δv} = \sqrt{{σ_{A}}^{2} + {σ_{B}}^{2}}

This formulation applies to:

Tangential velocity difference
Radial velocity difference
Total space velocity difference

Spectral Type Estimation

Estimates the stellar spectral class using the Gaia BP–RP color index based on empirical color boundaries:

    Spectral class | BP-RP Color-Index
    ---------------|---------------------
    O:             | bp_rp < -0,1
    B:             | -0,1  ≤ bp_rp ≤ 0,3
    A:             | 0,3   < bp_rp ≤ 0,5
    F:             | 0,5   < bp_rp ≤ 0,72
    G:             | 0,72  < bp_rp ≤ 1,14
    K:             | 1,14  < bp_rp ≤ 1,8
    M:             | bp_rp > 1,8
    -------------------------------------

The function returns one of the spectral classes: O, B, A, F, G, K, M.

This is an approximate classification intended for statistical and comparative analysis.

Mass–Luminosity Relation

Estimates stellar luminosity and mass using the absolute magnitude and an empirical mass–luminosity relation dependent on spectral type.

Luminosity:

\frac{L}{L_{☉}} = 10^{- 0.4 \cdot (M - M_{☉})}

Mass estimation:

\frac{M}{M_{☉}} = {\frac{L}{L_{☉}}}^{\frac{1}{α}}

where α is the mass–luminosity exponent corresponding to the estimated spectral type:

O => 3.5
B => 3.2
A => 3.0
F => 2.8
G => 2.6
K => 2.4
M => 2.2

Orbit Determination

The Double Star Calculator – Orbit Determination processes historical measurements of the separation and position angle of a double star and attempts to determine its orbital elements.

Optimization Summary

The Optimization Summary provides quantitative measures describing the quality, robustness, and statistical consistency of the orbit determination based on the supplied observations and measurement uncertainties. It should be interpreted as a whole. Individual metrics are most meaningful when considered together, particularly in relation to the number of observations, the orbital phase coverage, and the assumed measurement uncertainties.

N observations

Number of observations used in the orbit determination. Each observation consists of a measured separation (ρ) and position angle (θ), contributing two residuals (x and y) to the fit.

Time span

Time interval covered by the observations, given by the minimum and maximum observation epochs. A larger time span generally improves orbit determination, especially for long-period systems, as it increases orbital phase coverage.

χ² (Chi-squared)

The total chi-squared value of the fit, defined as the sum of squared, uncertainty-weighted residuals, where (x, y) are the Cartesian coordinates derived from the separation and position angle, and σ_x, σ_y are the corresponding propagated uncertainties. χ² measures the overall disagreement between the model and the observations. In this implementation, χ² is calculated from normalized Cartesian residuals.

χ² = Σ [ (x_model − x_obs) / σ_x ]² + [ (y_model − y_obs) / σ_y ]²

Degrees of freedom (dof)

Number of independent residuals available to evaluate the fit after accounting for the fitted model parameters. For N observations and seven fitted orbital elements (P, a, i, Ω, T, e, ω), the degrees of freedom are given by the equation below. A positive and sufficiently large number of degrees of freedom is required for a meaningful statistical interpretation of χ² and reduced χ².

dof = 2 · N − 7

Reduced χ²

The reduced χ² indicates how well the orbital model matches the observations relative to the assumed measurement uncertainties.

χ²_red = χ² / dof

χ²_red ≈ 1 : Residuals are consistent with the assumed uncertainties.
χ²_red ≪ 1 : Residuals are significantly smaller than the assumed uncertainties.
χ²_red ≫ 1 : Residuals are larger than expected, possibly indicating a poor fit or underestimated uncertainties.

RMS separation

Root-mean-square (RMS) residual of the separation (ρ), calculated from the differences between observed and modeled separations. This value quantifies the typical deviation of the observations from the model.

RMS (ρ) = \sqrt{[\frac{1}{N} \cdot \sum {(ρ_{model} - ρ_{obs})}^{2}]}

RMS position angle

Root-mean-square (RMS) residual of the position angle (θ), calculated from the differences between observed and modeled angles. The RMS position angle is given in degrees and reflects the typical angular discrepancy between the model and the observations.

RMS (θ) = \sqrt{[\frac{1}{N} \cdot \sum {(θ_{model} - θ_{obs})}^{2}]}

Max separation residual

Maximum absolute residual of the separation (ρ). This value highlights the largest individual deviation between an observed separation and the corresponding model prediction.

max | ρ_model − ρ_obs |

Max PA residual

Maximum absolute residual of the position angle (θ). This value highlights the largest individual deviation between an observed position angle and the corresponding model prediction.

max | θ_model − θ_obs |

Orbital phase coverage

Fraction of the orbital period covered by the observations, where P is the fitted orbital period. Values close to or exceeding one full period generally provide stronger constraints on the orbital elements, while smaller values indicate limited phase coverage and potential parameter degeneracies.

phase coverage = \frac{(t_{\max} - t_{\min})}{P}