Subroutines of the Wilson Devinney Program

Finis coronat opus (The end crowns the work)

The WD program is a practical expression of the WD model and there have been many more publications about the model than about the program. Even the best programs ordinarily have much shorter lifetimes than good ideas, so the model has been kept conceptually separate from its software implementation. Thus papers on the model essentially never mention the names of subroutines or main programs - only papers specifically about the program [e.g., Wilson (1993)] use those names. This appendix tries to provide some of the details which might improve the reader's understanding of the program and to connect the subroutines to the relevant features in the model.

The program's history has been one of various special-purpose versions that were developed for particular problems, followed by absorption of their capabilities into the general program. Overall, the idea has been to have one "direct problem" program (LC) and one "inverse problem" program (DC), each with multiple capabilities, rather than a library of special purpose programs. For example, the 1998 version computed light curves, radial velocity curves, spectral line profiles, and images, while versions that compute polarization curves and X-ray pulse arrival times now exist separately and may eventually be absorbed. The WD 2007 version incorporated improved atmosphere simulation, parameters of eclipses by clouds in extensive stellar atmospheres of hot stars, and third-body parameters. Generalizations that a user need not worry about are embedded invisibly wherever practical. For example, computational shortcuts for many special case situations speed execution without compromising more intricate cases.

The overall structure of the Wilson-Devinney program and its main routines LC and DC is shown in Figs. E.2 and E.3. The purpose of some subroutines is explained in Wilson (1993) and in Wilson's program manuals. Other subroutines lack a detailed description of the mathematics. Although the overall ideas of computing EB light curves usually can be presented in an elegant and systematic manner, the implementation of these ideas in practice requires a substantial amount of numerical analysis. An example is the problem of computing accurate derivatives and providing them with a derivative-based least-squares solver. This topic is related to the accuracy of representing the model in the computer (both the discretization of the surface of the components represented by a finite grid and the finiteness of the number representation in the machine).

A subroutine having an asterisk in its name is part of WD2 0 0 9, but not of the 2008 version of the Wilson-Devinney program.

E.1 ATM - Interfacing Stellar Model ATMospheres

This subroutine was the interface to a stellar model atmosphere. The original Wilson-Devinney program uses the Carbon-Gingerich model atmospheres. In WDx2 0 07, this subroutine is replaced by a routine written by C. Stagg (Milone et al. 1992b), implementing the Kurucz atmospheres. In both cases it returns the ratio between the flux based on the atmosphere and the blackbody law. The version of Jan. 23, 2004, and later versions incorporate the models of Kurucz (1993) in an use of external tables of metallicity, temperature, and log g to provide modern atmosphere simulations. The tables are automatically read in the 2007 version.

E.2 ATMx - Interfacing Stellar Model ATMospheres

The WD subroutine interfacing to Kurucz atmosphere is for the ATMx. It performs a four-point interpolation in log g and then an m-point Lagrangian interpolation. This routine exploits an a priori computed file that contains a block for each of the 19 compositions, each block listing the temperature limits and Legendre coefficients for every band, log g, and temperature sub-interval. With 11 log g's, 25 bands, 19 compositions, 4 temperature subintervals, and 10 Legendre coefficients with 2 temperature limits per subinterval, the data file contains 250,800 numbers.

Where grid elements require elements outside the grid, a tranfer is made to a Planckian through a "ramp" function that provides smooth continuity. That way Van Hamme & Wilson (2003) established an atmosphere to blackbody transition region in Teff and in log g and avoid any discontinuity. If the Teff, log g combination is outside the range of atmosphere applicability, the program smoothly connects atmosphere model intensities to bandpass blackbody intensities over built-in ranges in log g and Teff whose limits can easily be changed. Similarly to the files above, a Legendre blackbody file spans 500-500,000 K and is very small, as the only dimensions are bandpass and temperature sub-interval. Note that the files can be incorporated into other binary star programs to calculate model and blackbody intensities.

E.3 BBL - Basic BLock

Subroutine BBL is the so-called Basic BLock of the WD program (see Fig. E.1). It keeps the computation of numerical derivatives as simple as possible and at the same time tries to avoid redundant computations.

Fig. E.1 Structure of subroutine BBL. Courtesy R. E. Wilson

E.4 BinNum - A Search and Binning Utility

BinNum is a tool to find the bin in which a number is located. It is essentially similar to LOCATE in the Numerical Recipes by Press et al. (1992), which searches an ordered table by bisection.

E.5 BOLO - Bolometric Corrections

Subroutine BOLO uses Harris's (1963) calibration from T = 3, 970 K to 5, 800 K, Morton & Adams (1968) from T = 5, 800 to 37, 500 K, and the blackbody law (3.2.21) below T = 3, 970 K and above T = 37, 500 K to compute bolometric corrections. These corrections are needed to compute the ratio of bolometric luminosities involved in the reflection effect in Sect. 3.2.5.

E.6 CofPrep - Limb-Darkening Coefficient Preparation

Subroutine COFPREP computes the coefficients used in the limb-darkening interpolation scheme. It reads the numbers from the table selected based on input metal-licity into an array which is then available to LIMDARK.

E.7 CLOUD - Atmospheric Eclipse Parameters

Added in the 2002 version of the WD program, subroutine CLOUD computes atmospheric eclipse parameters as described in Sect. 3.4.4.5.

E.8 CONJPH - Conjunction Phases

This subroutine computes the phases of superior and inferior conjunctions. For eccentric orbits this becomes a relevant and subtle issue as described on page 88.

E.9 DGMPRD - Matrix-Vector Multiplication

Subroutine DGMPRD performs matrix multiplication and returns the resulting vector r = Ab, where matrix A ist stored in a one-dimensional chain and b is the input vector.

E.10 DMINV - Matrix Inversion

Subroutine DMINV computes the inverse of a n by n matrix A stored in a one-dimensional chain.

E.11 DURA - Constraint on X-Ray Eclipse Duration

This subroutine puts an explicit constraint on the size of a star based on the duration of an X-ray eclipse. Such constraints may be considered when X-ray eclipses of neutron stars, black-holes, or white dwarfs occur, as described in Wilson (1979) where the full mathematics is given. The basic observable is the semi-duration &e of an X-ray eclipse. For circular orbits and synchronous rotation a relation

has been derived by Chanan et al. (1977) which relates the inclination i, the Roche potential Q, the mass ratio q, and &e. The more general eccentric, nonsynchronous case is in Wilson (1979).

E.12 ELLONE - Lagrangian Points and Critical Potentials

This subroutine computes the x-coordinates of the equilibrium points L1 and Lp and the associated critical Roche potentials Q™1 and Q™1. The name equilibrium points is used here as a generalization of Lagrangian points, as used for the synchronous, circular case. The required input quantities are the mass ratio q = M2/M1, the ratio F1 = «i/«, and the distance d between the stars in units of the semi-major axis a of the relative orbit. The x-coordinate of the equilibrium points follows from the condition d^

As the equilibrium points are located on the line connecting the centers of the components in (6.3.2), the distances of a point x on that line take a special form. The star centers are at positions x1 = 0 and x2 = d. The coordinate xLp of the equilibrium point L1 fulfills the relation

so that its distances to the component centers are xLp and d - xLp. Therefore, for xL1p , the potential (6.3.2) takes the form

Telescopes Mastery

Telescopes Mastery

Through this ebook, you are going to learn what you will need to know all about the telescopes that can provide a fun and rewarding hobby for you and your family!

Get My Free Ebook


Post a comment