Machine solutions of light curves of eclipsing binary systems

Machine Solutions of Light Curves of Eclipsing Binary Systems I. J m ~ x r w c ~ Space Sciences Laboratory, General Electric Company and Adjust as indicated Flower and Cook Observatory, University of Pennsylvania Su~M~mY This paper considers various approaches to the solution of light curves by general-purpose digital computers. Two classes of solutions are defined. The direct methods attempt to arrive at all pertinent eclipse parameters simultaneously. The indirect methods, on the other hand, strive to isolate a specific parameter in terms of which the remaining elements of the eclipse can be determined. Advantages and weaknesses of various methods within these classes are discussed. The major part of the paper is devoted to the discussion of an example of implementing an actual computation. For this purpose the solution of a complete occultation eclipse in terms of the second method of Kopal is chosen. A detailed discussion of the computer programs and operator procedures is included. Typical computations are presented for the YZ Cassiopeiae system. The computer language employed is basic I~ORTRAlff. In conclusion it is argued that attempts to analyze complex binary systems by indirect methods hold little promise. The writer feels that such systems should be studied by direct methods implemented in terms of computer simulation techniques.

1. INTRODUCTION The advantages of applying general-purpose electronic digital machines to computations in physical sciences have been well established. Their further elaboration is hardly necessary. Suffice it to say t h a t these devices enabled investigators in various disciplines to work on problems unthinkable in terms of conventional computing techniques either because of the analytical complexity of theory involved or because of the amount of data required or employed in the computation. I n the case of eclipsing binaries both of these factors are present simultaneously. The basic computational aspects of the solution of eclipsing binary light curves have been established b y Russell and Shapley (1912) very early in their research on eclipsing binaries. Their ideas, undoubtedly influenced b y the computational means available at t h a t time, have set the tone for a greater p a r t of the work t h a t followed. To make the solution of light curves of eclipsing binaries tractable Russell and most subsequent investigators worked in the direction of isolating a specific parameter in terms of which elements of the eclipse could be determined. I t was eventually recognized t h a t this approach contains certain fundamental difficulties. These were discussed at length b y Kopal (a series of papers starting in 1941 and culminating in the monograph of 1959) and Tsessevich (1947). Their arguments indicated t h a t the only proper, although not necessarily the only useful, line of attack consisted in a t t e m p t s to determine all pertinent parameters simultaneously. Since such an approach was found 63

64

Machine Solutions of Light Curves of Eclipsing Binary Systems

to lead to iterative minimization procedures, the required computations were too extensive to become widely used as long as one was limited to manual computational aids. These conditions prevailed even for the simplest physical models of the binary systems observed at a very modest number of points. As the physical model becomes more complex, the conventional computing devices become increasingly inadequate. The value of automatic computers for such systems cannot be overemphasized. Another circumstance in which one seeks a recourse to such machines is connected with large amounts of observations usually available to the investigator. Since the size of the data sample presented an annoying problem for all conceivable approaches, it became necessary to form the normal points, "judiciously" censor the data, interpolate or extrapolate the curves over regions inadequately observed, etc., in order to obtain some information. All such operations modify, to a certain degree, the original data sample. Thus, the analyst begins to work not with the original, but with " m a n u f a c t u r e d " data. I n the process, certain information originally contained in the data m a y be obliterated or false information created which was never a p a r t of the original sample. Such difficulties are connected with one of the basic, and not easily answered, questions of information theory, namely the question of where in the data the information lies (e.g. Blackman and Tuckey, 1959). Yet, until the advent of electronic computers, there was no effective way to work with the large mass of original observations. Thus, if the binary system under study is both physically complex and extensively observed, the use of automatic machines becomes mandatory. Another typical problem in the solution of light curves which does not lend itself to an effective t r e a t m e n t b y manual computations occurs in the initial stages of a solution when parameters such as the depths of minima and the limb-darkening coefficients are either not known at all or known inadequately. The study of the response of the solution to changes in these parameters requires, short of a full-scale differential correction computation, a parametric study in which a repetitive re-computation of the light curve is undertaken with new initial conditions. For obvious reasons such numerical experimentation has not been widely used. However, one m a y expect that, just as in other areas, this approach m a y yield useful and, very often, unexpected results. For studies of the types just enumerated, an electronic digital machine is an indispensable tool.

2. APPROACHES TO THE AUTOMATIC SOLUTION OF LIGHT CURVES (a) General Remarks The analysis of light curves can follow numerous paths, each of which leads to specific computer requirements. However, before discussing these, consider the type of machine most suitable for the work contemplated here. The solution of light curves b y existing techniques involves various stages at which the investigator m a y desire to exercise his judgment before proceeding further. This is particularly i m p o r t a n t in studying light curves not previously analyzed or those light curves for which difficulties can be expected, but the results cannot be anticipated. I n such exploratory work a complete mechanization of the solution of the observed light curve is not too efficient a procedure. I t is more sensible to conduct a preliminary analysis on either a small machine or a remote station linked to a large machine operating in the time shared mode. I n both cases the computation is under the investigator's control in the sense t h a t the course of the computation can be

I. JURKEVICH

65

altered immediately depending on the nature of intermediate results observed b y the investigator. For problems of the type considered here a full.scale computation on a large machine is justified only when such a preliminary survey has been made. Depending on the computing installation available, the investigator's ability to conduct studies of a certain scope will depend on the compromise he may have to strike between factors such as the storage capacity, speed, cost, convenience, and flexibility of operation. These factors will also influence the degree of refinement necessary in his computer programs. Clearly, if no attention need be paid to such factors, his computer procedures may be unrefined. If, on the other hand, he cannot ignore them, his programs will have to be produced with great care. These considerations become more meaningful by discussing the computer requirements needed to implement the solution of a light curve in terms of specific methods, most of which can be assigned to one of the two fundamental groups which are defined subsequently. (b) Definition oI Classes o/Methods In any model of a binary system, however complex, the fractional loss of light during the eclipse phases of the orbit is functionally related to the appropriate parameters y~ and the time t by ~ ~-- ~(t, Yi). Each observed point on the light curve yields an equation of this type. The simultaneous solution of the resulting, usually highly overdetermined, system of equations for the parameters Yi will be referred to as a direct solution. There exists a number of ways in which the mechanics of the solution can be effected. However, regardless of the method, all direct approaches resort to a search for a set of parameters which, according to some criterion, best represent the observed light curve. In principle, it is immaterial whether the functions ~ are known in terms of explicit analytical equations or only numerically. For example, the analytical relations are available for the rectifiable Russell model and for the tidally and rotationally distorted model. For more complex models, such as those mentioned by Wood (1962), Grygar (1963), and Horhk (1966), the loss of fight would have to be evaluated numerically. Solutions based on the above approach constitute the group of direct methods. I t was pointed out b y Tsesevich (1947) and Kopal (1959) that the direct approaches permit a straightforward determination of the probable errors of the system parameters and do not require preliminary smoothing of observations. However, the practical application of the direct approaches usually leads to enormous computations even for the simplest physical models. This is particularly so when the solution must be carried out ab initio. The problem becomes more tractable only if a preliminary set of parameters is available, since then either a differential correction procedure or a parametric computation can be used. Thus, at the present time the direct determination of the preliminary values of the desired parameters is ineffective and ought to be resorted to only when all other approacheshave failed. For these reasons practically all methods extant are designed to produce, in an indirect way, a solution for simplified models of the binary systems. For example, for the Russell model the fractional loss of light is given by c¢ ~ ~(/c, p), where ]c and p have their usual meaning. The indirect procedure consists of inverting this relation to yield the function p ----p(¢¢, k) and constructing special functions which permit a relatively straightforward isolation of the parameter k, in terms of which the remaining parameters are obtained. Methods such as this comprise the group of indirect methods. Note that the above definitions differ from those of Kopal and are more in keeping with those of Tsessevich. Consider in greater detail the relative merits of direct and indirect approaches.

66

Machine Solutions of Light Curves of Eclipsing Binary Systems

Indirect methods The largest class of methods in existence is concerned with the problem of deducing parameters of the eclipsing system from its observed light curve in terms of the rectifiable Russell model (Russell and Merrill, 1952). The mechanization of such methods can be implemented in a number of ways. Historically, all methods within this class have become practical only because numerical values of functions appropriate to the given method were tabulated. Therefore, storing such tables in a computer's memory or in the peripheral mass storage area (e.g. magnetic tapes, disks or drums) offers the most direct way to the mechanization of the computation. The efforts of Wellmann {1962) and of Huffer and Collins (1963) are examples of the actual implementation of such solutions. There are several unattractive features in this approach. The storage of tables requires either a large memory or external storage units. This generally implies a computing installation of considerable size. Of greater importance is the fact that the use of tabular data implies the use of two-way interpolation to obtain the values of eclipse functions for the intermediate values of the ratio of radii k and the geometrical depth p. These are the two parameters usually employed in the tabulation of eclipse functions. The complexity of the interpolation depends on the tabulation interval and the precision desired. The smaller the interval or the precision required, the simpler the method. Conversely, for a reasonably high precision simple interpolation schemes require large storage capacity since the numerical values of the appropriate functions must be available at frequent intervals. Clearly, the use of any interpolation formula leads to limitations on the achievable precision. If, in addition, limb darkening is to be considered an adjustable parameter, one becomes concerned with three-way interpolation. With the presently available tables, the increment in limb darkening is generally too large to achieve the precision ultimately desired. The latter difficulty can, of course, be circumvented in some methods by working only with the functions for extreme values of limb darkening. The well-known interpolation formulae are then used to obtain functions for arbitrary values of the limb-darkening coefficient. A modification of the above scheme, which is in principle very attractive, is based on approximating the tables in question by an interpolation polynomial or a surface. Thus, for a given value of limb darkening x one would merely have to store the coefficients of the respective interpolation expressions. Since most functions to be approximated v a r y smoothly the number of such coefficients, commensurate with the required precision, is expected to be relatively small. Clearly, if an adequate representation could be found, a significant gain in storage area as well as in the computing speed would be realized. The difficulty with this approach is hidden in the selection of a suitable interpolation expression uniformly valid throughout the range of the appropriate independent variables. The author's experience with this problem was negative. For instance, for the fixed values of k, the method of least squares was used to fit various interpolation polynomials to Merrfll's p(c¢, k) functions. With one hundred tabular values and all computations carried out to eight significant figures, the best representation b y a polynomial without gaps (Hastings, 1955) required about five terms yielding residuals of the order of ~=0.005 throughout the range of k. A somewhat better representation was achieved by means of rational polynomial expressions. A cubic over a cubic was found to decrease the residuals b y almost a factor of 5. However, an attempt to extend the fit of/~(¢¢, k) in the k dimension resulted in a representation to no better than several units in the third decimal place. A significantly improved representation was achieved by subtracting from the tabular values of p(k, ~) a linear trend and then fitting a polynomial or a Fourier series to the function [p(k, ~)-linear trend]. With a five-term polynomial the residuals in the k direction

I. JU~K~VICH

67

were reduced to several units in the fourth decimal place. The representation in the direction was somewhat inferior. For a comparable precision the Fourier series representation required too many terms to be of practical value. The difficulty with the choice of suitable expressions could, perhaps, be alleviated if a reasonably simple analytical approximation to a(k, p) or T(~, k) could be found to serve as a guide to proper interpolating equations. However, the complexity of the defining expressions for ~(k, p) is well known, rendering attempts at series expansions and subsequent reversion of the series impractical. In addition to these analytical difficulties, the following was discovered. Many leastsquare approximations computed in this investigation indicated that the representation invariably deteriorated when polynomials of orders greater than five or six were used. Since in this instance one is not dealing with observational data, such behavior is not reasonable. An analysis which falls outside the main theme of this work indicated that for a polynomial of order n with no gaps, and an evenly spaced argument over the range of 0-1, the matrix of normal equations appears to be a multiple of a finite segment of the well-known Hflbert matrix which is extremely fll conditioned. Consequently, with eight digits employed in the computation the inverse of a matrix of the given rank may have a very low accuracy. I t appears that this critical loss of accuracy reveals itself in approximations employing the fifth- or sixth-order polynomials. An obvious way out of this difficulty is to use a higher machine precision, or, even better, an artifice of subtracting .known trends from the function in question and fitting polynomials to the remainders. Although this is an interesting problem in numerical analysis, it is only indirectly connected with the solution of the light curves and it was not pursued further. B y no means were all the possibilities in this area exhausted in the study mentioned. For instance, it was pointed out to this writer b y Dr. R. E. Wilson (1966) t h a t he succeeded in representing p(~, k) functions by polynomials in half integral powers of ~. The precision he achieved is approximately four decimal places, but an increase in the number of terms beyond a certain power of ~ led to no significant improvement. In the opinion of the present writer, the above approach in its present state of development is of a limited value in automatic computations where an accuracy greater than four decimal places is achievable by other methods. The difficulties characteristic of the foregoing approaches can be avoided by working with the defining equations of the eclipse functions appropriate to the given method of solution. In principle, this scheme will provide an unlimited precision in eclipse functions b y dispensing with the interpolation problems with respect to the ratio of radii k, the geometrical depth p and the coefficient of limb darkening x. In its demands for the machine storage capacity this method is intermediate between t h e methods discussed above. The required capacity is considerably larger than that which would be required by a set of constants describing approximating polynomials, but it is very much smaller than that required for the storage of tables. Examples of such approaches can be found in the works of Jurkevich (1964 and this paper) and Harris (1966). In summary, all indirect methods produce a solution by constructing special functions which permit one to isolate a critical parameter in terms of which the remaining parameters are determined. I t is only then that the parameters are used in a direct computation of a theoretical light curve which is compared with observations. If the representation is adequate, one can take a further direct step and differentially correct the preliminary parameters. Direct methods. The indirect methods of computing the eclipse parameters exist only for those models of binary system for which analytical relations between the light curve and the system parameters can be produced.

68

Machine Solutions of Light Curves of Eclipsing Binary Systems

I t is evident t h a t with the increasing complexity of the physical model such relations become so involved t h a t a n y a t t e m p t at an inverse solution becomes, for all practical purposes, impossible. Under such conditions, the direct methods seem to offer the only way of obtaining a solution. This hope is based on an observation t h a t for the known system parameters one can more readily compute the light curve even for a very complex physical model, whereas prospects for solving the inverse problem m a y be dim indeed. Thus, the direct methods attempt, within the framework of some theory, to establish from the beginning a set of parameters which satisfy either the observed light curve or some of its special properties. Clearly, this approach is equivalent to a parametric study which, as a rule, involves an immense amount of numerical experimentation and, as such, would be impossible without a large-capacity, high-speed, automatic machine. Some work along these lines has been reported b y Wood (1962) and Mauder (1962). For example, Wood (1962) assumes a number of eclipse parameters as well as certain physical properties of component stars, and from these he computes numerically the amount of light visible at any point in the orbit. Although he did not propose to do so, in principle a catalogue of light curves could be produced by parametric variation of all eclipse elements. I n its concept and application, such a catalogue would be equivalent to King's (1920) set of standard velocity curves for spectroscopic binaries. However, owing to the number of parameters necessary to describe the light curve of even a simple model, this would be an immense undertaking. Furthermore~ it is not clear t h a t the representation of light curves by a set of elements is unique. For these reasons this approach m a y be ineffective in obtaining a set of preliminary elements. However, within Wood's scheme of things, such elements are absolutely necessary, since he proposes to solve the observed light curves b y differentially correcting the appropriate preliminary elements. A partial way out of these potential difficulties is to utilize, whenever possible, a set of preliminary elements obtained b y one of the indirect methods before proceeding with the analysis of more complex physical models. The model itself can be extremely complex since within the present scheme the computational simplicity is of secondary consideration. Another approach along these lines has been followed by Mauder (1962). Starting with Kopal's (1961) suggestion t h a t a Fourier transform of the light curve would be a natural means of separating eclipse and proximity effects, Mauder proposes to compute as general a theoretical light curve as practical and subject it to a Fourier series analysis to produce its discrete amplitude spectrum. This is taken as the next best thing to the theoretical Fourier spectrum of the light curve. Once a catalogue of such spectra has been assembled, it can be used for comparison with the spectra of observed light curves. Since so far no results have been reported by Mauder, it remains to be seen how practical such an approach m a y be. One can, however, expect the occurrence of the same difficulties as in the method of Wood. I n its most general formulation Mauder's method would impose on the computing facility the requirements quite similar to those of Wood's scheme. Other, less ambitious, but immediately more practical examples of the indirect approaches are contained in papers of Grygar (1963) and ttor£k (1966). The latter author applies the direct method to a model consisting of an ellipsoidal star and a spherical star whose limb-darkening law is non-linear. Clearly, it is impractical to summarize in one article the computer implementation of all conceivable methods. At the present time insufficient experience exists with the direct methods. Consequently, it is too early to discuss intelligently their application. As far as indirect methods are concerned, it is useful to discuss only those methods which are naturally suitable for automatic treatment. Despite these restrictions, a complete discussion

I. JURX~.vIc~

69

of even the simplest model would require a prohibitive amount of space. For these as well as pedagogic reasons, it is proposed to devote the following pages to the discussion of the mechanization of the solution of a light curve by one of the existing methods. For this purpose we shall consider a total occultation eclipse arising in the rectifiable Russell model of a binary system. The computer implementation will be in terms of the basic FoRTRA~ language in order to exhibit the minute details of the problem. First, however, we shall explore the suitability of existing methods of solving light curves for automatic computation.

(c) Suitability ol the Existing Methods/or Automatic Computation Among the numerous indirect methods of deducing from the observed light curve the radii of the components of a binary system and the inclination of the orbital plane one can mention the following ones: (a) Method of Russell (1912), (b) Method of Sharbe (1924), (c) Method of Fetlaar (1923), (d) Method of K r a t (1934, 1935, 1936), (e) Method of Piotrowski (1937), (f) Method of Schneller (1949), (g) The first method of Kopal (1941), (h) The second method of Kopal (1946), (i) Method of Kitamura (1965) and various modifications of the above methods. In Section 2 (b) certain computer requirements for the mechanization of such methods have already been discussed. However, before making the final selection of a specific method, it is appropriate to discuss several fundamental concepts which will affect our decision. As mentioned earlier, beginning with Russell the main idea behind all older methods was to make a rather difficult computational problem tractable numerically. The intent of most methods was an isolation of k, the ratio of the radii. For instance, to achieve this, Russell adopted certain pivotal points which were used to construct special functions of k and ~ which enabled him to determine k in a straightforward way. Kopal discussed at length the underlying weakness of this and similar methods. Briefly, the principal points of Kopal's criticism are that the method of Russell assigns infinite weights to the pivotal points and that it does not employ the observed data. Some investigators subsequent to Russell (for instance, Fetlaar and Piotrowsky) must have recognized these difficulties at least intuitively, for they strove in various ways to minimize the above restrictions. Success of such attempts was, however, always limited as long as the method itself was, in principle, an extension of Russell's approach. Aside from the method of differential corrections discussed at length by Tsesevich (1947), the break with tradition came with the two methods of Kopal (1941, 1946), who abandoned all attempts at isolating the value of k. Although his methods employ the function p --~ p(~, k), which results from the inversion of co(k, p), a procedure characteristic of the indirect approaches, they are not, strictly speaking, of the indirect type since all pertinent parameters are determined simultaneously. Furthermore, Kopal succeeded in removing the restrictions which arise from the use of pivotal points on the :light curve. Unlike all previous methods, Kopal's approach required no special functions beyond c~

70

Machine Solutions of Light Curves of Eclipsing Binary Systems

and p. Unfortunately, the methods of Kopal had the common characteristic of all direct methods, namely, they were cumbersome computationally and thus have been little used as long as manual computations were involved. As already indicated, the computational difficulties are associated with the iterative least squares procedures employed. Since such calculations are readily mechanized, the two methods of Kopal, in addition to their conceptual merits, are superior to all others when high-speed digital machines are available. The choice between the first and the second method is a m a t t e r of convenience and, what is more important, of statistical considerations. The second method, discussed in great detail b y Kopal (1959), has the very attractive feature of being self starting. This means t h a t once the computation has been started with a suitably chosen initial value of k, the appropriate equations yield an improved value of k automatically as long as the process is a convergent one. The difficulty with the second method is t h a t it violates the Gauss-Markoff theorem on least squares. I t must be pointed out t h a t the least squares adjustment of data can be used to estimate the best numerical value of a quantity even though errors are not necessarily the observational ones. Generally, the Gauss-Markoff theorem is concerned with linear estimation of parameters appearing in linear equations. Assumptions which must be satisfied are : (a) Estimators of parameters of interest are unbiased linear combinations of the observed values of the drawn sample. (b) The " b e s t " unbiased linear estimator is t h a t one which minimizes the variance of the statistical variables (usually the observed quantity). The distribution law of the observed quantity is not restricted to a n y particular form. I f the expected value of a variable is E ( y i ) ~ a i l x l ~ ai~x~ -~ ... ~ atmXm

then the variance S is given by N

S ~- ~

(Yi - - a i l x l - - ai~x2 . . . . .

aimXm) 2 W i ,

i=1

where Wi are known as weights. Minimization of the variance S results in the famihar set of equations known as the normal equations. A very important fact to note is t h a t for the Gauss-Markoff theorem to hold, the coefficients a q must have known numerical values. I n fact, they m u s t be error free. I n practice, these coefficients are known from observations and thus contain the errors of observation. Therefore, the Gauss-Markoff theorem is not applicable to this situation. I t m u s t be pointed out t h a t the range of its application can be stretched if the coefficients are known with such accuracy t h a t they affect the expectation of the variable y to a smaller degree t h a n the standard deviation of a n y individually measured parameter. This situation exists in the first method of Kopal where the observed quantity (~ is related to time (or phase) b y the parameters to be estimated. Here, the precision of time measurement is several orders of magnitude better than the light intensity measurement, which yields the quantity ~. I n the second method of Kopal, time is related to a combination of the individual radii and the inclination b y means of coefficients which are functions of the geometrical depth. The latter quantity is not only unknown in the beginning, but also contains implicitly all the errors of observation. Furthermore, there is no guarantee t h a t the error

I. JuRKEVICII

71

distributions in the inclination i a n d the two radii ra a n d r b are the same. The current m a t h e m a t i c a l literature contains little discussion of practical value for these two situations. The violation of the Gauss-Markoff theorem b y the second m e t h o d of K o p a l has the i m p o r t a n t implication t h a t the best value of/c does n o t necessarily correspond to a minim u m value of the variance. Furthermore, the estimates of the desired parameters ~ l l be biased a n d the derived s t a n d a r d deviations m a y be questionable indicators of precision. Despite this serious drawback, the second m e t h o d of K o p a l is convenient in exploratory work because c o m p u t a t i o n a l l y it exhibits an excellent convergence rate even for poor first estimates of k. Also, /c is the ratio of certain coefficients occurring in the equations of this method, thus permitting one to work almost directly with k. These features are n o t available in the first method. I n subsequent discussion, the second m e t h o d of K o p a l will be used for the illustrative problem despite the fact t h a t the first m e t h o d is sounder f r o m the statistical point of view.

(d)

The Second Method o/Kopal

Summary o/ the equations required /or complete eclipses. Most equations necessary to solve a light curve b y the second m e t h o d of K o p a l can be f o u n d in K o p a l ' s recent monog r a p h (1959). Since the present discussion is concerned with a b i n a r y system consisting of two stars having similar elliptieity, the discussion will start with the so-called dynamical condition written for distorted systems. Thus, sin S 0 sin S i -k cos ~ i ---- r~(1 -- z cos 2 0) (1 -k

kP)2,

(1)

where 0 i rb ra

----phase angle, _--inclination of the orbital plane, ---- the radius of the larger component, ~- the radius of the smaller component, z --~ e2 sin 2 i, k = ra/rb,

p = geometrical depth at a given instant of eclipse. L e t 0" be the phase angle at the m o m e n t of internal tangency. T h e n at the latter instant Eq. (1) becomes sin S 0" sin S i -~ cos 2 i = r2b(1 - - Z COS~ 0") (1 --/C) ~.

(2)

This condition exists for complete eclipses of b o t h the occultation a n d transit types. Subtracting (2) from (1) a n d rearranging one obtains [(1 - - cos ~ 0) p 2 -- (1 - - z cos 2 0")] r~ cosec 2 i -k 2[(1 - - z cos 2 0) p -k (1 -- z cos ~ 0")]

x rarbcosec 2 i -k [(1 -- z cos ~ 0) -- (1 - - z cos 2 0")] r~ cosec ~ i -~ sin ~ 0" : sin S 0.

(3)

The second and the third terms on the left-hand side can be written as

zr~ cosec ~ i sin S 0 -~ (1 - - zr~ cosec ~ i) sin S 0". If, using K o p a l ' s notation, one defines E= E"

~---

1 - - z c o s ~0, 1 --ZCOS ~0'',

(4)

72

Machine Solutions of Light Curves of Eclipsing Binary Systems

t h e n t h e use of (4) in (3) yields 2 cosec s i

sin s O = ( E p s - - E " )

ra

-~- 2 ( E p ~-

Zr2b cosec s i

1 --

rar b cosec s i

E")

~ sin s

0"

,

(5)

1 - - zr~ cosec s i

where t h e factors c o n t a i n i n g t h e c o m b i n a t i o n s of r , , rb a n d i can for convenience be w r i t t e n as r~ cosec 4 i DI=

2 2 2 r a cosec 2 i - - Zrar b cosec 4 i

Ds _

r3rb cosec 4 i 2 cosecS i - - zrar 2 2b cosec 4 i ra

"

Also let sin s 0 " .

D 3 ~

W i t h this n o t a t i o n (5) becomes sin S 0 = ( E p s - - E " ) D x + 2 ( E p -4- E " ) D 2 -t- D s .

(6)

T h e i m p o r t a n t f a c t to n o t e is t h a t D 1 / D s = ra/rb = k. F o r s p h e r i c a l s t a r s K o p a l defines quantities 2 C 1 ~ r acOsec si, C s = rarb cosec s i ,

C3 = sin s 0 " .

(7)

T h e use of t h e s e q u a n t i t i e s c o n t i n u e s to be c o n v e n i e n t in t h e case of d i s t o r t e d c o m p o n e n t s . I t m a y be n o t e d t h a t D1 - - - - - - - - - 5 '

C 1 -- zC 2

D2 - D 3 ~

CIC2 - - - - - 2 '

C 1 -- zC 2 C s

f r o m w h i c h one o b t a i n s C 1 --

D2

C2 - -

D 1 + zD~ '

DID2 D 1 + zD 2

E q u a t i o n (6) will e x i s t for e v e r y o b s e r v e d p o i n t . Now, a c c o r d i n g to P i o t r o w s k i (1947) e v e r y p o i n t will h a v e a n i n t r i n s i c w e i g h t of t h e f o r m x/w=

1 --

;t 0o~

2D2~

Op

if t h e errors of o b s e r v a t i o n are c o n s t a n t on t h e i n t e n s i t y scale. I n t h e a b o v e 2 -~ f r a c t i o n a l l u m i n o s i t y a t t h e i n s t a n t of i n t e r n a l t a n g e n c y . F o r d i s t o r t e d s y s t e m s (~ m u s t be r e p l a c e d b y (1 - - z cos 2 0") (~ (e.g. B i n n e n d i j k , 1960). F u r t h e r m o r e , t h e discussion will be r e s t r i c t e d t o a single m i n i m u m a t a t i m e so t h a t t h e f a c t o r s (1 - - ~) a n d D 2 can be o m i t t e d since t h e y r e m a i n c o n s t a n t for t h a t p a r t i c u l a r m i n i m u m . Thus, t h e p r o p e r w e i g h t a s s o c i a t e d w i t h (6) is t h e n g i v e n b y

~

1 W ~

0~

--

2(1 - - z cos s 0") (1 ~- kp) ap

I. JUnKY.ricH

73

~/w sinZ 0 = (Ep' -- E") D,~/w + 2(Ep + E") D,x/w + Ds~/w.

(8)

Consequently Eq. (6) becomes

Equation (8) is fundamental to the solution of a light curve for complete eclipses. Before proceeding with the further discussion one may note from the second equation of (7) that ra= rb cosec 2 i "

Substitution in the first equation of (7) yields cl C 1 ~---

r2bcosec z i

from which one obtains =

(9)

C~ cosec' i

Also, from the first equations of (7)

C1

2_ ra

(I0)

cosec' ~ '

and from (2) cosec' i =

(Ca -- 1) C~ C1)' C22"

rbE (C, 2

rt

__

__

(11)

Employing (11) in (9) one obtains

(12)

=

(C, -- C,)*

E"

+

(I -- Cs) C,

Defining with Kopal (C, - - C1)' E " + (1 - - Cs) Cx = G,

it follows that C,

(13)

r b ~-

From (9) it immediately follows that cosec' i --

G C1

or

t

Finally, the use of (14) in (10) gives ra= 6

B-VI& VoL 12

C1

(14)

(15)

74

Machine Solutions of Light Curves of Eclipsing Binary Systems

A m o n g other quantities which m a y be of interest are the phase angles associated with the instants of internal and external tangencies. At the time of external t a n g e n c y p : 1 a n d Eq. (1) yields sin s 0e sin 2 i ~- cos 2 i ~ r~(1 -- z cos s 0e) (1 ~- k) 2. Solving this equation for sin 20e one has sin20e ~_ r2(1 + k) (1 -- z) -- cos s i sin s i -- r~(1 + k) 2 z

(16)

A t the i n s t a n t of internal t a n g e n c y p = --1 and a procedure similar to the one given above yields sin 20i = r2(1 -- k)2 ( 1 - - z ) - - c o s sin s i -- r~(1 -- k) 2 z

2 i

(17)

However, this q u a n t i t y is more directly obtainable from the coefficient D a. A t the i n s t a n t of mid-eclipse, 0 ~ 0 ° or 180 °, so t h a t cos 2 i ~-- r2(1 - - z) (1 ~- ]Cpo)S, from which one has p0 -

cos i

1

krh~/(1 -- z)

k

(18)

Once r,, rb and i have been determined, the dynamical condition (1) can be used to obtain the geometrical d e p t h at a n y point within the eclipse. Thus, 1 ~(1--sinSicos20~_

Finally, h a v i n g determined a value brightnesses from J, . . . J~

where

1

of k one can c o m p u t e the ratio of the m e a n surface L, 1 . . Lb k s

1--k,

,

(19)

~,]c2

L a ~ luminosity of the smaller component, Lb ~ luminosity of the larger component, 2a --~ fractional luminosity of the system during the total occultation eclipse.

Computational procedure /or complete eclipses. The solution of the light curve for a complete eclipse will commence with the c o m p u t a t i o n of the observed values of the fractional loss of light ~ given b y 1--1 ~(k,p)

1 --

where l --~ rectified intensity at a n y point, ~ fractional luminous intensity at the i n s t a n t of internal t a n g e n c y of a particular m i n i m u m being investigated. This q u a n t i t y ~ l l be estimated from the light curve. I t is now necessary to make a decision as to whether the m i n i m u m in question is an occultation or a transit. For complete eclipses when b o t h minima have been observed, this decision is n o t difficult to make from the appearance of the minima themselves.

I. JURKEWCH

75

The initial estimate of k necessary to start the iterative process using Eq. (8) can be obtained in two ways. If both minima are available, one can use the relation between the fractional intensities at the instants of internal tangency, namely kS

--

1

--

,~aY(k,

~b

(20)

--1)

where ),~ = fractional intensity at the moment of internal tangeney corresponding to the occultation eclipse, 2b ~ fractional intensity at the moment of internal tangency corresponding to the transit eclipse, and the function Y(k, --1) is defined by Kopal (1959) as r(k,

-1)

-

3(1 -- Xb)

+

3 -- xb

2xb

1(v)/¢

3 -- Xb

k2

-

-

--,

where xb -~ limb darkening of the larger component of the system. Equation (20) is solved by iteration and the program to do this has been discussed b y Jurkevich (1964) as part of the rectification computation. I n general, if the type of eclipse is not known for the minimum in question, Eq. (20) yields two solutions, obtained by setting the fractional luminosity at the internal tangency point successively equal to ~a and ~b. If only a single minimum is available, Eq. (1), written for the instants of internal and external tangency, yields sin 20isin ~ i + c o s 2i _ - - , { 1 - - k , ] z . sin z 02 sin ~ i + cos 2 i

\ 1+ k]

Solving for k one obtains k ----x/(sinZ 0~ -~ eotan 2 i) -- x/(sin~ 0i -- eotan ~ i) x/(sin ~ 0e ~- cotan ~ i) -~ ~/(sin ~ 0i -- eotan ~ i) If i is taken nearly equal to 90 ° and 0e and 0i are reasonably small, k - - Oe - - Oi

(21)

0 e -'~ 0 i

Equation (20) is generally referred to as the depth relation. This relation is effeetivin indicating the type of eclipse only if the two values of k obtained for 2a or ~b are signifie cantly different. As pointed out by Kopal, it is recommended to start the analysis by assuming that the primary minimum is due to a total eclipse if evidence to the contrary is not available. This assumption always leads to a physically realizable solution. I t is necessary to note that if the initial value of k is obtained from Eq. (20), the limb darkening of the larger star Xb has already been fixed. If, however, k is obtained from Eq. (21), the value of Xb as well as the type of eclipse remains open. Once a starting value of ]¢ has been determined, an initial value of p is assumed and theoretical values ¢¢(k, p) are computed until, for a given observation, theoretical and observed values of ¢¢(k, p) differ by less than a given amount. The value of p at such time is used to compute the coefficients in Eq. (8). This procedure is repeated for all observations. The resulting set of normal equations is solved yielding a value of k----D1/D ~. This ]¢ is taken as a new initial value and the entire procedure is repeated until successive values of/c change by less than a prescribed amount. 6*

76

Machine Solutions of Light Curves of Eclipsing Binary Systems

If the difference between the assumed k and the result of the second iteration fails to be substantially smaller than after the first iteration, the first inference is that one deals with the wrong type of eclipse. If both minima have been observed with sufficient precision, they can be subjected to analysis independently. Under such conditions one can at least in principle observe the following : (a) Tf ]~primary :~ kaecondary (within a reasonable tolerance), as obtained from the iteration process described above, and if ]¢ obtained from Eq. (20) is intermediate between the two, the indication is that limb darkening of both components was estimated incorrectly. (b) If k resulting from Eq. (20) agrees substantially with ~primary, but not with ksecondary, limb darkening of the smaller component xa is incorrect. (c) If ksecondary~ k obtained from Eq. (20), it follows that xa has been underestimated. (d) If kprimary ~ k obtained from Eq. (20), it follows that Xb has been overestimated.

Having obtained the final value of k, Eqs. (13)-(19) are used to obtain the desired parameters of the system. The representation of the light curve by these parameters is obtained by computing the theoretical values of p from r ~

(1-sin2ic°s~0)

1

1 -- z cos 2 0

k

P ----

(22)

and the corresponding values of ~(b, p). The latter quantities are then used in t =

1 -

(1 -

;t)~(k,

p)

(23)

to compute the theoretical values of intensities. It is appropriate to note that the value of inclination i obtained from Eq. (14) is the rectified inclination. The relation between the rectified inclination and the true inclination is given by sini,

~/\

1--z

'

COS i COS ir -- - ~/(1 -

z)

From these I

tani=tanir

~

1+ 1

Z

. -]1/2 cosec~ ~r[ •

(24)

l

The estimate appropriate for a triaxial ellipsoid model is tan i

b

tan

where b ~ c. The quantities b and c represent the smaller principal semiaxes of the ellipsoidal stars. Unfortunately, the present analysis is not capable of yielding the quantity c with any degree of certainty. The above discussion summarizes all the necessary expressions needed to analyze complete eclipses in the first approximation.

77

T. JURKEVICH 3. GENERATION OF NUMERICAL VALUES OF ECLIPSE FUNCTIONS

(a) Fundamental and Derived Functions It is evident from the foregoing discussion that in order to employ the second method of Kopal a source of values of p(a, k) must be available. However, explicit equations for this function in terms of the quantities k and ~ do not exist. For this reason, it will be necessary to derive the numerical values of p(~, k) from the fundamental functions a(k, p) by a suitable inversion technique. The functions ~(k, p) and p(k, ~) are the only functions required by the method of Kopal. All other indirect methods require, in addition to these, special functions constructed from ~(k, p) and p(k, ~). In Merrill's notation (1950) the equations for the fundamental functions, c~(k, p), are summarized below:

Zero limb-darkening o~_ 1 -- ~ [(2% -- sin 2%) + kZ(2~ -- sin 2q~,)], 2ztk ~

1 kpZ + 2p + k 2 l+kp

COS ~ s - -

(25)

sin • = k sin ~,.

Total limb.darkening Occultation

E'(O') V - E'(O') F(O, ~) + F'(O') F(O, ~) -1¢,oc = -1~ .I_T + F(O')U 6~+~/5 U:--258+(3k

2-10)

52+(8k ~-14)

V : 45 s -

av'(O')E(O,q~)],

5-(3k 4-9k 2+6),

(16k ~ -- 28) 5,

1/~

sino=lI(l+5+k)(l+5--k)] 5 sin q~ = p.

(26)

Transit-partial phas~ o~ eclipse lo~t r

-

-

1 1 r ~ -~ F'(O') W + E'(O') X 6 x/(k5) -- E'(O') F(O, of) + F'(O') F(O, qJ) -- F'(O') E(O, q~).,,

[2

W:--2k5

a+(3-10k

z) 5~ + ( 8 k -

14kz) 5 - - ( 3 - - 9 k

z+6k4),

X : 4k5 s - (16k - 28k 8) 5, sinO:

2 1 - - ' [ ' ( k + 5 - - 11) ( k/ - ~z - 5' +kl ) 5] sin~:

5--k.

(27)

78

Machine Solutions of Light Curves of Eclipsing Binary Systems

Transit-annular phases o/eclipse l~ann

__

2

[y(o')r_+y2o')_z

-

+ E'(O') F(O, q~) -- F'(O') F(O, el) + .F'(O') E(O, ~)],

7~lT L 6 X/(1 -- (k -- 6) s)

Y=64-Z:--6

s-k4),

(5 + 2 k s) 6 s + 6 6 - ( 2 - k

4~-2k6 a~- ( 5 - - 8k 2) 6~ - (8k-- 14k a ) ( ~ - ( 4 - - l l k s-~7k4), sin 0 : [ 11 ---- (k l / 2--~ ' ( k(~)2 6) $ ]

(28)

s i n g ~ = [ l1+/ks-' '-l~+]k + In the above expressions

6=1 +kp, IT

2 [3 arc sin x/k -- (3 + 2k -- 8k s) x/(k(1 -- k))], 3z 0' = ~/2 -

0,

and ~/2

F'(0') =

~/(1 -- cos 2 0 sin ~ ~v)

d~v,

0

.12

E'(O') = I x/( 1 -- c°sS 0 sin ~ ~v) d~o, 0 ¢

F(O, of) -~

~/(1 -- sin s 0 sin s ~v) 0

E(O, F) = I ~/(1 -- sin s 0 sin s y~) d~, 0

are the complete (F'(O'), E'(O')) and incomplete (F(O, q~), E(O, ~)) elliptic integrals of the first and second kind respectively. The normalized loss of light for intermediate values of limb darkening is given by the well-known expressions x~oc - - 3(1 -- x) 0~ + 3 --x ~t~

__ z ( 1 - - x ) oo~ + ~--X

_

2x _

1~ °c

3 --x (~ - - 1) x l~tr" g - - X

For the annular phases ~tr is to be replaced by ~ n . in the last equation.

I. JURK~VICH

79

The quantity u is given by 3k~ 3 k ~ _ 21v" Inspection of expressions (25)-(28) reveals that to compute the quantities 10~°c, l~tr and l~ann one must have a source of values of the necessary elliptic integrals. Furthermore, since it is planned to obtain p by a numerical inversion of ~(k, p) a suitable scheme to achieve this objective must be selected. The Newton-Raphson method of solving for the roots of an equation is chosen for reasons of its simplicity and the secondorder rate of convergence. To utilize the full advantage of this method, it is desirable that an analytical expression for Ox~/Op be obtained. Finally, the expressions for ~'s and Oa/Op become indeterminate for certain values of k and p. Since a machine cannot recognize such situations, the appropriate forms must be evaluated beforehand. These three problems are considered next in the order mentioned above.

(b)

Computation o/the elliptic integrals

General remarks. There exist infinite series for the complete elliptic integrals which converge with extreme rapidity throughout most of the range of the modulus. The discussion of typical cases was given by Wilson (1912} and many others. A somewhat different approach was taken by Hastings (1955) who developed very efficient approximations for the complete integrals based on asymptotic expansions. However, when it comes to the incomplete integrals, the situation is less satisfactory. Among many discussions on the subject one can mention Wilson (1912) and SehlSmilch (1879). Such series are found to converge either too slowly or if reasonable convergence rates are achieved, it is accomplished at the expense of computational simplicity. From the standpoint of machine computation it is desirable to search for a method which is iterative in nature. Under such circumstances it appears that the numerical evaluation of elliptic integrals can most profitably be based on Landen's transformation. A detailed discussion of this transformation is given by Hancock (1958). For the purposes of the following discussion it will be sufficient to state a few facts concerning Landen's transformation and write down the necessary equations. There are two forms of Landen's transformation. One form is characterized by the fact that the transformed modulus of the integral increases and the corresponding value of the amplitude decreases. In the notation of Hancock (1958) the transformation is given by (29)

F(k, ~) = F(k,,, %,) 1-] i=O

E(k, ~) = F(k, ~) [1

+ k 1 + - ~ 1 + klk----:-~- ..* + kl ... kn-1

I 2 -- k sin ~ ~- x/k sin ~1 ~- "'" -~

2 n-1

kl .U~n_i 2n

1

x/(kkl.." k2) sin %~-1 -- x/(kkl.." kn-1) sin ~n

(30)

80

Machine Solutions of Light Curves of Eclipsing Binary Systems

where km -- 2 ~/(km-1), 1 +

(31)

kin-1

sin (2q~m - - q~m-1) = k m - 1 sin qm-1, ko = k ,

(32)

q~o-----q~.

The other form decreases the modulus and increases the amplitude. The transformation equations are

F(k, ~) = P(k,~, %,))-I.= E(k,~):F(k,~)

_t_ k [ ~

1--~-

(33)

1-~-~-

--~-+...

sin ~ol % ~/(klk.) . . . sin1 ~2 +2

(34)

where km =

1 - - ~ / ( 1 - - ]¢m 2 --1) 1 + J ( 1 - km-1)

(35)

4( 1

(36)

sin (99m -- ~m-1) = ~/(1

-

-

--

k o -= k,

cfo = ¢p,

2 km-1) sin~m-1 2 kin-1 sins ~,,,-1)

m = 1, 2 . . . . . n .

The present case yields the complete elliptic integrals in the following form : 7~

E(k,:tl2)-----F(k,:tl2)

1--T

(i + kD, 1

2 +-7+...

.

The second form of the transformation is preferable when the modulus k is nearer zero than unity. The exact value of k which makes one of the forms superior to the other is relatively unimportant and is set at 0.8 in this discussion for reasons discussed by gurkevich (1964). In subsequent paragraphs the computer implementation of Eqs. (29)-(36) is considered.

D i s c u s s i o n o/ the elliptic integrals subroutine. In discussing this part of the computer program frequent reference will be made to the program listing (pp. 92-91, statement 500 to statement GO TO 401), and flow diagram shown in Fig. 1. Since this section of the program is designed to serve as a subroutine, certain quantities and parameters which appear in the discussion must be supplied from external sources. These will be mentioned in the appropriate places.

I. JURKEVICH

81

[ , . ~ "~,

FIG. 1. Computation of elliptic integrals. Flow chart. The subroutine is entered through the statement 500 which sends the computation through one of two available paths, depending on the value of CAY. The latter quantity denotes the starting value of the modulus of the integral being computed. If CAY is less t h a n or equal to 0.8, the program selects the proper Landen's transformation which decreases the moduli and increases the amplitudes. If CAY ~ 0.8 the path chosen increases the moduli and decreases the amplitudes. Consider the first of these alternatives. First, it is necessary to store the required values of the modulus CA¥ and the amplitude P H I at every iteration step. These quantities are therefore redefined as dimensioned variables A(I) and FEE(I). This is done starting at statement 1. I n addition, two parameters K and 1 are assigned specific values. The parameter K will be used to return the program from the quadrant testing subroutine. The parameter I will be employed in the iteration procedure to obtain successive values of A(I) and FEE(I). Experience has shown that even in the worst case approximately five iterations yielded the desired precision. However, throughout all programs storage for seven values of A(I) and F E E ( I ) is provided. The initial values of CAY and P H I are obtained from outside the subroutine. Five statements, starting with statement 20, are designed to compute the next value of the modulus A(I -b 1), and the sine and cosine functions of the difference of two s u c -

82

Machine Solutions of Light Curves of Eclipsing Binary Systems

cessive values of the amplitude. The latter quantities are assigned symbols SN and CS. At this point it is important to realize that since the present transformation increases the amplitude, F E E ( I ) may exceed a value of 2zt by a large factor. However, the evaluation of the inverse trigonometric functions by subroutines usually produces the angle whose value does not exceed 2z. Consequently, a scheme is needed which will yield the true value of the angle. Such a scheme is implemented in statements starting at AL ---- 1. and ending with GO TO 95. The idea is to reduce repetitively the value of FEE(I) by 27t until the difference is less than 2zt. When this condition is reached, statement 95 sends the computation to the quadrant testing subroutine which, using the previously assigned value of K, returns the computation to statement 40. The angle returned is designated as GAM and it must be restored to its true value by adding to it a multiple of 2zt by which the original angle FEE(I) has been reduced. This restoration as well as the evaluation of the transformed amplitude is accomplished in statement 40. At this point the transformed value of the modulus A(I Jr 1) is compared to a tolerance El)S, defined outside the subroutine. If A(I ~ 1) exceeds EPS, the program will increment I by 1 and return to statement 20 to obtain additional transformed values of the modulus and amplitude. If, however, A(I-4- 1) has been computed with the desired precision, the program proceeds to evaluate Eqs. (33) and (34). The evaluation starts with the statement 3 and ends with E = F*E -~ CAY*TT. The computation utilizes the loop D O 5 J - - ~ 1,57. Statement 3 defines the upper limit of the loop and the next five statements represent the starting values of the quantities FF, T, W, TT and WW. The quantity F F represents the product

needed to evaluate the integral of the first kind. The sum of the quantities T and W produces the term in the parentheses of Eq. (34). Finally, the combination of T T and WW produces the factor of k appearing in the second term of Eq. (34). The quantity E E represents the factor of F(k, ~) in the same equation. Finally F and E are the desired values of the elliptic integrals of the first and second kinds. I t remains to be decided whether these values are the complete or incomplete integrals. This is accomplished by statements starting with 59 and ending at, but not including, statement 2. The decision is made on the basis of a parameter L whose value is assigned outside the subroutine. If L ~ 1 the computed values F and E represent the incomplete integrals which for the subsequent use are designated as F 1 and E 1. The program exits from the subroutine in a manner which is somewhat different in each case. Thus, for the occultation case for complete eclipses, exit is by the statement GO TO 400. For the transit case of the complete eclipse, the statement is GO TO(400, 46), NN and the value of N N is assigned outside the subroutine. The latter case also applies to partial eclipses.

I. JURKEVICH

83

If L = 2 the computed values of F and E represent the complete elliptic integrals which are designated as F 2 and E2. Then E 2 and F 2 are transferred out of the subroutine as follows. I n the occultation case of a complete eclipse, the computation leaves the subroutine b y GO TO 401. I n the transit case as well as all partial eclipses, the subroutine is exited by a computed GO TO GO T0(39, 47), M. the value of M being assigned outside the subroutine. This completes the computation for the case when CAY < 0.8. I f CAY > 0.8 two slight complications must be accounted for. The first of these is t h a t lim F(1, ~) - , oo. ~0--~12

Depending on the computer there is a definite limit as to how large or small the numerical result of an operation m a y be. I t is then desirable to arrange the computation so t h a t the upper limit is not exceeded. The second complication is associated with the structure of Eq. (30). I t was found inconvenient to develop expressions which could be handled b y a single DO loop over all possible values of the upper index h r. The difficulty is associated with the case N = 0. For this reason it is handled b y a separate set of formulae. The present computation is entered through statement 2. The starting value of the amplitude is immediately compared to ~/2. If F E E ( l ) ----~/2, the statement 17 is used to examine the starting value of the modulus CAY. If CAY ---- 1 the case F(1, ~/2) has been encountered and the program is sent to statement 16. This and the following six statements are basically the same as the statement 59 and the six statements associated with it. The purpose is to identify the complete and incomplete integrals. I n the case of integrals of the first kind, infinity is replaced b y any arbitrary very large number such as 104s or 109e. The exit is by the same means as already discussed. If it is found t h a t F E E ( l ) < ~/2 or if F E E ( I ) ---- ~/2 but CAY =~ 1, the program branches to statement 15 which redefines CAY as a dimensioned variable A(1) and, as before, assigns specific values to parameters K and I for the very same purpose as before. Starting with the statement 51, the transformed value of the modulus and the sine and cosine of GAM ~ 2~0~ -- 9,,-1 are computed. The program is sent to the inverse trigonometric function and quadrant testing subroutine. The latter returns the computation to statement 41 which in turn evaluates the transformed value of the amplitude. Upon this, the departure of A(I -4- 1) from unity is compared with the preassigned tolerance. If the departure is too large, the parameter I is incremented b y unity and additional transformed values of A(I) and F E E ( I ) are computed. If, however, the discrepancy is within the assumed tolerance, the program branches to statement 11 which sets the upper limit of the DO loop defined b y DO6J=

1, N .

If only one iteration is needed to achieve the desired precision/V = 0, and according to Eq. (6) 2 F(k, ~) -- - F(k 1, ~1), l+]c 2

84

Machine Solutions of Light Curves of Eclipsing Binary Systems

Various terms occurring in these equations are implemented in three statements starting with statement 61. Following this, the program is sent to statement 63. This and the next statement are equivalent to

expressed, however, in terms of the available sine functions. The values of the desired elliptic integrals are computed b y F

----

FF*U,

E ----F*(l. ÷ CAY*T)

-- CAY*TT.

The program branches now to statement 59 and the statement immediately following the latter one. The function of these has already been discussed. I f the upper index of the DO loop in question is equal to 1 or greater, the computation is sent to statement 62. This statement and the four following ones set the initial values of various quantities needed to evaluate Eq. (30) b y means of DO 6 J = 1, h r. The quantities T and TT, which follow statement 6, represent the factor in parentheses of the first term of Eq. (30) and the factor in parentheses of the second t e r m respectively. Following this the program proceeds to statement 63. The remaining computation proceeds as indicated for the k ~-- 1, ~0 ----~/2 case. S t a t e m e n t 10 is either a STOP or a d u m m y statement. This completes the discussion of the elliptic integrals subroutine, The evaluation of the inverse trigonometric functions and the associated quadrant testing subroutine were discussed elsewhere (Jurkevich, 1964). The entire subroutine requires between 5012 and 5284 storage locations depending on the type of exit, if an I B M 1620 system is used.

(c) Partial derivatives o/the ]undamental /unctions

General remarks. I t was mentioned earlier in this discussion t h a t to obtain the geometrical depth p from ~(k, p) it is convenient to invert ~(k, p) numerically b y the N e w t o n - R a p h s o n method. This well-known method as applied to the present case is based on Taylor's expansion of a(k, p) in the neighborhood of p. Thus,

o~(k,p + Ap) = ~(k, p) + 0~(k, p) Ap +

. . .

~p The left-hand side of this equation will be known from observations. To start the numerical inversion, an estimate of the first term on the right-hand side is needed. Such an estimate can be obtained b y assuming a suitable value of p. I f the numerical process is started near the instant of internal or external tangency or near the time of minimum light, the required estimate of p is obtained b y noting t h a t in the case of complete eclipses p ~ --1 near the minimum light, and t h a t near external tangency p ~ ÷ 1. For partial eclipses the estimate p ~ 1 remains valid near external tangency. The estimate near minimum light is more difficult to make, although even in the worst case one a t t e m p t is sufficient to indicate what this estimate m a y be. I t is important to realize t h a t the estimate required need not be a close one. The method will converge as long as the initial estimate is within a few tenths of the true value of p.

I. JURKEWCH

85

Starting with an assumed k and the estimated value of p, p -=. ~o, o~(k, Po) appearing on the right-hand side is computed. This value will depart from the known value of ~(k, p) by an amount ~(k, p 4- Ap) -- a(k, Po). To find a closer estimate of p which will decrease this discrepancy, the computation of ~(k, p) is repeated with p given by

P = Po + A p , where Ap is given by

Ap = ~(k, Po + Ap) -- ~(k, Po)/ dta(k'l°u) ~p / The procedure is repeated until the discrepancy between the observed ~(k, l~) and the computed value is within a specified tolerance. To effect this computation most economically ao~]Opmust be evaluated from analytical expressions. Evaluation by numerical differentiation of ~(k, p) is lengthy, yields results of low accuracy, and leads to precision difficulties as [~(k, p + Ala) -- a(k, p)] decreases. The lengthy algebra involved in obtaining equations for O~o~/Ophas been considered by Jurkevich (1964). For the purpose of solving a light curve, one merely needs the final expressions. These are enumerated below.

expressions/or Oot/Op. For a uniform distribution of light over the stellar disk, regardless of the eclipse type, one has Op

zt

(37)

2(1 + kp)

For total limb darkening and an occultation eclipse O:°c

2{(1

1 -- k 2 + ~I } k z ~/~ E(x) ,

+ (~')-- k2

Op

2k' ~/t}

~=2\

F(u)

l+kp

(38)

]

For total limb darkening and a transit eclipse

~(k) °Y---~~ op

2 k, { (k + ~)' - 1 n 2 ~/(k~) F(~)

k2 -- 1 + ~ 1 ~/(k~) ~(~) '

1 [_l -- (k -- l -- kP)' ] i/2. 1 +kp

2

(39)

For annular phases of a transit eclipse

l~(k) _010~ _ ann Op

2

k , / ( 1 - (~ - ~),) ~. (k + ~), zt ~ 3 I =

k(1 + kp) l

-

--

x F(~)

Tn

~ - i - - kp),j

k~ - 1 + ~2 E(u)], 0

(40)

86

Machine Solutions of Light Curves of Eclipsing Binary Systems

The quantities F and E are the complete elliptic integrals. Note t h a t the incomplete integrals do not appear in the above expressions. Values of the partial derivatives for arbitrary values of limb darkening are generated from Ox~ °c 3(1 -- x) O°~ 2x al~ °c -

-

--

- -

ap

3 -- x

-

--

ap

Since for p

<

- -

~ -- x

-

-

,

3 -- x

~(1 -- x) 0%

axa tr -

~-

ap

(41)

ap

(~ -- 1) x ala tr +

(42)

ap

n -- x

ap

--1 o~ -- 1, ax~ ~"'

(~ -

ap

1)

x

~ -- x

alo¢ a " n

ap

I n equations (41--42) is as defined in section 3 (a). Aside from indeterminacy difficulties arising at certain points within the range 0 ~ k _~ 1 and - - 1 / k ~ p ~ 1, the equations given above are sufficient for the inversion of ~(k, p) with respect to p. R e s u l t s o/ evaluation o/ i n d e t e r m i n a t e /orms. Examination of the defining expressions for 0~, 1~oc, l~tr, l~ann and the corresponding partial derivatives with respect to p shows t h a t whenever k = O, as well as when k---- 1, p ---- --1, the defining equations become indeterminate. I n addition, for k = 0 the quantity u becomes indeterminate. These facts lead to difficulties whenever automatic machine computation is used. No m a t t e r what the actual value of the indeterminate expression m a y be the machine will interpret this situation as an overflow. Consequently, the equations leading to indeterminacy must either be rewritten in such a form t h a t indeterminacy does not occur, or, failing this, the computational code must provide for actual values of these expressions as evaluated beforehand. Numerical values of the various alpha functions at the troublesome points mentioned above are known from their definitions or from Tsesevich or Merrill's tables. I t is not known whether general expressions for these functions can be rewritten so as to remove indeterminacy at these points. The present analysis failed to produce positive results in this area. Special approximations are possible and have been produced by Linnell (1965). As far as the partial derivativs are concerned, neither the numerical values nor the analytical expressions at the points in question have been discussed in the literature in sufficient detail. For the computational needs treated here, only the final results are stated. I t is to be noted t h a t in this section the case k ----0 is not emphasized. Although it is of interest as a limiting case, it is of academic interest from the practical point of view. A s u m m a r y of values of ca, ic~oc and %tr at the critical points is given in Table 1 :

TABLE1. Summary of values of °a, la°e, 10¢tr, l•ann at critical points 1

~

~ O_
o~ 0

l~oc

--I

l~tr

l~ann

0

Not deILued

O~

l~oc

l~tr

1

1

l~ann

I. JURKEVICH

87

I n fact for arbitrary values of limb darkening x %d~(k, 1) = 0

%W(k, 1) = O,

xc~°c(k,--1) = 1

x0gtr(k, --1) = ].

regardless of the value of k. The above values follow, of course, from the definition of a. However, to show t h a t the defining expressions satisfy the above conditions, especially at the critical values, involves some lengthy algebra. Table 2 gives the results for the partial derivatives. TABLE2. Summary of values of O%/Op, Ola°¢/Op, 01atr/Op, 01aa~/Op at critical points --1

0

0

0

Not defined

0

0

0

defined

o

0

-f~

-~

0

0

0

2

Not

I t is also useful to note t h a t

O00¢Op k,p< -1 = O, olo~ann ]

is not defined,

[~=i,p>_i

Op

~v(k) 010~tr ]

--

Op [k,p = - l

4 ]¢2 % / [ k ( l -

k)],

:t

olo~ann lim Iv(k) - l ap

-- 0.

(43)

I n Table 2 the dot over c~'s denotes differentiation with respect to p.

Equations/or O~/Ok. Partial derivatives O~(p, k)/Ok are not needed for the computation of preliminary elements. However, they are listed below since they are potentially useful in the differential correction computations of the elements. 00~

2 --

ok 0%(0, p) _ Ok

~

(sin t g -- ~bg),

1 ~,/tll _p~)a) 3~z

(44)

0%(1, 4-1) ak

--0.

Machine Solutions of Light Curves of Eclipsing Binary Systems

88

al0¢oc ok

olo¢°e(0,p)

--

4

~~k 4

3

--

{2(1

+,u)

F

--

+,u)

E},

(45)

01~°c(1, + l )

(1 - - Io2) z

32

Ok

(3

- - O.

Ok

0ql0~tr

(1 + if) (3 + p) F - - 2(3 + ffp) E + 81~ tr

nlv(k)

Ok

a l ~ ( 0 , p) ok

5

-224

(1 + p) { - - ( 9

+ 8 p - - p~) F +

(46)

1 + kp ]J

(12 + 6 p - - 2 p ~) E } .

ol~tr(0,--1)

--0, ak

ala~(1, + l ) __ O. ak

o'~ °"" _

4 ((Ok")

ala ~" (0, p)

5

(p(1 + / x ) F - - (3 + lap) E} --

I(1--p'~

H2 ~/\--~--]

ak

16 t ~ , .

~/[/P(1 - - k)]

(47)

0 + p) {(3 + 2v - p') F -- (6 + 3p -- p') E}, oi~ann(o'--l)

=0,

Ok 01~ ~nn (1, - - 1) = O.

Ok

For central phases alaann/Ok decreases without bound with the decreasing values of k. For an arbitrary value of the limb darkening,

oi~xt~

--

Ok

4

X2

[8k al= ~/(1 - - k) -- ~tiv(k)] (lgtr __ Og) _~_

3~t k s

+

(~ - - 1) -- x

~(1

X)

g -- x x

8°c¢ Ok

010~ tr

ok

For all x ~ 0 or ~ 1 as k decreases the derivative Ox~tr/0k increases without bound. Knowing the general behavior of the above functions, one can proceed now with the computation of the fractional loss of light ~, the corresponding value of geometrical depth p, and finally the solution of light curves.

I. JUaXEWCH

89

4. SPECIFIC COMPUTERPROGRAMS (a) General Outline To implement the computations previously discussed on an IBM 1620 system, a minimum storage capacity of 40,000 decimal locations is required. With the computer memory of this size, it was found necessary to write separate programs for the occultation and transit cases of complete eclipses. A single program to cover partial eclipses of either type was found to be feasible. The simplest of these three cases is a complete occultation eclipse. Since it required the least amount of storage, this computation includes the case k = 0. Physically, it is equivalent to an occultation of a spherical star by a straight edge. Recall that handling of the case k = 0 requires special equations and hence additional storage capacity. The latter could not be made available in the programs for complete transit eclipses and the partial eclipses and consequently these two programs do not include this special case. In all three programs the computation is divided into the following main parts: 1. Computation of elliptic integrals. 2. Computation of the appropriate alpha functions and the corresponding derivatives with respect to p. 3. Formation of the coefficients of normal equations. 4. Inversion of the matrix of coefficients of normal equations. 5. Computation of r~, rb, 4, P0, etc. These blocks are connected by the appropriate logical decisions and some auxiliary computations. The details of the elliptic integrals subroutine have already been discussed. Those of a matrix inversion subroutine can be found elsewhere (Jurkevich, 1964). For reasons of clarity it is felt that before proceeding with the discussion of the entire program it is appropriate to discuss separately the computation of the alpha functions and their derivatives.

(b) Mechanization of x~oo and axaoo/ap The complete mechanization of the computation which yields x~o¢ and Ox~°c/Op for the case of a complete occultation eclipse is given in the F O R T R A N listing of the complete occultation program (p. 92). The flow chart of this subroutine is given in Fig. 2. The listing and the flow diagram will be continuously referred to in the subsequent discussion. The subroutine is entered at statement 700 where the auxiliary quantity (1 --~2) is computed. The next statement yields the apparent separation between the centers of the stellar disks, (~. Since the earlier analysis had failed to produce analytical expressions for x~o¢ and Ox~o¢/Op which would not exhibit indeterminacy at k = 1, p ~ --1, it was decided to store the actual values of these functions as either known from other analyses or as evaluated in the present work. Furthermore, it is to be recalled that the case k = 0 requires a different set of analytical expressions. The purpose of the next three statements is to establish what the values of ]c and are, and, depending on what is found, to send the computation through the appropriate path. 7

B-VIA. ¥oi. 12

90

Machine Solutions of Light Curves of Eclipsing Binary Systems

I-COMPUTEMOOULUS i ANO ~ ( N -

~¢ II9

*0 121

I~-,.p.-, ~L KT I'~

.o/4~.

"5 e (0

o

10

SlN÷m,~¢m

m

COMPUTE LAST THNEE I T[RI~ OF I II O0 L 403

I

~g(-,

-

¢..G~T[ I~00

V' }

UID ' l q

>O~ll?

I

F. (2~.-s,, 2~] SIN ~) g

(:OSqb9

1 U.V

35

,

.0

i

I

i

o"÷

1

FIG. 2. Computation of %, 1~oo, (ao~)/ap, (~]ao~)/a~ in the complete occultation program

Thus, if k (denoted in the program b y RATIO) is negative, which is physically impossible, the computation is stopped b y sending it to statement 10. I f k = 0 the program is branched to statement 120. Here cos ¢8 and sin Fs, as defined in Eq. (25), are computed and the program is sent to the arc tan ~0 subroutine at statem e n t 100. The parameter K, occurring in the computed GO TO instruction of this sub-

•. JURKEVICH

91

routine, is assigned an appropriate value before the subroutine is entered. In this case K ~ 1, and consequently the subroutine returns the computation to statement 71. This statement again examines the value of k. If it is found that k ~ 0 the computation is branched to statement 118. The four statements starting at this point represent 0~, lc~oc' Ooo¢/Op, and Oxa°c/Op for the k ~ 0 case. When these four quantities have been computed the program branches to statement 405. The statement 405 and the next statement yield x~oc and Ox~°c/Op computed from the interpolation formula given in Section 3 (a). The exit from the subroutine is effected by sending the computation to statement 410. Returning now to the original I F {RATIO) statement at location 10012, note that if k is positive the program proceeds to statement 119 which tests whether k ---- 1 or not. If it is found that k > 1 the program is stopped since this case is excluded by definition. If k is found to be less than unity, branching to statement 120 is effected. Computations between statements 120 and 71 have already been discussed. However, since in this case k is positive, statement 71 sends the program to statement 117. At this point K is set equal to 2, and the auxiliary quantity B ~ 1](2~k 2) is evaluated, together with sin Cg and cos ~g as defined in Eq. (25) of this discussion. The angle ~g is evaluated by the arc tan ~ subroutine and the computation is returned to statement 72. This and the following statement yield 0~ __ A L P H 0 and O°~/Op ---- DIAL 0. The next three statements compute the quantities U - ~ U U and V as defined by Eq. (26). Following this the sine and cosine of the amplitude of the appropriate incomplete elliptic integrals are computed. To prevent the computer from printing out the overflow message which would occur for the complete integrals for which cos ¢ = 0, the value of cos ~ is examined to see whether it is equal to zero or not. Since sin ~ ----p and cos ~ ----- ~/(1 -- Iv~) the quantity (1 -- p~) cannot be negative. This case is excluded b y stopping the machine should it occur. If it is found that cos ~ ~ 0, its value is changed to cos ~0 = 10-96 before sending the computation to statement 136. If cos ~ is found to be positive, the computation proceeds directly to statement 136 where the amplitude is evaluated by calling the arc tan ¢ subroutine. Note that the amplitude can be negative, and that its algebraic sign must be preserved. I t is for this reason t h a t the quadrant testing subroutine is not used since it would locate all angles in the positive angle domain. However, since the elliptic integrals subroutine can handle only positive amplitudes, statement 137 changes the algebraic sign of the amplitude should this quantity, as computed in statement 136, be negative. This change of sign will have to be accounted for later. Statement 129 and the following one yield the modulus and its complementary value from expressions defined in Eq. {26). Following this, the parameter 15 is assigned a value of unity and the program branches to the elliptic integrals subroutine. After the incomplete integrals have been computed the subroutine utilizes L to return the program to statement 400. Here the complementary modulus is employed along with L ~ 2 and the assigned value of the amplitude, ~/2, to enter the elliptic integrals subroutine in order to obtain the complete elliptic integrals. As before, L is used to return the computation to statement 401. Recalling now that the sign of the amplitude was determined by the sign of p and realizing t h a t F { k , - - ~ ) = --F(k, ~); E(k,--q~) = --E(k, qg), it is necessary to examine 7"

Machine Solutions of Light Curves of Eclipsing Binary Systems

92

t h e sign of p in o r d e r to assign t h e p r o p e r sign to t h e i n c o m p l e t e elliptic integrals. This f u n c t i o n is p e r f o r m e d b y t h e s t a t e m e n t s 401 a n d 402. I f t h e a m p l i t u d e is p o s i t i v e to s t a r t with, t h e p r o g r a m p r o c e e d s d i r e c t l y to s t a t e m e n t 403. I f not, t h e sign of t h e elliptic i n t e g r a l s is first i n v e r t e d . S t a t e m e n t 403 e v a l u a t e s t h e s u m of t h e t h i r d , fourth, a n d fifth t e r m s of Eq. (26). The n e x t t h r e e s t a t e m e n t s of t h e p r o g r a m y i e l d 1~oc a n d 01~x°c/Opa c c o r d i n g to t h e defining equations. T h e p r o g r a m t h e n p r o c e e d s to s t a t e m e n t 405 a n d t h e a s s o c i a t e d s t a t e m e n t s whose f u n c t i o n s h a v e a l r e a d y b e e n discussed. R e t u r n n o w to s t a t e m e n t 119. I f i t is f o u n d t h a t k ---- 1, i t is n e c e s s a r y to e x a m i n e w h a t t h e v a l u e of p is b e c a u s e a t t h e m o m e n t of i n t e r n a l t a n g e n c y t h e g e n e r a l e q u a t i o n s are n o t valid. I f b = 1, t h e p r o g r a m b r a n c h e s to s t a t e m e n t 121. I f i t is f o u n d t h a t p h a s a n y allowable v a l u e b u t u n i t y , t h e c o m p u t a t i o n is s e n t to s t a t e m e n t 120 a n d all t h e subs e q u e n t i n s t r u c t i o n s , a n d t h e c o m p u t a t i o n p r o c e e d s as d e s c r i b e d before. If, however, p ~ - - 1 , t h e i n d e t e r m i n a c y m e n t i o n e d earlier has been e n c o u n t e r e d . T h e c o m p u t a t i o n is s e n t to s t a t e m e n t 130 where t h e a p p r o p r i a t e values for o~ ( 1, - - 1 ), l a oc(1, - - 1 ), 0%(1, --1)/ap, a n d alaoc(1, --1)/ap are defined. T h e p r o g r a m t h e n b r a n c h e s in t h e usual m a n n e r to s t a t e m e n t 405 for e v a l u a t i o n of x~oc a n d Oxc~oc/Opa n d t h e n to s t a t e m e n t 410. This c o m p l e t e s t h e discussion of t h e s u b r o u t i n e in question. D e t a i l s of similar subr o u t i n e s for c o m p l e t e t r a n s i t eclipses a n d for p a r t i a l eclipses of b o t h k i n d s are g i v e n b y J u r k e v i c h (1964).

(c)

FORTRAN Listing o/ the Program /or Complete Occultation Eclipses and the Flow Diagram

COMPLETE ECLIPSES.COMPUTATION OF ELEMENTS. OCCULTATION. DIMENSION A(7),FEE(7),D(3,4),C(4),SE(3),TI(30),BR(30),P(30) DIMENSION XAL OC(30),DXAOC(30) 08300 READ 450,RATIO,PIN,EPS,X,PI,MI,TOL,Z,TDPR,FITI DO 114 1 = 1,MI O8432 08444 114 READ 452,TI(I),BR(I) KKK = 1 08564 P(1) = P I N 08588 08612 46 PRINT 569 08636 R = RATIO*RATIO 08672 RTIO1 = RATIO TRY -- 0. -~ 08696 08720 EEE = 0. : 08744 DO 7 1 = 1,3 08756 DO 7 J = 1,4 08768 7 D(I,J) = 0. 55 S1 = 0. 08924 08948 DO 600 J J = 1 , M / 08960 IF(SENSE SWITCH 3)56,831 08980 831 FRLL ---- (i. -- BR(JJ))/(I. -- FITI) 09088 I F ( K K K - - 1)10,855,561 O9156 855 I F ( J J - - 1)10,700,856 09224 856 P(JJ) = P ( J J -- 1) 09296 XALOC(JJ) = XALOC(JJ -- 1) 09368 DXAOC(JJ) -~ DXAOC(JJ -- 1) 09440 GO TO 561 09448 56 CS = (COS(TI(JJ)*PI/180.))**2 09544 SN = S Q R T ( ( 1 . - (SIN(GAMR * PI/180.)) * * 2 * CS)/(1. _ _ Z*CS)) 09712 P(JJ) = SN/RADA -- 1./RATIO 08300 C 08300 O83OO

I. JURK~,WCH 09832 09940 i0012 i0068 i0136 |O228 i0252 |0276 10334 10348 i0356 i0380 i0500 i0572 i0580 10636 108O4 i0960 i1044 i1092 11100 i1124 i1184 11340 11376 11484 i1492 i1624 i1852

12008

790 P P = 1. - - P ( J J ) * P ( J J ) D E L T A = 1. + R A T I O * P ( J J ) IF(RATIO)10,120,119 119 I F ( R A T I O - - 1.)129,121,19 121 I F ( P ( J J ) + 1.)120,139,129 139 A L P H 9 = 1. ALOC1 = 1. D I A L 0 = --(2./PI) DAOC1 = 0. GO TO 405 120 K = 1 CS = (P(JJ) + (P(JJ) + R A T I O ) / D E L T A ) / 2 . SN ---=(1. - - CS*CS)**.5 GO TO 10O 71 IF(RATIO)10,118,117 I18 A L P H 0 = GAM/PI - - ( P ( J J ) / P I ) * P P * * . 5 ALOCI = . 5 - - . 7 5 * P ( J J ) + .25*P(JJ)**3 D I A L 0 ---- (2./PI)*SQRT(PP) DAOC1 = - - . 7 5 * P P GO TO 405 117 K -~ 2 B = I./(2.*PI*R) A L P H 0 = B*R*(2.*GAM - - SIN(2.*GAM)) SN =- RATIO*SIN(GAM) CS = 1. - - R * P P / ( 2 . * D E L T A ) GO TO 100 72 A L P H 0 -~ A L P H 0 + B*(2.*GAM SIN(2.*GAM)) D I A L 0 = --(2./PI)*SQRT(1. - - ( ( D E L T A * * 2 - - 1. + R)/(2.*RATIO*DELTA))**2) T E M P ------ - 2 . * D E L T A * * 3 + (3.*R - - 10.)*(DELTA**2) U U = T E M P + (8.*R - - 14.)*DELTA - - 3.*(R*(R - - 3.) + 2.) V = 4.*DELTA**3 - - (16.*R - - 28.)*DELTA SN = P ( J J ) CS = P P * . 5 IF(PP)10,135,136 135 CS = 1.E - - 96 136 P H I = ATAN(SN/CS) IF(PHI)137,139,129 137 P H I ---- - - P H I 129 CAY = .5*SQRT(((1. + DELTA)**2 - - R ) / D E L T A ) C A Y P R = .5*SQRT((R -- ( D E L T A - - 1,)**2)/DELTA) -

93

-

-

-

12176 i2308 13356 13393 12448 12472 12520 12576 12612 12708 12816 12840 GO TO 500 12848 400 CAY = C A Y P R 12872 P H I = PI/2. 13908 L=2 i2932 GO TO 500 12940 401 IF(P(JJ))402,403,403 1302O 402 F1 = - - F 1 E1 = - - E l 13056 13O92 403 T E M P = - - E 2 * F 1 + F2*F1 - - F2*E1 ALOC1 = (PI/2. + ( F 2 * U U + E2*V)/(6.*(RATIO**3)*(DELTA**.5)) + T E M P ) / P I 13336 13500 DAOC1 -~ ((((1. + DELTA)**2 -- R)/2.)*F2 - - (1. - - R + D E L T A * * 2 ) * E 2 ) / S Q R T ( D E L T A ) 1374O DAOC1 = --(2./(PI*R))*DAOC1 i3812 405 X A L O C ( J J ) = 3.*(1. - - X)*ALPH0/(3. - - X) + 2.*X*ALOC1/(3. - - X) 14053 D X A O C ( J J ) = 3.*(1. - - X)*DIAL0/(3. - - X) + 2.*X*DAOC1/(3. - - X) GO TO 410 14292 i4300 500 IF(CAY - - .8)1,1,2 14368 1 A(1) = CAY 14393 FEE(l) ~ PHI 14416 K=3

94 14440 i4464 i4584 i4692 14932 i5004 i5064 i5o88 i5172 i5228 i5264 i5324 i5332 i5488 i5580 i5616 i5624 i5648 i5672 i5696 i5720 i5744 i5768 i5780 i5864 i5936 i5972 i6056 i6164 i6248 i6308 i6392 i6460 i6484 i6508 i6516 i6540 i6564 i6572 i6596 i6688 i6756 X6780 i6804 i6828 |6996 17236 17320 i7344 i7352 i7448 i7552 |7588 i7596 i7632 17700 i7748 i7808 17940

Machine Solutions of L i g h t Curves of Eclipsing B i n a r y Systems I=l 20 S = (1. -- A(I)*A(I))**.5 A(I + 1) = (1. -- S)/(1. + S) D E N = (1. -- A(I)*A(I)*SIN(FEE(I))*SIN(FEE(I)))**.5 SN = S * S I N ( F E E ( I ) ) / D E N CS = C O S ( F E E ( I ) ) / D E N A L = 1. A N G = F E E ( I ) -- 2.*PI 95 IF(ANG)100,100,96 96 A L = AL + 1. A N G = A N G - - 2.*PI GO TO 95 40 F E E ( I + 1) = F E E ( I ) + GAM + (AL -- 1.)*2.*PI I F ( A ( I + 1) -- EPS)3,3,4 4I=I+l GO TO 20 3N=I F F = 1. T----I. T T = 0. W----1. WW=I.

DO 5 J = 1,N F F ----FF*(I. + A(J + 1))/2.

W = W * A ( J + 1 )12. T=T+W W W = W W * ( A ( J + 1)**.5)/2. 5 TT = TT + WW*SIN(FEE(J + I)) E E ----I. -- C A Y * C A Y * T / 2 . F ----F F * F E E ( N + I) E = F*EE + CAY*TT

59 I F ( L -- 1)10,138,139 138 F 1 = F E1 = E GO TO 400 139 F2 = F E2 = E GO TO 401 2 FEE(l) = PHI I F ( F E E ( I ) -- PI/2.)15,17,10 17 IF(CAY -- 1.)15,16,10 15 A(1) = CAY K=4 I=l 51 A(I + 1) = 2.*A(I)**.5/(1. + A(I)) D E N = SQRT(1. -- A(I)*A(I)*SIN(FEE(I))*SIN(FEE(I))) SN = A ( I ) * S I N ( F E E ( I ) ) CS = D E N GO TO 100 41 F E E ( I + 1) = ( F E E ( I ) + GAM)/2. IF(1. - - A(I + 1) -- E P S ) l l , l l , 1 2 121=1+1 GO TO 51 IIN=I--I IF(I -- I)10,61,62

61 F F = 2./(1. + A(1)) T = 1. -- 2./A(2) TT = S I N ( F E E ( l ) ) -- 2.*SIN(FEE(2))/(A(1)**.5) GO TO 63

I. JURKEVICH 17948 17996 is020 18044 i8068 is092 I8104 18212 18284 28320 18428 IS536 18596 is788 is932 18968 is004 19124 19132 i9200 i9224 i9248 I9256 19280 I9304 I9312 i9368 19392 I9440 I9496 19532 i9540 I9596 I9644 i9728 19748 19844 i9904 I9952 20060 20128 20162 20160 20220 90276 20312 20380 20476 20560 20568 $0660 $0728 20752 20760 20844 $0864 50960 21008 21080

62 FF = 2./(1. + A(1)) T=l. TT = SIN(FEE(1)) W=l. WW= 1. DO6J=l,N FF = FF*2./(1. + A(J + 1)) W = W*2./A(J + 1) T=T+W WW = WW*2./(A(J)**.5) 6 TT = TT $ WW*SIN(FEE(J + 1)) T = T - 2.*W TT = TT - 2.*WW*SIN(FEE(N + 2))/(A(N + 1)**.5) 63 U = (1. + SIN(FEE(N + 2)))/(1. - SIN(FEE(N + 2))) U = .5*LOG(U) F = FF*U E = F*(l. + CAY*T) - CAY*TT GO TO 59 16 IF(L - 1)10,18,81 lSFl=l.E+48 El=l. GO TO 400 SIF2=1.E+48 E2=1. GO TO 401 100 IF(CS)21,22,21 22 CS = 1.E - 96 21 GAM = ATAN(SN/CS) IF(CS)31,32,32 3lGAM=GAM+PI GOT034 32 IF(SN)33,34,34 33 GAM = GAM + 2.*PI 34 GO T0(71,72,40,41),K 410 IF(SENSE SWITCH 3)58,60 58 FRLL = 1. - (1. - FITI)*XAL,OC(JJ) DEL = BR(JJ) - FRLL Sl = Sl + DEL*DEL PRINT 65,TI(JJ),BR(JJ),FRLL,DEL IF(JJ - MI)600,66,10 66 PRINT 43,Sl GO TO 46 60 DIFF = FRLL - XALOC(JJ) IF(DIFF)550,551,551 550 DIFF = -DIFF 551 IF(DIFF - TOL)560,560,561 561 DELP = (FRLL - XALOC(JJ))/DXAOC(JJ) P(JJ) = P(JJ) + DELP GO TO 700 560 IF(P(JJ) + 1.)10,406,563 406 IF(RATI0 - 1.)563,407,10 407 RTW = 0. GO TO 564 563 RTW = -DXAOC(JJ)/(2.*DELTA) 564 IF(SENSE SWITCH 2)570,111 570 PRINT 453,TI(JJ),P(JJ),RTW 111 TDPR = TDPR*PI/l80. TT = TI(JJ)*PI/180. IF(SENSE SWITCH 1)19,24

95

96 21100 21124 21132 21240 21348 21492 21564 ~1612 21636 ~1708 21756 21780 21792 21804 22056 22128 22164 22176 22188 22404 ~2440 22464 22500 22524 22536 22584 22596 22848 ~3028 23124 23220 23232 23364 23508 ~3520 23820 23888 23924 23932 23944 23992 ~4004 24220 24352 24400 24448 24508 24568 24628 24640 24808 24832 24016 24976 25~00 25036 25096 25120 25168

Machine Solutions of Light Curves of Eclipsing Binary Systems 19 R T W = 1. GO TO 27 24 I~TW = I~TW/(1. - - Z*COS(TT)*COS(TT)) 27 AA ----I~TW*(1. - - Z*COS(TDPR)*COS(TDPR)) B ~ RTW*(1. - - Z*COS(TT)*COS(TT))*P(JJ) C(1) = B * P ( J J ) - - AA c(2) = 2.*(AA + B) C(3) = R T W C(4) = RTW*SIN(TT)*SIN(TT) E E E = E E E + C(4)*C(4) L-~I DO 13 J ~ 1,3 DO 35 K ~ L,4 35 D(J,K) = D(J,K) + C(J)*C(K) 13L=L+I 600 CONTINUE DO 36 1 ~ 1,3 DO 36 J ~ 1,3 36 D(J,I) ----D(I,J) D 0 , I ) = 1./D(1,1) K=2 37 K K = K - - 1 P P = 0. DO 38 1 ~- 1 , K K c(I) = o. D o 69 J = 1 , K K 69 C(I) = C(I) + D(I,J)*D(ff,K) 38 P P = P P - - D(I,K)*C(I) P P = P P + D(K,K) D(K,K) = 1./PP DO 8 I = 1 , K K D(I,K) = - - C ( I ) / P P D(K,I) = D(I,K) DO 8 J = 1 , K K 8 D(I,J) = D(I,J) + C(I)*C(J)/PP I F ( K -- 3)9,39,39 9K:K+I GO TO 37 39 DO 67 I = 1,3 c ( i ) = o. D o 14 j = 1,3 14 C(I) = C(I) + D(I,J)*D(J,4) 67 TRY = T R Y + C(I)*D(I,4) P P --~ MI - - 4 $2 ~ ( E E E - - TI~Y)/PP R12 = D(1,2)/SQRT(D(1,1)*D(2,2)) R13 = D(1,3)/SQRT(D(1,1)*D(3,3)) R23 = D(2,3)/SQRT(D(2.2)*D(3,3)) DO 42 I = 1,3 42 SE(I) = SQRT(S2*D(1,1)) P R I N T 568 P R I N T 43,C(1),SE(1),C(2),SE(2),C(3),SE(3) P R I N T 47,R12,R13,R23,S2 PRINT n0,RTI01 RATIO = C(1)/C(2) D E N = C(1) + Z*C(2)*C(2) SN ~ SQRT(C(3)) CS ~ SQRT(1. - - C(3)) T D P R = ATAN(SN/CS)*180./PI

I. JURK~WCB 25240 25288 25336 ~5540 25576 ~5612 25660 25728 25752 25776 ~5848 ~5904 25928 25976 26060 26120 26216 26312 ~6516 26588 26660 26708 26756 26816 26876 26912 26932 26980 27000 27024 27032 27100 27198 27226 27258 27296 27324 27410 27552 27676 27802 27946 27988 ~8036

(d)

97

CC1 = C(1)*C(1)/DEN CC2 = C(1)*C(2)/DEN G = (CC2 -- CC1)*'2"(1. -- Z*COS(TDPI~)*COS(TDPR)) + CCI*(1. -- C(3)) RADA = CC1/SQRT(G) RADB = CC2/SQRT(G) SN = SQRT(CC1/G) IF(SN -- 1.)52,52,50 59 PRINT 43,SN SN = 1. 52 CS = SQRT(1. -- SN*SN) IF(CS)10,45,44 45 CS = 1.E -- 96 44 GAM = ATAN(SN/CS) SB = SQRT(1. -- Z/(SN*SN)) ANG = ATAN(SN/(CS*SB)) RSBR = (1. -- FITI)/(FITI*RATIO*RATIO) AA = RADB*RADB*(1. -t- RATIO)**2 SN = SQRT((AA*(1. -- Z) -- COS(GAM)**2)/(SIN(GAM)**2 -- AA*Z)) CS = SQRT(1. -- SN*SN) TH E T E = ATAN(SN/CS)*lSO./PI GAMR = GAM*180./PI ANG = ANG*180./PI PRINT 48, RATIO,RADA,RADB,GAMR PRINT 49,RSBR,TDPR,THETE,ANG KKK = KKK + 1 IF(SENSE SWITCH 4)26,25 26 RATIO = (RATIO + RTIO1)/2. 25 IF(SENSE SWITCH 3)830,46 830 PRINT 94 GO TO 55 450 FORMAT(F8.O,F4.O,F8.O,F4.O,F10.0,I4,F8.0,F6.0,F6.O,F5.0) 569 FORMAT(/6H TI,11H POC,10H RTW) 452 FORMAT(F9.4,F15.8) 453 FORMAT(F.4,F10.5,F.4) 65 FORMAT(FS.3,F9.4,F9.4,F12.8) 43 FORMAT(E14.8,E17.8) 568 FORM_AT(8H C(I),19H SE(I)) 47 FORM_AT(5H R12 = FS.3,3X 5H R13 = F5.3,3X 5H R23 ----F5.3,3X 4H $2 = E14.8) 48 FORMAT(6HRATIO = F8.5,7H RADA = F7.5,8H RAI)B = F7.5,8H GAMR = F7.3) 49 FORMAT(6HJA/JB = F7.5,8H TDPR = F7.3,8H THETE = F7.3,8H ANG = F7.3) 94 FORMAT(3 X 2HTI,8 X 4HOBFL,5 X 4HCOFL,6 X 3HDEL) 119 FORMAT(6HRTIO1 = F8.5) 10 STOP END

Program/or the Solution o/Light Curves o~ Complete Occultation Ecllpses

T h e specific p r o g r a m for c o m p l e t e o c c u l t a t i o n eclipses is g i v e n in t h e F O R T R A N listing w h i c h was e m p l o y e d in t h e discussion of t h e m e c h a n i z a t i o n of x~oc a n d Ox~°c/Op. T h e flow d i a g r a m of t h e p r e s e n t c o m p u t a t i o n is shown in Fig. 3.

98

Machine Solutions of Light Curves of Eclipsing Binary Systems START

TOPR. F I T ;

(-i

5

~s3

..j

40T

~ C O M P URi~lO T[ O~S~AO~em~,~T~ ANO 0

OFF 24 C~ J~

L,, I

I nt,ls N(coca Fo* co(,~',c*cu[:~s,J°~. . . . ~lm*Oi

OF~, i

>~

t~-~o~t rH

NO

Z5 OFF

25

r SET~

COMPUTE TH(O~[T~CAL p ~,VENoe~,.~o va~m OF ~ s [ J~*~E 8 c~mcotl oa~o ( ~

C ....' ~. . . . (, .... w u o~ pKv~o~

t

n, ~ m . ~ umN0mm KZl , |

)

:.,=~..)

N

FIG. 3. Complete occultation eclipse. Solution of light curve. The program sets aside sufficient room in the memory to store the following quantities: A (7) :FEE (7) D(3, 4) C(4)

SE(3)

----- Seven values of the modulus of the elliptic integrals. ~ Seven values of the amplitude of the incomplete elliptic integrals. = Matrix of normal equations and the column matrix of the right-hand sides of the normal equations. ~ Four auxiliary quantities used to form the coefficients of the normal equations. Locations set aside for three of these quantities are also used to store the u n k n o w n s resulting from the solution of the normal equations. ~ Standard errors associated with the three quantities C1, C2, Ca--the solutions of the normal equations.

I. JURKEWC~ T I (30) B R (30) P (30) XALOC (30) DXAOC (30) The input RATIO PIN EPS X PI MI TOL Z TDPR

FITI

99

---- Thirty values of phase angle. ---- Thirty values of rectified intensity. = Thirty values of geometrical depth. ~ Thirty values of x~°c(k, p). = Thirty values of 0~°¢(k, p)/Op. record consists of: = Initial estimate of the ratio of radii k. ~ The starting value of the geometrical depth associated with a specific observed value of ~o¢. = Tolerance employed in computing the elliptic integrals. ~ Limb darkening of the smaller component of the binary system. ~ 3.1415926. = Number of observations used. -~ Tolerance employed in computing geometrical depth p. ~ 2 e sin 2 i resulting from the rectification computation. This quantity enters the expression for photometric ellipticity. ~ An initial estimate of the phase angle associated with the instant of internal tangency. Equations of condition used by this program are written in terms of this angle rather than the phase at the moment of external tangency. = An estimate of the intensity at the moment of internal tangency. This quantity is read off the graph of the light curve.

The two instructions following the input record store the observed phase angles and rectified intensities in the memory. Next, the parameter K K K is assigned the value of unity. This parameter is needed for the following reason. During the first iteration in computing the ratio of radii, k, all values of p and x~oc are not defined. Therefore in computing each value of p and x~oc for a given phase the program starts with p and ~oc which have been computed for the previous observation. However, during the second and succeeding iterations the appropriate locations store the values of p and ~oc computed during the previous iteration. I t is these values which are now used as the starting values. The control over which values are to be used is effected by the value of the parameter KKK. Due to the fact that no value of p is defined initially, it is necessary to define one. This is accomplished by the instruction P(1) = PIN.

The next instruction causes the following heading to be printed TI

POC

RTW.

Here T I = observed phase angle. POC = the value of geometrical depth computed for the given value of k and corresponding to the given observed value of x~oc. RTW =

Ox~°c/aP. This quantity is proportional to the value of the intrinsic

2~ weight of the given point as given by =

-

2~(1

--

z cos 2 0')

100

Machine Solutions of Light Curves of Eclipsing Binary Systems

If desired, the corresponding quantities will be printed in columns under these headings. Printing can be suppressed b y employing SENSE S W I T C H 2 in the OFF position. However, printing of the heading cannot be suppressed. The next four instructions define the following quantities: R ~ k 2. R T I O 1 = The input value of k used for the given iteration. TRY ~ 0. This is the initial value of the quantity which is used to form the sum of the products of the unknowns and the right hand sides of the normal equations. This sum is used to form the residual variance. EEE -----O. The initial value of the quantity which eventually results in ~ w sin 40i. I t is also employed in computing the residual variance, i I t should be recalled t h a t this method of forming sums is standard in computer practice. The next four instructions set initially all terms of the matrix of the normal equations and the quantity S 1 equal to zero. The latter quantity represents ~ (l observed-/computed)y, i where the /'s are the rectified intensities. I t must be pointed out t h a t this quantity is of interest since the residual variance obtained from the solution of normal equations is a measure of the quantity (sin 20ob s -- sin 2 0 comp)~, which is not a measure of the theoretical fit of the light curve. Reasons for this have been mentioned in the remarks on the choice of the method for solving light curves (Section 2 (c)). The next instruction DO 600 J J ~ 1, MI is the first instruction of the overall DO loop which is used to compute the geometrical depth p, intrinsic weights x/w and the coefficients of normal equations. The function of SENSE S W I T C H 3 is to cause the program to compute either a theoretical light curve or to solve the observed light curve. I n the former case this switc h is in the ON position and the program branches to statement 56. During the first iteration the O N position must not be used, since at this time no elements are available. Consider now the second case; namely, SENSE S W I T C H 3 is in the O F F position. The program proceeds to statement 831 which computes the observed value of xc~oc from x oc ¢~obs

1 --

-

-

-

1jj

-

1 --

~a

for the J J - t h observation. The quantity X~obsOCis in the program represented b y the symbol FRLL. Following this, the value of K K K is tested. Since initially it is equal to 1, the program branches to statement 855 where the value of J J is tested. For the first observation J J = 1 and the program proceeds to statement 700 which is the first statement of the subroutine for the computation of ~ o c and Oxo~°c/Op. This subroutine has been discussed in detail in the previous sections. I t m u s t be pointed out t h a t neither K K K nor g g can assume values smaller t h a n unity. The initial input to the above subroutine is P I N which, hopefully, corresponds reasonably well to the first value of . .C%bs, .. and the quantity RATIO. The computation leaves the subroutine through statement 410, and since the SENSE S W I T C H 3 is OFF, the program branches to statement 60. This statement takes the difference between the observed ~ o c and the computed value ~o¢°c(k, P I N ) and compares its absolute value to the preset tolerance

I. JURKEVICH

101

T O L . This c o m p a r i s o n is p e r f o r m e d i n s t a t e m e n t 551. I f t h e d i s c r e p a n c y is larger t h a n t h e tolerance, s t a t e m e n t 561 p r o d u c e s a c o r r e c t i o n i n c r e m e n t in p f r o m X

OC

Ap -~ D E L P - - aobs

X

OC

O~com O ~¢oJap X

-

-

O¢

c o m p u t e s a new v a l u e of p f r o m P ( J J ) ---- P ( J J ) + D E L P , a n d r e t u r n s t h e p r o g r a m t o t h e x~oc, oxo~oc/ap s u b r o u t i n e ( s t a t e m e n t 700). This p r o c e d u r e is r e p e a t e d u n t i l t h e difference (x~obs oc _ x~com) oc is r e d u c e d to a v a l u e e q u a l t o or less t h a n TOL. W h e n this occurs, t h e p r o g r a m proceeds to s t a t e m e n t 560. S t a t e m e n t s 560 a n d 406 h a v e t h e p u r p o s e of d e t e c t i n g w h e t h e r k = 1 a n d p _----1 s i m u l t a n e o u s l y . S h o u l d t h i s be t h e case, t h e p r o g r a m proceeds to s t a t e m e n t 407 which defines t h e w e i g h t of t h e p o i n t to be zero. The r e a s o n for t h i s is t h e f a c t t h a t t h e weight a t i n t e r n a l t a n g e n c y grows w i t h o u t limit, a n d t h e r e f o r e t h e entire solution w o u l d be d e t e r m i n e d b y a single o b s e r v a t i o n (Jurkevich, 1964). F o r a n y o t h e r c o m b i n a t i o n of k a n d p t h e p r o g r a m b r a n c h e s to s t a t e m e n t 563 which e v a l u a t e s t h e q u a n t i t y R T W . H o w e v e r , w h a t e v e r t h e source of R T W , t h e p r o g r a m t h e n proceeds to s t a t e m e n t 564 which allows t h e o p e r a t o r , a t t h e c o n t r o l of S E N S E S W I T C H 2, to p r i n t t h e r e s u l t i n g v a l u e s of P ( J J ) a n d R T W or to b y p a s s t h e p r i n t e d o u t p u t b y sending t h e c o m p u t a t i o n to s t a t e m e n t 111. I t is a t this p o i n t t h a t t h e c o m p u t a t i o n of t h e e x p a n d e d m a t r i x of n o r m a l e q u a t i o n s begins. To clarify t h e m e c h a n i z a t i o n of this c o m p u t a t i o n , i t is useful to digress f r o m t h e m a i n discussion a n d to i n d i c a t e t h e expressions t h a t are being p r o g r a m m e d . R e c a l l t h a t for t h e c o m p l e t e eclipses t h e e q u a t i o n of c o n d i t i o n is given b y x/w [(1 - - z cos 2 0) p~ - - (1 - - z cos 2 0")] D 1 ~- 2 x/w[(1 - - z cos 2 0) p + (1 - - z cos 2 O")] D 2

+ x/w" D2 -----x/w sin 2 0, which w i t h o b v i o u s definitions can be w r i t t e n as

aiD1 -}- a,z~D2"-k aaDa ~- ~/w sin 2 0. T h e n o r m a l e q u a t i o n s are t h e n given b y

A3A 1 AaA 2

A~]

D3

x/w" A 3sin ~ 0

T h e coefficients are given b y A~

-= ~ w[(1 - - z cos 2 0) p~ - - (1 - - z cos 2 0")] 2,

A~A 2 ----~ 2w[(1 - - z cos ~ 0) p~ - - (1 - - z cos 2 0 , ) ] [(1 - - z cos ~ 0) T + (1 - - z cos 2 0")], A1A a = ~ w[(1 - - z cos ~ O) p9 __ (1 - - z cos 2 0")], A22

= ~ 4w[(1 - - z cos 2 0) p + (1 - - z cos 2 0")] 2,

=Ew, A2Aa ---- ~ 2w[(1 - - z cos 2 0) p + (1 - - z cos 2 0")] 2.

102

Machine Solutions of Light Curves of Eclipsing Binary

Systems

Nineteen statements starting with statement 111 generate the triangular matrix

and the column matrix

The first two instructions in this sequence convert 6” = TDPR and the observed phase angle to radians. The next three instructions are included only in the program treating complete occultation eclipses. At the control of SENSE SWITCH 1, the operator can set the weight of observation equal to 1. In the OFF position this switch causes the program to proceed to statement 24, which produces the actual value of the intrinsic weight of an observation. The expression involved is given by Jw

=

a”~79

_

2S(l -

2 COGf3) *

The quantity AA computed by statement 27 is equal to Jw(1 -

2 COG,“)

and the quantity B to Jw( 1 -

2 cos2 e) p .

Both of these quantities occur in every coefficient of the normal equations except AX. It should be clear that C(1) represents A,, C(2) is equivalent to A,, C(3) represents A,, and C(4) gives the right-hand side of the equation of condition. The variable EEE will eventually represent c w sin 40i. The last six instructions in the sequence produce partial sums which eventu&ly yield the coefficients themselves. The last instruction 600 CONTINUE causes the program to return to the first instruction of the overall DO loop DO600JJ=

l,MI

and to start processing the next observation. As before the value of %$& is computed. Since KKK continues to be equal to 1, the program proceeds to statement 855. However, this time JJ = 2 and the program branches to statement 856. The three instructions starting at this point define the starting values of p, %E,C,r, and ax&,,,/ap to be used in obtaining I, for the second point. Note that the starting values are those which have been computed for the first point. These values are sufficiently close to those for the second point provided the interval between observations is not too large. Once these starting values have been defined the program proceeds to statement 561 rather than 766 because all quantities needed to compute Ap are available. Furthermore, one may as well do this without testing for the departure of %z& from %~&, , because the departure is almost certainly bound to be larger than TOL. Thus, for the second and succeeding observations the %oc, iYVc/i3~ subroutine is entered only after evaluation of

I. JURK~WCH

103

/lp and p ---- Po ~- A p . The computations which follow differ in no way from those already described. After the last observation has been processed, the quantities D ( J , K ) represent the triangular coefficient matrix of normal equations and the column matrix of the right-hand sides of these equations. The quantity E E E is now equal to ~ w sin 40 i . Instructions i DO 36 1 = 1, 3 36 D(J, I) = D(I, J) expand the triangular matrix to a square one by setting the appropriate off diagonal elements equal to each other. The next nineteen instructions accomplish the inversion of the matrix of normal equations according to the method of successive supplementation of rows and columns of a known inverse matrix. Following the inversion, the unknown coefficients, now designated by C(1), C(2), C(3) and the quantity [C(1) ~ x/w" A1 sin2 0 + C(2) ~ ~/w. A~ sin2 0 ~- C(3) ~ x/w. A a sin~ 0l are computed. This is accomplished by instructions 39 DO 67 1 = 1, 3 ,

•

•

67 TRY ---- TRY ~- C(I)*D(I, 4). The computation of C(I)'s is accomplished by multiplying each row of the inverted matrix by the column matrix representing the right-hand sides of the normal equations. Seven instructions following this computation produce the residual variance $2, correlation coefficients R 12, R 13, R23, and the standard errors SE(1), SE(2), SE(3) associated with the quantities C(1), C(2), and C(3). At this point the following printed format is produced, C(I)

SE(I)

c(1) SE(1) C(2) SE(2)

c(3) SE(3) R12

R13 R23

$2

RTIO1 The last quantity is the value of b which was used as the input at this iteration stage. Immediately after this, what is hopefully, an improved value of k is computed from RATIO = C(1)/C(2). The improved value of 0", following directly from the coefficient C(3), is produced by 25168 TDPR -----ATAN(SN/CS)* 180./PI. To avoid confusion, it is appropriate to point out that the variables C(1), C(2), C(3) employed at this stage in the program are equivalent to the coefficients D1, D2, D a employed in the section summarizing the required equations (Section 2 (d)). The quantities CC1

104

Machine Solutions of Light Curves of Eclipsing Binary Systems

and CC2 of the program are equivalent to C 1 and C 2 of t h a t summary. The latter, computed in terms of C(I)'s, are given b y CC 1 : C(1)* C(1)/DEN, CC 2 = C(1)* C(2)/DEN. The next instruction produces the quantity G used in Eqs. (13), (14) and (15) of the summary. Having obtained CC1, CC2 and G, one can compute ra and rb which in the present program are designated by R A D A and R A D B respectively. The above quantities are also sufficient to compute the improved rectified inclination. This computation is carried out b y instructions SN : SQRT(CC l/G), 44 GAM : ATAN(SN/CS). Since for the total eclipses i m a y be very close to 90 ° and since we are dealing with a least squares solution, there is no guarantee t h a t sin ir will not exceed unity. A test for this condition is included, and should sin ir exceed unity, the resulting value is printed. However, before the next instruction is executed sin ir is set equal to unity. The next two instructions yield an estimate of the true value of inclination designated b y the symbol ANG. The quantity R S B R yields the ratio of mean surface brightnesses according to Eq. (19) of the summary. The next four instructions produce an improved estimate of the phase angle at the external tangency, designated b y T H E T E . At this point the following quantities are printed RATIO RSBR

RADA TDPR

RADB THETE

GAMR ANG

where GAMR is the computed rectified inclination. Following this output, the value of the parameter K K K is increased b y unity. Presently, it must be pointed out t h a t had the original choice of k = R T I O 1 been quite close to the true value, R A T I O as computed above would differ from R T I O 1 b y a small amount. A still better approximation m a y then be obtained b y repeating the entire computation with R A T I O as the new value of k. The procedure will converge quite rapidly. Should, however, the initial estimate of k be relatively far from the true value, the above procedure will converge slowly. There exist numerous methods to speed up the convergence. However, most of these require the computation of several additional values of /c before some kind of interpolation can be performed. A simpler procedure, and quite an efficient one, is based on the fact t h a t for most initial estimates of k, as long as the process is a convergent one, the successive approximations to k will oscillate around some stable value. Thus, if the initial estimate of ]c is relatively poor, a much better approximation is obtained b y merely taking a mean value of two successive approximations of k as the starting value for the next iteration. The present program provides the operator with a means of using either each successive value of k or the mean of two successive values of k as the starting point for the next iteration. This is accomplished under control of SENSE S W I T C H 4. If this switch is in the ON position, the mean value is taken; otherwise, the computed value of k is taken directly. I n either case the program proceeds to statement 25 whose function remains the same as before.

I. JUI~KEVIOH

105

Assuming t h a t it is too early to examine a theoretical fit of the light curve, this switch will return the computation to statement 46 which initiates the next iteration. The only difference between the second and succeeding iterations, and the first one is the fact t h a t the instruction I F ( K K K -- l) 10, 855, 561 will now branch the computation to statement 561 before sending it to the zoo°%O~°¢/Op subroutine. The program does not terminate the computation automatically. This decision is left to the operator who can manually interrupt the computation when the two successive values of k differ b y less t h a n an amount decided upon b y the operator. When this occurs the operator can, upon the last complete printed output, set the SENSE S W I T C H 3 to the ON position. Under these conditions statement 25 will branch the computation to statement 830 which causes the machine to print the following heading TI where T I OBFL COFL DEL

OBFL

COFL

DEL

---- phase angle ~ observed fractional luminosity ----computed fractional luminosity = OBFL -- COFL.

Following this the machine branches to statement 55. Since SENSE S W I T C H 3 is now in the ON position the program proceeds to statement 56. The purpose of this and the next two instructions is to compute the theoretical value of p. Using the last value of k and the theoretical value of p, the program proceeds to compute ~aoc. Statement 410 sends the computation to statement 58. This instruction computes tim theoretical value of the intensity. The next two instructions yield the difference between the observed and computed values of the intensity and the partial sum of the squares of these differences. At this point the phase angle, the observed and computed values of the intensity, and the difference between the two are printed. The above computations are repeated for every observation. When the last observation has been processed, instruction 66 causes the program to print the quantity S1 = ~{/obs -- l¢omp)~" The program then returns to statem e n t 46 ready for the next iteration should the operator so desire. This completes the discussion of the complete occultation program. Details of similar programs for the complete transit eclipses as well as for partial eclipses can be found in the work of Jurkevich (1964).

(e) Program Summary The meaning of various symbols appearing in the input record have been defined in the previous section. Columns on the card reserved for these quantities are as follows :

Card 1:RATIO(1 -- 8), PIN(9 -- 12), EPS!13 -- 20), X(21 -- 24), PI(25 -- 34), MI(35 -- 38), TOL(39 -- 46), Z(47 -- 52), TDPR(53 - - 58), FITI(59 -- 63). All quantities except MI are in the floating point form without an exponent. The quantity MI is a fixed point number. The next MI records contain observations in the form T I (1-9). B R (10-24), one pair to a card. Both quantities are in the floating point form without an exponent. 8

B-VIA Vol. 12

106

Machine Solutions of Light Curves of Eclipsing Binary Systems

O b s e r v a t i o n s m u s t be o r d e r e d according to t h e v a l u e s of phase. H o w e v e r , i t does n o t m a t t e r w h e t h e r one s t a r t s a t a p o i n t n e a r t h e i n t e r n a l t a n g e n c y of one b r a n c h , p r o c e e d i n g u p t h i s b r a n c h , a n d t h e n d o w n t h e o t h e r b r a n c h or w h e t h e r one s t a r t s n e a r t h e e x t e r n a l t a n g e n c y of one b r a n c h , p r o c e e d i n g t h r o u g h t h e m i n i m u m to t h e o t h e r p o i n t of e x t e r n a l t a n g e n c y . I t is also permissible t o a r r a n g e cards so t h a t t h e y describe t h e reflected minim u m . I n a n y of these cases one m u s t choose P I N such t h a t i t r e a s o n a b l y c o r r e s p o n d s to a v a l u e which m a y be e x p e c t e d to e x i s t a t t h e first d a t u m p o i n t . PROGRAM SWITCH settings: SWITCH

ON

OFF Theoretical intrinsic weights used in the computation Printout by-passed

All weights set equal to unity The following quantities are printed: TI(JJ) = observed phase angle, P(JJ) = computed geometrical depth, RTW = quantity proportional to weight Theoretical light curve is computed. The following quantities are printed: TI(JJ), BR(JJ), F R L L = computed intensity, DEL = BR(JJ) -- F R L L The value of k employed for the next iteration is given by

Computer processes data solve the light curve

The value of k used for the next iteration is equal to the last value computed

k(m_l) • k(m-1) km= 2 where m = order of iteration

to

'

Output Solution o] light curve: U n d e r c o n t r o l of S E N S E S W I T C H 2 t h e following t a b l e is p r o duced TI POC RTW ,

°

•

,

,

,

°

TI(JJ) P(JJ)

R e g a r d l e s s of t h e p o s i t i o n of S E N S E errors are p r e s e n t e d as follows:

SWITCH

C(I) c(i) C(2)

SE(I) SE0) SE(2)

C(3)

SE(3)

°

,

Ox~°c

1

Op

26

2 the unknowns and the associated standard

S i m u l t a n e o u s l y , t h e c o r r e l a t i o n coefficients a n d t h e t r i a l v a l u e of k are p r e s e n t e d in t h e following f o r m : R 12 - - (-4-).XXX S 2 -= + . X X X X X X X X

R 13 = ( ± ) . X X X E(±) YY

R 23 ---- ( - V ) . X X X

R T I O 1 ---- A - X . X X X X X

I. JURKEVIeH

107

If, due to the least squares adjustment of data, sin ir ~ 1, the value of sin ir is printed as +.XXXXXXXX E(±)YY. The final output appears as RATIO---- - ~ . X X X X X RADA ~ + . X X X X X JA/JB

= q-.XXXXX

TDPR

~ q- X X . X X X

RADB = - ~ . X X X X X THETE

= q- X X . X X X

GAMR

=

q-

XX.XXX

ANG

=

q- X X . X X X

The above output is repeated in its entirety for every iteration. Theoretical light curve If, following the printing of ANG, SENSE SWITCH 3 is set in the ON position, the following tabulation is produced

OBFL

TI .

.

.

.

TI(JJ)

.

.

.

.

COFL

DEL

FRLL(com-

BR(JJ)-FRLL

.

BR(JJ)

,

° ,

puted intensity) .

,

.

.

.

°

,

°

•

°

°

°

Operating Procedure. For the first iteration in the solution of a light curve keep SENSE SWITCH 3 in the OFF position. Positions of SENSE SWITCHES 1 and 2 are set at the operator's convenience according to the PROGRAM SWITCH summary. The position of SENSE SWITCH 4 is immaterial for the first iteration. Following the printing of ANG the operator can interrupt the computation by pressing the INSTANT STOP key on the computer console and examine the discrepancy between RTIO1 and RATIO. If the discrepancy is too large it will pay to set the SENSE SWITCH 4 in the ON position and then restart the computation by pressing the START key on the console. Following the second printout the operator could repeat the above procedure and either keep the SENSE SWITCH 4 ON or if the discrepancy is reasonably small, turn it OFF. In the latter case one can leave the machine running without further interruption until such time when the operator feels that the difference (RTIO1-RATIO) has reached a sufficiently small value. At such time the machine can be stopped by pressing the INSTANT STOP key, SENSE SWITCH 3 turned ON in order to compute a theoretical light curve, and the computation of the latter initiated by pressing the START key. Upon printing the sum of the squares of residuals, the program is ready for the next iteration. This is signaled by the printing of the heading TI

POC

RTW.

If one desires an additional iteration, SENSE SWITCH 3 must be turned OFF during printing of the above heading. The above procedure is somewhat awkward. The entire sequence of events describe dabove could be made automatic if more storage capacity were available. However, the user, depending on his computing facility, can effect such modifications without difficulty. (f) Typical Computations The application of the program discussed in the previous sections is illustrated by applying it to the ¥ Z Cassiopeiae system for which high-quality observations as well as a good solution are available (Kron, 1939). The primary minimum of this system is due to a transit eclipse. Independent solutions for both the primary and secondary eclipses are 8*

108

Machine Solutions of Light Curves of Eclipsing B i n a r y Systems TABLE 3.

-,~0

t . O OP~ . 5 7 4 0 0 . 0 0 0 0 0 1

~[~Ydi~fle

YZ Cassiopeiae (Min I) Sample C o m p u t a t i o n .50 3,1"415926,00010

m. I . ~ 1 I . x n ~ J m n m a . q U - ], a l U . q - x ~ l a ~ l q e a

,00000 3.200

.9357

a q a ~s' qmu n q ~ ~ n qm a ~ a a q a

x ~,~ n x . ~ . x - nJ

-1,,05 -S, 2 7 .7080 L ! ~ 4 1 1 OI e ~ 1 e " ~ a ~ u ~ i 1 ~ t ~ I n ~ p ~ n m ~ q m m ~ q ~ x a q ~ M ~ e a ~ q ~ s ~ q u ~ s q ~ n ~ a a ~ a q a ~ x ~ n a ~

TI

PTR

RTW

+3.2640 +3.7640 +4.2150 +4.4890 +4.9960 +5.3670 +5.7540 +6.2450 -b6'6160 %7.0510 -~7.4140 ~-8.0340 -~8.4370 ~-8.8240 +9.2920 +9.5900 ~-10.1300 ~-10-6630 ~-11.1050 ~-11.5320 ~-11.9190 %21.4750 +2.6430 ~-2.1440 +1.5310 +.9590 +.3300

--.97951 --.88560 --.80685 --.73051 --.65076 --.55870 --.46691 --.37175 --.29640 --.20105 --.11315 +.00414 +.09470 ~-.17786 +.27888 ~-.33441 +.44401 ~-.57778 +.65711 +.74947 +.84436 +.92616 --1.06153 --1.14867 --1-20270 --1.24765 --1"30788

+.4121 +.5311 +.5452 ~-.5361 -~.5159 +.4865 ~-4548 ~.4214 ~-3954 %.3634 -~.3349 +.2986 +.2717 ~-.2479 -~.2199 ~.2049 -~.1759 +.1411 +.1205 ~.0960 +.0694 +.0433 +-2465 ~-.2283 +.2208 +.2159 -~.2106

c(1) +.59520510E -- 02 ~-.11221750E -- 01 +.29331100E -- 02 R12 = + . 8 1 9 R A T I O ---- +.53040 J A / J B ~ +.24426 P0 = --1.3458

First Iteration

SE(I) +.15314873E -- 03 +.74357132E -- 04 +.58355787E - - 04 R13 = + . 4 3 8 R A D A = +.07708 TDPR ~ +3.104 RTIO1 = +.57400

R23 = - - . 0 2 0 R A D B = +.14532 T H E T E - - ~ +12.631

$2 : +.37460869E -- 98 GAMI~ : +87.616 ANG---- +87.616

Fourth Iteration

c(I) + ' 5 9 6 6 4 1 4 0 E -- 02 +.11212708E - - 01 +.28991870E - - 02 R12 = + . 8 0 6 R A T I O ~ +.53211 J A / J B = +.24269 P 0 = -- 1.3432

SE(I) +.16677224E -- 03 +.77844952E -- 04 +.66906128E -- 04 R13 = + . 4 7 0 R23 = --.011 $2 = +.41460869E - - 9 8 R A D A = +.07717 R A D B = +.14503 GAMR ~ +87-629 TDPR +3.086 T H E T E = +12.621 ANG : +87.629 R T I 0 1 = +.53196

I. JuRxEvm~ Theoretical Light Curve

TI OBFL +3.2640 +.70990 +3.7640 +.72240 +4.2150 +.73520 +4.4890 +.74870 +4.9060 +-76350 +5.3670 +.78110 +5.7540 +.79890 +6.2450 +.81740 +6.6160 +.83190 +7.0510 +.85000 +7.4140 +.86630 +8.0340 +.88730 +8.4370 +.90280 +8.8240 +.91640 +9.2920 +.93200 +9.5900 +.94010 +10.1300 +.95500 +10.6630 +.97100 +11.1050 +.97920 +11.5320 +.98740 +11.9190 +.99410 +12.4750 +.99820 +2.6430 +.7fl420 +2.1440 +.69980 +1.5310 +.69750 +.9590 +.69580 +.3300 +.69380 +.63454703E -- 04

COFL +.7111 +.7239 +-7382 +.7478 +.7667 +.7812 +.7966 +.8165 +.8316 +.8491 +'8635 +.8873 +.9021 +.91~8 +.9315 +.9409 +.9567 +.9704 +.9803 +.9884 +.9942 +-9994 +.7032 +.6991 +.6954 +'6932 +.6919

DEL --.00127210 --.00156080 --.00309480 +.00084080 --.00323460 --.00011440 +.00221610 +.00~82830 +.00027810 +.00087400 +.00279250 --.000~4370 +.00~60490 +.00052580 +.00046620 --.00083780 --.00171520 +.00051100 --.00116930 --.00100430 --.00011720 --.00129960 +.00096180 +.00~67570 +.00204530 +.00257990 +.00182690

TABLE4. YZ Ca~iopeiae(MinII) Sample Computation ,50000

+ , 9 ;O0000t ,50 ~ , 1 4 1 8 9 2 6 0020 .00001 ,00000 3 , 4 8 4 , 9 3 5 7

TI +192.1600 +191.7000 +191-1360 +190.5560 +190.1290 +189.5970 +189.1540 +188.6380 +188.1380 +187.7030 +187.3160 +186.9300 +186:5100 +186-0750 +185.6880 +185.2770 +184.9630 +184.5680 +184.1570 +183.7220

POC +.78201 -~.72400 +.61834 +.48206 +.41356 +.30996 +-20312 +.14329 +.00653 --.02391 --.13498 ---23610 --.29149 --.44013 --.47181 --.54454 --.64348 --.71995 --'77878 --.83832

I~TW +.1072 +.1265 +.1598 +.2004 +.2205 +.2498 +.2791 +.2950 +.3297 +.3371 +'3625 +.3830 +'3929 +.4120 +.4143 +.4161 +.4fl86 +.3912 +'3678 +-3309

First Iteration

109

110

Machine Solutions of Light Curves ofEdipsing Binary Systems

c(I)

SE(I)

+.79030400E--02 ~-.12130253E--01 +.31774800E--02 R12=+.644 R T I O I = +.50000 +.10004584E+01 RATIO = -}-.65151 JA/JB=+.16189

+.66148755E--03 +.18676123E--03 -~.40441624_E--03 R13=+.827

R23=-~.156 RADB = +.1~51 THETE = +13.029

RADA = +.08~3 TDPR = +3.231

$2=+.25513750E--07 GAMR=+89.999 ANG=+89.999

Eighth Iteration C(I)

SE(I)

+.73380200E--02 +.12115433E--01 +.30690100E--02 R12 = +.722 RTIO1 = +.60621 RATIO = +.60567 JA/JB = +.18732

+.66818388E--03 +.20779252E--03 +.35931085E--03 R13 = ~-.815

TI

RADA----+.08566 TDPR = +3.175

OBFL

+192.160 +.9974 +191.700 +-9962 +191.136 +.9936 +190.556 +.9896 +190.129 +.9873 +189.597 +.9836 -~189.154 +.9795 +188.638 +.9771 +188.138 +.9714 +187.703 +.9701 +187.316 +.9653 +186.930 +.9609 +186.510 +.9585 +186.075 +.9522 +185.688 +.9509 -~185.277 4.9480 -}-184.963 ~-.9443 -~184.568 +.9417 +184.157 +.9399 -~183.722 +.9383 ~-.12618932E--04

COFL

R23 = -~.246 RADB = +.14142 THETE-----+13.120

82 = +.28203750E -- 07 GAM~ = +89.631 A.NG= +89.631

DEL

+.9979 +.9961 +.9934 +.9901 +.9873 +.9835 +.9802 +.9761 +.9719 +.9682 +.9648 +.9615 +-9578 +.9540 +.9508 +.9474 +'9450 +.9422 +.9395 +.9372

--.00057870 +.00000950 +.00014460 --.00050770 --.00004260 +.00000960 --.00074160 ~-.00097210 --.00056740 +.00185450 +.00041800 --.00060400 -~.00066610 --.00187810 +.00008520 +.00052370 --.00074910 --.00050130 -~.00033560 -~-00101500

Theoretwal Light Curve

shown in an abbreviated form in Tables 3 and 4 (on pp. 108-110). Changes in k with successive iterations are shown in Fig. 4. Note that the program given in this paper is applicable only to the secondary minimum. The procedure followed in computations is as follows. For complete eclipses the depth relation was solved for k over a range of assumed values of X L - - t h e coefficient of limb darkening of the larger component. The resulting values of k can be plotted against XL. This plot will be useful in subsequent discussion. The depth relation together with an estimate for k taken from k-

0ext - - 0int 0ex t -]- 0in t

I. JUI~EVIOH

111

.7,574 ".574 .534

| I |

~.s33

i

>-

.,t .53t

i/v

.6-

• 530

.65 .64

/11

.63 >.61

z o

I r

.60 .59

.4

.58 5O

.50 I START I

I 2

I i I I 3 4 5 6 ITERATION NO,

FIG. 4

I 7

i 8

I

0

I

I

I

I

I

I

I

I

J

.I .2 .3 . 4 . . 5 .6 .7 .8 .9 1.0 XL FIG..5

FIG. 4. Stabilization of k with successive iterations for YZ Cassiopeiae. FIG. 5. Variations of k from the depth and shape relations as functions of XL for YZ Cassiopeiae yields a reasonable starting value of k and determines the type of eclipse at hand. Taking now a starting value of k indicated by the depth relation, the complete occulation or transit program is entered as appropriate for the eclipse under study, and an improved value of k is obtained. I t can be noted t h a t the corresponding value of XL need not be exactly the same as t h a t indicated b y the depth relation. I t is of interest to repeat the main computation of k for various value of XL and then to plot k as a function of XL. This plot gives the values of k obtained from the shape relation. The functions k(XL) obtained from the depth and shape relations are expected to intersect at a point corresponding to the common values of XL and k. This interesection then represents the actual solution. Examples of this are shown in Fig. 5. An independent solution of the secondary minimum can be carried out in a similar manner. However, in this case, due to a number of reasons mentioned below, the solution m a y not be meaningful. I t is only too well known t h a t the secondary minimum is rarely adequately observed. Even in cases when there are a sufficient number of observations their precision is inferior to t h a t of observations in the primary minimum. I n fact, the shallower the secondary minimum the more pronounced becomes the difference in quality of observations in the two minima. This m a y lead to values of k widely different from those obtained from the primary minimum. Also recall t h a t the depth relation is a function of XL but not of Xs. This circumstance is fortunate because, in connection with the shape relation of the minimum corresponding to the eclipse of the larger component, it permits one to estimate the value of XL. However, no corresponding relations exist for the determination of Xs. If the secondary minimum has

112

Machine Solutions of Light Curves of Eclipsing Binary Systems

been adequately observed, one can formally accept k obtained from the primary minimum as the nearly correct one and search for those values of Xs t h a t will yield this value of k from the secondary minimum. However, there is no certainty t h a t this procedure will lead to reasonable values of X s. The coefficient of limb darkening of the small component obtained in this manner can either fall outside the allowable range 0 ~ Xs g 1 or be at variance with the spectral type of the smaller component. Indeed, the data used for the illustrative problem fail, for all possible values of X s , to yield a k consistent with t h a t obtained from the primary minimum. I n such cases one can always a t t e m p t to improve the solution by changing other parameters entering the problem, such as the intensities at the instant of internal tangency. The programs at hand permit such experimentation and, in fact, one can conduct a limited parametric study b y varying the appropriate parameters over a suitable range. However, such study falls short of an overall differential correction computation for which the present programs are not designed. As far as the practical application of the programs for complete eclipses is concerned, the following remarks are appropriate. To evaluate the fractional loss of light ~, it is necessary to have a value of the light intensity at the instant of internal tangency for the eclipse under consideration. I n the case of complete transit eclipses some care is needed in selecting this parameter. The reason for this is as follows. Recall t h a t for a given k the limiting value of p for a central eclipse is given b y Pc ~ -- l/k, for which, there is a definite limiting value of a. However, since both the light intensity at a given phase and at the phase of internal tangency are subject to observational and selection errors, the quantity 1--1 1--~ representing the observed fractional loss of light may, in the neighborhood of the minimum light, exceed the theoretical limit. Should this be the case, the procedure used in evaluating p corresponding to the observed value of intensity 1 will fail. The machine will signal this condition b y printing an error message. Most likely the error code encountered will be F 6 which corresponds to a negative argument in floating point square root or floating point exponent subroutines. This condition will occur in the a routines used in the present programs whenever an a t t e m p t is made to take p out of its allowable theoretical range. Therefore, it is recommended t h a t before the main program is used, one or two observations in the immediate neighborhood of the minimum be examined for the above condition. This is done easily by manipulating quantities MI, MM and F I T I T in the input to the complete transit program. An analogous situation exists in the complete occultation case. Here li and $ must never be such as to yield o¢ > 1. Difficulties associated with p being outside its allowable range can also arise in the neighborhood of the external tangency. I n this case it m a y be advisable to remove the offending observation from the computation. A similar problem can be encountered in computing a theoretical light curve. This computation employs the elements ra, k and i corresponding to a given iteration step, to produce theoretical values of p from r~(1 P=

-- sin~ i cos~ 0 ) 1 -- z cos ~ 0

1 k

I. JURKEVICH

113

The uncertainties in the computed elements m a y be sufficient to take p outside its permissible range or to produce a negative argument under the radical sign. If at all, this difficulty should occur only in the neighborhood of the extrema of a given minimum. All of the above difficulties are of the same nature and can be removed b y an adjustment of the intensity at the instant of internal tangcncy and b y removing troublesome points in the neighborhood of external tangency. Such adjustments are m a n d a t o r y with the present forms of the programs since there are no provisions for proceeding to the next observation if the indicated difficulties are encountered. I n the latter case, the machine, unless manually stopped, will print several F 6 and perhaps other error messages before proceeding either to a machine stop, a C H E C K STOP or R E A D E R NO F E E D condition. This situation can, of course, be avoided, provided additional storage capacity is available to include suitable diagnostic statements. There is little one can add to the above as far as the operation of the partial eclipses program is concerned. Aside from obvious changes such as one's inability to determine Xs and XL even in principle, it must be realized that, at least initially, the computation is an undisguised guessing game in which there are few, if any, helpful rules to go by. Success or failure of the computation depends almost completely on the convergence of successive iterations in k. However, even this does not insure t h a t the computations have been carried out correctly.

5. CONCLUSIONS AND THE ~)IEECTION OF FUTURE WORK

The power of automatic computers in the study of binary systems has so far not been sufficiently exploited. I t s appreciation, if any, is based mostly on viewing the automatic computer as an oversize slide rule or a desk calculator rather than a research tool which opens totally new possibilities. This state of affairs is not unlike t h a t which existed a few years ago in celestial mechanics and other disciplines. Such attitudes will change only with an increased application of machines to practical problems. At the present time, the practical experience in studying binaries b y means of automatic programmed machines is being accumulated sporadically. The work of the past few years resulted in a number of operational computer programs based, with very few exceptions, on the indirect methods. Exceptional cases, such as exist, are of the differential correction type. The existing programs are not fully automatic. An ab initio solution of a light curve is obtained in stages, with the computation being interrupted at the end of each stage. This is done to minimize the memory requirements and the program complexity, as well as to permit the study of intermediate results. Such a mode of operation is practically inevitable in the initial stages of the study of a system and, in fact, is not as detrimental as it m a y appear. Among the indirect methods there are m a n y t h a t are not worth mechanizing since they cannot utilize the most desirable features of the automatic computation. There are other existing methods whose practical use ought to be pursued as far as possible. The foremost example of the latter is the perturbation theory of Kopal which contains a very detailed dynamical account of binary systems. Yet, this theory has so far found little practical application. The future course of mechanization and effective application of the indirect approaches will inevitably be limited to those physical models for which analytical relations between

114

Machine Solutions of Light Curves of Eclipsing Binary Systems

the loss of light and system parameters can be not only found, but also inverted to permit isolation of some critical parameter. Theoretical studies of the past two decades, almost exclusively due to Kopal, indicate that attempts to obtain the indirect solutions of more complex systems hold little promise, at least for the foreseeable future. Such solutions would be particularly difficult to achieve for systems in which the gravitational dynamics may not be able to account for all observed complications. Thus, in recent years, it has become customary to attribute many unexplained effects in the observed light curves to gas streams surrounding the components, long-lived dark spots encompassing large fractions of stellar surfaces, and other ad hoc hypothesis which so far have not been subjected to any credible computational tests. The problem resides not so much in lack or difficulty of the required theory as with the magnitude of the computational work necessary. I t is this writer's conviction that such tests can be carried out only in terms of direct methods implemented by computer simulation techniques. The efforts of Wood, Mauder and ttor~k represent but specialized examples of this approach. Although detailed discussion of simulation techniques is outside the scope of this paper, a few general observations are in order. Computer simulation refers to numerical experimentation with a mathematical model of a physical system which is too complicated for analytical treatment. A particularly valuable feature of computer simulation is its ability to yield information on the relative importance of variables as well as on the complex interaction which may exist between them. Any model is an abstraction which replaces the physical system under study by a model of similar but simpler structure. Clearly, the mathematical model must contend with two conflicting requirements--realism and simplicity. The more realistic the model, the more complex is its structure and the more difficult it becomes to understand and operate with. On the other hand, simple models are seldom realistic. In the field of binary stars the available simple models leave too many observed features unaccounted for, whereas the more realistic models have not really been formulated. The lack of a physical theory, however, should not preclude computer simulation since the unknown relations describing the interaction between certain variables can be replaced, for the time being, by formal mathematical constructs. Subsequent experimentation with the model will reveal whether the assumed construct affects the overall model adversely and whether suitable changes can be made to improve its operation. Such studies may, in fact, lead to a discovery of the underlying physical connections between the variables under consideration. The above approach need not excessively disturb the binary star investigator since the use of purely mathematical models is not a totally unfamiliar concept to him. To illustrate, consider the use of limb darkening laws for stars of early spectral type. As the model atmospheres of these stars were developed, it became clear that the limb darkening law was distinctly non linear. Unfortunately, the model atmosphere work does not yield an explicit analytical expression for this quantity. Subsequently, various investigators represented the law by second and third degree polynomials. However, such expressions are merely formal mathematical constructs adopted for convenience and having, as yet, little basis in physical theory. This is but one example of the use of formal mathematical models in binary star astronomy. I t must be emphasized that computer simulation has its limitations. First, one may fail to find a mathematical model which yields a satisfactory approximation to the physical model. Second, a satisfactory model may yield a solution at an exorbitant price. The latter is closely connected with the complexity of the model. A complex model involving many variables may impose unreasonable requirements on the computer memory and lead to excessive computing time. I t is highly important to formulate a model which accurately describes the behavior of the system under study, but at the same time minimizes the com-

I. JURK~.VIC]t

115

puting and programming time. Unfortunately, these factors are totally interdependent. The question of computing efficiency--that is the use of efficient numerical and analytical techniques to perform the requisite c o m p u t a t i o n s ~ i s vital in minimizing the computing time. I n treating an observed fight curve b y simulation techniques one searches for a set of parameters which best represent the observations. This, of necessity, leads us into area of statistics concerned either with the iterative variance reduction techniques or the exploration of the so-called response surfaces via regression analysis. I n both cases the amount of computation required is very large, particularly if a n y of the variables are statistical in nature. For example, in simulating a binary system embedded in a gas envelope, the velocity of injection of gas particles into the circumstellar space must be considered a random variable. I t will also be necessary to consider the question of the validity of the model. A great deal of uncertainty will exist in this area, particularly in those cases for which the output variable (e.g. the loss of fight) is an insensitive function of some input variables. I n such cases, the solution m a y not be unique. Similarly, not all variables are equally significant, b u t at the outset of the computation their relative significance m a y not be known. The variables excluded from consideration at any particular stage of the study could then play the role of the so-called " l u r k i n g " variables about which one m a y know very little. The cumulative effect of such variables could be very significant and yet very difficult to estimate. I n conclusion, recall t h a t the amount of realism built into the model must be adequate to describe not only the present state of the system, but its past and future behavior as well. I f these objectives are not achieved, the model must be considered a failure. Finally, the investigator must make every effort to guard against making actual simulation runs with an inadequate model.

REFERENCES BLA.CK~_A~,R. B. and TUCKI~Y,J. W. (1959) The Measurement o/Power Spectra, Dover Publications, New York. GRYGAI~,J. (1963) Applieation of the non-linear laws of limb darkening on eclipsing binaries, Bull. Astr. Inst. Czech. (BAC) 14, No. 4,127. F~.TL~R, J. (1923) A contribution to the theory of eclipsing binaries, Reef. Astr. Obs. Utrecht 9, 1. HA~COCK,H. (1958) Elliptic Integrals, Dover Publications, New York. H~u~m, A. L. (1966) Photometric Studies o/the Eclipsing Binary Systems B V 412, B Y 332, and V 836 Cygni. A dissertation presented to the University of Pennsylvania (Philadelphia, Pa.). HASTinGS, C., Jr. (1955) Approximations/or Digital Computers, Princeton University Press, Princeton. HOR~K, T. 1966, New elements of eclipsing variable UX Herculis, B.A.C. 177 (1) 27. HUFF~R, C. M. and COLLI~S, G.W. (1963) Computation of elements of eclipsing binary stars by high speed computing machines, Ap. J. Suppl. Set. 77 No. 71. JURKEWCH,I. (1964) Machine Reduction, Computation and Analysis o] Light Curves o/Eclipsing Binary Systems~ Technical Information Series Report R 64 SD 8, Space Sciences Laboratory, General Electric Co., Valley Forge STC, Pa. KI~(~, E. S. (1920) Harvard Ann. 81, 231. KING, L. V. (1924) On the Direct Numerical Calculation o] Elliptic Functions and Integrals, Cambridge University Press, Cambridge. KITXMU~A,M. (1965) Determination of the elements of eclil3singvariables from Fourier transforms of their light curves, Advances in Astronomy and Astrophysics, Vol. 3, Academic Press, New York and London.

116

Machine Solutions of Light Curves of Eclipsing Binary Systems

KOPAL, Z. (1959) Close Binary Systems, John Wiley, New York. KRAT, V. (1934) Astr. Zh. 11, 407. KRAT, V. (1935) Astr. Zh. 12, 21. KRAT, V. (1936) Astr. Zh. 18, 521. KRO~, G. E. (1939) Lick Obs. Bull. No. 499, pp. 59-71. LII~NELL, A. P. (1965) The calculation of direct eclipse functions, L A p . J. Suppl. Series, Vol. XI (99) 185; II. Ap. J. Suppl. Series, XII, (109) 263. LI~NELL, A. P. (1966) III-V. Ap. J. Suppl. Series, Vol. XIII, (120) 413. MXUDER, H. (1962) Bahnbestimmung yon Bedeckungsver~nderlichen mit elektronischen Rechenmaschinen, 2. Ver~inderlichen-Colloquium, Kleine Ver6~entlichung.en der Remeis-Sternwarte Bamberg Nr. 34, 30. M),UDER, H. (1966) Die Berechnung der Elemente yon Bedeckungsveriinderlichen durch Fouriertransformation der Gesamtlichtkurve, Kleine Ver6~entlichungen der Remeis-Sternwarte Ba~berg Nr. 38. Tables of required functions appear in the same series as publication Nr. 39. MERRILL,J. E. (1950) Tables for solution of light curves of eclipsing binaries. Coefficient of limb-darkening x ~ 0.0, Contr. Princeton Univ. Obs., No. 23. MERRILL, J. E. (1953) Limb-darkening coefficient x =: 1.0. Contr. Princeton Univ. Obs., No. 23. PIOTROWSKI,S. L. (1947) Ap. J. 106~ 472. PmTROWS~L S. L. (1937) Acta Astr. (a) 4~ 1. I:~USSELL,H. ~T. (1912) Ap. J. 85, 340. RUSSELL, H. N. (1912) Ap. J. 86, 54. RUSSELL, H. N. and MERRILL, J. E. (1952) The determination of the elements of eclipsing binaries, Contrib. Princeton Univ. Ob8., No. 26. RUSSELL, H. N. and S~PLEY, H. (1912) Ap. J. 86~ 239. RUSSELL, H. 5[. and SHXeLEY, H. (1912) Ap. J. 86, 385. ScnL6~[CH, O. (1879) Compendium der HSheren Analysis, Zweiter Band, Friedrich Vieweg, Braunschweig. SCH~ELLER, H. (1949) Eine einfache Methede zur Bestimmung der Systemkonstanten bei Bedeckungsver~nderlichen, Ver6~entlichungen der Sternwarte zu Sonneberg, Bd. 1, 4. SHARnE, S. B. (1924) Izv. Glav. Astr. Obs. 10~ No. 1. TSESEVICH, V. P. (1939) Tablitsy Dlia Opredelenia Elementov Orbit Zatmennykh Tipa Algolia, BiuL Astr. Inst., No. 45, Leningrad. TSES~WCH, V. P. (1940) Tablitsy Fotometricheskikh Faz Zatmenii Peremennykh Zviozd Tipa Algolia, Biul. Astr. Inst., No. 50. TSESEVIC]~, V.P. (1947) Metody Izuchenia Peremennykh Zviod, Vol. 3, Gosudarstvennoe Izdatel'stvo Tekhniko-Teoreticheskoi Literatury (OGIZ), Moscow-Leningrad. WELL~X~, P. (1962) Solution of light curves using the IBM 650, Prec. X I General Assembly o/the I A U~ Berkeley, 1961, edited by D. H. SXDL~R, Academic Press, London and l~ew York. WILSOn, E. B. (1912) Advanced Calculus, Ginn, Boston. WrLso)r, R. E. (1966) University of South Florida, Tampa, private communication. WOOD, D. B. (1962) Machine method for calculation of orbital parameters, Prec. X I General Assembly o/the I A U , p. 368, Berkeley, 1961, Academic Press, London and New York.