1 """
2 Roche models of binary or multiple stars.
3
4 You'd better not use this module, but talk to Pieter Degroote.
5
6 The premisse of the Roche potential is that the stellar mass can be represented
7 by a point source.
8
9 Currently, three types of stellar shapes are implemented:
10 1. asynchronously rotating eccentric binary (L{get_binary_roche_radius})
11 2. fast rotating star (L{get_fastrot_roche_radius})
12 3. differentially (fast) rotating star (L{get_diffrot_roche_surface_gravity})
13
14 This module can be used to calculate the following information:
15 1. the distorted shape of the star due to a Roche potential
16 2. the local surface gravity assuming Von Zeipel gravity darkening
17 3. the L{local effective temperature<local_temperature>}
18 4. the local velocity vector due to rotation
19 5. the radial velocity
20 6. the total distorted surface area
21 7. the total distorted luminosity
22 8. the L{projected intensity<projected_intensity>} in the line of sight
23 9. the mean radial velocity (e.g. Rossiter effect in binaries)
24 10. a L{synthetic spectrum<spectral_synthesis>} starting from a library of spectra
25 11. a L{synthetic light curve<binary_light_curve_synthesis>} from the projected intensities
26
27 Section 1. Binary star
28 ======================
29
30 As an example we take the binary star SX Aurigae. This is the minimum set of
31 parameters.
32
33 >>> P = 1.2100802 # period in days
34 >>> e = 0. # eccentricity
35 >>> q = 0.54369 # mass ratio M2/M1
36 >>> asini = 11.9*constants.Rsol/constants.au # semi-major axis
37 >>> i = 81.27 # inclination angle
38 >>> F = 1. # synchronicity parameter
39 >>> d = 1. # separation in semi-major axis units
40 >>> Phi1 = 3. # gravitational potential primary
41 >>> Phi2 = 5.05 # gravitational potential secondary
42 >>> T_pole1 = 25000. # polar temperature primary
43 >>> T_pole2 = 18850. # polar temperature secondary
44 >>> M1 = 10.3 # primary mass
45 >>> M2 = 5.6 # secondary mass
46
47 Note that we actually do not need the masses, since we can derive them from
48 the third kepler law with the other parameters. We do so for convenience.
49
50 Everything is calculated in units of semi-major axis. Thus, to convert to SI
51 units we need the following conversion factor:
52
53 >>> a = asini * sin(i/180.*pi)
54 >>> to_SI = a*constants.au
55
56 We need some constants for the calculations: the polar radii and angular velocity
57
58 >>> r_pole1 = newton(binary_roche_potential,1e-5,args=(0,0,Phi1, q,d,F))
59 >>> r_pole2 = newton(binary_roche_potential,1e-5,args=(0,0,Phi2,1/q,d,F))
60 >>> omega_rot = F * 2*pi/(P*24*3600) * 1/d**2 * sqrt( (1+e)*(1-e))
61 >>> omega_rot_vec = np.array([0.,0.,-omega_rot])
62
63 Derive the shape of the two stars
64
65 >>> radius1 = np.array([get_binary_roche_radius(itheta,iphi,Phi=Phi1,q= q,d=d,F=F,r_pole=r_pole1) for itheta,iphi in zip(thetas,phis)]).reshape(theta.shape)
66 >>> radius2 = np.array([get_binary_roche_radius(itheta,iphi,Phi=Phi2,q=1/q,d=d,F=F,r_pole=r_pole2) for itheta,iphi in zip(thetas,phis)]).reshape(theta.shape)
67
68 We focus on the primary, then repeat everything for the secondary: The local
69 surface gravity can only be calculated if we have Cartesian coordinates.
70
71 >>> x1,y1,z1 = vectors.spher2cart_coord(radius1,phi,theta)
72 >>> g_pole1 = binary_roche_surface_gravity(0,0,r_pole1*to_SI,d*to_SI,omega_rot,M1*constants.Msol,M2*constants.Msol,norm=True)
73 >>> Gamma_pole1 = binary_roche_potential_gradient(0,0,r_pole1,q,d,F,norm=True)
74 >>> zeta1 = g_pole1 / Gamma_pole1
75 >>> dOmega1 = binary_roche_potential_gradient(x1,y1,z1,q,d,F,norm=False)
76 >>> grav_local1 = dOmega1*zeta1
77 >>> grav_local1 = np.array([i.reshape(theta.shape) for i in grav_local1])
78 >>> grav1 = vectors.norm(grav_local1)
79
80 Now we can, as before, calculate the other local quantities:
81
82 >>> areas_local1,cos_gamma1 = surface_elements((radius1,theta,phi),-grav_local1)
83 >>> teff_local1 = local_temperature(grav1,g_pole1,T_pole1,beta=1.)
84 >>> ints_local1 = local_intensity(teff_local1,grav1,np.ones_like(cos_gamma1),photband='OPEN.BOL')
85 >>> velo_local1 = np.cross(np.array([x1,y1,z1]).T*to_SI,omega_rot_vec).T
86
87
88 Now make some plots showing the local quantities:
89
90 >>> quantities = areas_local1,np.log10(grav1*100),teff_local1,ints_local1
91 >>> names = 'Area','log g', 'Teff', 'Flux'
92 >>> p = pl.figure()
93 >>> rows,cols = 2,2
94 >>> for i,(quantity,name) in enumerate(zip(quantities,names)):
95 ... p = pl.subplot(rows,cols,i+1)
96 ... p = pl.title(name)
97 ... q = quantity.ravel()
98 ... vmin,vmax = q[-np.isnan(q)].min(),q[-np.isnan(q)].max()
99 ... if vmin==vmax: vmin,vmax = q[-np.isnan(q)].mean()-0.01*q[-np.isnan(q)].mean(),q[-np.isnan(q)].mean()+0.01*q[-np.isnan(q)].mean()
100 ... p = pl.scatter(phis/pi*180,thetas/pi*180,c=q,edgecolors='none',vmin=vmin,vmax=vmax)
101 ... p = pl.colorbar()
102 ... p = pl.xlim(0,360)
103 ... p = pl.ylim(180,0)
104
105 ]include figure]]ivs_binary_rochepotential_binarystar.png]
106
107 >>> p = pl.figure()
108 >>> p = pl.subplot(121,aspect='equal')
109 >>> p = pl.title('top view')
110 >>> p = pl.scatter(x1,y1,c=teff_local1.ravel(),edgecolors='none')
111 >>> p = pl.subplot(122,aspect='equal')
112 >>> p = pl.title('side view')
113 >>> p = pl.scatter(x1,z1,c=teff_local1.ravel(),edgecolors='none')
114
115 ]include figure]]ivs_binary_rochepotential_binarystar_shape.png]
116
117 Section 5. Synthesizing spectra
118 ===============================
119
120 In this case study, we calculate the projected radial velocity of a uniformly
121 rotating star, but in the limit of no deformation due to the rotation. We do this,
122 so that we can easily compare rotational broadening calculated with the ROTIN
123 package, and numerically. So let us build a star with these parameters:
124
125 >>> omega = 0.0 # ratio to critical velocity
126 >>> T_pole = 13000. # K
127 >>> r_pole = 2.0 # solar radii
128 >>> M = 1.5 # solar mass
129 >>> view_angle = pi/2 # radians
130 >>> theta,phi = get_grid(20,100,full=True,gtype='spher')
131 >>> thetas,phis = np.ravel(theta),np.ravel(phi)
132
133 Then calculate the shape of this star
134
135 >>> radius = np.array([get_fastrot_roche_radius(itheta,r_pole,omega) for itheta in thetas]).reshape(theta.shape)
136 >>> grav_local = np.array([fastrot_roche_surface_gravity(iradius,itheta,iphi,r_pole,omega,M) for iradius,itheta,iphi in zip(radius.ravel(),thetas,phis)]).T
137 >>> grav_local = np.array([i.reshape(theta.shape) for i in grav_local])
138 >>> g_pole = fastrot_roche_surface_gravity(r_pole,0,0,r_pole,omega,M)[-1]
139 >>> grav = vectors.norm(grav_local)
140
141 and the local quantities
142
143 >>> areas_local,cos_gamma = surface_elements((radius,theta,phi),-grav_local)
144 >>> teff_local = local_temperature(vectors.norm(grav_local),g_pole,T_pole,beta=1.)
145 >>> ints_local = local_intensity(teff_local,grav,photband='OPEN.BOL')
146 >>> x,y,z = vectors.spher2cart_coord(radius.ravel(),phis,thetas)
147
148 Assume, with a shape of a non-rotating star, that we have a velocity component
149 on the surface of the star:
150
151 >>> myomega = 0.5
152 >>> velo_local = diffrot_velocity((phi,theta,radius*constants.Rsol),myomega,myomega,r_pole,M)
153
154 Collect all the necessary information in one record array.
155
156 >>> starlist = [x,y,z] + [i.ravel() for i in velo_local] + [i.ravel() for i in grav_local] + [teff_local.ravel(),areas_local.ravel()]
157 >>> starnames = ['x','y','z','vx','vy','vz','gravx','gravy','gravz','teff','areas']
158 >>> star = np.rec.array(starlist,names=starnames)
159
160 Project the star in some line-of-sight. The velocity component in the X-direction
161 is the radial velocity.
162
163 >>> view_angle = pi/2 # edge on
164 >>> mystar = project(star,view_long=(0,0,0),view_lat=(view_angle,0,0),photband='OPEN.BOL',only_visible=True,plot_sort=True)
165
166 We can calculate the synthetic spectra for all surface elements between 7055 and
167 7075 angstrom:
168
169 >>> loggs = np.log10(np.sqrt(mystar['gravx']**2 + mystar['gravy']**2 + mystar['gravz']**2)*100)
170 >>> iterator = zip(mystar['teff'],loggs,mystar['vx']/1000.)
171 >>> wave_spec = np.linspace(7055,7075,750)
172 >>> spectra = np.array([spectra_model.get_table(teff=iteff,logg=ilogg,vrad=ivrad,wave=wave_spec)[1] for iteff,ilogg,ivrad in iterator])
173
174 The total observed spectrum is then simply the weighted sum with the local projected
175 intensities:
176
177 >>> average_spectrum = np.average(spectra,weights=mystar['projflux'],axis=0)
178
179 We compare with the ROTIN package, which assumes a linear limbdarkening law
180
181 >>> original = spectra_model.get_table(teff=mystar['teff'][0],logg=loggs[0],wave=wave_spec)[1]
182 >>> rotin1 = spectra_model.get_table(teff=mystar['teff'][0],logg=loggs[0],vrot=vmax/1000.*sin(view_angle),wave=wave_spec,fwhm=0,epsilon=0.6,stepr=-1)[1]
183
184 And a plot can be made via
185
186 >>> colors = mystar['eyeflux']/mystar['eyeflux'].max()
187 >>> p = pl.figure(figsize=(7,7))
188 >>> ax = pl.axes([0,0.5,1,0.5])
189 >>> p = ax.set_aspect('equal')
190 >>> p = ax.set_axis_bgcolor('k')
191 >>> p = pl.box(on=False)
192 >>> p = pl.xticks([]);p = pl.yticks([])
193 >>> p = pl.scatter(mystar['y'],mystar['z'],c=colors,edgecolors='none',cmap=pl.cm.gray,vmin=0,vmax=1)
194 >>> p = pl.scatter(mystar['y']+1.02*mystar['y'].ptp(),mystar['z'],c=mystar['vx'],edgecolors='none',vmin=vmin,vmax=vmax,cmap=pl.cm.RdBu)
195 >>> ax = pl.axes([0.13,0.1,0.78,0.4])
196 >>> p = pl.plot(wave_spec,average_spectrum,'k-',lw=2,label='Numerical')
197 >>> p = pl.plot(wave_spec,original,'r-',label='No rotation')
198 >>> p = pl.plot(wave_spec,rotin1,'b-',label='Rotin v=%.1f km/s'%(vmax/1000.),lw=2)
199 >>> p = pl.ylim(0.95,1.0)
200 >>> p = pl.legend(loc='lower right',prop=dict(size='small'))
201 >>> p = pl.xlabel('Wavelength [angstrom]',color='1')
202 >>> p = pl.ylabel('Normalised flux',color='1')
203
204 ]include figure]]ivs_binary_rochepotential_norot_spectrum_b.png]
205
206 ]include figure]]ivs_binary_rochepotential_norot_spectrum_a.png]
207
208
209 """
210 import logging
211 import os
212 import pylab as pl
213 import numpy as np
214 from numpy import pi,cos,sin,sqrt,nan
215 from scipy.optimize import newton
216 from scipy.spatial import KDTree
217 try:
218 from scipy.spatial import Delaunay
219 except ImportError:
220 print "import Error Delaunay"
221 import time
222 from ivs.timeseries import keplerorbit
223 from ivs.units import constants
224 from ivs.units import conversions
225 from ivs.coordinates import vectors
226 from ivs.sed import model as sed_model
227 from ivs.sed import limbdark
228 from ivs.spectra import model as spectra_model
229 from ivs.roche import local
230 from ivs.aux import loggers
231 from ivs.inout import ascii
232 from ivs.inout import fits
233
234 logger = logging.getLogger("BIN.ROCHE")
235
236
237
239 """
240 Unitless eccentric asynchronous Roche potential in spherical coordinates.
241
242 See Wilson, 1979.
243
244 The synchronicity parameter F is 1 for synchronised circular orbits. For
245 pseudo-synchronous eccentrical orbits, it is equal to (Hut, 1981)
246
247 F = sqrt( (1+e)/ (1-e)^3)
248
249 Periastron is reached when d = 1-e.
250
251 @param r: radius of Roche volume at potential Phi (in units of semi-major axis)
252 @type r: float
253 @param theta: colatitude (0 at the pole, pi/2 at the equator)
254 @type theta: float
255 @param phi: longitude (0 in direction of COM)
256 @type phi: float
257 @param Phi: Roche potential value (unitless)
258 @type Phi: float
259 @param q: mass ratio
260 @type q: float
261 @param d: separation (in units of semi-major axis)
262 @type d: float
263 @param F: synchronicity parameter
264 @type F: float
265 @return: residu between Phi and roche potential
266 @rtype: float
267 """
268 lam,nu = cos(phi)*sin(theta),cos(theta)
269 term1 = 1. / r
270 term2 = q * ( 1./sqrt(d**2 - 2*lam*d*r + r**2) - lam*r/d**2)
271 term3 = 0.5 * F**2 * (q+1) * r**2 * (1-nu**2)
272 return (Phi - (term1 + term2 + term3))
273
275 """
276 Gradient of eccenctric asynchronous Roche potential in cartesian coordinates.
277
278 See Phoebe scientific reference, http://phoebe.fiz.uni-lj.si/docs/phoebe_science.ps.gz
279
280 x,y,z,d in real units! (otherwise you have to scale it yourself)
281
282 @param x: x-axis
283 @type x: float'
284 @param y: y-axis
285 @type y: float'
286 @param z: z-axis
287 @type z: float'
288 @param q: mass ratio
289 @type q: float
290 @param d: separation (in units of semi-major axis)
291 @type d: float
292 @param F: synchronicity parameter
293 @type F: float
294 @param norm: flag to return magnitude (True) or vector form (False)
295 @type norm: boolean
296 @return: Roche potential gradient
297 @rtype: ndarray or float
298 """
299 r = sqrt(x**2 + y**2 + z**2)
300 r_= sqrt((d-x)**2 + y**2 + z**2)
301 dOmega_dx = - x / r**3 + q * (d-x) / r_**3 + F**2 * (1+q)*x - q/d**2
302 dOmega_dy = - y / r**3 - q * y / r_**3 + F**2 * (1+q)*y
303 dOmega_dz = - z / r**3 - q * z / r_**3
304
305 dOmega = np.array([dOmega_dx,dOmega_dy,dOmega_dz])
306 if norm:
307 return vectors.norm(dOmega)
308 else:
309 return dOmega
310
311
313 """
314 Calculate surface gravity in an eccentric asynchronous binary roche potential.
315 """
316 q = M2/M1
317 x_com = q*d/(1+q)
318
319 r = np.array([x,y,z])
320 d_cf = np.array([d-x_com,0,0])
321 d = np.array([d,0,0])
322 h = d - r
323
324 term1 = - constants.GG*M1/vectors.norm(r)**3*r
325 term2 = - constants.GG*M2/vectors.norm(h)**3*h
326 term3 = - omega**2 * d_cf
327
328 g_pole = term1 + term2 + term3
329 if norm:
330 return vectors.norm(g_pole)
331 else:
332 return g_pole
333
334
336 """
337 Calculate the eccentric asynchronous binary Roche radius in spherical coordinates.
338
339 This is done via the Newton-Raphson secant method. If r_pole is not given
340 as a starting value, it will be calculated here (slowing down the function).
341
342 If no radius can be calculated for the given coordinates, 'nan' is returned.
343
344 @param theta: colatitude (0 at the pole, pi/2 at the equator)
345 @type theta: float
346 @param phi: longitude (0 in direction of COM)
347 @type phi: float
348 @param Phi: Roche potential value (unitless)
349 @type Phi: float
350 @param q: mass ratio
351 @type q: float
352 @param d: separation (in units of semi-major axis)
353 @type d: float
354 @param F: synchronicity parameter
355 @type F: float
356 @param r_pole: polar radius (serves as starting value for NR method)
357 @type r_pole: float
358 @return r: radius of Roche volume at potential Phi (in units of semi-major axis)
359 @rtype r: float
360 """
361 if r_pole is None:
362 r_pole = newton(binary_roche_potential,1e-5,args=(0,0,Phi,q,ds.min(),F))
363 try:
364 r = newton(binary_roche_potential,r_pole,args=(theta,phi,Phi,q,d,F))
365 if r<0 or r>d:
366 r = nan
367 except RuntimeError:
368 r = nan
369 return r
370
371
372
373
374
376
377
378 reflection_iter = 0
379 while (reflection_iter<max_iter):
380 R1,R2 = np.ones(len(primary['teff'])/4),np.ones(len(primary['teff'])/4)
381 for i in xrange(len(R1)):
382 if i%100==0: print i,len(R1)
383
384 s12 = np.array([ primary['x'][i]-secondary['x'],
385 primary['y'][i]-secondary['y'],
386 primary['z'][i]-secondary['z']])
387 psi2 = vectors.angle(+s12,-np.array([secondary['gravx'],secondary['gravy'],secondary['gravz']]))
388 psi1 = vectors.angle(-s12,-np.array([primary['gravx'][i:i+1],primary['gravy'][i:i+1],primary['gravz'][i:i+1]]))
389 s12 = s12.ravel()
390
391 keep = (psi2<pi/2.) & (psi1<pi/2.)
392 Lambda_2_bol = np.array([limbdark.get_itable(teff=iteff,logg=np.log10(igrav*100),theta=ipsi2,photbands=['OPEN.BOL'],absolute=False)[0]\
393 for iteff,igrav,ipsi2 in zip(secondary['teff'][keep],secondary['grav'][keep],psi2[keep])])
394
395 s = vectors.norm(s12[keep])
396 J_21_entrant = A1 * R2[i] * np.sum(secondary['flux'][keep] * cos(psi1[keep])*cos(psi2[keep])*Lambda_2_bol*secondary['areas'][keep]/s**2)
397 if np.isnan(J_21_entrant).any(): raise ValueError
398 J_1_tot = primary['flux'][i]
399 R1_new = 1+ J_21_entrant/J_1_tot
400
401
402 s12 = np.array([primary['x']-secondary['x'][i],
403 primary['y']-secondary['y'][i],
404 primary['z']-secondary['z'][i]])
405 psi2 = vectors.angle(+s12,-np.array([primary['gravx'],primary['gravy'],primary['gravz']]))
406 psi1 = vectors.angle(-s12,-np.array([secondary['gravx'][i:i+1],secondary['gravy'][i:i+1],secondary['gravz'][i:i+1]]))
407 s12 = s12.ravel()
408
409 keep = (psi2<pi/2.) & (psi1<pi/2.)
410 Lambda_2_bol = np.array([limbdark.get_itable(teff=iteff,logg=np.log10(igrav*100),theta=ipsi2,photbands=['OPEN.BOL'],absolute=False)[0]\
411 for iteff,igrav,ipsi2 in zip(primary['teff'][keep],primary['grav'][keep],psi2[keep])])
412
413 s = vectors.norm(s12)
414 J_12_entrant = A2 * R1[i] * np.nansum(primary['flux'][keep] * cos(psi1[keep])*cos(psi2[keep])*Lambda_2_bol*primary['areas'][keep]/s**2)
415 J_2_tot = secondary['flux'][i]
416 R2_new = 1+ J_12_entrant/J_2_tot
417
418 R1[i] = R1_new
419 R2[i] = R2_new
420
421
422
423
424
425
426
427
428
429
430
431
432 break_out = True
433 trash,trash2,R1,R2 = stitch_grid(theta,phi,R1.reshape(theta.shape),R2.reshape(theta.shape))
434 del trash,trash2
435 R1 = R1.ravel()
436 R2 = R2.ravel()
437
438 if (R1[-np.isnan(R1)]>1.05).any():
439 print "Significant reflection effect on primary (max %.3f%%)"%((R1.max()**0.25-1)*100)
440 primary['teff']*= R1**0.25
441 primary['flux'] = local_intensity(primary['teff'],primary['grav'],np.ones_like(primary['teff']),photbands=['OPEN.BOL'])
442 break_out = False
443 else:
444 print 'Maximum reflection effect on primary: %.3f%%'%((R1.max()**0.25-1)*100)
445
446 if (R2[-np.isnan(R2)]>1.05).any():
447 print "Significant reflection effect on secondary (max %.3g%%)"%((R2.max()**0.25-1)*100)
448 secondary['teff']*= R2**0.25
449 secondary['flux'] = local_intensity(secondary['teff'],secondary['grav'],np.ones_like(secondary['teff']),photbands=['OPEN.BOL'])
450 break_out = False
451 else:
452 print 'Maximum reflection effect on secondary: %.3g%%'%((R1.max()**0.25-1)*100)
453
454 if break_out:
455 break
456
457
458 reflection_iter += 1
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480 return primary,secondary
481
482
483
485 """
486 Generate a synthetic spectrum of one or more stars.
487
488 If you give more than one star, you get more than one synthesized spectrum
489 back. To compute the total spectrum, just add them up.
490
491 WARNING: the spectra are scaled with the value of the projected intensity.
492 This is usually calculated within a specific photometric passband. The total
493 spectrum will thus be dependent on the photometric passband, unless you
494 specified 'OPEN.BOL' (which is the bolometric open filter). If you want to
495 know the RV for a specific line, you might want to look into calculating the
496 projected intensities at specific wavelength intervals.
497
498 @param stars: any number of (projected) star record arrays
499 @type stars: tuple of star record arrays
500 @keyword wave: wavelength template to compute the synthetic spectrum on
501 @type wave: numpy array
502 @return: tuple of synthetic spectra, scaled to total intensity
503 @rtype: tuple of numpy arrays
504 """
505 wave = kwargs.get('wave',np.linspace(7050,7080,1000))
506
507 spectra = []
508 get_spectrum = spectra_model.get_table
509 for mystar in stars:
510 loggs = np.log10(np.sqrt(mystar['gravx']**2 + mystar['gravy']**2 + mystar['gravz']**2)*100)
511 iterator = zip(mystar['teff'],loggs,mystar['vx']/1000.)
512 print "Temperature range:",mystar['teff'].min(),mystar['teff'].max()
513 print "logg range:",loggs.min(),loggs.max()
514
515
516 spectra_elements = np.array([get_spectrum(teff=iteff,logg=ilogg,vrad=ivrad,wave=wave)[1] for iteff,ilogg,ivrad in iterator])
517
518 average_spectrum = np.sum(spectra_elements*mystar['projflux'],axis=0)
519 spectra.append(average_spectrum)
520 return spectra
521
522
524 """
525 Generate a synthetic light curve of a binary system.
526
527 @keyword name: name of the binary system, used for output
528 @type name: string
529 @keyword direc: directory to put output files
530 @type direc: string (if empty string, then current directory will be taken,
531 if None, then not output will be written)
532 @parameter gres: grid resolution for primary and secondary. If integer, the
533 resolution of both components and both longitude and latitude will be the
534 same. A 2-tuple will be interpreted as (latitude resolution, longitude resolution)
535 @type gres: integer, 2-tuple or 4-tuple
536 @parameter tres: number of phase steps to comptue the light curve on
537 @type tres: integer
538 """
539
540
541 name = parameters.setdefault('name','mybinary')
542 direc = parameters.setdefault('direc','')
543
544 res = parameters.setdefault('gres',20)
545 gtype = parameters.setdefault('gtype','spher')
546 tres= parameters.pop('tres',125)
547 photband = parameters.setdefault('photband','JOHNSON.V')
548 max_iter_reflection = parameters.setdefault('ref_iter',1)
549
550 gamma = parameters.setdefault('gamma',0.)
551 incl = parameters.setdefault('incl',90.)
552 q = parameters.setdefault('q',1.)
553 e = parameters.setdefault('e',0.)
554 omega = parameters.setdefault('omega',0.)
555
556 F = parameters.setdefault('F1',sqrt((1+e)/(1-e)**3))
557 F2= parameters.setdefault('F2',sqrt((1+e)/(1-e)**3))
558 A1= parameters.setdefault('A1',1.)
559 A2= parameters.setdefault('A2',1.)
560 beta1 = parameters.setdefault('beta1',1.0)
561 beta2 = parameters.setdefault('beta2',1.0)
562
563
564 T_pole = parameters['Tpole1']
565 T_pole2 = parameters['Tpole2']
566 P = parameters['P']
567 asini = parameters['asini']
568 Phi = parameters['Phi1']
569 Phi2 = parameters['Phi2']
570
571
572 a = asini/sin(incl/180.*pi)
573 a1 = a / (1.+1./q)
574 a2 = a / (1.+ q)
575 q2 = 1./q
576 view_angle = incl/180.*pi
577
578 a_ = conversions.convert('au','m',a)
579 P_ = conversions.convert('d','s',P)
580 M = (4*pi**2 * a_**3 / P_**2 / constants.GG ) / constants.Msol
581 M1 = M / (1.+ q)
582 M2 = M / (1.+1./q)
583
584 mu = M2 / (M1+M2)
585 z = (mu/3.)**(1./3.)
586 L1 = z- 1/3.*z**2 - 1/9.*z**3 + 58./81*z**4
587 L2 = z+ 1/3.*z**2 - 1/9.*z**3 + 50./81*z**4
588 parameters['M1'] = M1
589 parameters['M2'] = M2
590 parameters['L1'] = L1*a*constants.au/constants.Rsol
591 parameters['L2'] = L2*a*constants.au/constants.Rsol
592
593
594 if not hasattr(tres,'__iter__'):
595 times = np.linspace(0.5*P,1.5*P,tres)
596 else:
597 times = tres
598 r1,theta1 = keplerorbit.orbit_in_plane(times,[P,e,a1,gamma],component='primary')
599 r2,theta2 = keplerorbit.orbit_in_plane(times,[P,e,a2,gamma],component='secondary')
600 RV1 = keplerorbit.radial_velocity([P,e,a1,pi/2,view_angle,gamma],theta=theta1,itermax=8)
601 RV2 = keplerorbit.radial_velocity([P,e,a2,pi/2,view_angle,gamma],theta=theta2,itermax=8)
602
603 r1 = r1/constants.au
604 r2 = r2/constants.au
605
606 x1o,y1o = r1*cos(theta1)/a,r1*sin(theta1)/a
607 x2o,y2o = r2*cos(theta2)/a,r2*sin(theta2)/a
608
609 rot_i = -(pi/2 - view_angle)
610 x1o_,z1o = vectors.rotate(x1o,np.zeros_like(y1o),rot_i)
611 x2o_,z2o = vectors.rotate(x2o,np.zeros_like(y2o),rot_i)
612
613
614 ds = np.sqrt( (x1o-x2o)**2 + (y1o-y2o)**2)
615
616
617 r_pole = newton(binary_roche_potential,1e-5,args=(0,0,Phi,q,ds.min(),F))
618 r_pole2 = newton(binary_roche_potential,1e-5,args=(0,0,Phi2,q2,ds.min(),F2))
619 parameters['Rp1'] = r_pole*a*constants.au/constants.Rsol
620 parameters['Rp2'] = r_pole2*a*constants.au/constants.Rsol
621
622
623 phi1_crit = binary_roche_potential(L1,0,0,Phi,q,ds.min(),F)
624 phi2_crit = binary_roche_potential(L1,0,0,Phi2,q2,ds.min(),F2)
625 R_roche1 = a * 0.49 * q2**(2./3.) / (0.6*q2**(2./3.) + np.log(1+q2**(1./3.)))
626 R_roche2 = a * 0.49 * q**(2./3.) / (0.6*q**(2./3.) + np.log(1+q**(1./3.)))
627
628
629 tdyn1 = np.sqrt(2*(r_pole*constants.au)**3/(constants.GG*M1*constants.Msol))
630 tdyn2 = np.sqrt(2*(r_pole2*constants.au)**3/(constants.GG*M2*constants.Msol))
631 lumt1 = (M1<2.) and M1**4 or (1.5*M1**3.5)
632 lumt2 = (M2<2.) and M2**4 or (1.5*M2**3.5)
633 tthr1 = constants.GG * (M1*constants.Msol)**2 / (r_pole*constants.au*lumt1*constants.Lsol)
634 tthr2 = constants.GG * (M2*constants.Msol)**2 / (r_pole2*constants.au*lumt2*constants.Lsol)
635 tnuc1 = 7e9 * M1 / lumt1
636 tnuc2 = 7e9 * M2 / lumt2
637
638 logger.info('=================================================')
639 logger.info('GENERAL SYSTEM AND COMPONENT PROPERTIES')
640 logger.info('=================================================')
641 logger.info("Period P = %.3g d"%(P))
642 logger.info("Inclination angle i = %.3g deg"%(incl))
643 logger.info("Systemic velocity gamma = %.3g km/s"%(gamma))
644 logger.info("Eccentricity e = %.3g"%(e))
645 logger.info("Mass ratio q = %.3g"%(q))
646 logger.info("Mass primary M1 = %.3g Msun"%(M1))
647 logger.info("Mass secondary M2 = %.3g Msun"%(M2))
648 logger.info("Polar Radius primary R1 = %.3g Rsun (Roche = %.3g)"%(r_pole*a*constants.au/constants.Rsol,R_roche1*constants.au/constants.Rsol))
649 logger.info("Polar Radius secondary R2 = %.3g Rsun (Roche = %.3g)"%(r_pole2*a*constants.au/constants.Rsol,R_roche2*constants.au/constants.Rsol))
650 logger.info("Center-of-mass = %.3g Rsun"%(q/(1+q)*a*constants.au/constants.Rsol))
651 logger.info('Lagrange L1 = %.3g Rsun'%(L1*a*constants.au/constants.Rsol))
652 logger.info('Lagrange L2 = %.3g Rsun'%(L2*a*constants.au/constants.Rsol))
653 logger.info('Phi primary = %.3g (critical = %.3g)'%(Phi,phi1_crit))
654 logger.info('Phi secondary = %.3g (critical = %.3g)'%(Phi2,phi2_crit))
655 logger.info('Luminosity primary L1 (MS) = %.3g Lsun'%(lumt1))
656 logger.info('Luminosity secondary L2 (MS) = %.3g Lsun'%(lumt2))
657 logger.info('System semi-major axis a = %.3g Rsun'%(a*constants.au/constants.Rsol))
658 logger.info('------')
659 logger.info('TDYN primary = %.3g hr'%(tdyn1/3600.))
660 logger.info('TDYN secondary = %.3g hr'%(tdyn2/3600.))
661 logger.info('TNUC primary = %.3g yr'%(tnuc1))
662 logger.info('TNUC secondary = %.3g yr'%(tnuc2))
663 logger.info('TTHERM primary = %.3g yr'%(tthr1/(24*3600*365)))
664 logger.info('TTHERM secondary = %.3g yr'%(tthr2/(24*3600*365)))
665 logger.info('=================================================')
666
667
668 if hasattr(res,'__iter__'):
669 mygrid = local.get_grid(res[0],res[1],gtype=gtype)
670 theta,phi = mygrid[:2]
671 else:
672 mygrid = local.get_grid(res,gtype=gtype)
673 theta,phi = mygrid[:2]
674 thetas,phis = np.ravel(theta),np.ravel(phi)
675
676 light_curve = np.zeros_like(times)
677 RV1_corr = np.zeros_like(times)
678 RV2_corr = np.zeros_like(times)
679 to_SI = a*constants.au
680 to_CGS = a*constants.au*100.
681 scale_factor = a*constants.au/constants.Rsol
682
683 fitsfile = os.path.join(direc,'%s.fits'%(name))
684 if direc is not None and os.path.isfile(fitsfile):
685 os.remove(fitsfile)
686 logger.warning("Removed existing file %s"%(fitsfile))
687 if direc is not None:
688 parameters['scalefac'] = a*constants.au/constants.Rsol
689 parameters.pop('gres')
690 outputfile_prim = os.path.join(direc,'%s_primary.fits'%(name))
691 outputfile_secn = os.path.join(direc,'%s_secondary.fits'%(name))
692 outputfile_prim = fits.write_primary(outputfile_prim,header_dict=parameters)
693 outputfile_secn = fits.write_primary(outputfile_secn,header_dict=parameters)
694
695 ext_dict = {}
696 for di,d in enumerate(ds):
697 report = "STEP %04d"%(di)
698
699
700
701 omega_rot = F * 2*pi/P_ * 1/d**2 * sqrt( (1+e)*(1-e))
702 omega_rot_vec = np.array([0.,0.,-omega_rot])
703
704 if e>0 or di==0:
705
706 out = [[get_binary_roche_radius(itheta,iphi,Phi=Phi,q=q,d=d,F=F,r_pole=r_pole),
707 get_binary_roche_radius(itheta,iphi,Phi=Phi2,q=q2,d=d,F=F2,r_pole=r_pole2)] for itheta,iphi in zip(thetas,phis)]
708 rprim,rsec = np.array(out).T
709
710
711
712 radius = rprim.reshape(theta.shape)
713 this_r_pole = get_binary_roche_radius(0,0,Phi=Phi,q=q,d=d,F=F,r_pole=r_pole)
714 x,y,z = vectors.spher2cart_coord(radius,phi,theta)
715 g_pole = binary_roche_surface_gravity(0,0,this_r_pole*to_SI,d*to_SI,omega_rot,M1*constants.Msol,M2*constants.Msol,norm=True)
716 Gamma_pole = binary_roche_potential_gradient(0,0,this_r_pole,q,d,F,norm=True)
717 zeta = g_pole / Gamma_pole
718 dOmega = binary_roche_potential_gradient(x,y,z,q,d,F,norm=False)
719 grav_local = dOmega*zeta
720
721
722
723 grav_local = np.array([i.reshape(theta.shape) for i in grav_local])
724 grav = vectors.norm(grav_local)
725 areas_local,cos_gamma = local.surface_elements((radius,mygrid),-grav_local,gtype=gtype)
726 teff_local = local.temperature(grav,g_pole,T_pole,beta=beta1)
727 ints_local = local.intensity(teff_local,grav,np.ones_like(cos_gamma),photband='OPEN.BOL')
728 velo_local = np.cross(np.array([x,y,z]).T*to_SI,omega_rot_vec).T
729
730
731
732 lumi_prim = 4*pi*(ints_local*areas_local*to_CGS**2).sum()/constants.Lsol_cgs
733 area_prim = 4*areas_local.sum()*to_CGS**2/(4*pi*constants.Rsol_cgs**2)
734 logger.info('----PRIMARY DERIVED PROPERTIES')
735 logger.info('Polar Radius primary = %.3g Rsun'%(this_r_pole*a*constants.au/constants.Rsol))
736 logger.info("Polar logg primary = %.3g dex"%(np.log10(g_pole*100)))
737 logger.info("Luminosity primary = %.3g Lsun"%(lumi_prim))
738 logger.info("Surface area primary = %.3g Asun"%(area_prim))
739 logger.info("Mean Temp primary = %.3g K"%(np.average(teff_local,weights=areas_local)))
740 ext_dict['Rp1'] = this_r_pole*a*constants.au/constants.Rsol
741 ext_dict['loggp1'] = np.log10(g_pole*100)
742 ext_dict['LUMI1'] = lumi_prim
743 ext_dict['SURF1'] = area_prim
744
745
746
747 radius2 = rsec.reshape(theta.shape)
748 this_r_pole2 = get_binary_roche_radius(itheta,iphi,Phi=Phi2,q=q2,d=d,F=F2,r_pole=r_pole2)
749 x2,y2,z2 = vectors.spher2cart_coord(radius2,phi,theta)
750 g_pole2 = binary_roche_surface_gravity(0,0,this_r_pole2*to_SI,d*to_SI,omega_rot,M2*constants.Msol,M1*constants.Msol,norm=True)
751 Gamma_pole2 = binary_roche_potential_gradient(0,0,this_r_pole2,q2,d,F2,norm=True)
752 zeta2 = g_pole2 / Gamma_pole2
753 dOmega2 = binary_roche_potential_gradient(x2,y2,z2,q2,d,F2,norm=False)
754 grav_local2 = dOmega2*zeta2
755
756
757
758 grav_local2 = np.array([i.reshape(theta.shape) for i in grav_local2])
759 grav2 = vectors.norm(grav_local2)
760 areas_local2,cos_gamma2 = local.surface_elements((radius2,mygrid),-grav_local2,gtype=gtype)
761 teff_local2 = local.temperature(grav2,g_pole2,T_pole2,beta=beta2)
762 ints_local2 = local.intensity(teff_local2,grav2,np.ones_like(cos_gamma2),photband='OPEN.BOL')
763 velo_local2 = np.cross(np.array([x2,y2,z2]).T*to_SI,omega_rot_vec).T
764
765
766
767 lumi_sec = 4*pi*(ints_local2*areas_local2*to_CGS**2).sum()/constants.Lsol_cgs
768 area_sec = 4*areas_local2.sum()*to_CGS**2/(4*pi*constants.Rsol_cgs**2)
769 logger.info('----SECONDARY DERIVED PROPERTIES')
770 logger.info('Polar Radius secondary = %.3g Rsun'%(this_r_pole2*a*constants.au/constants.Rsol))
771 logger.info("Polar logg secondary = %.3g dex"%(np.log10(g_pole2*100)))
772 logger.info("Luminosity secondary = %.3g Lsun"%(lumi_sec))
773 logger.info("Surface area secondary = %.3g Asun"%(area_sec))
774 logger.info("Mean Temp secondary = %.3g K"%(np.average(teff_local2,weights=areas_local2)))
775 ext_dict['Rp2'] = this_r_pole2*a*constants.au/constants.Rsol
776 ext_dict['loggp2'] = np.log10(g_pole2*100)
777 ext_dict['LUMI2'] = lumi_sec
778 ext_dict['SURF2'] = area_sec
779
780
781
782
783
784
785
786
787 theta_,phi_,radius,gravx,gravy,gravz,grav,areas,teff,ints,vx,vy,vz = \
788 local.stitch_grid(theta,phi,radius,grav_local[0],grav_local[1],grav_local[2],
789 grav,areas_local,teff_local,ints_local,velo_local[0],velo_local[1],velo_local[2],
790 seamless=False,gtype=gtype,
791 vtype=['scalar','x','y','z','scalar','scalar','scalar','scalar','vx','vy','vz'])
792
793 theta2_,phi2_,radius2,gravx2,gravy2,gravz2,grav2,areas2,teff2,ints2,vx2,vy2,vz2 = \
794 local.stitch_grid(theta,phi,radius2,grav_local2[0],grav_local2[1],grav_local2[2],
795 grav2,areas_local2,teff_local2,ints_local2,velo_local2[0],velo_local2[1],velo_local2[2],
796 seamless=False,gtype=gtype,
797 vtype=['scalar','x','y','z','scalar','scalar','scalar','scalar','vx','vy','vz'])
798
799
800 x_of,y_of,z_of = vectors.spher2cart_coord(radius.ravel(),phi_.ravel(),theta_.ravel())
801 x2_of,y2_of,z2_of = vectors.spher2cart_coord(radius2.ravel(),phi2_.ravel(),theta2_.ravel())
802 x2_of = -x2_of
803
804 primary = np.rec.fromarrays([theta_.ravel(),phi_.ravel(),radius.ravel(),
805 x_of,y_of,z_of,
806 vx.ravel(),vy.ravel(),vz.ravel(),
807 gravx.ravel(),gravy.ravel(),gravz.ravel(),grav.ravel(),
808 areas.ravel(),teff.ravel(),ints.ravel()],
809 names=['theta','phi','r',
810 'x','y','z',
811 'vx','vy','vz',
812 'gravx','gravy','gravz','grav',
813 'areas','teff','flux'])
814
815 secondary = np.rec.fromarrays([theta2_.ravel(),phi2_.ravel(),radius2.ravel(),
816 x2_of,y2_of,z2_of,
817 vx2.ravel(),-vy2.ravel(),vz2.ravel(),
818 -gravx2.ravel(),gravy2.ravel(),gravz2.ravel(),grav2.ravel(),
819 areas2.ravel(),teff2.ravel(),ints2.ravel()],
820 names=['theta','phi','r',
821 'x','y','z',
822 'vx','vy','vz',
823 'gravx','gravy','gravz','grav',
824 'areas','teff','flux'])
825
826
827 primary,secondary = reflection_effect(primary,secondary,theta,phi,
828 A1=A1,A2=A2,max_iter=max_iter_reflection)
829
830
831
832 rot_theta = np.arctan2(y1o[di],x1o[di])
833
834
835 if direc is not None:
836 prim = local.project(primary,view_long=(rot_theta,x1o[di],y1o[di]),
837 view_lat=(view_angle,0,0),photband=photband,
838 only_visible=False,plot_sort=False,scale_factor=scale_factor)
839 secn = local.project(secondary,view_long=(rot_theta,x2o[di],y2o[di]),
840 view_lat=(view_angle,0,0),photband=photband,
841 only_visible=False,plot_sort=False,scale_factor=scale_factor)
842
843 com_x = (x1o[di] + q*x2o[di]) / (1.0+q)
844 com_y = (y1o[di] + q*y2o[di]) / (1.0+q)
845 com = np.array([com_x,com_y,0.])
846 rot_i = -(pi/2 - view_angle)
847 com[0],com[1] = vectors.rotate(com[0],com[1],rot_theta,x0=x1o[di],y0=y1o[di])
848 com[0],com[2] = vectors.rotate(com[0],com[2],rot_i)
849 prim_header = dict(x0=x1o[di],y0=y1o[di],i=view_angle,comx=com[0],comy=com[1],com_z=com[2],nr=di,time=times[di])
850 secn_header = dict(x0=x2o[di],y0=y2o[di],i=view_angle,comx=com[0],comy=com[1],com_z=com[2],nr=di,time=times[di])
851
852 if direc is not None and (di%20==0): close = True
853 else: close = False
854
855 outputfile_prim = fits.write_recarray(prim,outputfile_prim,close=close,header_dict=prim_header)
856 outputfile_secn = fits.write_recarray(secn,outputfile_secn,close=close,header_dict=secn_header)
857
858 prim = local.project(primary,view_long=(rot_theta,x1o[di],y1o[di]),
859 view_lat=(view_angle,0,0),photband=photband,
860 only_visible=True,plot_sort=True)
861 secn = local.project(secondary,view_long=(rot_theta,x2o[di],y2o[di]),
862 view_lat=(view_angle,0,0),photband=photband,
863 only_visible=True,plot_sort=True)
864 prim['vx'] = -prim['vx'] + RV1[di]*1000.
865 secn['vx'] = -secn['vx'] + RV2[di]*1000.
866
867
868
869
870
871
872
873 if secn['x'].min()<prim['x'].min():
874 front,back = prim,secn
875 front_component = 1
876 report += ' Primary in front'
877 else:
878 front,back = secn,prim
879 front_component = 2
880 report += ' Secondary in front'
881 coords_front = np.column_stack([front['y'],front['z']])
882 coords_back = np.column_stack([back['y'],back['z']])
883
884 if gtype!='delaunay':
885
886
887 tree = KDTree(coords_front)
888 distance,order = tree.query(coords_back)
889
890
891
892
893 in_eclipse = distance < np.sqrt(front['areas'][order])
894 else:
895
896
897 eclipse_detection = Delaunay(coords_front)
898 in_eclipse = eclipse_detection.find_simplex(coords_back)>=0
899 if np.sum(in_eclipse)>0:
900 report += ' during eclipse'
901 else:
902 report += ' outside eclipse'
903
904
905
906 total_intensity = front['projflux'].sum() + back['projflux'][-in_eclipse].sum()
907 light_curve[di] = total_intensity
908 report += "---> Total intensity: %g "%(total_intensity)
909
910 if di==0:
911 ylim_lc = (0.95*min(prim['projflux'].sum(),secn['projflux'].sum()),1.2*(prim['projflux'].sum()+secn['projflux'].sum()))
912 back['projflux'][in_eclipse] = 0
913 back['eyeflux'][in_eclipse] = 0
914 back['vx'][in_eclipse] = 0
915 back['vy'][in_eclipse] = 0
916 back['vz'][in_eclipse] = 0
917
918
919 if front_component==1:
920 RV1_corr[di] = np.average(front['vx']/1000.,weights=front['projflux'])
921 RV2_corr[di] = np.average(back['vx'][-in_eclipse]/1000.,weights=back['projflux'][-in_eclipse])
922 front_cmap = pl.cm.hot
923 back_cmap = pl.cm.cool_r
924 else:
925 RV2_corr[di] = np.average(front['vx']/1000.,weights=front['projflux'])
926 RV1_corr[di] = np.average(back['vx'][-in_eclipse]/1000.,weights=back['projflux'][-in_eclipse])
927 front_cmap = pl.cm.cool_r
928 back_cmap = pl.cm.hot
929 report += 'RV1=%.3f, RV2=%.3f'%(RV1_corr[di],RV2_corr[di])
930 logger.info(report)
931
932
933 if direc is not None:
934
935 if di==0:
936 size_x = 1.2*(max(prim['y'].ptp(),secn['y'].ptp())/2. + max(ds))
937 size_y = 1.2*(max(prim['z'].ptp(),secn['z'].ptp())/2. + max(ds) * cos(view_angle))
938 vmin_image,vmax_image = 0,max([front['eyeflux'].max(),back['eyeflux'].max()])
939 size_top = 1.2*(max(prim['x'].ptp(),secn['x'].ptp())/2. + max(ds))
940
941 pl.figure(figsize=(16,11))
942 pl.subplot(221,aspect='equal');pl.title('line of sight intensity')
943 pl.scatter(back['y'],back['z'],c=back['eyeflux'],edgecolors='none',cmap=back_cmap)
944 pl.scatter(front['y'],front['z'],c=front['eyeflux'],edgecolors='none',cmap=front_cmap)
945 pl.xlim(-size_x,size_x)
946 pl.ylim(-size_y,size_y)
947 pl.xlabel('X [semi-major axis]')
948 pl.ylabel('Z [semi-major axis]')
949
950
951 pl.subplot(222,aspect='equal');pl.title('line of sight velocity')
952 if di==0:
953 vmin_rv = min((prim['vx'].min()/1000.+RV1.min()),(prim['vx'].min()/1000.+RV2.min()))
954 vmax_rv = max((secn['vx'].max()/1000.+RV2.max()),(secn['vx'].max()/1000.+RV2.max()))
955 pl.scatter(front['y'],front['z'],c=front['vx']/1000.,edgecolors='none',cmap=pl.cm.RdBu_r,vmin=vmin_rv,vmax=vmax_rv)
956 pl.scatter(back['y'],back['z'],c=back['vx']/1000.,edgecolors='none',cmap=pl.cm.RdBu_r,vmin=vmin_rv,vmax=vmax_rv)
957 cbar = pl.colorbar()
958 cbar.set_label('Radial velocity [km/s]')
959 pl.xlim(-size_x,size_x)
960 pl.ylim(-size_y,size_y)
961 pl.xlabel('X [semi-major axis]')
962 pl.ylabel('Z [semi-major axis]')
963
964
965 pl.subplot(223,aspect='equal');pl.title('Top view')
966 pl.scatter(prim['x'],prim['y'],c=prim['eyeflux'],edgecolors='none',cmap=pl.cm.hot)
967 pl.scatter(secn['x'],secn['y'],c=secn['eyeflux'],edgecolors='none',cmap=pl.cm.cool_r)
968 pl.xlim(-size_top,size_top)
969 pl.ylim(-size_top,size_top)
970 pl.xlabel('X [semi-major axis]')
971 pl.ylabel('Y [semi-major axis]')
972
973
974 pl.subplot(224);pl.title('light curve and RV curve')
975 pl.plot(times[:di+1],2.5*np.log10(light_curve[:di+1]),'k-',label=photband)
976 pl.plot(times[di],2.5*np.log10(light_curve[di]),'ko',ms=10)
977 pl.xlim(times.min(),times.max())
978 pl.ylim(2.5*np.log10(ylim_lc[0]),2.5*np.log10(ylim_lc[1]))
979 pl.ylabel('Flux [erg/s/cm2/A/sr]')
980 pl.legend(loc='lower left',prop=dict(size='small'))
981
982 pl.twinx(pl.gca())
983
984 pl.plot(times[:di+1],RV1_corr[:di+1],'b-',lw=2,label='Numeric 1')
985 pl.plot(times[:di+1],RV1[:di+1]-gamma,'g--',lw=2,label='Kepler 1')
986 pl.plot(times[di],RV1[di]-gamma,'go',ms=10)
987 pl.plot(times[di],RV1_corr[di],'bo',ms=10)
988
989 pl.plot(times[:di+1],RV2_corr[:di+1],'r-',lw=2,label='Numeric 2')
990 pl.plot(times[:di+1],RV2[:di+1]-gamma,'c--',lw=2,label='Kepler 2')
991 pl.plot(times[di],RV2[di]-gamma,'cs',ms=10)
992 pl.plot(times[di],RV2_corr[di],'rs',ms=10)
993 pl.ylim(1.2*min(min(RV1-gamma),min(RV2-gamma)),+1.2*max(max(RV1-gamma),max(RV2-gamma)))
994 pl.xlim(times.min(),times.max())
995 pl.ylabel('Radial velocity [km/s]')
996 pl.xlabel('Time [d]')
997 pl.legend(loc='upper left',prop=dict(size='small'))
998 pl.savefig(os.path.join(direc,'%s_los_%04d'%(name,di)),facecolor='0.75')
999 pl.close()
1000
1001
1002 pl.figure(figsize=(7,size_y/size_x*7))
1003 ax = pl.axes([0,0,1,1])
1004 ax.set_aspect('equal')
1005 ax.set_axis_bgcolor('k')
1006 pl.xticks([]);pl.yticks([])
1007 if front_component==1: sfront,sback = 16,9
1008 else: sfront,sback = 9,16
1009 pl.scatter(back['y'][-in_eclipse],back['z'][-in_eclipse],s=sback,c=back['eyeflux'][-in_eclipse],edgecolors='none',cmap=pl.cm.gray,vmin=vmin_image,vmax=vmax_image)
1010 pl.scatter(front['y'],front['z'],s=sfront,c=front['eyeflux'],edgecolors='none',cmap=pl.cm.gray,vmin=vmin_image,vmax=vmax_image)
1011 pl.xlim(-size_x,size_x);pl.ylim(-size_y,size_y)
1012 pl.savefig(os.path.join(direc,'%s_image_%04d'%(name,di)),facecolor='k')
1013 pl.close()
1014
1015
1016
1017 if direc is not None:
1018 outputfile_prim.close()
1019 outputfile_secn.close()
1020 return times, light_curve, RV1_corr, RV2_corr
1021
1022 if __name__=="__main__":
1023 import doctest
1024 doctest.testmod()
1025 pl.show()
1026
1027 logger = loggers.get_basic_logger("")
1028 import sys
1029