非接觸型與半接觸型食雙星光變曲線合成與分析

全文

(1)國立台灣師範大學地球科學研究所天文組碩士論文. 非接觸型與半接觸型食雙星光變曲線合成與分析 Light Curve Synthesis and Analysis of Semi-Detached and Detached Eclipsing Binaries. 研究生 : 羅卓延駿 Yen-Chun Luo Cho 指導教授 : 傅學海教授 Prof. Hsieh-Hai Fu. 中華民國一○五年八月.

(2) 中文摘要食雙星的光變曲線合成可以幫助我們得到雙星系統的物理參數包括了恆星質量、半徑等資訊，有助於瞭解雙星系統的結構及演化過程。針對非接觸型雙星系統，我們模型使用簡單球體模型搭配三角幾何運算並加入臨邊昏暗效應的修正來合成光變曲線；而對於半接觸型雙星系統，則基於 Kopal 的方法計算系統的洛希重力等位面並以線性關係式對臨邊昏暗效應做修正。我們選擇 XY Cet 和 DD Mon 這兩個分別是非接觸與半接觸型雙星系統來做模型的測試與驗證，它們的星等觀測資料均來自 All Sky Automated Survey (ASAS)，物理參數則參考 Gazeas、Southworth 等人的研究。從光變曲線合成的結果來看，我們使用的方法的確可以成功建立食雙星的光變曲線模型，同時也可以藉由對觀測資料的擬合得出系統的物理參數。我們也分析了 CQ Ser 這個食雙星系統的光變曲線，它的星等觀測資料部份來自於 ASAS，另一部份則使用我們自己以台師大 RC16 望遠鏡觀測所得到的資料，而擬合過程所需的模型初始參數則是直接從光變曲線估計得出。擬合的結果成功給出了 CQ Ser 的物理參數同時也顯示了這個雙星系統的狀態可能是兩個成員星都幾乎脹滿了洛希瓣的半接觸型雙星系統。. 關鍵字:. 食雙星, 雙星, 光變曲線合成, 光變曲線擬合, Kopal 模型,. 洛希模型, XY Cet, DD Mon, CQ Ser. i.

(3) Abstract Light curve synthesis of eclipsing binaries can help us build a physical model of the binary systems and get the physical parameters of them. For detached binaries, we use the simple spherical model with the concepts of trigonometry and linear limb-darkening law to synthesize the light curve; on the other hand, for semi-detached systems, the model we used is base on Kopal’s methods that assuming the systems have circular orbits and using Roche model with considering the limb-darkening effect. We synthesized and analyzed the light curve of XY Cet and DD Mon which are detached and semi-detached system respectively with the photometric data taken from All Sky Automated Survey (ASAS) and the parameters used in our model were obtained from Gazeas (2010) and Southworth (2011). The results of synthetic light curves show that the methods we used have ability to construct a physical model of eclipsing binaries with certain reliability and could derive the system parameters by light curve fitting process. Then we analyzed the light curve of CQ Ser using photometric data combining with ASAS and our observation which observed with RC16 telescope at FRO of NTNU. The initinal parameters used in the fitting process were estimated from the light curve directly. The result of light curve fitting shows some system parameters and indicates that the system may be a semi-detached system with the both components almost filling the critical Roche lobe.. Keyword:. eclipsing binaries, light curve synthesis, light curve fitting,. binaries, Kopal model, Roche model, XY Cet, DD Mon, CQ Ser ii.

(4) Contents 中文摘要 ................................................................................................... i Abstract..................................................................................................... ii Contents ................................................................................................... iii List of Tables............................................................................................. v List of Figures .......................................................................................... vi. I Introduction............................................................................. 1 II The Model of Light curve Synthesis ..................................... 4 2.1 Simple Spherical Model for Detached System................................4 2.2 Roche Equipotential Model for Semi-Detached System .................7. III Observation and Photometry Data ................................... 15 3.1 Introduction of three targets: XY Cet, DD Mon and CQ Ser.........15 3.2 Online Photometric data .............................................................18 3.3 Observation of CQ Ser .................................................................20. IV Light curve synthesis and analysis .................................... 23 4.1 The result of light curve synthesis ...............................................23 4.2 Light curve analysis......................................................................25. V Discussion and Conclusions................................................. 32 References ................................................................................ 36 Appendixes ............................................................................... 38. iii.

(5) A. The PyRAF code of images data reduction ....................................38 B. The source code of light curve synthesis-XY Cet...........................42 C. The source code of light curve synthesis-DD Mon.........................45 D. The source code of CQ Ser light curve fitting ................................52 E. Table of our photometric data of CQ Ser ........................................64. iv.

(6) List of Tables Table 1. Photometric Data of Three Targets from ASAS (2001~2009) .... 18 Table 2. Observing Log of CQ Ser .......................................................... 21 Table 3. The parameters using in our model of XY Cet and DD Mon...... 23 Table 4. The result of light curve analysis of XY Cet and DD Mon ......... 27 Table 5. The result of light curve analysis of CQ Ser............................... 31. v.

(7) List of Figures Fig. 1.1 The morphological classification of eclipsing binaries. ................ 2 Fig. 2.1 The detached systems geometry in center of mass coordinates. .... 6 Fig. 2.2 The summation of the arcuate region of detached systems. .......... 6 Fig. 2.3 The semi-detached systems in Cartesian coordinates.................... 7 Fig. 2.4 Surface elements distortion are increased with angle β............... 10 Fig. 2.5 The coordinate of surface elements was projected to plane of sky. ................................................................................................................ 13 Fig. 2.6 The initial boundary elements for visible region......................... 14 Fig. 3.1 ASAS Observation light curves of XY Cet, DD Mon and CQ Ser. ................................................................................................................ 19 Fig. 3.2 Image of our observed field including CQ Ser and reference stars. ................................................................................................................ 21 Fig. 3.3 Differential magnitude light curve of CQ Ser and reference stars. . . ................................................................................................................ 22 Fig. 3.4 Our Observation light curve of CQ Ser..................................... . 22 Fig. 4.1 The synthetic light curves of XY Cet and DD Mon. ................... 24 Fig. 4.2 The fitting light curve of XY Cet and DD Mon. ......................... 28 Fig. 4.3 The fitting light curve of CQ Ser. ............................................... 31 Fig. 5.1 The configuration of XY Cet and DD Mon at phase 0.75. .......... 34 Fig. 5.2 The configuration of CQ Ser at phase 0.75. ................................ 35 Fig. 5.3 The configuration of CQ Ser in 3D view at phase 0.16. ............ 35. vi.

(8) I. Introduction We try to make a computer program (code source seen on appendixes. A, B, C, and D) of light curve synthesis of eclipsing binaries using the physical model with some physical parameters found from the observed light curve. Except for stellar masses, the information, orbital period; orbital eccentricity; orbital inclination; the radii of both components in units of their separation radii; the relation between individual surface temperatures of both components and even the stellar spots could be obtain from the light curve directly. Those parameters are necessary for studing of stellar structure and evolution. In this chapter, we will introduce the classification of eclipsing binaries; outline the history of modeling light curve of eclipsing binaries and the research purposes of this paper. The classical classification of eclipsing binaries is corresponding to the shape of light curves such as Algol, β Lyrae, and W UMa light curves and also known as EA, EB, and EW type, respectively. The EA light curves are typically almost flat topped with a large difference between the depths of the two minima and there may be an increase in light near the phase of secondary minimum due to the reflection effect. The EB light curves are continuously variable because of tidally distorted components and the large difference in depths of minima indicating components with quite different surface brightness. The EW light curves are also continuously variable but with only a small difference in the depths of the minima. The variation outside the eclipse indicates that EW systems are over-contact. The morphological classification of eclipsing binaries is according to the shape of Roche surface of the system. The equipotential surfaces also 1.

(9) called Roche surfaces that the sum of gravitational and rotational energy per unit mass is constant on these level surfaces. The Roche surface passing through Lagrangian point L1 consists of two ovoid surfaces called the Roche lobes of the components. The two ovoids touch at L1 and have three situations according to the shape of components, as shown in Fig. 1.1 , which called detached, semi-detached and over-contact systems. The detached systems are neither component fills its Roche lobe; semi-detached systems are one component fills its Roche lobe and the other does not; and over-contact systems are both components exceed their Roche lobes. Roughly specking, the EA type light curves may be produced by detached systems and the EB type light curves may be produced by semi-detached systems.. Fig 1.1 The morphological classification of eclipsing binaries. Picture from Bradley W. Carroll, Dale A. Ostlie . An introduction to modern astrophysics (2nd ed.) (pp. 659).. 2.

(10) Russell (1912) provided the first light curve synthesis methods which assuming the stars are spherical with a circular orbit and also considering the limb-darkening effect. Russell & Merrill (1952) modified the Russell model by the correction of ellipticity of stars and adding the correction of gravity darkening and reflection effect. Nelson & Davis (1972) presented the model which close to Russell & Merrill (1952) but using the more efficiency method of calculating the limb-darkening effect and improving the accuracy of the results. They also introduced the eccentric orbits into the Russell model first. Generally, this model is quite suitable analysis of detached binary systems with minimal shape distortion. The first using Roche geometry to synthesize light curve was presented by Kopal (1959) assuming the stars have the circular orbit around the center of mass and considering the limb-darkening effect. Wilson & Devinney (1971) introduced the eccentric orbits with asynchronous rotation and some other effects like gravity darkening and reflection effect. This model also can work on over-contact systems. The W-D model is the most widely used of the light curve analysis method today. In this paper, we focus on building our own light curve synthesis program for detach and semi-detach systems base on the above-mentioned methods. We choose two eclipsing binaries systems as test targets and put their physical parameters which come from literature into the program to confirm it have considerable credibility or not. Finally, we using our program analyze the light curve of selected binaries to get the light curve solution of the systems.. 3.

(11) II Light curve Synthesis Method 2.1 Simple Spherical Model for Detached System In this model we assume the shape of stars is spherical and the stars have circular orbits with synchronous rotation. The system geometry was illustrated in Fig. 2.1 and the stars coordinate can be estiblished by following equations: X = r sinθ Y = r cosθ cos i Z = r cosθ sin i. x1= - X / [1+(1/q)] ,x2 = X / (1+q) y1= - Y / [1+(1/q)] ,y2 = Y / (1+q) z1 = - Z / [1+(1/q)] ,z2 = Z / (1+q). (2.1.1). where r is distance between the stars centers, i is orbital inclination, q is mass ratio and θ is phase angle. We consider the limb-darkening effect and assume the emergent intensity, I0, which emitted from stellar dick center as the blackbody radiation that is described by using Stefan–Boltzmann law, so the flux emitted from each stellar dick element could be written as: 4. I   Te (1  x  x cos  ). (2.1.2). where x is linear limb-darkening coefficient, γ is an angle between the radius vector and the line-of-sight vector and could be written as:.   sin 1 (. dR ) R. (2.1.3). where R is radius of stellar dick and dR is a distance changing from disk center to disk edge. As shown in Fig. 2.2, we calculate the area of arcuate region (show in blue and pink) which was considered having the same limb-darkening effect in a specific thickness dr and radius: dr1-1/2*dr, then the area multiply (2.1.2) to get the specific flux. Summation the flux in 4.

(12) various dr1 which from 0 to R1 to determine how much flux was blocked by star2 when eclipse is happening. Using Cosine rule to get the angle θa and θb then the flux from star1 blocked by star2 could be written as: R1. Flux1blocked .  I . a. ( dr12  ( dr1  dr ) 2 )   (dr1  dr )( b   a ) dr . dr10. (2.1.4) , so the total flux emmited from system when star1 was blocked by star2 is given by:. Flux  Flux1  Flux2  Flux1blocked. (2.1.5). ,where the Flux1 and Flux2 are the flux emitted from star1 and star2 respectively when eclipse is not occurring.. 5.

(13) Fig. 2.1 The detached system geometry in center of mass coordinates. θ is phase angle and i is inclination angle.. Fig. 2.2 The summation of the arcuate region (blue and pink region) of detached systems with the various dr1 and a specific thickness dr determine the area was blocked by star2 when eclipse is happening. Using Cosine rule to find the θa and θb . 6.

(14) 2.2 Roche Equipotential Model for Semi-Detached System 2.2.1 System Geometry and Coordinates The main concepts are calculating the Roche equipotential surface (Ω) in spherical coordinate determines the shape of stars and assuming the system have circular orbits with synchronous rotation which based on Kopal (1959) . Fig. 2.3 illustrate the geometry of the system in Cartesian coordinate. M1 and M2 are masses of individual stars; A is the distance between the two; m is the small element of unit mass in the atmosphere of M1; CM is the center of mass; d1 and d2 are the distances of m from M1 and M2 respectively. y-axis m d1 M1. d2 M2. CM A. x-axis. Fig. 2.3 The semi-detached systems in Cartesian coordinates. The origin is in the M1. The shape of both stars can be determined by the forces on the m within the system:. FG  FP  m 2 ~ r. (2.2.1). where ω is the angular velocity of the system, ~ r is the vector from CM to m, FG is gravitational force and FP is the pressure force may be written as :. FP  . 7. m p . (2.2.2).

(15) where ρ is the mass density of volume. The gradient of a potential Ψ can be expressed as:. FG  m 2 ~ r . GM 1m ~ GM 2 m ~ d1  d 2  m 2 ~ r   m 2 2 d1 d2 p   . (2.2.3) (2.2.4). Since p and  are parallel and the curl of above equation:.   p    (  )      -     0. (2.2.5). shows that.     0. (2.2.6). , this suggests that surfaces of potential, pressure and density coincide. In other words, we could find the possible shapes of the stars by calculating the surface potential Ψ:. . GM 1 GM 2 1 2 2    r d1 d2 2. (2.2.7). The d1, d2, r and ω can be expressed in Cartesian coordinates of stars surface:. d1 . x2  y2  z2. d2 . d 1  2 xA  A 2. r.  . 2. 2. x CM  2 x CM x  x 2  y 2 G M 1  M A3. 2. (2.2.8). . where xCM is position of CM lies on the x-axis (q = M2 / M1):. xCM . q A 1 q. (2.2.9). The x, y, z coordinate can be transformed to spherical coordinate system.: 8.

(16) x  d 1 sin cos  d 1  y  d 1 sin sin  d 1 . (2.2.10). z  d 1 cos  d 1 where θ = altitude angle, 0 ~ π, from +z to -z φ = azimuth angle, 0 ~ 2π, from x-axis, CCW Now the surface potential Ψ can be rewritten as: 2 GM1  A  d1  d1 A q2  2   q 1 1    (r,, )   q   A d1  d12  2d1 A  A2 A  2A2 21  q    . (2.2.11) According to a modified Roche potential, Ω, defined by Kopal (1959):. q2 A    GM1 2 1  q. (2.2.12). and introducing a variable R = d1 /A , the Ω can be expand to :  .  1 1  1 2 2      q   R  q  1 1   R  R  1  2R  R 2  2. (2.2.13). It should be note that the q term in Ω with respect to M2 star is using 1 / q. We could find the R for each λ and υ in given Ωpole. The Ωpole is the potential at the pole of M1 star:.  pole. so,. R.   1 1    q 2 R pole  1  R pole    1 q. 1 2 2  pole  qR    q  1   1    R 2 1  2R  R 2. (2.2.14). (2.2.15). Finally, solving this equation iteratively for each λ, υ and using Rpole as an initial guess for R determines the surface (also known shape) of star. 9.

(17) 2.2.2 Surface Elements and Radiative Properties Because of the distortion of the stellar surface depends on the value of Roche equipotential, Ω, the area of each surface elements will be changed by an angle β which between surface element normal, n~ , and unity radius vector ~ ri , as shown in Fig. 2.4 The surface element normal may be written as Linnell (1984) presented:.  n~   . (2.2.16). and. cos   n~  ~ ri  n x   n y   n z . (2.2.17). ,so the area of differential surface element could be written as:. d . 1 R 2 sin dd cos . (2.2.18). z. ~ Surface element normal, n β Unity radius vector, ~ ri γ. y. Line-of-sight vector, ~ s. Fig. 2.4 Surface elements (eg. the red region) distortion are increased with angle β. The unity radius vector is a function of λ, μ, ν which be described in previous section. The line of sight vector parallel to y-axis. We assume each surface element has the same emergent intensity, I0, which emitted from star center as the blackbody radiation that is described 10.

(18) by the Planck function. The emergent intensity was assumed by using Stefan–Boltzmann law:.  Te I0  . 4. (2.2.19). where σ is Stefan–Boltzmann constant, so the flux could be defined by the I0 on all wavelengths emitted from ‘effective’ unity stellar surface area in unity time interval, the effective means that need to considering the line of sight vector, and flux could be written as:. f 0   cos I 0. (2.2.20). where γ is an angle between the unity radius vector and the line-of-sight vector. The I0 was done the limb-darkening correction by using linear limb-darkening law, getting the local intensity as:. I  I 0 (1  x  x cos  ). (2.2.21). where x is linear limb-darkening coefficient, and the f0 becomes to:. f 0   cos I. (2.2.22). Now we could get the local flux, f, which defined by local intensity emitted from the surface element:. f   f 0  f 0 d. (2.2.23). Finally, the total flux of ‘visible’ surface region on the one of components is computed by the following expression: 2. F  0. . T 0. e. 4. (1  x  x cos  ). cos  2 R sin dd cos . (2.2.24). and the total flux of system is given by summation of flux of the two components. 11.

(19) 2.2.3. The Eclipse Effect and Visible Surface Elements In order to determine which surface elements could be visible by observer, we transform the stellar surface coordinates which are generated by Ω to the plane of sky coordinates with simple rotation for the given phase and inclination. The origin of the coordinate of surface elements is reset to the center of mass at first, then the coordinate of surface elements were rotated around x-axis with an angle ic which is a complementary angle of inclination in CCW, then rotated around new z-axis with a phase angle θ in CCW. The rotation matrix is show below:.  x  cos   y    sin      z   0. - sin cos  0. 0  1 0 0 0 cos ic  1  0 sin ic.   xori   sin ic   yori  (2.2.25)   cos ic   zori  0. The x and z values in [ x y z ]T are the coordinate of surface elements in x-z plane which are the original coordinate of surface elements were projected to plane of sky, as shown in Fig. 2.5. Now we have the coordinates of surface element in specific phase and inclination, so it’s time to verify which elements could be visible. We set the line of sight vector is tower to +y-axis and the δ is an angle between line of sight vector and surface element normal, as shown in Fig. 2.4. The surface elements which have cos δ＜0 (δ >90°), means that could be visible by observer and just dependent on eclipsing or not. The eclipse occurs when the distance between the centers of components is less than sum of maximum radius of components. To test which star is front, we defined star2 is front when phase≦0.25 and phase≧0.75, otherwise the star1 is front. We choose the elements which have -0.1＜cos δ＜0 as the 12.

(20) initial boundary elements of stars when eclipses occur. Then put its coordinates into the special algorithm to determine which surface elements are not eclipsed by the other star. The main concept of this algorithm is that the model evaluates the distance, dtest, between front star’s center, m1, and back star’s coordinates of surface element, b1, then compare with the distance, db, between m1 and c1. c1 is a cross point between the line of m1 to b1 and front star’s boundary element, see Fig. 2.6. If the dtest＞db , then the surface element of b1 is visible, the result of this special algorithm was shown in Fig. 2.5.. Fig. 2.5 The coordinate of surface elements were projected to plane of sky, x-z plane. This figure also show the result of visible elements determined from the special algorithm when eclipse occurs.. 13.

(21) m1. db. c1 dtest. ●. b1. Fig. 2.6 The initial boundary elements for visible region shown in blue and red dot. The green dots represent the boundary elements of front star which inside the region of back star.. 14.

(22) III Observation and Photometry Data Three targets, XY Cet, DD Mon and CQ Ser, are chosen to correspond to detached and semi-detached system respectively. The photometric data of XY Cet and DD Mon with related physical parameters from literatures will be used to verify the feasibility of our model because that has been studied on many occasions and the physical properties were firmly established. However, CQ Ser is suitable as the target because there are fewer related literatures and almost not yet have in-depth discussion.. 3.1 Introduction of XY Cet, DD Mon and CQ Ser Morrison & Morrison (1968) did the first photometric study of XY Cet. They found the orbital period was 2.780712 day and gave the preliminary result for the elements of the system using the method of Russell & Merrill (1952), finding that the system was detached system and two stars were very similar. Popper (1971) did the first spectroscopic study of XY Cet. He combined his spectroscopic results with the photometric elements of Morrison & Morrison (1968) to make the first determination of the masses and radii of the component stars of XY Cet. He also gave a spectral classification of A2m + F0m. Srivastava (1987) analyzed the photometric data from Srivastava & Padalia (1975) and used the spectroscopic elements given by Popper (1971). He obtained result for the elements of the system using the method of Kopal (1952), finding that values of radii which were quite different from Popper (1971). Southworth J. et al. (2011) presented a new radial velocity result and modeled the WASP light curve using the JKTEBOP code which based on Nelson & Davis (1972), giving the comprehensive physical properties of the XY Cet 15.

(23) system to accuracies of 1 - 2 percent. The result basically agrees with Srivastava (1987) and also indicated that both stars show enhanced abundances of iron-group elements confirming their classification as Am stars and XY Cet is a prime candidate for detailed spectroscopic analyses of metallic-lined stars. The object DD Mon was done the first photoelectric photometric with BV filter and the spectroscopic observation by Yamasaki et al. (1990). They assigned the spectral type of the primary component as F5V and indicated that it is a detached near contact binary system with the primary component almost filling the critical Roche lobe determined from their own light curve synthesis model. Qian et al. (1997) also obtained BV light curves, but they suggested that the system is a semi-detached near-contact binary with the secondary component filling the critical Roche lobe, that was slightly difference between previous result. They got the photometric solution based on spectral type F5V and Wilson-Devinney program showing that DD Mon was a post-mass-transfer semi-detached system. Qian et al. (2009) analyzed the variations of the orbital period, finding the period increase at a rate of 1.40*10-7 day/yr. They said the reason could be a mass transfer from the less-massive component to the more-massive one. Recently, the more accurate radial velocities of DD Mon were obtained by Pribulla et al. (2009) which show that the value 133.32±3.06(km/s) and 198.85±3.18(km/s) correspond to primary and secondary component respectively. Gazeas et al.,(2010) analyzed the light curve of DD Mon and combined with radial velocities from Pribulla et al. (2009), giving the absolute elements of DD Mon by using Wilson-Devinney program. Brancewicz & Dworak (1980) did the first light curve analysis of CQ 16.

(24) Ser. They used the simplified model based on Kopal (1952) and assumed that the components of system behave as non-rotating stars and also ignored the effect of limb darkening and reflection. They obtained the variety of physical parameters based on spectral type A8 but intrinsic accuracy is questionable. Shaw (1994) referred clearly that the system was near-contact system but he did not give any physical parameters of the system. So far there is no radial velocity analysis result of CQ Ser.. 17.

(25) 3.2 Online Photometric Data The photometric data of these three targets we were taken from All Sky Automated Survey (ASAS) online database (Pojmanski, G. 2002) from 2001 to 2009 and taken our observations only for CQ Ser that will be described in the next section. The ASAS systems is located at Las Campanas Observatory and equipped with 200mm/f2.8 Minolta telephoto lenses and 2K x 2K CCD camera, covers 8.8 x 8.8 degree of the sky through a standard V-band filter. ASAS data have five different aperture magnitudes for each individual observation, and the aperture we used according to which aperture have the smallest magnitude dispersion. The detail of ASAS photometric data we used are given in Table 1. Because of the models we mentioned in Chapter 2 are calculating the normalized flux of the system, the photometric data entering the models must be the same as the flux. Therefore the differential magnitudes with respect to the reference magnitudes, see Table 1, were converted to normalized flux as. flux  10. 0.4 ( m  m0 ). where Δm is differential. magnitudes and Δm0 is the median of Δm at the quarter phase. The light curves of each target were shown in Fig. 3.1 Table 1. Photometric Data of Three Targets from ASAS (2001~2009) Points. Reference mag.. Reference Epoch(HJD). Period(day). XY Cet. 345. 8.713. 2438372.9455*. 2.78071392*. DD Mon. 419. 10.747. 2450100.2938**. 0.56801844**. CQ Ser. 674. 10.872. 2456075.8392***. 0.760060***. Reference: *:Southworth(2011); **:Qian(2009); ***:IBVS 6029. 18.

(26) Fig. 3.1 ASAS Observation light curves of three targets. Upper panel: XY Cet; middle panel: DD Mon; bottom panel: CQ Ser. The photometric data all tacken from ASAS (Pojmanski, G. 2002). 19.

(27) 3.3 The observation of CQ Ser The photometric observation of CQ Ser was carried out on May 2016 and July 2016 at the FRO located on campus of NTNU in Taipei City. Using RC16 telescope and ST-10XME CCD camera with Johnson V filter got 469 effective images with the exposure time of 70 second in each frame. The ephemeris used for phase calculations and the log of photometric observations are given in Table 2. The dark and flat-field correction of images were done by using MaxIm DL 5, then the images were processed with the XYXYMATCH and the DAOPHOT package in IRAF to perform aperture photometry and get the instrumental magnitude of stars. The aperture for photometry was set to 1.69 times the FWHM of the stars in all images. The photometric errors for individual observations were 0.002 – 0.04 mag. Differential magnitudes were determined and four stars were chosen as the reference stars, the Fig. 3.2 shows its relative position of CQ Ser. The differential magnitudes of all measured stars were calculated with respect to the average magnitude of four reference stars, and its light curve was shown in Fig. 3.3. The phase in that figure was set the reference of minima epoch to HJD 2456075.8392 and period 0.760060 day from IBVS 6029. Because of the model demand, we transform the differential magnitudes into normalization of flux by using formula which shows in Section 3.2. The combination of our observations and ASAS photometric data will become the input data of the model. The flux light curves of our observations combine with ASAS data are shown in Fig. 3.4.. 20.

(28) Fig. 3.2 Image of our observed field including CQ Ser ( Ra = 18h08m53s , Dec = -14°15’55’’ ) and other reference stars. The FOV is 15’x10’. Table 2. Observing Log of CQ Ser Date. Start Time. End Time. Phase*. Frames. (HJD-2400000). (HJD-2400000). 2016.05.19. 57527.6935. 57527.7595. 0.84 - 0.92. 10. 2016.05.29. 57537.6793. 57537.7753. 0.98 - 0.10. 46. 2016.06.03. 57542.6170. 57542.8511. 0.47 - 0.78. 73. 2016.06.04. 57543.5951. 57543.8571. 0.76 - 0.18. 92. 2016.06.15. 57555.1244. 57555.3372. 0.27 - 0.55. 115. 2016.06.18. 57558.1115. 57558.3014. 0.20 - 0.45. 83. 2016.07.16. 57586.1936. 57586.2113. 0.15 - 0.17. 9. 2016.07.21. 57591.0868. 57591.1796. 0.59 - 0.71. 41. (V-filter). *The phase was calculated by setting minima epoch to HJD 2456075.8392 and period 0.76006 day from IBVS 6029.. 21.

(29) Fig. 3.3 Differential magnitude light curves of CQ Ser (right panel) and other reference stars (left panel).. Fig. 3.4 Our Observation light curve of CQ Ser which shown in red dots The blue plus markers indicate the data from ASAS. 22.

(30) IV Light curve synthesis and analysis. 4.1 The result of light curve synthesis As chapter III mention that the parameters of XY Cet and DD Mon will be used to verify the feasibility of our model. The parameters we used were chosen from Southworth J. et al.(2011) and Gazeas et al.,(2010) respectively because they both done the radial velocity analysis and got the absolute physical parameters of the systems. The details of parameters that we input to our model of two targets are presented in Table 3. The synthetic light curves were shown in Fig. 4.1.. Table 3. The parameters using in our model of XY Cet and DD Mon XY Cet* Parameter. Star 1. Star 2. DD Mon** Star 1. Star 2. Mass (M⊙). 1.773±0.016 1.615±0.016 1.39±0.07 0.93±0.06. Radius (R⊙). 1.873±0.035 1.773±0.029 1.44±0.03 1.29±0.03. Teff (K). 7870±115. 7620±125. 6250. 5195±4. linear LD coefficient. 0.49±0.10. 0.40±0.12. 0.547. 0.690. Orbital separation (R⊙) Orbital inclination (°). 12.493±0.035. 4.00±0.08***. 87.66±0.12. 77.2±0.1. *: Southworth J. et al.(2011) ; **: Gazeas et al.,(2010) and fixed mass ratio ***: Radius divided by Radiuspole (0.3598±0.0005) of star1. 23.

(31) Fig. 4.1 The synthetic light curves of XY Cet (upper panel) and DD Mon (bottom panel). The black dots show the photometric data from ASAS and the red line is the computed light curve.. 24.

(32) 4.2 Light curve analysis We use the least-squares method on the light curve fitting process, the chi-square cost function could be written as : n. 2 .  y. i.  f ( x i ,  ) 2. i 1. n - np. (4.2.1). where yi is the observed flux , f (xi,θ) is a computed flux at the certain unknown parameter sets, θ, and np is a number of adjusted parameters. We choose the Levenberg - Marquard algorithm to minimize this function and use a Python package LMFIT version 0.9.5 to do the fitting process. This algorithm was first applied to eclipsing binaries by Djurasevic (1992). The partial derivatives of the synthetic light curve at a given phase and model parameters in the Jacobian matrix using in the algorithm are calculated in two ways: the analytical form for limb-darkening coefficients and the temperature of the second component; the numerical derivatives approximated by using symmetric finite difference for inclination, mass ratio, and the radii of two components. We can take partial derivatives of equation (2.2.24) with respect to the parameters then get the analytical form: 2  F cos  2 3  0 0 4 Te,2 (1  x2  x2 cos  ) R sin dd T2 cos . (4.2.2). 2  F cos  2 4  0 0  Te,1 (cos   1) R sin dd x1 cos . (4.2.3). 2  F cos  2 4  0 0  Te,2 (cos   1) R sin dd x2 cos . (4.2.4). The numerical derivatives approximated by symmetric finite difference of 25.

(33) other 4 parameters are following the equation by Djurasevic (1992):. F F ( p  x / 2)  F ( p  x / 2)  p x. (4.2.5). where p is the model parameter and Δx is increment of the parameter. The size of Δx was setted to 0.01 for three parameters except inclination which was setted to 0.1 since the significant figures of the initial value of the parameters. In the analysis of XY Cet light curve, the mass ratio of the system fixed in q = 0.911 according to the radial velocity obtained by Southworth J., et al.(2011). The temperature of the two components was fixed in 7950K and 7600K respectively because of their spectral type A8+F0. The orbital separation was also fixed in 12.493 R⊙ obtained by Southworth J., et al.(2011). The adjustable parameters are the orbital inclination, radii of two stars and linear limb-darkening coefficients. The initinal guess value of these parameters are using the result of Southworth J. et al.(2011) which shown in Table 3. The fitting results of light curve were shown in Table 4 and the fitting curve presented in Fig. 4.2 bottom panel. In the analysis of DD Mon light curve, the mass ratio of the system was fixed in q = 0.67 according to the spectroscopic mass ratio obtained by Pribulla et al. (2009). The temperature of the primary component was fixed in 6250K due to the spectral type F5V. The following parameters were adjusted during the fitting iterations: orbital inclination, radii of two stars, orbital separation and linear limb-darkening coefficients. The initinal value of these parameters are using the result of Gazeas et al.,(2010) which shown in Table 3. The fitting results were shown in Table 4 and the fitting light curve was shown on the upper panel of Fig. 4.2. 26.

(34) Table 4 The result of light curve analysis of XY Cet and DD Mon Parameter. XY Cet. DD Mon. q. 0.911*. 0.67*. T1eff (K). 7950*. 6250*. T2eff (K). 7600*. 5265± 181. R1 (R⊙). 1.879±0.037. 1.468±0.027. R2 (R⊙). 1.816±0.030. 1.302±0.001. LD coefficient, x1. 0.56±0.25. 0.30±0.18. LD coefficient, x2. 0.58±0.24. 0.63±0.25. Orbital separation (R⊙). 12.493*. 3.982±0.001. Orbital inclination (°). 87.74±0.37. 75.82±1.51. chi-square. 0.00013. 0.00054. *: fixed. 27.

(35) Fig. 4.2 The fitting light curve of XY Cet (upper panel) and DD Mon (bottom panel). The black dots show the data from data from ASAS and the red line show the fitted light curve. 28.

(36) In the analysis of CQ Ser light curve, we encountered some difficulties since the relevant literature lacks the radial velocity required to directly determine the mass ratio, worst of all, there may be only one paper, Brancewicz & Dworak (1980), have done the light curve analysis of CQ Ser in the past four decades, so the initinal parameters need to be estimated from light curve directly in this case. There are two fixed parameters that one is the temperature of the primary component which was fixed in 7350K of the spectral type A8 and the other is the orbital separation which was fixed in 1 because that the radii of stars are in unit of the orbital separation. The adjustable parameters are the orbital inclination, mass ratio, the temperature of secondary component, radii of two stars and linear limb-darkening coefficients. The ways of estimating initinal guess from the CQ Ser light curve of these parameters are shown below: (a) The temperature : The ratio of primary to secondary eclipse depths is equal to the corresponding ratio of eclipsed star surface brightnesses so the ratio of temperature could be approximately:. 1  flux primary_eclipse 1  fluxsecondary_ eclipse. T    1   T2 . 4. (4.2.6). where we set the primary component having the higher temperature than secondary component. According to the light curve, fluxprimary_eclipse = 0.501, fluxsecondary_eclipse = 0.770, the temperature of secondary component could be setted to 6050K or lower value. (b) The radii : Here we assume the orbital inclination angle is very close to 90 degree. The ratio of eclipsed duration to the 0.25 phase may equal to the 2 times of 29.

(37) sum of radii of two components in the unit of orbital separation approximately, the 0.25 phase means the phase duration between the minimum and maximum brightness, so that could be approximately:.  2( R1  R2 ) Eclipsed duration     0.25 O rbital seperation  . (4.2.7). According to the light curve, the eclipsed duration is about 0.304 phase so the lower limit of the sum of radii in the unit of orbital separation could be 0.608. (c) The orbital inclination : Since there is no reliable way to know orbital inclination from the light curve, we could select the value by plotting the series of synthetic light curves keeping all parameters except inclination and observing the changes. (d) The mass ratio : We have been set an important constraint of star’s radii by caculating the critcal surface potential and Lagrangian point with mass ratio to avoid the over-contect situation in the program. The initinal value of mass ratio will be adjusted manually according to the initinal value of the radii of two components to support this constraint. (e) The linear limb-darkening coefficients : We set the initinal value to 0.5 both for two components because of the range of the coefficients is 0 to 1, and we select the half-point of this range. We repeated manual adjustment of these 7 parameters and the best parameter sets which have the relative minimum of cost function will be the initial parameters. The initial parameters and fitting results were shown in Table 5 and the fitting light curve in Fig. 4.3. 30.

(38) Fig. 4.3 The fitting light curve of CQ Ser. The red plus markers show the observation data and the blue line show the best fitting light curve. Table 5 The result of the light curve analysis of CQ Ser Parameter. Initinal value. Fitting results. q. 0.900. 0.899± 0.014. T1eff (K). 7350*. 7350*. T2eff (K). 5400. 5400± 189. R1/A. 0.366. 0.366±0.001. R2/A. 0.349. 0.348±0.005. LD coefficient, x1. 0.500. 0.499±0.077. LD coefficient, x2. 0.500. 0.499±0.369. Orbital separation, A. 1*. 1*. Orbital inclination (deg). 74.70. 74.70 ±0.15. chi-square. 0.0002673. *: fixed 31.

(39) V Discussion and Conclusions As a result from the synthetic and fitting light curves this work, XY Cet shows a detached system with the two components are absolutely not filling the critical Roche lobe, as shown in Fig. 5.1 upper panel. On the other hand, DD Mon shows that is a detached system but close to semi-detached system with the secondary component almost filling the critical Roche lobe, as shown in Fig. 5.1 middle panel. Comparing with the Fig. 5.1 lower panel, the degree of filling up the Roche lobe in our fitting result is lower than the result from Gazeas et al.,(2010) by using our model. This slightly difference shows that the model we used might underestimate the radii of stars. The fitting result of CQ Ser suggests that the system may be a semi-detached system with the both component almost filling the critical Roche lobe, as shown in the Fig. 5.2. The fitting light curve of CQ Ser is not perfect fit to the data. We can see that the curve have a bump around the 0.13, 0.4, 0.6 and 0.87 phase. These stages are exactly the beginning and the ending of the eclipsing events. We think that may be related to the limb-darkening effect since that the higher error of the limb-darkening coefficient of component 2 may be caused by the linear form of limb-darkening correction and the radii of component 2 underestimating. This underestimating could be improved by increasing the initial value of the fitting process. From the above results, the models we used have the ability to construct a physical model of the detached and semi-detached eclipsing binaries systems with certain reliability and also have ability to analyze 32.

(40) light curve and derive the system parameters. If the model could consider more physical property of the system such as gravity darkening, reflection effect and the property of non-blackbody radiation, we think that will be improve the applicability of the model.. 33.

(41) Fig 5.1 The configuration of XY Cet (upper panel) and DD Mon(middle and bottom panel) at phase 0.75. The upper panel show the shape of stars derived from our fitting results and the bottom panel is based on result from Gazeas et al.,(2010) .The red dots region indicate the area of Roche lobe which through the L1 point. 34.

(42) Fig. 5.2 The configuration of CQ Ser at phase 0.75. The red squares region indicate the area of Roche lobe which through the L1 point. It shows that both components are almost filling up with this region.. Fig. 5.3 The configuration of CQ Ser in 3D view at phase 0.16 with parameters obtained by our fitting result. The red plus markers indicate the component 2 and the blue dots show the component 1. 35.

(43) References Bradley W. Carroll, Dale A. Ostlie. An introduction to modern astrophysics 2nd ed. Brancewicz, H. K.,Dworak, T. Z., 1980,AcA, 30, 501 Djurasevic, G., 1992, Astrophys. Space. Sci. 197, 17 Gazeas, K., et al., 2010, PASP, 435 Josef Kallrath, Eugene F. Milone.,2009, Eclipsing Binary Stars: Modeling and Analysis Kopal, Z., 1959, Close Binary Systems Linnell, A.P., 1984, ApJ Suppl, 54, 17 Linnell, A.P., 1989, ApJ, 342, 449 Morrison D., Morrison N. D., 1968, AJ, 73, 777 Nelson B., Davis W. D., 1972, ApJ, 174, 617 Pojmanski, G. 2002, Acta Astronomica, 52, 397. Popper D. M., 1971, ApJ, 169, 549 Qian, S., et al., 1997, A&AS, 125,475 Qian, S., et al., 2009, PASJ,61, 331 Russell H. N., Merrill J. E., 1952, The Determination of the Elements of Eclipsing Binaries. Shaw J. S., 1994, MmSAI, 65, 95 Srivastava R. K., 1987, Ap&SS, 138, 197 Southworth J., et al., 2011, MNRAS, 414, 3740 Van Hamme, W., 1993, AJ 106, 2096 Wilson, R.E. & Devinney, E.J., 1971, ApJ, 166, 605 Yamasaki, A., et al., 1990, AJ, 99, 1218 36.

(44) Eclipsing Binary Stars, http://www.midnightkite.com/index.aspx?URL=Binary. 37.

(45) Appendix – A The PyRAF code of images data reduction # -*- coding: utf-8 -*import numpy as np from pyraf import iraf from astropy.io import fits def set_daopars(expadu,fwhm): """ Set all the iraf parameter files for daophot """ iraf.daophot.phot.verbose="no" iraf.daophot.phot.interactive="no" iraf.daophot.phot.verify="no" iraf.datapars.scale=1. iraf.datapars.fwhmpsf=fwhm iraf.datapars.emiss="yes" iraf.datapars.sigma=5. iraf.datapars.datamin=-0.001 iraf.datapars.datamax="INDEF" iraf.datapars.noise="poisson" iraf.datapars.epadu=expadu iraf.datapars.airmass="AIRMASS" iraf.datapars.filter="FILTER" iraf.datapars.itime=1. iraf.centerpars.calgo="centroid" iraf.centerpars.cbox=12.41 #1.7*fwhm iraf.centerpars.minsn=1. iraf.centerpars.cmaxiter=10 iraf.centerpars.maxshift=1. iraf.centerpars.clean="no" iraf.centerpars.mkcenter="no" iraf.daophot.daopars.function="auto" iraf.daophot.daopars.varorder=1 iraf.daophot.daopars.sat="no" iraf.daophot.daopars.matchrad=fwhm iraf.daophot.daopars.psfrad=2*fwhm iraf.daophot.daopars.fitrad=fwhm iraf.daophot.daopars.recenter="yes" iraf.daophot.daopars.fitsky="yes" iraf.daophot.daopars.groups="yes" iraf.daophot.daopars.sannu=fwhm*0.4246*6 iraf.daophot.daopars.wsann=fwhm*0.4246*9 iraf.daophot.daopars.proferr=2.5 iraf.daophot.daopars.maxiter=50 iraf.daophot.daopars.mergerad="INDEF" iraf.daophot.daopars.text="yes" def find_objects(inputImage,fwhm): """ Find the objects in the field use DAOfind """ output_locations= inputImage + ".stars" #find sky sigma state = iraf.imstat(inputImage,Stdout=1)[1].split() sigma = float(state[3]) print '%s avg_sigma: %.3f '%(inputImage,sigma) #set up some finding parameters iraf.daophot.datapars.sigma=sigma+7 iraf.daophot.datapars.fwhmpsf=fwhm iraf.daophot.findpars.threshold=3 #3sigma detections only iraf.daophot.findpars.mkdetections="no" #Do Daofind iraf.daofind(image=inputImage,output=output_locations,interactive="no",verify="no",verbose="no") print "Saved DaoFind output: %s \n"%(output_locations) return output_locations #return the name of the saved file def do_phot(inputImage, coord_list, aper, sky_annulus, width_sky):. 38.

(46) """ aperture photmoetry on the input image at the specified locations """ print "Phot aperture: %s pixels"%(aper) print "Sky Annulus: %i -> %i pixels"%(sky_annulus,sky_annulus+width_sky) mag_output_locations = inputImage+ ".phot" iraf.photpars.aperture=aper iraf.photpars.weighting="constant" iraf.photpars.mkapert="no" iraf.fitskypars.salgo="centroid" iraf.fitskypars.annu=sky_annulus iraf.fitskypars.dannu=width_sky iraf.fitskypars.smaxi=10 iraf.fitskypars.binsi=0.1 iraf.daophot.phot.interactive="no" iraf.daophot.phot.radplots="no" iraf.daophot.phot.verify="no" iraf.daophot.phot.verbose="no" #results on the screen or not iraf.phot.coords=coord_list iraf.phot.output=mag_output_locations iraf.daophot.phot(inputImage,coords=coord_list,output=mag_output_locations,verbose="no",verify=" no",interactive="no") print "Saved Phot output: %s \n"%(mag_output_locations) return mag_output_locations ########################################################################### img_file = open('name_list.txt','r') imglist = img_file.read().split() img_file.close() #FIRST: Find star's coordinate for each image using DAOfind coord_list = [] for img in imglist: #get info from header hdulist = fits.open(img) exptime = hdulist[0].header['EXPOSURE'] #read header name airmass = hdulist[0].header['AIRMASS'] #read header name egain = hdulist[0].header['EGAIN'] expadu = egain*exptime #set fwhm for psf and apphot size fwhm = 7.3 #do process set_daopars(expadu,fwhm) coord_list.append(find_objects(img,fwhm)) #filename #SECOND: Match star's coordinate for each image and # transform image to new image which have same coordinate iraf.images(_doprint=0) iraf.immatch(_doprint=0) #for xyxymatch, geomap, geotran of coords firstimg = '0519_007.fit' counting = 0 for img in imglist: #Runing xyxymatch and geomap and geotran to create new img if counting >= 0: # >= 0 is do all Nmatch=30 iraf.xyxymatch.unlearn()#defult match = iraf.xyxymatch(input=img+".stars",reference=firstimg+".stars", output=img+"_xymatch.out", toler=2, matching="triangles", nmatch=Nmatch,Stdout=1) match_number = int(match[12].split()[0]) print 'Number of stars Matched= ',match_number if match_number < 6: # Making sure atleast a few stars were matched. otherwise geoxytran will exit with error. print 'WRONG!!!' print("Failed to find the coordinates for "+img) print("We need to find the transformation interactively,you should keep ds9 opened") print("IMPORTANT: First select 3 stars in Frame 1 of ds9. Second image to select is in Frame 2") iraf.display(firstimg,1) iraf.display(img,2) iraf.xyxymatch.unlearn() match = iraf.xyxymatch(input=img+".stars",reference=firstimg+".stars", output=img+"_xymatch.out", toler=5, matching="tolerance", nmatch=10,interactive="yes",Stdout=1). 39.

(47) match_number = int(match[21].split()[0]) print 'Number of 3 stars Matched= ',match_number iraf.geomap(input=img+"_xymatch.out", database=img+"_tran.db", xmin=1, xmax=2184, ymin=1, ymax=1472, interactive=0) #Do image transform iraf.geotran(input=img, output=img+"_Trans.fit",database=img+"_tran.db",transforms=img+"_xymatch.out") counting = counting+1 print ('\nTransformed DONE !!\n') #THIRD: photometry by apphot #Choice which stars need to Phot from reference image and creat its coordinate and keep ds9 opened refimage = firstimg+'_Trans.fit' iraf.display(refimage,1,Stdout=1) print ('Choice which stars needed to do PHOT\n') print ('Press _a_ over stars which you want') print ('Press _q_ to leave when you choose finished.') print ('Note that the first star must be your target and be remember the order of your selection.\n') imx = iraf.imexam(Stdout=1) print 'Number of stars you selected:', (len(imx)-2)/2 selected_star = open('Selected_phot_stars.coo',"w") selected_star.write('# Coordinates of Selected Phot Stars from %s\n'%(refimage)) selected_star.write('# x y\n') for s in range((len(imx)-2)/2): #creat stars coordinate x = imx[2+s*2].split()[0] #coord at idx 2 4 6 8... y = imx[2+s*2].split()[1] selected_star.write('%s %s\n'%(x,y)) print '%s %s %i'%(x,y,s+1) selected_star.close() #Creat image(transformed image) list from original image list imglist_phot = '_Trans.fit '.join(imglist).split() imglist_phot[-1] = imglist_phot[-1]+'_Trans.fit' #last element needed #Do Apphot mag_list = []#store the result .mag file list for img in imglist_phot: hdulist = fits.open(img) exptime = hdulist[0].header['EXPOSURE'] airmass = hdulist[0].header['AIRMASS'] egain = hdulist[0].header['EGAIN'] expadu = egain*exptime #set fwhm for psf and apphot size fwhm = 7.3 aper = str(fwhm*0.4246*4) #4 sigma sky = fwhm*0.4246*6 #6 sigma width = fwhm*0.4246*3 #3 sigma set_daopars(expadu,fwhm) #Do Phot, creat .phot file for orignal photometry result phot_result = do_phot(img,'Selected_phot_stars.coo', aper=aper, sky_annulus=int(sky), width_sky=int(width)) #catch xcenter,ycenter,mag,merr from .phot file iraf.txdump(phot_result,fields='xcenter,ycenter,mag,merr',expr='yes',headers='yes',parameters='n o',Stdout='%s.mag'%img) mag_list.append('%s.mag'%img) #phot_result .mag list print 'Phot End!!!' #FOURTH: Creat final .txt output file, combine all mag from each .mag files of selected stars # alos combine observe_time from .fit header final_result = open('final_result.txt',"w") final_result.write('# Result of AppPhot of Star: CQ Ser\n') count_img = 0 for mag_file in mag_list: #do each .mag file hdulist = fits.open(imglist_phot[count_img]) date = hdulist[0].header['DATE-OBS'] #read header Time mag_data = np.loadtxt(mag_file,comments='#',usecols=(2,3),dtype=str) #read .mag file stars_number = mag_data[:,0].size if count_img == 0: line_label = '# Date' for n in range(stars_number): temp1 = ' '+'mag%i'%(n+1) temp2 = ' '+'%i_err'%(n+1) line_label = line_label+temp1+temp2 print line_label final_result.write('%s\n'%(line_label)) line_data = date for n in range(stars_number):. 40.

(48) temp1 = ' '+'%s'%(mag_data[n][0]) #mag temp2 = ' '+'%s'%(mag_data[n][1]) #err line_data = line_data+temp1+temp2 final_result.write('%s\n'%(line_data)) count_img = count_img +1 print line_data final_result.close() print 'final_result.txt has been created success!!' print '\nAll work Done!!'. 41.

(49) Appendix – B The source code of light curve synthesis-XY Cet # -*- coding: utf-8 -*""" For GCVS XY Cet (EA/DM) inc = 87.66 #inclination Rh = 1.873 #hot , solar unit Rc = 1.773 #cold, solar unit Th = 7870.0 #hot K Tc = 7620.0 #cold K Mh = 1.773 #hot sun mass Mc = 1.615 #cold sun mass sma = 12.493 #semimajor Axis , solar unit period = 2.780712 #days xh = 0.49 #LD coefficient hot xc = 0.40 """ from astropy import constants as const from astropy import units as u import numpy as np import matplotlib.pyplot as plt #constant sigma = const.sigma_sb # W*m^-2*k^-4 ,l=A*sigma*T^4 pi = np.pi ##parameter(h=hot,c=cold) totalphase = 1.2 sampling = 2 # n times per degree inc = 87.66 #orbital inclination 87.66 70~90 Th = 7970.0 #temperature per Sun or real temperature Tc = 7620.0 Mh = 1.773# mass per Sun Mc = 1.615 Rh = 1.873 #radius per Sun 0.132 Rc = 1.773 #0.108 sma = 12.493/Rh period = 2.78071392 #days fluxh = sigma*(Th*u.K)**4 #hot flux in real fluxc = sigma*(Tc*u.K)**4 #cold flux in real xh = 0.49 #LD coefficient hot xc = 0.4 #check primary always larger, p=primary,s=secondary if Rh > Rc: # primary is hot for XY Cet binary mp = Mh/Mh ms = Mc/Mh rp = Rh/Rh rs = Rc/Rh fp = fluxh/fluxh fs = fluxc/fluxh phi = 0 #starting angle for orbit, sec. infront of pri. (hot is behind) else: mp = Mc/Mc ms = Mh/Mc rp = Rc/Rc rs = Rh/Rc fp = fluxc/fluxc fs = fluxh/fluxh phi = pi #sec. behind of pri. ,in radian. #obsdata obs1 = np.loadtxt('xycetflux_asas.txt',dtype= float, ndmin= 2) p = 2.78071392 zerodate = 2438372.9455#By John Southworth (2011) obsphase1 = obs1[:, 0] obsflux = obs1[:, 1] #no eclips drp1 = np.linspace(0, rp , 151) drs1 = np.linspace(0, rs , 151) ldangp = np.arcsin((drp1-(drp1[1]-drp1[0])/2)/rp) ldangp[0] = 0.0. 42.

(50) ldangs = np.arcsin((drs1-(drs1[1]-drs1[0])/2)/rs) ldangs[0] = 0.0 ldp = fp*(1-xh*(1-np.cos(ldangp))) lds = fs*(1-xc*(1-np.cos(ldangs))) ldip = ldp*(drp1-(drp1[1]-drp1[0])/2)*(drp1[1]-drp1[0])*2*pi ldis = lds*(drs1-(drp1[1]-drp1[0])/2)*(drs1[1]-drs1[0])*2*pi itp = np.sum(ldip) its = np.sum(ldis) print "total intensity p = ", itp print "total intensity s = ", its print (fp*pi*rp**2)*(1-(xh/3)) print (fs*pi*rs**2)*(1-(xc/3)) ######################################################## i = np.deg2rad(inc) #change inclination to radians q = ms / mp #mass ratial ldiset = np.array([] , dtype=float) ldataset = np.array([] , dtype=float) #y-axis for luminusity phaseset = np.array([] , dtype=float) #x-axis for phase phaserange = np.linspace(0,360,totalphase*sampling*360 + 1) ##compute star positions and light curve (in circle orbital) drp = rp/150.0 drs = rs/150.0 for t in phaserange: phase = t/360. phaseangle = phi + 2*pi*phase X = sma*np.sin(phaseangle) #compute the coord. projected to x-y plane Y = sma*np.cos(phaseangle)*np.cos(i) Z = sma*np.cos(phaseangle)*np.sin(i) xp = -X/(1+(1/q)) #pri.star x-position yp = -Y/(1+(1/q)) zp = -Z/(1+(1/q)) xs = X/(1+q) ys = Y/(1+q) zs = Z/(1+q) dist = np.hypot(xp-xs,yp-ys) ldiset = ldiset*0.0 itotal = 0.0 anglast = 0.0 if dist > rp+rs: #no eclipses itotal = itp+its elif dist < rp-rs and zp > zs: # full eclipses and p 在前 itotal = itp ############################################################## if zs > zp and dist < rp+rs: # s 在前且有蝕 r_up = rp r_lim = 0.0 while r_up-drp/2 >= r_lim: #當 r up 仍然大於下限時則繼續 if r_up > abs(rs-dist): thetap = np.arccos((r_up**2-rs**2+dist**2)/(2*dist*r_up)) #cosine rule else: # pi or 0 if dist-rs > 0: thetap = 0.0 else: thetap = pi ldang = np.arcsin((r_up-drp/2)/rp)#取 dr 的中間為 ld 的半徑 ld = fp*(1-xh*(1-np.cos(ldang))) ldplus = 0.0###### if thetap < pi:###### deltaang = thetap-anglast###### ldplus = ld*drp*r_up*deltaang####### ldi = ld*(r_up-drp/2)*drp*2*thetap+ldplus # ldf*r*dr*ang = 該半徑之 arc 之總強度 ldiset = np.append(ldiset,ldi) #將 ldi 寫入 ldi array anglast = thetap####### r_up -= drp #由邊界往後縮 dr,直到 r up < rlim 後終止 loop eclips = np.sum(ldiset) itotal = itp+its-eclips ############################################################# if zp > zs and rp-rs < dist < rp+rs: # p 在前且有蝕 r_up = rs r_lim = 0.0 while r_up-drs/2 >= r_lim: if r_up > abs(rp-dist): thetap =np.arccos((r_up**2-rp**2+dist**2)/(2*dist*r_up)). 43.

(51) else: # pi or 0 if dist-rp > 0: thetap = 0.0 else: thetap = pi ldang = np.arcsin((r_up-drs/2)/rs) ld = fs*(1-xc*(1-np.cos(ldang))) ldplus = 0.0###### if thetap < pi: deltaang = thetap-anglast ldplus = ld*drs*r_up*deltaang ldi = ld*(r_up-drs/2)*drs*2*thetap+ldplus # ldf*r*dr*ang = 該半徑之 arc 之總強度 ldiset = np.append(ldiset,ldi) #將 ldi 寫入 ldi array anglast = thetap r_up -= drs eclips = np.sum(ldiset) itotal = itp+its-eclips ############################################################# ldataset = np.append(ldataset,itotal) phaseset = np.append(phaseset,phase) obsphase1 = np.append(obsphase1,obsphase1-1) obsflux = np.append(obsflux,obsflux) fmax = np.max(ldataset) fratio = ldataset/fmax fratio = np.append(fratio,fratio) phaseset = np.append(phaseset,phaseset-1) idx = np.argsort(phaseset) phaseset = phaseset[idx] fratio = fratio[idx] plt.figure().set_size_inches(10,6.3) plt.ylim(0.4 , 1.1) plt.xlim(-0.2 , 1.0) plt.plot(phaseset, fratio, 'r-',label='synthesis') plt.plot(obsphase1, obsflux,'k.',label='ASAS Data') plt.title('XY Cet Light Curve') plt.ylabel('Flux') plt.xlabel('Phase') plt.legend(loc='lower right') plt.minorticks_on() plt.savefig('xycet_synthesis_asas.png'). 44.

(52) Appendix – C The source code of light curve synthesis-DD Mon # -*- coding: utf-8 -*""" ################################################################### DD Mon (EA/SD) K. Gazeas(2010) inc = 78.2(qph),77.2(qsp) Rh = 1.62(qph),1.44(qsp) Rc = 1.37(qph),1.29(qsp) Mh = 1.94(qph),1.39(qsp) Mc = 1.12(qph),0.93(qsp) t1 = 6250 t2 = 5202 5195(qsp) rpole1 = 0.3702(qph),0.3598(qsp) rpole2 = 0.3123(qph),0.3231(qsp) x1=0.547(qph)0.547(qsp)v x2=0.689(qph)0.690(qsp)v omega1=3.230(qph),3.41(qsp) omega2=3.037(qph),3.19(qsp) sma = r1/r1pole or r2/r2pole(rpole=r/sma) 1.44/0.3598=4.0022 p = 0.568019day mag=11.1 - 11.8 p #################################################################### RV from Pribulla(2009) k1=133.32 km/s k2=198.85 km/s p=0.568020 d q=k1/k2=0.67=m2/m1 (m1+m2)sin^3(i)=2.157 solar mass a*sin(i)=(p/2pi)(k1+k2) = 2594525.2888104985 km a*sin(i)=3.73040 solar radius p0=52500.1745 + 0.568020 #################################################################### """ import numpy as np import matplotlib.pyplot as plt from astropy import constants as con import astropy.units as u sigma_bb = con.sigma_sb # W*m^-2*k^-4 ,l=A*sigma*T^4 lsun = con.L_sun #3.846e+26 W msun = con.M_sun #1.9891e+30 kg rsun = con.R_sun #695508000.0 m tsun = 5778.0 pi = np.pi def transform(xo,yo,zo,delta_i,phaseangle,sight,cmd1,cmd2): #XYplan 逆時針轉 phaseangle 為正(+z 軸), #2=center mass rotation for m2 #YZplan 逆時針轉 delta_i 為正(+x 軸) #1=center mass rotation for m1 if sight == 1 : xo = xo-cmd1 if sight == 2 : xo = xo+cmd2#+sma-cmd1# if sight == 3 : #for gradient sight vector iii = np.matrix([[1,0,0], [0,np.cos(delta_i),-np.sin(delta_i)], [0,np.sin(delta_i),np.cos(delta_i)]]) fi = np.matrix([[np.cos(phaseangle),-np.sin(phaseangle),0], [np.sin(phaseangle),np.cos(phaseangle),0], [0,0,1]]) ori = np.matrix([xo,yo,zo]).reshape(3,1) xyz = np.asarray(fi*iii*ori) x = xyz[0] y = xyz[1] z = xyz[2] else: ori = np.matrix([xo,yo,zo]) iii = np.matrix([[1,0,0], [0,np.cos(delta_i),np.sin(delta_i)], [0,-np.sin(delta_i),np.cos(delta_i)]]) fi = np.matrix([[np.cos(phaseangle),np.sin(phaseangle),0], [-np.sin(phaseangle),np.cos(phaseangle),0],. 45.

(53) [0,0,1]]) xyz x = y = z = return. = np.asarray(iii*fi*ori) xyz[0] xyz[1] xyz[2] x, y, z. def gradient(xset,yset,zset,q,delta_i,phaseangle,cmd1,cmd2): gdpx = (-xset*(xset**2+yset**2+zset**2)**-1.5) \ +(q*(1-xset)*((1-xset)**2+yset**2+zset**2)**-1.5)-q \ +xset*(q+1) gdpy = (-yset*(xset**2+yset**2+zset**2)**-1.5) \ -(q*(yset)*((1-xset)**2+yset**2+zset**2)**-1.5) \ +yset*(q+1) gdpz = (-zset*(xset**2+yset**2+zset**2)**-1.5) \ -(q*(zset)*((1-xset)**2+yset**2+zset**2)**-1.5) mag = (gdpx**2+gdpy**2+gdpz**2)**0.5 surface_nx = -gdpx/mag #outward sureface no negtive surface_ny = -gdpy/mag surface_nz = -gdpz/mag sight_nx, sight_ny, sight_nz = transform(0,1,0,delta_i,phaseangle,3,cmd1,cmd2) cosr = surface_nx*sight_nx + surface_ny*sight_ny + surface_nz*sight_nz cosridx = np.where(cosr <= 0) surface_sight_nx=surface_nx[cosridx] surface_sight_ny=surface_ny[cosridx] surface_sight_nz=surface_nz[cosridx] r_nx = xset[cosridx] r_ny = yset[cosridx] r_nz = zset[cosridx] magr = (r_nx**2+r_ny**2+r_nz**2)**0.5 cosb = (r_nx/magr)*surface_sight_nx + (r_ny/magr)*surface_sight_ny + (r_nz/magr)*surface_sight_nz cosridxall = np.array(range(cosr.size),dtype=int) surface_all_nx=surface_nx[cosridxall] surface_all_ny=surface_ny[cosridxall] surface_all_nz=surface_nz[cosridxall] r_nx = xset[cosridxall] r_ny = yset[cosridxall] r_nz = zset[cosridxall] magr = (r_nx**2+r_ny**2+r_nz**2)**0.5 all_cosb = (r_nx/magr)*surface_all_nx + (r_ny/magr)*surface_all_ny + (r_nz/magr)*surface_all_nz return cosr,cosb,all_cosb def bdoriidx(cosr,component): if component == 1: lim = -0.1 #-0.047 elif component == 2: lim = -0.1 # cosr_u = np.where(cosr<=0) cosr_l = np.where(cosr>=lim) idx = np.intersect1d(cosr_u,cosr_l) #取交集 return idx def bdsort(bx,by,comp): if comp==2: xin = (np.max(bx)+np.min(bx))/2. yin = (np.max(by)+np.min(by))/2. bx = bx-xin # 圖型平移使中心點為(0,0) by = by-yin bdtheta = np.arccos(bx/np.hypot(bx,by))#只有 0~pi bdtheta[np.where(by<0)] = bdtheta[np.where(by<0)]+pi#y 為負的角度再加 pi bdsortidx1 = np.argsort(bdtheta) #(要的)重新排好序的 idx 值,即原角度的 idx 值按照新的次序(0~2pi)排列 bdtheta = bdtheta[bdsortidx1] # 排好序的角度值 but 超過 pi 的 idx 需再反轉 bx = bx[bdsortidx1]#+xin### 已照次序排好的邊界點群 by = by[bdsortidx1]#+yin### bx[np.where(bdtheta>=pi)] = bx[np.where(bdtheta>=pi)][::-1] by[np.where(bdtheta>=pi)] = by[np.where(bdtheta>=pi)][::-1] if comp==2: bx = bx+xin by = by+yin return bx,by. 46.

(54) def eclips(bx1,by1,bx2,by2,sightx,sighty,phase):#bx1=前星邊界 ex.bdx1ori #sightx=後星可視點 ex.sight2x xc = (np.max(bx1)+np.min(bx1))/2. #前星取中心點 yc = (np.max(by1)+np.min(by1))/2. if phase > 0.75 or phase <0.25: bx1,by1=bdsort(bx1,by1,2) else: bx1,by1=bdsort(bx1,by1,1) distmin = np.array([],dtype=float)# for bb in range(bx1.size): # 計算"前星邊界"與"後星邊界"的間距 bd1tobd2 = np.hypot(bx1[bb]-bx2,by1[bb]-by2) dmin=np.min(bd1tobd2) distmin = np.append(distmin,dmin)#記錄每一個前星邊界對應的最小距離 idx = np.argsort(distmin) #兩邊界相距由小到大 distmin = distmin[idx] bxtop10 = bx1[idx][0:20] bytop10 = by1[idx][0:20] ctotop10 = np.hypot(xc-bxtop10,yc-bytop10)#min2max if np.abs(ctotop10[0]-ctotop10[1]) < 10e-6 :#x 軸對稱 bxmin=bxtop10[0] bymin=bytop10[0] bxmax=bxtop10[1] bymax=bytop10[1] else: bxmin=bxtop10[np.where(ctotop10==np.min(ctotop10))[0][0]] bymin=bytop10[np.where(ctotop10==np.min(ctotop10))[0][0]] bxmax=bxtop10[np.where(ctotop10==np.max(ctotop10))[0][0]] bymax=bytop10[np.where(ctotop10==np.max(ctotop10))[0][0]] #print 'bd_crosspoint= \n', bxmin,bymin,'| \n',bxmax,bymax for cc in range(bx1.size): #idx if bx1[cc]==bxmin and by1[cc]==bymin :p1idx=cc if bx1[cc]==bxmax and by1[cc]==bymax :p2idx=cc if p1idx == 0: p1idx = p1idx+1 if p2idx == 0: p2idx = p2idx+1 if p1idx == p2idx : p2idx = p2idx+1 if p1idx < p2idx:#起始點角度小於終點角(idx 大=角大) if phase >= 0.25 and phase <=0.75: if ((bxmin*bxmax>0 and bymin*bymax<0) or (bxmin*bxmax<0 and bymin*bymax<0)) and phase>0.5:#1-4, 1-3, 2-3, bdpartx=bx1[0:p1idx] bdparty=by1[0:p1idx] bdpartx=np.append(bdpartx,bx1[p2idx:bx1.size-1]) bdparty=np.append(bdparty,by1[p2idx:bx1.size-1]) else: bdpartx=bx1[p1idx:p2idx] bdparty=by1[p1idx:p2idx] if phase > 0.75 or phase <0.25: if ((bxmin*bxmax<0 and bymin*bymax<0) or (bxmin*bxmax>0 and bymin*bymax<0)) and phase>0.75 :# bdpartx=bx1[p1idx:p2idx] bdparty=by1[p1idx:p2idx] elif ((bxmin*bxmax<0 and bymin*bymax<0) or (bxmin*bxmax>0 and bymin*bymax<0)) : #idx1=0 idx2=41 bdpartx=bx1[0:p1idx] bdparty=by1[0:p1idx] bdpartx=np.append(bdpartx,bx1[p2idx:bx1.size-1]) bdparty=np.append(bdparty,by1[p2idx:bx1.size-1]) else: bdpartx=bx1[p1idx:p2idx] bdparty=by1[p1idx:p2idx] elif p1idx > p2idx : if phase >= 0.25 and phase <=0.75: if ((bxmin*bxmax>0 and bymin*bymax<0) or (bxmin*bxmax<0 and bymin*bymax<0)) and phase>0.5:#1-4 or 1-3 bdpartx=bx1[p1idx:bx1.size-1] bdparty=by1[p1idx:bx1.size-1] bdpartx=np.append(bdpartx,bx1[0:p2idx]) bdparty=np.append(bdparty,by1[0:p2idx]) else: bdpartx=bx1[p2idx:p1idx] bdparty=by1[p2idx:p1idx] if phase > 0.75 or phase <0.25: if ((bxmin*bxmax<0 and bymin*bymax<0) or (bxmin*bxmax>0 and bymin*bymax<0)) and phase>0.75 :# bdpartx=bx1[p2idx:p1idx]. 47.

(55) bdparty=by1[p2idx:p1idx] elif ((bxmin*bxmax<0 and bymin*bymax<0) or (bxmin*bxmax>0 and bymin*bymax<0)) :# bdpartx=bx1[p1idx:bx1.size-1] bdparty=by1[p1idx:bx1.size-1] bdpartx=np.append(bdpartx,bx1[0:p2idx]) bdparty=np.append(bdparty,by1[0:p2idx]) else: bdpartx=bx1[p2idx:p1idx] bdparty=by1[p2idx:p1idx] ## 計算 xsight = np.array([],dtype=float)#back star's visible point ysight = np.array([],dtype=float)# idxsight = np.array([],dtype=int) for ss in range(sightx.size): if xc-sightx[ss] == 0.: ms = 0. else: ms = (yc-sighty[ss])/(xc-sightx[ss]) b = yc-ms*xc #ms*x-y+b=0 bptoline = np.abs((ms*bdpartx-bdparty+b))/((ms**2+1.)**0.5) if bptoline.size == 0: print bptoline bptolinex = bdpartx[np.where(bptoline == np.min(bptoline))][0] bptoliney = bdparty[np.where(bptoline == np.min(bptoline))][0] ctosight = np.hypot(xc-sightx[ss],yc-sighty[ss]) ctobp = np.hypot(xc-bptolinex,yc-bptoliney) #print ctobp if ctosight >= ctobp: xsight= np.append(xsight,sightx[ss]) ysight= np.append(ysight,sighty[ss]) idxsight = np.append(idxsight,ss) return bdpartx,bdparty,xsight,ysight,idxsight def lagrangian1(delta_i,phaseangle,cmd1,cmd2): rini = 0.5 q = q1 l1 = rini for t in range(10): g = (-1/l1**2)-q*(((l1-1)/(l1**2-2*l1+1)**1.5)+1)+(1+q)*l1 gd1 = (2/(l1**3))+(2*q/(l1**2-2*l1+1)**1.5)+(1+q) l1 = l1-g/gd1 l1xo = np.array([l1*sma], dtype=float) # for l1x>>(l1*sma)-cmd1 l1yo = np.array([0.], dtype=float) l1zo = np.array([0.], dtype=float) return l1xo, l1yo, l1zo def lagrangian2(delta_i,phaseangle,cmd1,cmd2): rini = 0.5 q = 1/q1 l1 = rini for t in range(10): g = (-1/l1**2)-q*(((l1-1)/(l1**2-2*l1+1)**1.5)+1)+(1+q)*l1 gd1 = (2/(l1**3))+(2*q/(l1**2-2*l1+1)**1.5)+(1+q) l1 = l1-g/gd1 l2xo = np.array([l1*sma], dtype=float) l2yo = np.array([0.], dtype=float) l2zo = np.array([0.], dtype=float) return l2xo, l2yo, l2zo #Read_data############################################# obsdata=np.loadtxt('ddmonflux_asas.txt',comments='#',usecols=(0,1),dtype='float') period = 0.56801844 #qian 2009 time_min = 2450100.2938 #qian 2009 phase_obs = obsdata[:,0] fratio_obs = obsdata[:,1] ############################################################################## surface_part = 73# 10d=37 3d=121 5d=73 2.5d=145 2d=181 phaseset = phase_obs phi0 = pi/2#need to plus pi/2 XYplan 逆時針轉為正角 phaseangleset = phi0 + 2*pi*phaseset theta = np.deg2rad(np.linspace(0, 180, surface_part/2+1)) # 0~pi , 90deg=at 赤道 phiset = np.deg2rad(np.linspace(0, 360, surface_part)) # 0 = +x, 180 = -x, 90 = +y, -90 = -y dtheta = pi/((surface_part-1)/2.)#wei do dphi = (2.*pi)/(surface_part-1) nu = np.cos(theta) #array inc = 77.2 delta_i = (pi/2)-(inc*pi/180). 48.

(56) q1 = 0.93/1.39 q2 = 1/q1 t1 = 6250.#sp t2 = 5195 r1 = 1.44 # sp r2 = 1.29 # sp sma = r1/0.3598 #+1.5 limb_x1= 0.547#sp limb_x2=0.690#sp cmd1 = sma*(q1/(q1+1))#center mass dist. from 1 cmd2 = sma*(q2/(q2+1)) flux1 = (sigma_bb *(t1*u.K)**4)/pi flux2 = (sigma_bb *(t2*u.K)**4)/pi ftotal_all = np.array([] , dtype=float) ############################################################################## xset = np.array([] , dtype=float) yset = np.array([] , dtype=float) zset = np.array([] , dtype=float) rr_1 = np.array([] , dtype=float) theta1 = np.array([] , dtype=float) phi1 = np.array([] , dtype=float) r1pole = r1/sma const = 1/r1pole + q1/(r1pole**2+1)**0.5 rr = r1pole for t in range(20): ans = 1/(const-q1/(rr**2+1)**0.5) rr = ans #array (whole wei do at specific gen do) omega1 = const for phi in phiset: lamda = np.sin(theta)*np.cos(phi) #theta(wei do) is array at specific gen do for t in range(20): ans = 1/(const-q1/(rr**2-2*rr*lamda+1)**0.5+q1*rr*lamda-0.5*(1+q1)*(1-nu**2)*rr**2) rr = ans #array (whole wei do at specific gen do) if phi == 0: r1xface = (max(rr)*sma) #哪一個緯度的 rr 最大,face to face(phi=0) if 3.1 < phi <3.23: r1xback = (max(rr)*sma)#back(phi=pi) if 1.52 < phi <1.61: r1yside = (max(rr)*sma)#side(phi=pi/2) x1o = rr*np.sin(theta)*np.cos(phi) y1o = rr*np.sin(theta)*np.sin(phi) z1o = rr*np.cos(theta) xset = np.append(xset,x1o)#save x-axis value at this specific gen do yset = np.append(yset,y1o)# zset = np.append(zset,z1o)# rr_1 = np.append(rr_1,rr*sma)#save these point's radious at this specific gen do theta1 = np.append(theta1,theta)#save these point's theta (wei do) phi1 = np.append(phi1,phi)#save these point's phi (gen do) (same loop smae value) ############################################################################## x2set = np.array([] , dtype=float) y2set = np.array([] , dtype=float) z2set = np.array([] , dtype=float) rr_2 = np.array([] , dtype=float) theta2 = np.array([] , dtype=float) phi2 = np.array([] , dtype=float) r2pole = r2/sma const = 1/r2pole + q2/(r2pole**2+1)**0.5 omega2 = const rr = r2pole for phi in phiset: lamda = np.sin(theta)*np.cos(phi) for t in range(20): ans = 1/(const-q2/(rr**2-2*rr*lamda+1)**0.5+q2*rr*lamda-0.5*(1+q2)*(1-nu**2)*rr**2) rr = ans if phi == 0: r2xface = (max(rr)*sma) #哪一個緯度的 rr 最大,face to face if 3.1 < phi <3.23: r2xback = (max(rr)*sma)#back if 1.52 < phi <1.61: r2yside = (max(rr)*sma)#side x2o = -rr*np.sin(theta)*np.cos(phi) #-x+sma y2o = rr*np.sin(theta)*np.sin(phi) z2o = rr*np.cos(theta) x2set = np.append(x2set,x2o) y2set = np.append(y2set,y2o) z2set = np.append(z2set,z2o) rr_2 = np.append(rr_2,rr*sma) theta2 = np.append(theta2,theta) phi2 = np.append(phi2,phi) ##############################################################################. 49.

(57) l1x, l1y, l1z = lagrangian1(delta_i,0.,cmd1,cmd2)#center is 0 l2x, l2y, l2z = lagrangian2(delta_i,0.,cmd1,cmd2)#center is 0 print 'L1 position_1=', l1x[0] if r1xface > l1x[0] or r2xface > l2x[0] : print 'OVER Contect' ##################################################################################### for ph in range(phaseset.size): phase = phaseset[ph] phaseangle = phaseangleset[ph] x1, y1, z1 = transform(xset*sma,yset*sma,zset*sma,delta_i,phaseangle,1,cmd1,cmd2)#center mass rotation for m1 cosr1,sight_cosb1,all_cosb1 = gradient(xset,yset,zset,q1,delta_i,phaseangle,cmd1,cmd2)#用質點中心旋轉後的座標 cosridx1 = np.where(cosr1 <= 0)# return index #cosrvalue = np.extract(cosr<0, cosr) #return value sight_cosr1 = np.cos(pi-np.arccos(cosr1[cosridx1])) sight1x = x1[cosridx1]# sight1y = y1[cosridx1] sight1z = z1[cosridx1]# ##################################################### x2, y2, z2 = transform(x2set*sma,y2set*sma,z2set*sma,delta_i,phaseangle,2,cmd1,cmd2)#center mass rotation for m2 cosr2,sight_cosb2,all_cosb2 = gradient(x2set,y2set,z2set,q2,delta_i,phaseangle,cmd1,cmd2) cosridx2 = np.where(cosr2 <= 0) sight_cosr2 = np.cos(pi-np.arccos(cosr2[cosridx2])) sight2x = x2[cosridx2] sight2y = y2[cosridx2] sight2z = z2[cosridx2] ##################################################### cm1x, cm1y, cm1z = transform(np.zeros(1),np.zeros(1),np.zeros(1),delta_i,phaseangle,1,cmd1,cmd2) cm2x, cm2y, cm2z = transform(np.zeros(1),np.zeros(1),np.zeros(1),delta_i,phaseangle,2,cmd1,cmd2) dist = np.hypot(cm1x-cm2x,cm1z-cm2z)[0] distr1 = np.max(np.hypot(cm1x-sight1x,cm1z-sight1z)) distr2 = np.max(np.hypot(cm2x-sight2x,cm2z-sight2z)) if dist <= distr1+distr2: #兩星中心相距小於 r1max+r2max # print 'eclips, phase=',phase bdidx1 = bdoriidx(cosr1,1) bdx1ori = x1[bdidx1] bdy1ori = y1[bdidx1] bdz1ori = z1[bdidx1] bdidx2 = bdoriidx(cosr2,2) bdx2ori = x2[bdidx2] bdy2ori = y2[bdidx2] bdz2ori = z2[bdidx2] if phase >= 0.25 and phase <=0.75: #star1 前 bdpartx,bdparty,xsight,ysight,idxsight = eclips(bdx1ori,bdz1ori,bdx2ori,bdz2ori,sight2x,sight2z,phase) sight2x = sight2x[idxsight] sight2y = sight2y[idxsight] sight2z = sight2z[idxsight] sight_cosr2 = sight_cosr2[idxsight] sight_cosb2 = sight_cosb2[idxsight] rr_2_new = rr_2[cosridx2][idxsight] theta2_new = theta2[cosridx2][idxsight] rr_1_new = rr_1[cosridx1] theta1_new = theta1[cosridx1] if phase > 0.75 or phase <0.25: #star2 前 bdpartx,bdparty,xsight,ysight,idxsight = eclips(bdx2ori,bdz2ori,bdx1ori,bdz1ori,sight1x,sight1z,phase) sight1x = sight1x[idxsight] sight1y = sight1y[idxsight] sight1z = sight1z[idxsight] sight_cosr1 = sight_cosr1[idxsight] sight_cosb1 = sight_cosb1[idxsight] rr_1_new = rr_1[cosridx1][idxsight] theta1_new = theta1[cosridx1][idxsight] rr_2_new = rr_2[cosridx2] theta2_new = theta2[cosridx2] else: rr_1_new = rr_1[cosridx1] theta1_new = theta1[cosridx1] rr_2_new = rr_2[cosridx2] theta2_new = theta2[cosridx2]. 50.