☰Double Star Calculator - Technical Documentation

Version : 1.1 (Draft)
Author  : H.Wichmann
Date    : 06.04.2026
License : CC BY 4.0

$$\mathrm{\Delta \delta }={\delta }_{\mathrm{B}}-{\delta }_{\mathrm{A}}$$
$$\mathrm{\Delta \alpha }={\alpha }_{\mathrm{B}}-{\alpha }_{\mathrm{A}}$$
Calculation of Separation ρ :

The final angular separation is derived using the atan2(y, x) function. The atan2 implementation is preferred because it is better conditioned for all angles and avoids potential domain errors due to floating-point rounding.

Error Propagation

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

---

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.

Arguments

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

Results

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

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.

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:

---

Separation and position angle as a function of Epoch

The Double Star Calculator extrapolates the position of two stars, as well as their separation and position angle, from the Gaia DR3 reference epoch J2016.0 to user-defined target epochs (Δt = target-epoch - J2016.0). The calculation assumes a linear motion model over time.

Coordinates are projected using a linear motion model. The calculation accounts for the spherical geometry of the celestial coordinate system. The total uncertainty is calculated using Gaussian error propagation, combining the initial position error with the cumulative error of the proper motion over the elapsed time.

Right Ascension for target epoch

Declination for target epoch

Error Propagation

Where i ∈ {α, δ}

Parameters

• α_ref : Right Ascension at reference epoch J2016.0
• δ_ref : Declination at reference epoch J2016.0
• μ_α* : Proper motion in Right Ascension (μ_α ⋅ cos(δ) (mas/year))
• μ_δ : Proper motion in Declination (mas/year)
• σ_i : Standard uncertainty of the respective parameter

---

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.

Arguments

• d_A : Distance of star A
• d_B : Distance of star B
• θ : Angular separation in radians

Results

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

Spatial separation

Error propagation

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

Parallax Significance Indicator

This function calculates the Parallax Significance Indicator, referred to as the Reliability Index (RI) of the 3D separation. It quantifies whether the observed difference in parallaxes between two stars is statistically significant or merely a result of measurement uncertainties.

A high significance indicates that the stars are physically located at different distances, whereas a low significance suggests that the stars could be at the same distance within their measurement errors. The indicator is calculated by dividing the absolute difference of the parallaxes by the combined standard deviation (quadratic sum of the individual errors):

Classification Thresholds:

The Parallax Significance Indicator is mapped to qualitative labels to simplify the interpretation of the 3D data:

Parallax Significance (σ)	Label	Reliability	Description
0 ≤ σ < 3	LOW	Stat. unresolved	The parallax difference is within the 3-sigma measurement uncertainty; 3D separation is statistically unresolved.
3 ≤ σ < 5	MED	Moderate	Difference is likely real, but the 3D distance carries significant uncertainty.
σ ≥ 5	HIGH	Significant	Highly significant difference. Stars are clearly located at different depths.

Example:

Spatial Separation                 :       38.97     pc
  1-Sigma Range                    : [     36.55,    41.39] pc
  Significance                     :       13.09     σ 
  Note: High-confidence 3D 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.

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 position angle (σ_θ) is derived via Gaussian error propagation of the atan2 function. By using partial derivatives with respect to μ_α and μ_δ, the following equation provides an error estimation that accounts for the relative magnitudes of the proper motion components:

---

Total Proper Motion

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

Uncertainty of Total Proper Motion

---

Tangential Velocity

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

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

---

Space Velocity

By default, the Radial Velocity is retrieved from the Gaia DR3 record. If this data is missing, the Double Star Calculator automatically queries the SIMBAD database as a fallback. In this case, the field is labeled Radial Velocity (Simbad) and includes the corresponding Bibcode to identify the original scientific publication and its precision. If a Radial Velocity measurement is available, the total space velocity is computed as:

Uncertainty of Space Velocity

---

Results

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

---

Kinematic and Dynamical Linkage Analysis

This chapter describes algorithms to analyze the association between the two components of a double star. It details the calculation of the traditional Harshaw rPM (2D) algorithm and evaluates potential associations based on the Uncertainty Weighted Likelihood Estimator (UWLE). Furthermore, the equations describing orbital and binding dynamics are presented.

---

Harshaw 2D Methodology (rPM)

The Harshaw method (JDSO Vol. 12 No. 4 April 22, 2016, Richard W. Harshaw, CCD Measurements of 141 Proper Motion Stars: The Autumn 2015 Observing Program at the Brilliant Sky Observatory, Part 3) focuses on the ratio of Proper Motion (rPM), a dimensionless quantity comparing the angular velocity difference to the maximum magnitude. It is robust against distance uncertainties but ignores depth (parallax) and radial motion.

Right ascension: α_A, α_B (degrees)

Declination: δ_A, δ_B (degrees)

Harshaw Classification Levels

Indicator Type	Threshold	Classification
Harshaw rPM	< 0.2	CPM (Common Proper Motion)
Harshaw rPM	0.2 - 0.6	SPM (Similar Proper Motion)
Harshaw rPM	> 0.6	DPM (Distinct Proper Motion)

---

Double Star Calculator 3D Kinematics

This approach converts observed angular motion and radial velocity into a unified 3D velocity vector in km/s. By utilizing the parallax (ϖ), we account for the physical scale of the system.

Space Velocity Transformation
Relative Velocity Magnitude

This equation calculates the 3D relative velocity magnitude including the error between two stars by combining their tangential (dv_ra, dv_dec) and radial (dv_r) velocity differences. The weighting factor (w_r) acts as a binary toggle (1.0 or 0.0), allowing the formula to seamlessly switch between a full 3D space velocity and a 2D transverse velocity depending on data availability. For details refer to chapter Relative Velocity (3D) Difference.

---

Uncertainty-Weighted Likelihood Estimator

The core algorithm uses a variance-summation model. The Evidence Indicator "E" is defined by a Gaussian kernel multiplied by a normalization pre-factor that accounts for the relative weight of the measurement uncertainty.

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

Results

• 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°.

$$\mathrm{\Delta \theta }=|\theta {}_B-\theta {}_A|$$

Uncertainty:

$${\sigma }_{\mathrm{\Delta \theta }}=\sqrt{{{\sigma }_A}^2+{{\sigma }_B}^2}$$

---

Total Proper Motion Difference

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

$$\mathrm{\Delta \mu }=|\mu {}_A-\mu {}_B|$$

Uncertainty:

$${\sigma }_{\mathrm{\Delta \mu }}=\sqrt{{{\sigma }_A}^2+{{\sigma }_B}^2}$$

---

Velocity Differences

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

Generic Equation:

$$\mathrm{\Delta v}=|v{}_A-v{}_B|$$

Uncertainty:

$${\sigma }_{\mathrm{\Delta v}}=\sqrt{{{\sigma }_A}^2+{{\sigma }_B}^2}$$

This formulation applies to:

• Tangential velocity difference
• Radial velocity difference
• (Scalar) Speed difference

---

Relative Velocity (3D) Difference

The magnitude of the vectorial difference between the stars' 3D velocity vectors. Unlike the scalar version, this value accounts for the stars' actual directions of motion. It represents the total relative speed between the two components. The following function computes the 3D velocity difference.

Velocity Component Calculation:

$$v_{\alpha ,i}=k·\frac{{\mu }_{\alpha *,i}}{{\varpi }_i}\text{,}{v}_{\delta ,i}=k·\frac{{\mu }_{\delta ,i}}{{\varpi }_i}$$

Error Propagation of Individual Components:

$${\sigma }_{v_{\alpha ,i}}=\sqrt{{\left(\left(,\frac{k}{{\varpi }_i}{\sigma }_{{\mu }_{\alpha *,i}},\right)\right)}^2+{\left(\left(,\frac{k·{\mu }_{\alpha *,i}}{{\varpi }_i^2}{\sigma }_{{\varpi }_i},\right)\right)}^2}$$

Total Velocity Difference (Δv_rel):

$${\mathrm{\Delta v}}_{\mathrm{rel}}=\sqrt{{\mathrm{\Delta v}}_{\alpha }^2+{\mathrm{\Delta v}}_{\delta }^2+{(w_r·{\mathrm{\Delta v}}_r)}^2}$$

Final Error Propagation (σ_Δv):

$${\sigma }_{{\mathrm{\Delta v}}_{\mathrm{rel}}}=\sqrt{{\left(\left(,\frac{{\mathrm{\Delta v}}_{\alpha }·{\sigma }_{{\mathrm{\Delta v}}_{\alpha }}}{{\mathrm{\Delta v}}_{\mathrm{rel}}},\right)\right)}^2+{\left(\left(,\frac{{\mathrm{\Delta v}}_{\delta }·{\sigma }_{{\mathrm{\Delta v}}_{\delta }}}{{\mathrm{\Delta v}}_{\mathrm{rel}}},\right)\right)}^2+{\left(\left(,\frac{w_r·{\mathrm{\Delta v}}_r·{\sigma }_{{\mathrm{\Delta v}}_r}}{{\mathrm{\Delta v}}_{\mathrm{rel}}},\right)\right)}^2}$$

Nomenclature & Constants

Symbol	Description	Value / Unit / Definition
$k$	AU to km/s conversion factor	4.74057
$w_r$	Radial velocity weighting factor	1.0 (enabled) or 0.0 (disabled)
$\varpi$	Parallax	mas
${\mu }_{\alpha *}$	Proper Motion in Right Ascension (μα · cos δ)	mas/yr
${\mu }_{\delta }$	Proper Motion in Declination	mas/yr
$\sigma$	Standard uncertainty (Sigma)	Associated error of a variable
$v_{\tan }$	Tangential velocity component	km/s
$v_r$	Radial velocity component	km/s
${\mathrm{\Delta v}}_{\mathrm{rel}}$	3D Space Velocity Difference	km/s (Vectorial Difference)

---

Spectral Type Estimation

Estimates the stellar spectral class and temperature using the Gaia BP–RP color index based on empirical color boundaries. The temperatures are derived from data provided by the University of Northern Iowa / https://sites.uni.edu/morgans/astro/course/Notes/section2/spectraltemps.html.

    Spectral class |  BP-RP  Color-Index    |  Temperature 
    ---------------|--------------------------------------
    O:             |         bp_rp < -0,1   |    50000
    B:             | -0,1  ≤ bp_rp ≤  0,3   |    15000
    A:             |  0,3  < bp_rp ≤  0,5   |     8750
    F:             |  0,5  < bp_rp ≤  0,72  |     6700
    G:             |  0,72 < bp_rp ≤  1,14  |     5800
    K:             |  1,14 < bp_rp ≤  1,8   |     4800
    M:             |         bp_rp >  1,8   |     3300
    ------------------------------------------------------

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

This chapter describes the estimation of stellar luminosity and mass for main-sequence stars based on their absolute magnitude and spectral type. The implemented method follows the segmented mass–luminosity relations presented by Eker et al. (2018), adapted here to broad spectral-type intervals for practical application. The luminosity is derived from the absolute magnitude relative to the Sun. The stellar mass is subsequently computed by inverting the logarithmic mass–luminosity relation.

Luminosity:

The stellar luminosity expressed in units of the solar luminosity where $M$ is the absolute magnitude of the star and $M_{☉}$ is the solar absolute magnitude, adopted here as $M_{☉}=\; 4.74$. A Bolometric correction factor acc. to Pecaut and Mamajek (2013) is applied:

	"O" => -3.5
	"B" => -1.4
	"A" => -0.05
	"F" => -0.05
	"G" => -0.15
	"K" => -0.70
	"M" => -3.30
	

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

Given an uncertainty ${\sigma }_M$ in the absolute magnitude, lower and upper luminosity limits are computed as:

$$L_{\min }={10}^{-0.4\cdot (M+{\sigma }_M-M_{☉})}$$ $$L_{\max }={10}^{-0.4\cdot (M-{\sigma }_M-M_{☉})}$$

Mass–Luminosity Relation and Mass estimation:

Eker et al. (2018) describe the mass–luminosity relation for main-sequence stars as a segmented power law:

$${\log }_{10}\left(L\right)=\alpha \cdot {\log }_{10}\left(M\right)+C$$

Here, $\alpha$ is the slope and $C$ the intercept of the logarithmic relation, both depending on the stellar mass regime. For mass estimation, the relation is inverted to yield:

$$M={10}^{\frac{{\log }_{10}\left(L\right)-C}{\alpha }}$$ Notes and Assumptions

• The coefficients $\alpha$ and $C$ are assigned according to broad spectral-type intervals, mapped from the segmented relations given by Eker et al. (2018). This approach provides a practical approximation for main-sequence stars when detailed stellar parameters are unavailable
  

      "O" => { α: 3.96, C:  0.00 }, # High Mass
      "B" => { α: 3.96, C:  0.00 }, # High Mass
      "A" => { α: 4.67, C: -0.15 }, # Intermediate Mass
      "F" => { α: 4.67, C: -0.15 }, # Intermediate Mass
      "G" => { α: 5.56, C: -0.03 }, # Low Mass
      "K" => { α: 4.20, C: -0.20 }, # Ultra-Low Mass
      "M" => { α: 2.27, C:  0.60 }  # Ultra-Low Mass

• The method is valid for main-sequence stars only
• The mass estimate is computed from the nominal luminosity value
• Luminosity limits reflect the uncertainty in absolute magnitude
• Spectral-type boundaries are approximate and intended for practical use

---

Hertzsprung-Russell Diagram

The positions of both components are plotted in a simplified Hertzsprung–Russell diagram (HRD). The positions are derived from their calculated stellar luminosities, expressed in solar luminosities and converted to logarithmic scale, and from their effective temperatures (T_eff) as provided by Gaia (plotted as solid circles). The effective temperatures are taken from the Gaia GSP-Phot Aeneas solution based on BP/RP spectra (teff_gspphot) and are likewise converted to logarithmic scale. In case the effective temperature is not available, the temperature is estimated from the color index (bp_rp) measured by Gaia (plotted as open circles). Upper and lower bounds of both luminosity and effective temperature are propagated into logarithmic space to derive symmetric uncertainties for the HRD coordinates.

In addition to the stellar positions, the HRD includes schematic regions indicating the approximate locations of the main sequence, giant stars, supergiants, and white dwarfs. These regions are intended for qualitative orientation within the HR diagram and are based on standard tabulated relations between spectral type and effective temperature, such as those provided by the University of Northern Iowa / https://sites.uni.edu/morgans/astro/course/Notes/section2/spectraltemps.html.

---

---

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.

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 χ².

Reduced χ²

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

• χ²_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 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.

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 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.

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.

---

History

Version 1.0 : https://doi.org/10.5281/zenodo.19160070

---

Disclaimer

The algorithms used in the Double Star Calculator were developed and implemented by me to the best of my knowledge and belief. I have thoroughly tested the tool, and it provides plausible results based on GAIA data. However, I do not warrant the accuracy, completeness, or correctness of the results. This tool is provided "as is" without any guarantees regarding its suitability for a specific purpose or non-infringement of third-party rights.

The use of this tool is entirely at your own risk. I assume no liability for any damages, including but not limited to direct, indirect, incidental, or consequential damages, arising from the use of this tool or any errors or inaccuracies in the results.

This tool is intended for informational and research purposes only and is not a substitute for professional scientific or astronomical analysis. I reserve the right to modify, update, or discontinue the tool at any time without prior notice.

If the Double Star Calculator has been helpful in your work, a citation or acknowledgment would be greatly appreciated.

Privacy Notice: This tool does not collect, store, or process any personal data from users.

Copyright (C) 2026. The documentation is licensed under CC BY 4.0. Unless otherwise noted, all output generated by the Double Star Calculator including orbit diagrams, ephemerides, and analytical reports is licensed under CC BY 4.0.

Acknowledgements

The Double Star Calculator makes use of data from the European Space Agency (ESA) mission Gaia, processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Additionally, it incorporates data from the Washington Double Star Catalog, maintained at the U.S. Naval Observatory, makes use of the Aladin Sky Atlas developed at CDS, Strasbourg Observatory, France, and makes use of the SIMBAD database, operated at CDS, Strasbourg, France.

---

www.stella-vega.de