1
2 """
3 Fit and evaluate binary (Keplerian) orbits.
4
5 Be careful with the inclination angle: if you only want an image projected onto
6 the sky, it is not case-sensitive. Only for the radial velocities this becomes
7 important (the primary and secondary can be mixed).
8
9 Example usage: compute the absolute orbits of the two components of
10 mu Cassiopeia, and compute the relative orbit of the second component.
11
12 See ps file on http://www.chara.gsu.edu/~gudehus/binary.html and Drummond et al.,
13 1995 (the latter's orbital parameters are used).
14
15 Neccesary imports:
16
17 >>> from ivs.units.constants import *
18
19 Orbital elements, Euler angles and masses of the system
20
21 >>> P = 21.753*365 # period in days
22 >>> e = 0.561 # ecc
23 >>> omega,Omega,i = 332.7, 47.3, 106.8 # degrees
24 >>> T0 = 1975.738*365 # year in days
25 >>> plx = 133.2 # mas
26 >>> M1,M2 = 0.719,0.168 # Msun
27 >>> times = np.linspace(T0,T0+0.95*P,100)
28
29 Derive semi-major axis of absolute orbits (semi-major axis of relative orbit is a+a2)
30
31 >>> fA = M2**3/ (M1+M2)**2 * Msol
32 >>> a1 = ((P*24*3600)**2 * GG * fA / (4*np.pi**2))**(1/3.)
33 >>> a1 = a1/au
34 >>> a2 = a1*M1/M2
35
36 Convert the angles to radians
37
38 >>> omega,Omega,i = omega/180*np.pi,Omega/180*np.pi,i/180*np.pi
39
40 Compute the orbits of the primary and secondary in the plane of the sky
41
42 >>> x1,y1,z1 = orbit_on_sky(times,(P,e,a1,T0,omega,Omega,i),component='prim',distance=(1000./133.,'pc'))
43 >>> x2,y2,z2 = orbit_on_sky(times,(P,e,a2,T0,omega,Omega,i),component='sec',distance=(1000./133.,'pc'))
44
45 Plot it:
46
47 >>> import pylab as pl
48 >>> colors = times/365.
49 >>> p = pl.title('mu Cassiopeia')
50 >>> p = pl.scatter(x1,y1,c=colors,edgecolors='none',cmap=pl.cm.spectral,label='_nolegend_')
51 >>> p = pl.scatter(x2,y2,c=colors,edgecolors='none',cmap=pl.cm.spectral,label='_nolegend_')
52 >>> p = pl.scatter(x2-x1,y2-y1,c=colors,edgecolors='none',cmap=pl.cm.spectral,label='_nolegend_')
53 >>> p = pl.plot(x1,y1,'b-',label='Primary')
54 >>> p = pl.plot(x2,y2,'g-',label='Secondary')
55 >>> p = pl.plot(x2-x1,y2-y1,'r-',label='Secondary relative')
56 >>> p = pl.colorbar()
57 >>> p = pl.legend(loc='lower right')
58 >>> p = pl.grid()
59 >>> p = pl.xlabel(r'Angular size (arcsec) --> N')
60 >>> p = pl.ylabel(r'Angular size (arcsec) --> E')
61
62 ]include figure]]ivs_binary_keplerorbit_muCas.png]
63
64 Example usage 2: compute the absolute orbits of the two components of
65 Capella, and compute the relative orbit of the primary component.
66
67 See Hummel et al. 1994 (?).
68
69 Neccesary imports:
70
71 >>> from ivs.units.constants import *
72
73 Orbital elements, Euler angles and masses of the system
74
75 >>> P = 104.022 # period in days
76 >>> e = 0.000 #ecc
77 >>> i = 137.18 # degree
78 >>> T0 = 2447528.45 # year
79 >>> omega = 0 # degrees
80 >>> Omega = 40.8 # degrees
81 >>> plx = 76.20 # mas
82 >>> M1 = 2.69 #Msun
83 >>> M2 = 2.56 #Msun
84 >>> RV0 = 0.
85 >>> times = np.linspace(T0,T0+0.95*P,100)
86
87 Derive semi-major axis of absolute orbits (semi-major axis of relative orbit is a+a2)
88
89 >>> fA = M2**3/ (M1+M2)**2 * Msol
90 >>> a1 = ((P*24*3600)**2 * GG * fA / (4*np.pi**2))**(1/3.)
91 >>> a1 = a1/au
92 >>> a2 = a1*M1/M2
93
94 Convert the angles to radians
95
96 >>> omega,Omega,i = omega/180*np.pi,Omega/180*np.pi,i/180*np.pi
97
98 Compute the orbits of the primary and secondary in the plane of the sky
99
100 >>> x1,y1,z1 = orbit_on_sky(times,(P,e,a1,T0,omega,Omega,i),component='prim',distance=(1000./plx,'pc'))
101 >>> x2,y2,z2 = orbit_on_sky(times,(P,e,a2,T0,omega,Omega,i),component='sec',distance=(1000./plx,'pc'))
102
103 Plot it:
104
105 >>> import pylab as pl
106 >>> colors = times/365.
107 >>> p = pl.figure()
108 >>> p = pl.title('Capella')
109 >>> p = pl.scatter(x1,y1,c=colors,edgecolors='none',cmap=pl.cm.spectral,label='_nolegend_')
110 >>> p = pl.scatter(x2,y2,c=colors,edgecolors='none',cmap=pl.cm.spectral,label='_nolegend_')
111 >>> p = pl.scatter(x1-x2,y1-y2,c=colors,edgecolors='none',cmap=pl.cm.spectral,label='_nolegend_')
112 >>> p = pl.plot(x1,y1,'b-',label='Primary')
113 >>> p = pl.plot(x2,y2,'g-',label='Secondary')
114 >>> p = pl.plot(x1-x2,y1-y2,'r-',label='Secondary relative')
115 >>> p = pl.colorbar()
116 >>> p = pl.legend(loc='lower right')
117 >>> p = pl.grid()
118 >>> p = pl.xlabel(r'Angular size DEC (arcsec) --> N')
119 >>> p = pl.ylabel(r'Angular size RA (arcsec) --> E')
120
121 ]include figure]]ivs_binary_keplerorbit_capella.png]
122
123 """
124 import logging
125 import numpy as np
126 from scipy import optimize
127 from ivs.units import conversions
128 from ivs.units.constants import *
129 from ivs.coordinates import vectors
130
131 logger = logging.getLogger('IVS.KEPLER')
132
133
134
136 """
137 Evaluate Radial velocities due to kepler orbit.
138
139 These parameters define the Keplerian orbit if you give times points (C{times}, days):
140 1. Period of the system (days)
141 2. time of periastron passage T0 (not x0!) (HJD)
142 3. eccentricity
143 4. omega, the longitude of periastron gives the orientation
144 of the orbit within its own plane (radians)
145 5. the semiamplitude of the velocity curve (km/s)
146 6. systemic velocity RV0 (RV of centre of mass of system) (km/s)
147
148 These parameters define the Keplerian orbit if you give angles (C{theta}, radians):
149 1. Period of the system (days)
150 2. eccentricity
151 3. a, the semi-major axis of the absolute orbit of the component (au)
152 4. omega, the longitude of periastron gives the orientation of the
153 orbit within in its own plane (radians)
154 5. inclination of the orbit (radians)
155 6. systemic velocity RV0 (RV of centre of mass of system) (km/s)
156
157 The periastron passage T0 can be derived via x0 by calculating
158
159 T0 = x0/(2pi*Freq) + times[0]
160
161 See e.g. p 41,42 of Hilditch, 'An Introduction To Close Binary Stars'
162
163 @parameter parameters: parameters of Keplerian orbit (dependent on input)
164 @type parameters: iterable
165 @parameter times: observation times (days)
166 @type times: numpy array
167 @parameter theta: orbital angles (radians)
168 @type theta: numpy array
169 @param itermax: number of iterations to find true anomaly
170 @type itermax: integer
171 @return: fitted radial velocities (km/s)
172 @rtype: ndarray
173 """
174
175
176 if times is not None:
177
178 P,T0,e,omega,K,RV0 = parameters
179 freq = 1./P
180 x0 = T0*2*np.pi*freq
181
182 E,true_an = true_anomaly(times*2*np.pi*freq-x0,e,itermax=itermax)
183
184 RVfit = RV0 + K*(e*np.cos(omega) + np.cos(true_an+omega))
185
186 elif theta is not None:
187 P,e,a,omega,i,RV0 = parameters
188 P = conversions.convert('d','s',P)
189 a = conversions.convert('au','m',a)
190 K = 2*np.pi*a*np.sin(i)/ (P*np.sqrt(1.-e**2))
191 RVfit = RV0 + K*(e*np.cos(omega) + np.cos(theta+omega))/1000.
192
193 return RVfit
194
195 -def orbit_in_plane(times,parameters,component='primary',coordinate_frame='polar'):
196 """
197 Construct an orbit in the orbital plane.
198
199 Give times in days
200
201 Parameters contains:
202 1. period (days)
203 2. eccentricity
204 3. semi-major axis (au)
205 4. time of periastron passage T0 (not x0!) (HJD)
206
207 Return r (m) and theta (radians)
208
209 @param times: times of observations (days)
210 @type times: array
211 @param parameters: list of parameters (P,e,a,T0)
212 @type parameters: list
213 @param component: component to calculate the orbit of. If it's the secondary,
214 the angles will be shifted by 180 degrees
215 @type component: string, one of 'primary' or 'secondary'
216 @param coordinate_frame: type of coordinates
217 @type coordinate_frame: str, ('polar' or 'cartesian')
218 @return: coord1,coord2
219 @rtype: 2xarray
220 """
221 P,e,a,T0 = parameters
222 P = conversions.convert('d','s',P)
223 a = conversions.convert('au','m',a)
224 T0 = conversions.convert('d','s',T0)
225 times = conversions.convert('d','s',times)
226
227 n = 2*np.pi/P
228 ma = n*(times-T0)
229 E,theta = true_anomaly(ma,e)
230 r = a*(1-e*np.cos(E))
231 PR = r*np.sin(theta)
232
233
234
235
236 if 'sec' in component.lower():
237 theta += np.pi
238
239 if coordinate_frame=='polar':
240 return r,theta
241 elif coordinate_frame=='cartesian':
242 return r*np.cos(theta),r*np.sin(theta)
243
244 -def velocity_in_plane(times,parameters,component='primary',coordinate_frame='polar'):
245 """
246 Calculate the velocity in the orbital plane.
247
248 @param times: times of observations (days)
249 @type times: array
250 @param parameters: list of parameters (P,e,a,T0)
251 @type parameters: list
252 @param component: component to calculate the orbit of. If it's the secondary,
253 the angles will be shifted by 180 degrees
254 @type component: string, one of 'primary' or 'secondary'
255 @param coordinate_frame: type of coordinates
256 @type coordinate_frame: str, ('polar' or 'cartesian')
257 @return: coord1,coord2
258 @rtype: 2xarray
259 """
260
261 r,theta = orbit_in_plane(times,parameters,component=component,coordinate_frame='polar')
262 P,e,a,T0 = parameters
263 P = conversions.convert('d','s',P)
264 a = conversions.convert('au','m',a)
265
266
267 l = r*(1+e*np.cos(theta))
268 L = 2*np.pi*a**2/P*np.sqrt(1-e**2)
269 rdot = L/l*e*np.sin(theta)
270 thetadot = L/r**2
271
272
273 if coordinate_frame=='polar':
274 return rdot,thetadot
275 elif coordinate_frame=='cartesian':
276 vx,vy,vz = vectors.spher2cart((r,theta,np.pi/2.),(rdot,r*thetadot,0))
277 return vx,vy
278
279
281 """
282 Project an orbit onto the plane of the sky.
283
284 Parameters contains the Euler angles:
285 1. omega: the longitude of periastron gives the orientation of the
286 orbit within in its own plane (radians)
287 2. Omega: PA of ascending node (radians)
288 3. i: inclination (radians), i=pi/2 is edge on
289
290 Returns x,y (orbit in plane of the sky) and z
291
292 See Hilditch p41 for a sketch of the coordinates. The difference with this
293 is that all angles are inverted. In this approach, North is in the positive
294 X direction, East is in the negative Y direction.
295 """
296 omega,Omega,i = parameters
297 x = r*(np.cos(Omega)*np.cos(theta+omega) - np.sin(+Omega)*np.sin(theta+omega)*np.cos(i))
298 y = r*(np.sin(Omega)*np.cos(theta+omega) + np.cos(+Omega)*np.sin(theta+omega)*np.cos(i))
299 z = r*(np.sin(theta+omega)*np.sin(i))
300 return x,y,z
301
302
303
304 -def orbit_on_sky(times,parameters,distance=None,component='primary'):
305 """
306 Construct an orbit projected on the sky.
307
308 Parameters contains:
309 1. period (days)
310 2. eccentricity
311 3. semi-major axis (au)
312 4. time of periastron passage T0 (not x0!) (HJD)
313 5. omega: the longitude of periastron gives the orientation of the
314 orbit within in its own plane (radians)
315 6. Omega: PA of ascending node (radians)
316 7. i: inclination (radians), i=pi/2 is edge on
317
318 You can give an extra parameter 'distance' as a tuple (value,'unit'). This
319 will be used to convert the distances to angular scale (arcsec).
320
321 Else, this function returns the distances in AU.
322
323 See Hilditch p41 for a sketch of the coordinates. The difference with this
324 is that all angles are inverted. In this approach, North is in the positive
325 X direction, East is in the negative Y direction.
326 """
327
328 pars_in_plane,euler_angles = parameters[:4],parameters[4:]
329
330 r,theta = orbit_in_plane(times,pars_in_plane,component=component)
331
332 x,y,z = project_orbit(r,theta,euler_angles)
333
334
335 if distance is not None:
336 d = conversions.convert(distance[1],'m',distance[0])
337 x = 2*np.arctan(x/(2*d))/np.pi*180*3600
338 y = 2*np.arctan(y/(2*d))/np.pi*180*3600
339 z = 2*np.arctan(z/(2*d))/np.pi*180*3600
340 return x,y,z
341 else:
342 return x/au,y/au,z/au
343
344
345
346
348 """
349 Calculation of true and eccentric anomaly in Kepler orbits.
350
351 M is the phase of the star, e is the eccentricity
352
353 See p.39 of Hilditch, 'An Introduction To Close Binary Stars'
354
355 @parameter M: phase
356 @type M: float
357 @parameter e: eccentricity
358 @type e: float
359 @keyword itermax: maximum number of iterations
360 @type itermax: integer
361 @return: eccentric anomaly (E), true anomaly (theta)
362 @rtype: float,float
363 """
364
365 Fn = M + e*np.sin(M) + e**2/2.*np.sin(2*M)
366
367 for i in range(itermax):
368 F = Fn
369 Mn = F-e*np.sin(F)
370 Fn = F+(M-Mn)/(1.-e*np.cos(F))
371 keep = F!=0
372 if hasattr(F,'__iter__'):
373 if np.all(abs((Fn-F)[keep]/F[keep])<0.00001):
374 break
375 elif (abs((Fn-F)/F)<0.00001):
376 break
377
378
379 true_an = 2.*np.arctan(np.sqrt((1.+e)/(1.-e))*np.tan(Fn/2.))
380 return Fn,true_an
381
382
383
384
386 """
387 Compute orbital phase from true anomaly T
388
389 @parameter T: true anomaly
390 @type T: float
391 @parameter omega: argument of periastron (radians)
392 @type omega: float
393 @parameter e: eccentricity
394 @type e: float
395 @parameter pshift: phase shift
396 @type pshift: float
397 @return: phase of superior conjunction, phase of periastron passage
398 @rtype: float,float
399 """
400 E = 2.0*np.arctan(np.sqrt((1-e)/(1+e)) * np.tan(T/2.0))
401 M = E - e*np.sin(E)
402 return (M+omega)/(2.0*np.pi) - 0.25 + pshift
403
404
405
407 """
408 Compute phase of superior conjunction and periastron passage.
409
410 Example usage:
411 >>> omega = np.pi/4.0
412 >>> e = 0.3
413 >>> print calculate_critical_phases(omega,e)
414 (-0.125, -0.057644612788576133, -0.42054512757020118, -0.19235538721142384, 0.17054512757020118)
415
416 @parameter omega: argument of periastron (radians)
417 @type omega: float
418 @parameter e: eccentricity
419 @type e: float
420 @parameter pshift: phase shift
421 @type pshift: float
422 @return: phase of superior conjunction, phase of periastron passage
423 @rtype: float,float
424 """
425
426 Phi_omega = (omega - np.pi/2.0)/(2.0*np.pi) + pshift
427
428 Phi_conj = calculate_phase(np.pi/2.0-omega,e,omega,pshift)
429 Phi_inf = calculate_phase(3.0*np.pi/2.0-omega,e,omega,pshift)
430 Phi_asc = calculate_phase(-omega,e,omega,pshift)
431 Phi_desc = calculate_phase(np.pi-omega,e,omega,pshift)
432 return Phi_omega,Phi_conj,Phi_inf,Phi_asc,Phi_desc
433
434
436 """
437 Calculate the eclipse separation between primary and secondary in a light curve.
438
439 Minimum separation at omega=pi
440 Maximum spearation at omega=0
441
442 @parameter e: eccentricity
443 @type e: float
444 @parameter omega: argument of periastron (radians)
445 @type omega: float
446 @return: separation in phase units (0.5 is half)
447 @rtype: float
448 """
449 radians = np.pi+2*np.arctan(e*np.cos(omega)/np.sqrt(1-e**2)) + 2*e*np.cos(omega)*np.sqrt(1-e**2)/(1-e**2*np.sin(omega)**2)
450 return radians/(2*np.pi)
451
452
454 """
455 Caculate the argument of periastron from the eclipse separation and eccentricity.
456
457 separation in phase units.
458
459 @parameter separation: separation in phase units (0.5 is half)
460 @type separation: float
461 @parameter e: eccentricity
462 @type e: float
463 @return: omega (longitude of periastron) in radians
464 @rtype: float
465 """
466 minsep_omega = np.pi
467 maxsep_omega = 0.
468 minsep = eclipse_separation(e,minsep_omega)
469 maxsep = eclipse_separation(e,maxsep_omega)
470 if separation<minsep or maxsep<separation:
471 logger.warning('Phase separation must be between %.3g and %.3g when e=%.3g'%(minsep,maxsep,e))
472 return np.nan
473 else:
474 omega = optimize.bisect(lambda x:separation-eclipse_separation(e,x),maxsep_omega,minsep_omega)
475 return omega
476
477
478
479
480
482 """
483 Kepler's third law.
484
485 Give two quantities, derived the third.
486
487 M = total mass system (solar units)
488 a = semi-major axis (au)
489 P = period (d)
490
491 >>> print third_law(M=1.,a=1.)
492 365.256891359
493 >>> print third_law(a=1.,P=365.25)
494 1.00003773538
495 >>> print third_law(M=1.,P=365.25)
496 0.999987421856
497 """
498 if a is not None:
499 a *= au
500 if P is not None:
501 P *= (24*3600.)
502 if M is not None:
503 M *= Msol
504
505 if M is None:
506 return 4*np.pi**2*a**3/P**2/GG/Msol
507 if a is None:
508 return (GG*M*P**2/(4*np.pi**2))**(1./3.)/au
509 if P is None:
510 return np.sqrt(4*np.pi**2*a**3/(GG*M))/(24*3600.)
511
512
513 if __name__=="__main__":
514 import doctest
515 import pylab as pl
516 doctest.testmod()
517 pl.show()
518