\documentstyle[twoside,epsf,hyperref]{article}
\input ibvs2.sty

\begin{document}
\IBVSheadDOI{63}{6224}{11 December 2017}

\IBVStitle{V500 C\lowercase{yg} -- A Classical Algol}


\IBVSauth{Nelson, Robert H.$^{1,2,3}$}

\IBVSinst{Mountain Ash Observatory, 1393 Garvin Street, Prince George, BC, Canada, V2M 3Z1, bob.nelson@shaw.ca}
\IBVSinst{Guest investigator, Dominion Astrophysical Observatory, Herzberg Institute of Astrophysics, National Research Council of Canada}
\IBVSinst{Desert Blooms Observatory, Benson AZ, 31\deg 56.454\arcm N,  110\deg 15.450\arcm W}

\SIMBADobj{V500 Cyg}


\begintext

The discoverer of the variability of V500 Cyg (AN 1939.0081; TYC 2693--139--1) appears to be undocumented.  The first available reference (in the GCVS and SIMBAD) is Whitney (1959) who provided revised elements, three new eclipse timings, and notes regarding a companion separated by 0.3\arcm.  Since then, there have been numerous eclipse timings published, but no light curve or analysis.

In order to rectify this lack, the author first secured, in the autumns of 2010, 2013, 2014, and 2015,  a total of eight medium resolution (R $\sim$ 10000 on average) spectra of V500 Cyg at the Dominion Astrophysical Observatory (DAO) in Victoria, British Columbia, Canada using the Cassegrain spectrograph attached to the 1.85 m Plaskett Telescope.  He used the 21181 configuration and a grating with 1800 lines/mm, blazed at 5000 \AA, and giving a reciprocal linear dispersion of 10 \AA/mm in the first order.  The wavelengths ranged from 5000 to 5260 \AA, approximately.  A log of observations is given in Table 1 and an eclipse timing diagram, in Figure 9 later in the paper.  The latter was used to derive the following elements, used for both radial velocity (RV) and photometric phasing:

\begin{equation}
{\rm JD~(Hel)~Min~I} = 24 57914.8640 (49) + 0.9242233 (2) E
\end{equation}

\noindent where the quantities in brackets are the standard errors of the preceding quantities in units of the last digit.

\begin{table}
\centering
\caption{Log of DAO observations}
\begin{tabular}{lccccc}
\hline
DAO&Mid Time&Exposure&Phase at&V1&V2\\
Image \#&(HJD$-$2400000)&(sec)&Mid$-$exp&(km/s)&(km/s)\\ 
\hline
10-17392&55474.7097&3600&0.778&77.4 (2.8)&$-$215.3 (14.8)\\
13-09641&56544.8987&3600&0.712&74.1 (4.2)&$-$225.2 (10.8)\\
12-24533&56912.6665&3600&0.633&42.3 (1.3)&$-$196.3 (0.9)\\
15-13142&57295.8492&3600&0.232&$-$123.5 (4.8)&159.9 (10.7)\\
15-13144&57295.8926&3600&0.279&$-$126.1 (5.0)&174.0 (16.5)\\
15-13176&57296.8290&3600&0.292&$-$120.4 (4.6)&163.3 (13.5)\\
15-13238&57298.7427&3600&0.363&$-$94.6 (2.6)&104.2 (7.0)\\
15-13265&57299.6278&3600&0.321&$-$113.8 (4.8)&134.6 (16.3)\\ 
\hline
\end{tabular}
\end{table}

Frame reduction was performed by software {\sc RaVeRe} (Nelson 2013).  See Nelson (2010) and Nelson et al.\ (2014) for further details.  The normalized spectra are reproduced in Fig.~1, sorted by phase (the vertical scale is arbitrary).  Note towards the right the strong neutral iron lines (at 5167.487 and 5171.595 \AA) and the strong neutral magnesium triplet (at 5167.33, 5172.68, and 5183.61 \AA). 

\IBVSfig{8cm}{6224-f1.eps}{V500 Cyg spectra at phases 0.232, 0.279, 0.292, 0.321, 0.363, 0.633, 0.712, 0.778 (from top to bottom).  Each has been shifted vertically for clarity.  The vertical scale is arbitrary.}
\IBVSfigKey{6224-f1.eps}{V500 Cyg}{spectrum}

Radial velocities were determined using the Rucinski broadening functions (Rucinski, 2004, Nelson, 2010) as implemented in software Broad25 (Nelson, 2013).  See Nelson et al.\ (2014) for further details. An Excel worksheet with built-in macros (written by him) was used to do the necessary radial velocity conversions to geocentric and back to heliocentric values (Nelson 2014).  The resulting RV determinations are also presented in Table 1 (along with standard errors in units of the last digits, enclosed in brackets).  The mean rms errors for RV$_1$ and RV$_2$ are 3.8 and 11.3 km/s, respectively, and the overall rms deviation from the (sinusoidal) curves of best fit is 9.7 km/s. The best fit yielded the values $K_1 = 98.6 (2.7)$ km/s, $K_2 = 196.8 (4.9)$ km/s and $V_{\gamma} = -129.1 (2.2)$ km/s, and thus a mass ratio $q_{\rm sp} = K_1/K_2 = M_2/M_1 = 0.50 (1)$.

Representative broadening functions, at phases 0.232 and 0.778 are depicted in Figs.~2 and 3, respectively (the vertical scale is arbitrary).  Smoothing by a Gaussian filter is routinely done in order to centroid the peak values for determining the radial velocities.

\IBVSfig{5cm}{6224-f2.eps}{Broadening functions at phase 0.232--smoothed and unsmoothed.}
\IBVSfigKey{6224-f2.eps}{V500 Cyg}{radial velocity}

\IBVSfig{6cm}{6224-f3.eps}{Broadening functions at phase 0.778--smoothed and unsmoothed.}
\IBVSfigKey{6224-f.eps}{V500 Cyg}{radial velocity}

During twelve nights in 2017, May 24 -June 14, the author took a total of 198 frames in $V$, 197 in $R_{\rm C}$ (Cousins) and 199 in the $I_{\rm C}$ (Cousins) band at the newly-opened Desert Blooms Observatory, jointly owned by the author and Dr.~Kevin B.\ Alton.  Hosted at the San Pedro Observatory complex located near Benson, Arizona, the telescope is operated remotely. It consists of a Software Bisque Taurus 400 equatorial fork mount, a Meade LX-200 40 cm Schmidt-Cassegrain optical assembly operating at f/7, a SBIG STT-1603 XME CCD camera (with a field of view $11\arcm\times18\arcm$), and a filter wheel with the usual $B$, $V$, $R_{\rm C}$, and $I_{\rm C}$ filters.  For unattended operation, automatic focusing is required owing to the large temperature changes throughout the night (typically +35 \deg C to +10 \deg C in late spring).

Standard reductions were then applied (see Nelson et al.\ 2014 for more details).  The variable, comparison, and check stars are listed in Table 2.  The coordinates and magnitudes for V500 Cyg, the comparison, and check 2 are from the Tycho Catalogue, Hog et al.\ (2000), with magnitudes converted to standard Johnson values using relations due to Henden (2001).  For check 1, the $V$ magnitude is from the GSC catalogue and the approximate $B-V$ value is from our photometry.  Quantities in brackets are standard errors, in units of the last digit.

\begin{table}
\centering
\caption{Details of variable, comparison and check stars.}
\begin{tabular}{llllll}
\hline
Object&GSC&RA (J2000)&Dec (J2000)&$V$ (mag)&$B-V$ (mag)\\
\hline
Variable&2693-0139&20\hr 24\mm 40\fsec379&+34\deg 57\arcm 05\farcs40&11.91 (16)&+0.25 (21)\\
Comparison&2693-0828&20\hr 24\mm 39\sec&+34\deg 56\arcm 59\arcs&11.20 (7)&0.22 (9)\\
Check 1&2693-1630&20\hr 24\mm 28\sec&+34\deg 55\arcm 45\arcs&12.1&0.34\\
Check 2&2693-1230&20\hr 24\mm 16\fsec9528&+34\deg 58\arcm 39\farcs642&10.91 (7)&+0.73 (14)\\
\hline
\end{tabular}
\end{table}

The author used the 2003 version of the Wilson-Devinney (WD) light curve and radial velocity analysis program with Kurucz atmospheres (Wilson \& Devinney, 1971, Kurucz, 1979, Wilson, 1990, Kallrath \& Milone, 1998, Wilson, 1998) as implemented in the Windows front-end software WDwint (Nelson, 2013) to analyze the data.  To get started, the spectral type F4--5 (taken from SIMBAD, no reference given; main sequence assumed) was adopted.  Interpolated tables from Flower (1996) gave a temperature $T_1 = 6610 \pm 134$ K ($T_1$ is the mean of the two sub-classes)  and $\log g = 4.348 \pm 0.014$.  (The quoted errors refer to one spectral sub-class.)  An interpolation program by Terrell (1994, available from Nelson 2013) gave the Van Hamme (1993) limb darkening values; and finally, a logarithmic (LD = 2) law for the limb darkening coefficients was selected, appropriate for temperatures  $<$ 8500 K (ibid.).  The limb darkening coefficients are listed in Table 3.  (The values for the second star are based on the later-determined temperature of 4584 K and assumed spectral type of K5.)  Convective envelopes for both stars were used, appropriate for cooler stars (hence values gravity exponent $g = 0.32$ and albedo $A = 0.500$ were used for each).


\begin{table}
\centering
\caption{Limb darkening values from Van Hamme (1993).}
\begin{tabular}{lcccc}
\hline
Band&$x_1$&$x_2$&$y_1$&$y_2$\\
\hline
Bol&0.640&0.628&0.242&0.150\\
$V$&0.707&0.797&0.278&0.015\\
$R_{\rm C}$&0.634&0.753&0.286&0.104\\
$I_{\rm C}$&0.550&0.667&0.276&0.150\\
\hline
\end{tabular}
\end{table}

From the GCVS 4 designation (EA/SD) and from the shape of the light curve, mode 5 (classical Algol) mode was used.  Later on, mode 2 (detached) was tried but DC adjustments required decreases in potential 2 below the critical value; consequently mode 2 was abandoned.

Convergence using differential corrections (DC) and the method of multiple subsets was reached in a small number of iterations.  (See Wilson \& Devinney, 1971 and Kallrath \& Milone 1998 for an explanation of the method.)  The subsets were: ($a$, $V_{\gamma}$, $i$, $L_1$), ($T_2$, $q$), and ($T_2$, $\Omega_1$). However, the visual fit was poor in that the calculated depth of the secondary minimum was too deep.  Therefore, in LC mode temperature $T_2$ was lowered until the fit was satisfactory.  Then, switching back to DC mode, temperature $T_2$ was held constant while all other parameters allowed to vary.  Once convergence was obtained, $T_2$ was again allowed to vary with only small changes indicated.

Detailed reflections were tried, with the number of reflections, $n_{\rm ref} = 3$, but there was little--if any--difference in the fit from the simple treatment.

The model is presented in Table 4 (for an explanation of column 3, see below).  For the most part, the error estimates are those provided by the WD routines and are known to be under-estimated; however, it is a common practice to quote these values and we do so here.  Also, estimating the uncertainties in temperatures $T_1$ and $T_2$ is somewhat problematic.  A common practice is to quote the temperature difference over--say--one spectral sub-class (assuming that the classification is good to one spectral sub-class, the precision being unknown in this case).  In addition, various different calibrations have been made (Cox, 2000, page 388--390 and references therein, and Flower, 1996), and the variations between the various calibrations can be significant.  If the classification is $\pm$ one sub-class, an uncertainty of $\pm$ 200 K to the absolute temperatures of each, would be reasonable.  The modelling error in temperature $T_2$, relative to $T_1$, is indicated by the WD output to be much smaller, around 9 K.)

\begin{table}
\centering
\caption{Wilson--Devinney parameters.}
\begin{tabular}{lcccc}
\hline
WD Quantity&Value&Revised values&error&Unit\\
\hline
Temperature, $T_1$&6610&6610&[fixed]&K\\
Temperature, $T_2$&4584&4594&200&K\\
$q = m_2/m_1$&0.557&0.554&0.005&---\\
Potential, $\Omega_1$&3.703&3.690&0.015&---\\
Potential, $\Omega_2$&2.984&2.978&[fixed]&\\
Inclination, $i$&83.06&83.38&0.10&degrees\\
Semi-major axis $a$&5.38&5.38&0.12&solar radii\\
$V_\gamma$&$-$25.3&$-$25.3&2.6&km/s\\
Fill-out, $f_1$&$-$2.186&$-$2.177&0.001&\\
$L_1/(L_1+L_2)$ (V)&0.8664&0.8664&0.0003&---\\
$L_1/(L_1+L_2)$ (R$_{\rm C}$&0.8245&0.8245&0.0004&---\\
$L_1/(L_1+L_2)$ (I$_{\rm C}$)&0.7866&0.7866&0.0006&---\\
$r_1$ (pole)&0.3153&0.3153&0.0015&orbital radii\\
$r_1$ (point)&0.3377&0.3377&0.0022&orbital radii\\
$r_1$ (side)&0.3234&0.3234&0.0017&orbital radii\\
$r_1$ (back)&0.3317&0.3317&0.0019&orbital radii\\
$r_2$ (pole)&0.3083&0.3083&0.0007&orbital radii\\
$r_2$ (point)&0.4402&0.4402&0.0027&orbital radii\\
$r_2$ (side)&0.3220&0.3220&0.0007&orbital radii\\
$r_2$ (back)&0.3544&0.3544&0.0007&orbital radii\\
Phase shift&0.0011&0.0016&0.0001&---\\
$\Sigma\omega_{\rm res}^2$&0.06012&0.03943&---&---\\
\hline
\end{tabular}
\end{table}

The light curve data and the fitted curves are depicted in Figures 4--6.  The residuals (in the sense observed-calculated) are also plotted, shifted upwards by 0.25 units.  


It is not clear why, in all three light curves, a few points near phase 0.03 (and all from the same night) are deviant, other than possibly due to a passing cloud which could have differentially affected the flux from one of the stars (variable, comparison) compared to the other.  In response to a referee's concerns about these errant points, new modelling trials were undertaken with these points deleted.  The result was slight differences in the resultant parameters at convergence; these are reported in column 3.  The reader will note that, for the most part, these lie inside the estimated (one sigma) confidence intervals and are therefore not significantly different.


\IBVSfig{6cm}{6224-f4.eps}{$V$ light curves for V500 Cyg  -- data, WD fit, and residuals.}
\IBVSfigKey{6224-f4.eps}{V500 Cyg}{light curve}

\IBVSfig{6cm}{6224-f5.eps}{$R$ light curves for V500 Cyg  -- data, WD fit, and residuals.}
\IBVSfigKey{6224-f5.eps}{V500 Cyg}{light curve}

\IBVSfig{6cm}{6224-f6.eps}{$I$ light curves for V500 Cyg  -- data, WD fit, and residuals.}
\IBVSfigKey{6224-f6.eps}{V500 Cyg}{light curve}

The radial velocities are shown in Fig.~7.  A three-dimensional representation from Binary Maker 3 (Bradstreet, 1993) is shown in Fig.~8.  (The crosses are the centres of mass of the individual stars and of the system as a whole.  The ellipses are of the respective centres of mass.)



The WD output fundamental parameters and errors are listed in Table 5.  Most of the errors are output or derived estimates from the WD routines.  From Kallrath \& Milone (1998), the fill-out factor is f = ($\Omega_{\rm I}-\Omega)/ (\Omega_{\rm I}-\Omega_{\rm O})$, where $\Omega$ is the modified Kopal potential of the system, $\Omega_{\rm I}$ is that of the inner Lagrangian surface, and $\Omega_{\rm O}$, that of the outer Lagrangian surface, was also calculated.


To determine the distance, the analysis proceeded as follows: first the WD routine gave the absolute bolometric magnitudes of each component; these were then converted to the absolute visual ($V$) magnitudes of both, $M_{\rm V,1}$ and $M_{\rm V,2}$, using the bolometric corrections BC = $-0.135$ and $-0.72$ for stars 1 and 2 respectively.  The latter were taken from interpolated tables constructed from Cox (2000).  The absolute $V$ magnitude was then computed in the usual way, getting $M_V = 2.63 \pm 0.06$ magnitudes.  The apparent magnitude in the $V$ passband was $V = 11.93 \pm 0.02$, taken from the Tycho values (Hog et al.\ 2000) and converted to the Johnson magnitude $11.91 \pm 0.02$ using relations due to Henden (2001).


\begin{table}
\centering
\caption{Fundamental parameters.}
\begin{tabular}{lccc}
\hline
Quantity&Value&Error&unit\\
\hline
Temperature, $T_1$&6610&200&K\\
Temperature, $T_2$&4584&200&K\\
Mass, $m_1$&1.58&0.10&M0\\
Mass, $m_2$&0.88&0.04&M0\\
Radius, $R_1$&1.74&0.01&R0\\
Radius, $R_2$&1.77&0.01&R0\\
$M_{\rm bol,1}$&3.00&0.02&mag\\
$M_{\rm bol,2}$&4.55&0.02&mag\\
$\log g_1$&4.15&0.01&cgs\\
$\log g_2$&3.88&0.01&cgs\\
Luminosity, $L_1$&5.20&0.10&L0\\
Luminosity, $L_2$&1.25&0.02&L0\\
Fill-out factor 1&$-$2.219&0.010&---\\
Fill-out factor 2&0&[fixed]&\\
Distance, $r$&602&27&pc\\
\hline
\end{tabular}
\end{table}

Ignoring interstellar absorption, we calculated a preliminary value for the distance $r = 717$ pc from the standard relation:	

\begin{equation}
r = 10^{0.2(V - M_V -  A_V +  5)}	{\rm~parsecs}
\end{equation}

Galactic extinction was obtained from a model by Am\^{o}res \& L\'epine (2005).  The code (available in IDL and converted by the author to a Visual Basic routine) assumes that the interstellar dust is well mixed with the gas, that the Galaxy is axisymmetric, that the gas density in the disk is a function of the Galactic radius and of the distance from the Galactic plane, and that extinction is proportional to the column density of the gas,  Using Galactic coordinates of $l = 74.0787\deg$ and $b = -1.5709\deg$ (SIMBAD), and the initial distance estimate of $d = 0.717$ kpc, a value of $A_V = 0.451$ mag was determined,  Further iteration of several steps resulted in final values of $A_V = 0.382$ mag and $r = 602$ pc.

The errors were assigned as follows: $\delta M_{\rm bol,1} = \delta M_{\rm bol,2} = 0.02$, $\delta {\rm BC}_1$ = $\delta {\rm BC}_2 = 0.09$ (the variation of 1 spectral sub-class), $\delta V = 0.02,$ $\delta A_V = 0.02$, all in magnitudes.  Combining the errors rigorously (i.e., by adding the variances) yielded an estimated error in $r$ of 27 pc.

\IBVSfig{7cm}{6224-f7.eps}{Radial velocity curves for V500 Cyg -- data and WD fit.}
\IBVSfigKey{6224-f7.eps}{V500 Cyg}{radial velocity}

\IBVSfig{4.5cm}{6224-f8.eps}{Binary Maker 3 representation of the system -- at phases 0.48 and 0.75.}
\IBVSfigKey{6224-f8.eps}{V500 Cyg}{other}

Four new times of minima emerged from the observations; these are reported in Table~6.  Each is the mean of three values (one for each filter).  Four methods of minimum determination, as implemented in software Minima23 (Nelson 2013), were used: the digital tracing paper method, sliding integrations (Ghedini 1982), curve fitting using five Fourier terms, and Kwee and van Woerden (Kwee \& Woerden 1956, Ghedini 1982).  Because, in the literature, many (or perhaps most) error estimates can be shown to be low (sometimes unrealistically so), the estimated errors were taken as double the standard deviations of the various determinations.

\begin{table}
\centering
\caption{New times of minima for V500 Cyg obtained in this study.}
\begin{tabular}{lcc}
\hline
Min (Hel)$-$2400000&Type&Error (days)\\
\hline
57901.9264&I&0.0002\\
57908.8590&II&0.0006\\
57913.9397&I&0.0004\\
57914.8639&I&0.0009\\
\hline
\end{tabular}
\end{table}

Some comments regarding the period variation are in order.   An eclipse timing difference (O--C) plot using timings from 1988 is depicted in Fig.~9.  Although there is considerable scatter, a linear relation over the data collection interval (cycles 28800 to 30770 for the RVs and cycles 31420 to 31440 for the light curve data) is assumed.  This yielded a weighted best-fit linear solution and ephemeris of Equation (1) above.  (Standard weighting was used:  pg = 0.2, vis = 0.1, and PE = CCD = 1.  Two nearly identical points lying more than three standard deviations from the curve of best fit were rejected.)

\IBVSfig{8cm}{6224-f9.eps}{V500 Cyg -- eclipse timing (O--C) diagram with linear (solid blue) and quadratic (dashed red) fits for points after cycle 20000 (see equation 1). (Note:  pg = photographic; vis = visual; PE = photoelectric; and CCD = charge coupled device.}
\IBVSfigKey{6224-f9.eps}{V500 Cyg}{O-C diagram}

\IBVSfig{8cm}{6224-f10.eps}{V500 Cyg -- eclipse timing (O--C) diagram with a quadratic fit for all available points.}
\IBVSfigKey{6224-f10.eps}{V500 Cyg}{O-C diagram}

Also, all the available timing data since the earliest in 1935 (available online at Nelson 2016) are plotted in Fig.~10. There may well be a quadratic relation;  the relevant parameters for which are given in Equation 3.

\begin{equation}
{\rm JD~(Hel)~Min~I} = 24 57914.8651 (29) + 0.9242105 (5) + 2.1 (2) \times 10^{-10}  E^{2}
\end{equation}

However, the quadratic relation does not fit the data since cycle 20000 particularly well (see Fig.~9) and was not used in the analysis.  The period behaviour might perhaps be better explained by the light time effect (LiTE; Irwin 1952, 1959) due to a third star.  However, due to the obvious scatter in the early photographic data near cycle 0, (due to Wachmann, cited in the O--C Gateway with only the ambiguous reference of AAAN 11.5.43), a LiTE analysis does not appear to be justified at this time.  High quality data over the coming decades will be required to settle the matter.  The reader is referred to Nelson et al.\ (2014, 2015, 2016) for further discussions on this difficulties encountered in period analysis.


\noindent{\bf Acknowledgements:}
It is a pleasure to thank the staff members at the DAO (Dmitry Monin, David Bohlender, and the late Les Saddlmyer) for their usual splendid help and assistance.  Many thanks are also due to the San Pedro Observatory resident astronomer/technician Dean Salman for his tireless help.  Much use was made of the SIMBAD database during this research.

\references

Am\^{o}res, E.B., L\'epine, J.R.D., 2005, {\it AJ}, {\bf 130}, 659 \DOI{10.1086/430957}

Bradstreet, D.H., 1993, {\it IAUCB}, {\bf 21}, 151 \BIBCODE{1993IAUCB..21..151B} \DOI{10.1007/978-1-4612-2742-7}

Cox, A.N., ed, 2000, Allen's Astrophysical Quantities, 4th ed., (Springer, New York, NY) \DOI{10.1007/978-1-4612-1186-0}

Flower, P.J., 1996, {\it ApJ}, {\bf 469}, 355 \DOI{10.1086/177785}

Ghedini, S., 1982, Software for Photometric Astronomy (Willmann-Bell, Inc.) \BIBCODE{1982spa..book.....G}

Henden, A., 2001, {\tt http://www.tass-survey.org/tass/catalogs/tycho.old.html}

Hog, E., et al., 2000, {\it A\&A}, {\bf 355}, L27

Irwin, J.B., 1952, {\it ApJ}, {\bf 116}, 211 \DOI{10.1086/145604}

Irwin, J.B., 1959, {\it AJ}, {\bf 64}, 149 \DOI{10.1086/107913}

Kallrath, J., Milone, E.F., 1998, Eclipsing Binary Stars--Modeling and Analysis (Springer-Verlag). \DOI{10.1007/978-1-4419-0699-1}

Kurucz, R.L., 1979, {\it ApJS}, {\bf 40}, 1 \DOI{10.1086/190589}

Kwee, K.K. and Woerden, H., 1956, {\it BAN}, {\bf 12}, 327

Nelson, R.H., 2010, ``Spectroscopy for Eclipsing Binary Analysis'' in The Alt-Az Initiative, Telescope Mirror \& Instrument Developments (Collins Foundation Press, Santa Margarita, CA), R.M. Genet, J.M. Johnson and V. Wallen (eds) 

Nelson, R.H., 2013, Software by Bob Nelson, \\ {\tt https://www.variablestarssouth.org/bob-nelson/}

Nelson, R.H., 2014, Spreadsheets, by Bob Nelson, \\ {\tt https://www.variablestarssouth.org/bob-nelson/}

Nelson, R.H., \c{S}enavci, H.V., Ba\c{s}t\"{u}rk, \"{O}, and Bahar, E., 2014, {\it NewA}, {\bf 29}, 57 \DOI{ 10.1016/j.newast.2013.11.006}

Nelson, R.H., Terrell, D., Milone, E.F., 2014, {\it NewAR}, {\bf 59}, 1 \DOI{10.1016/j.newar.2014.01.001}

Nelson, R.H., Terrell, D., Milone, E.F., 2015, {\it NewAR}, {\bf 69}, 1 \DOI{10.1016/j.newar.2015.07.001}

Nelson, R.H., Terrell, D., Milone, E.F., 2016, {\it NewAR}, {\bf 70}, 1 \DOI{10.1016/j.newar.2016.02.001}

Nelson, R.H., 2016, Bob Nelson's O--C Files, \\ {\tt http://www.aavso.org/bob-nelsons-o-c-files}

O--C Gateway, Paschke, A., {\tt http://var2.astro.cz/ocgate/}

Rucinski, S. M., 2004, {\it IAUS}, {\bf 215}, 17

Terrell, D., 1994, Van Hamme Limb Darkening Tables, vers. 1.1.

Van Hamme, W., 1993, {\it AJ}, {\bf 106}, 2096 \DOI{10.1086/116788}

Whitney, B.S., 1959, {\it AJ}, {\bf 64}, 258 \DOI{10.1086/107930}

Wilson, R.E., and Devinney, E.J., 1971, {\it ApJ}, {\bf 166}, 605 \DOI{10.1086/150986}

Wilson, R.E., 1990, {\it ApJ}, {\bf 356}, 613 \DOI{10.1086/168867}

Wilson, R.E., 1998, Documentation of Eclipsing Binary Computer Model (available from the author)

\endreferences 

\end{document}
