Module binary

source code

Roche models of binary or multiple stars.

You'd better not use this module, but talk to Pieter Degroote.

The premisse of the Roche potential is that the stellar mass can be represented by a point source.

Currently, three types of stellar shapes are implemented:

1. asynchronously rotating eccentric binary (get_binary_roche_radius)
2. fast rotating star (get_fastrot_roche_radius)
3. differentially (fast) rotating star (get_diffrot_roche_surface_gravity)

This module can be used to calculate the following information:

1. the distorted shape of the star due to a Roche potential
2. the local surface gravity assuming Von Zeipel gravity darkening
3. the local effective temperature
4. the local velocity vector due to rotation
5. the radial velocity
6. the total distorted surface area
7. the total distorted luminosity
8. the projected intensity in the line of sight
9. the mean radial velocity (e.g. Rossiter effect in binaries)
10. a synthetic spectrum starting from a library of spectra
11. a synthetic light curve from the projected intensities

Section 1. Binary star

As an example we take the binary star SX Aurigae. This is the minimum set of parameters.

>>> P = 1.2100802        # period in days
>>> e = 0.               # eccentricity
>>> q = 0.54369          # mass ratio M2/M1
>>> asini = 11.9*constants.Rsol/constants.au # semi-major axis
>>> i = 81.27            # inclination angle
>>> F = 1.               # synchronicity parameter
>>> d = 1.               # separation in semi-major axis units
>>> Phi1 = 3.            # gravitational potential primary
>>> Phi2 = 5.05          # gravitational potential secondary
>>> T_pole1 = 25000.     # polar temperature primary
>>> T_pole2 = 18850.     # polar temperature secondary
>>> M1 = 10.3            # primary mass
>>> M2 = 5.6             # secondary mass

Note that we actually do not need the masses, since we can derive them from the third kepler law with the other parameters. We do so for convenience.

Everything is calculated in units of semi-major axis. Thus, to convert to SI units we need the following conversion factor:

>>> a = asini * sin(i/180.*pi)
>>> to_SI = a*constants.au

We need some constants for the calculations: the polar radii and angular velocity

>>> r_pole1 = newton(binary_roche_potential,1e-5,args=(0,0,Phi1,  q,d,F))
>>> r_pole2 = newton(binary_roche_potential,1e-5,args=(0,0,Phi2,1/q,d,F))
>>> omega_rot = F * 2*pi/(P*24*3600) * 1/d**2 * sqrt( (1+e)*(1-e))
>>> omega_rot_vec = np.array([0.,0.,-omega_rot])

Derive the shape of the two stars

>>> 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)
>>> 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)

We focus on the primary, then repeat everything for the secondary: The local surface gravity can only be calculated if we have Cartesian coordinates.

>>> x1,y1,z1 = vectors.spher2cart_coord(radius1,phi,theta)
>>> 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)
>>> Gamma_pole1 = binary_roche_potential_gradient(0,0,r_pole1,q,d,F,norm=True)
>>> zeta1 = g_pole1 / Gamma_pole1
>>> dOmega1 = binary_roche_potential_gradient(x1,y1,z1,q,d,F,norm=False)
>>> grav_local1 = dOmega1*zeta1
>>> grav_local1 = np.array([i.reshape(theta.shape) for i in grav_local1])
>>> grav1 = vectors.norm(grav_local1)

Now we can, as before, calculate the other local quantities:

>>> areas_local1,cos_gamma1 = surface_elements((radius1,theta,phi),-grav_local1)
>>> teff_local1 = local_temperature(grav1,g_pole1,T_pole1,beta=1.)
>>> ints_local1 = local_intensity(teff_local1,grav1,np.ones_like(cos_gamma1),photband='OPEN.BOL')
>>> velo_local1 = np.cross(np.array([x1,y1,z1]).T*to_SI,omega_rot_vec).T

Now make some plots showing the local quantities:

>>> quantities = areas_local1,np.log10(grav1*100),teff_local1,ints_local1
>>> names = 'Area','log g', 'Teff', 'Flux'
>>> p = pl.figure()
>>> rows,cols = 2,2    
>>> for i,(quantity,name) in enumerate(zip(quantities,names)):
...    p = pl.subplot(rows,cols,i+1)
...    p = pl.title(name)
...    q = quantity.ravel()
...    vmin,vmax = q[-np.isnan(q)].min(),q[-np.isnan(q)].max()
...    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()
...    p = pl.scatter(phis/pi*180,thetas/pi*180,c=q,edgecolors='none',vmin=vmin,vmax=vmax)
...    p = pl.colorbar()
...    p = pl.xlim(0,360)
...    p = pl.ylim(180,0)

]include figure]]ivs_binary_rochepotential_binarystar.png]

>>> p = pl.figure()
>>> p = pl.subplot(121,aspect='equal')
>>> p = pl.title('top view')
>>> p = pl.scatter(x1,y1,c=teff_local1.ravel(),edgecolors='none')
>>> p = pl.subplot(122,aspect='equal')
>>> p = pl.title('side view')
>>> p = pl.scatter(x1,z1,c=teff_local1.ravel(),edgecolors='none')

]include figure]]ivs_binary_rochepotential_binarystar_shape.png]

Section 5. Synthesizing spectra

In this case study, we calculate the projected radial velocity of a uniformly rotating star, but in the limit of no deformation due to the rotation. We do this, so that we can easily compare rotational broadening calculated with the ROTIN package, and numerically. So let us build a star with these parameters:

>>> omega = 0.0        # ratio to critical velocity
>>> T_pole = 13000.    # K
>>> r_pole = 2.0       # solar radii
>>> M = 1.5            # solar mass
>>> view_angle = pi/2  # radians
>>> theta,phi = get_grid(20,100,full=True,gtype='spher')
>>> thetas,phis = np.ravel(theta),np.ravel(phi)

Then calculate the shape of this star

>>> radius = np.array([get_fastrot_roche_radius(itheta,r_pole,omega) for itheta in thetas]).reshape(theta.shape)
>>> 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
>>> grav_local = np.array([i.reshape(theta.shape) for i in grav_local])
>>> g_pole = fastrot_roche_surface_gravity(r_pole,0,0,r_pole,omega,M)[-1]
>>> grav = vectors.norm(grav_local)

and the local quantities

>>> areas_local,cos_gamma = surface_elements((radius,theta,phi),-grav_local)
>>> teff_local = local_temperature(vectors.norm(grav_local),g_pole,T_pole,beta=1.)
>>> ints_local = local_intensity(teff_local,grav,photband='OPEN.BOL')
>>> x,y,z = vectors.spher2cart_coord(radius.ravel(),phis,thetas)

Assume, with a shape of a non-rotating star, that we have a velocity component on the surface of the star:

>>> myomega = 0.5
>>> velo_local = diffrot_velocity((phi,theta,radius*constants.Rsol),myomega,myomega,r_pole,M)

Collect all the necessary information in one record array.

>>> starlist = [x,y,z] + [i.ravel() for i in velo_local] + [i.ravel() for i in grav_local] + [teff_local.ravel(),areas_local.ravel()]
>>> starnames = ['x','y','z','vx','vy','vz','gravx','gravy','gravz','teff','areas']
>>> star = np.rec.array(starlist,names=starnames)

Project the star in some line-of-sight. The velocity component in the X-direction is the radial velocity.

>>> view_angle = pi/2 # edge on
>>> mystar = project(star,view_long=(0,0,0),view_lat=(view_angle,0,0),photband='OPEN.BOL',only_visible=True,plot_sort=True)

We can calculate the synthetic spectra for all surface elements between 7055 and 7075 angstrom:

>>> loggs = np.log10(np.sqrt(mystar['gravx']**2 + mystar['gravy']**2 + mystar['gravz']**2)*100)
>>> iterator = zip(mystar['teff'],loggs,mystar['vx']/1000.)
>>> wave_spec = np.linspace(7055,7075,750)
>>> spectra = np.array([spectra_model.get_table(teff=iteff,logg=ilogg,vrad=ivrad,wave=wave_spec)[1] for iteff,ilogg,ivrad in iterator])

The total observed spectrum is then simply the weighted sum with the local projected intensities:

>>> average_spectrum = np.average(spectra,weights=mystar['projflux'],axis=0)

We compare with the ROTIN package, which assumes a linear limbdarkening law

>>> original = spectra_model.get_table(teff=mystar['teff'][0],logg=loggs[0],wave=wave_spec)[1]
>>> 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]

And a plot can be made via

>>> colors = mystar['eyeflux']/mystar['eyeflux'].max()
>>> p = pl.figure(figsize=(7,7))
>>> ax = pl.axes([0,0.5,1,0.5])
>>> p = ax.set_aspect('equal')
>>> p = ax.set_axis_bgcolor('k')
>>> p = pl.box(on=False)
>>> p = pl.xticks([]);p = pl.yticks([])
>>> p = pl.scatter(mystar['y'],mystar['z'],c=colors,edgecolors='none',cmap=pl.cm.gray,vmin=0,vmax=1)
>>> 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)
>>> ax = pl.axes([0.13,0.1,0.78,0.4])
>>> p = pl.plot(wave_spec,average_spectrum,'k-',lw=2,label='Numerical')
>>> p = pl.plot(wave_spec,original,'r-',label='No rotation')
>>> p = pl.plot(wave_spec,rotin1,'b-',label='Rotin v=%.1f km/s'%(vmax/1000.),lw=2)
>>> p = pl.ylim(0.95,1.0)
>>> p = pl.legend(loc='lower right',prop=dict(size='small'))
>>> p = pl.xlabel('Wavelength [angstrom]',color='1')
>>> p = pl.ylabel('Normalised flux',color='1')

]include figure]]ivs_binary_rochepotential_norot_spectrum_b.png]

]include figure]]ivs_binary_rochepotential_norot_spectrum_a.png]