Module rotation

source code

Roche models of fast and differential rotation

Section 1. Non-rotating spherically symmetric single star

The most trivial usage of this module would be to construct a non-rotating, spherically symmetric single star. All quantities should be constant over the stellar surface, except of course the projected ones. We compute this via the zero-rotation limit of the fast rotating model:

First set some parameters: the rotation rate, effective temperature at the pole, polar radius, mass and viewing angle.

>>> omega = 0.00       # ratio to critical velocity
>>> T_pole = 5500.     # K
>>> r_pole = 1.        # solar radii
>>> M = 1.             # solar mass
>>> view_angle = pi/2  # radians

Construct a coordinate grid with a resolution of 50 points, covering the whole star. Unravel, so that we have access to 1D coordinate arrays

>>> theta,phi = local.get_grid(50,full=True)
>>> thetas,phis = np.ravel(theta),np.ravel(phi)

Calculate the shape of the stellar surface and the local gravity for every point in the grid. This is a 3D vector in Cartesian coordinates, so reshape all components to match the grid shape. As a reference, also explicitly calculate the polar surface gravity, which is the z-component of the gravity vector.

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

Calculate the value of the surface gravity in SI-units:

>>> grav = vectors.norm(grav_local)

Compute the size of all surface elements, and the angle with respect to the radius vector (should be all zero). The normal to the surface is the negative of the surface gravity vector, which is directed inwards:

>>> areas_local,cos_gamma = local.surface_elements((radius,(theta,phi)),-grav_local)

Now we can calculate the local effective temperature (assuming radiative atmosphere in Von Zeipel law), and bolometric intensities:

>>> teff_local = local.temperature(vectors.norm(grav_local),g_pole,T_pole,beta=1.)
>>> ints_local = local.intensity(teff_local,grav,photband='OPEN.BOL')

Define the line-of-sight vector:

>>> line_of_sight = np.zeros_like(grav_local)
>>> line_of_sight[0] = -sin(view_angle)
>>> line_of_sight[2] = -cos(view_angle)

Compute the angles between surface normals and line-of-sight. We do this in 1D:

>>> gravity = np.array([igrav.ravel() for igrav in grav_local])
>>> line_of_sight = np.array([ilos.ravel() for ilos in line_of_sight])
>>> angles = vectors.angle(-gravity,line_of_sight)
>>> mus = cos(angles)

Now compute the bolometric projected intensity, taking care of limbdarkening effects:

>>> intens_proj = local.intensity(teff_local.ravel(),grav.ravel(),mus,photband='OPEN.BOL').reshape(theta.shape)

The total and projected luminosity can then be calculated the following way:

>>> print pi*(ints_local*areas_local*constants.Rsol_cgs**2).sum()/constants.Lsol_cgs
0.992471247895
>>> print pi*np.nansum(intens_proj*areas_local*constants.Rsol_cgs**2)/constants.Lsol_cgs
0.360380373413

Now make some plots showing the local quantities:

>>> quantities = areas_local,np.log10(grav*100),teff_local,ints_local,intens_proj,angles/pi*180,radius
>>> names = 'Area','log g', 'Teff', 'Flux', 'Proj. flux', 'Angle'
>>> p = pl.figure()
>>> rows,cols = 2,3    
>>> 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_nonrotstar.png]

Section 2. Fast and uniformly rotating star

Repeating the same example as above, but now with parameters

>>> omega = 0.98

Gives the figure below:

]include figure]]ivs_binary_rochepotential_fastrotstar.png]

Changing the viewing angle to 80 degrees

>>> view_angle = 80/180.*pi

gives: Note that the maximum projected intensity is higher than in the previous example: this is because we view directly at higher latitudes now, which are hotter than the equator. The least flux is coming from the equatorial zone.

]include figure]]ivs_binary_rochepotential_fastrotstar_80incl.png]

To generate an image of the star, it is wise to convert the coordinates to Cartesian coordinates:

>>> x,y,z = vectors.spher2cart_coord(radius.ravel(),phis,thetas)
>>> p = pl.figure()
>>> p = pl.subplot(121,aspect='equal')
>>> p = pl.title('top view')
>>> p = pl.scatter(x,y,c=teff_local.ravel(),edgecolors='none')
>>> p = pl.subplot(122,aspect='equal')
>>> p = pl.title('side view')
>>> p = pl.scatter(y,z,c=teff_local.ravel(),edgecolors='none')

The movie below is obtained by calculating the projected intensity in the line of sight

>>> proj_intens,mus = local.projected_intensity(teff_local,grav_local,areas_local,np.array([0,-sin(view_angle),-cos(view_angle)]),photband='OPEN.BOL')

for different lines of sight, and then for each line of sight computing the Cartesian coordinates, and finally rotating in Y-Z plane.

>>> y,z = vectors.rotate(y,z,-(pi/2-view_angle))

]include figure]]ivs_binary_rochepotential_fastrotstar_shape.png]

]include figure]]ivs_binary_rochepotential_fastrot.gif]

Section 3. Differentially rotating star

The differential rotation rate has to follow the stereotypical law

Omega(theta) = Omega_pole + b * sin(theta)

The critical velocity can now easily be exceeded at the equator, if the rotation rate increases towards the pole.

First set some parameters: the rotation rate of the equator and pole, effective temperature at the pole, polar radius, mass and viewing angle.

>>> omega_eq = 0.1       # ratio to critical rotation rate (equatorial)
>>> omega_pl = 0.9       # ratio to critical rotation rate (polar)
>>> T_pole = 5500.       # K
>>> r_pole = 1.0         # solar radii
>>> M = 1.               # solar mass
>>> view_angle = pi/2    # radians

We now do very similar stuff as in Section 1, except for the different Roche potential. (We can skip making the grid now)

>>> radius = np.array([get_diffrot_roche_radius(itheta,r_pole,M,omega_eq,omega_pl) for itheta in thetas]).reshape(theta.shape)
>>> grav_local = np.array([diffrot_roche_surface_gravity(iradius,itheta,iphi,r_pole,M,omega_eq,omega_pl) 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 = diffrot_roche_surface_gravity(r_pole,0,0,r_pole,M,omega_eq,omega_pl)[-1]

Calculate the value of the surface gravity in SI-units:

>>> grav = vectors.norm(grav_local)

Compute the size of all surface elements, and the angle with respect to the radius vector (should be all zero). The normal to the surface is the negative of the surface gravity vector, which is directed inwards:

>>> areas_local,cos_gamma = local.surface_elements((radius,(theta,phi)),-grav_local)

Now we can calculate the local effective temperature (assuming radiative atmosphere in Von Zeipel law), and bolometric intensities:

>>> teff_local = local.temperature(vectors.norm(grav_local),g_pole,T_pole,beta=1.)
>>> ints_local = local.intensity(teff_local,grav,photband='OPEN.BOL')

Compute the angles between surface normals and line-of-sight. We do this in 1D:

>>> gravity = np.array([igrav.ravel() for igrav in grav_local])
>>> line_of_sight = np.array([ilos.ravel() for ilos in line_of_sight])
>>> angles = vectors.angle(-gravity,line_of_sight)
>>> mus = cos(angles)

Now compute the bolometric projected intensity, taking care of limbdarkening effects:

>>> intens_proj = local.intensity(teff_local.ravel(),grav.ravel(),mus,photband='OPEN.BOL').reshape(theta.shape)

Now make some plots showing the local quantities:

>>> quantities = areas_local,np.log10(grav*100),teff_local,ints_local,intens_proj,angles/pi*180,radius
>>> names = 'Area','log g', 'Teff', 'Flux', 'Proj. flux', 'Angle'
>>> p = pl.figure()
>>> rows,cols = 2,3    
>>> 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_diffrotstar.png]

To generate an image of the star, it is wise to convert the coordinates to Cartesian coordinates:

>>> x,y,z = vectors.spher2cart_coord(radius.ravel(),phis,thetas)
>>> p = pl.figure()
>>> p = pl.subplot(121,aspect='equal')
>>> p = pl.title('top view')
>>> p = pl.scatter(x,y,c=teff_local.ravel(),edgecolors='none')
>>> p = pl.subplot(122,aspect='equal')
>>> p = pl.title('side view')
>>> p = pl.scatter(y,z,c=teff_local.ravel(),edgecolors='none')

]include figure]]ivs_binary_rochepotential_diffrotstar_shape.png]

>>> Omega = critical_angular_velocity(1.,1.)
>>> P = 2*pi/Omega
>>> P = conversions.convert('s','d',P)
>>> print 'Critical rotation period of the Sun: %.3f days'%(P)
Critical rotation period of the Sun: 0.213 days

>>> critical_velocity(1.,1.,definition=1)
356.71131858379499
>>> critical_velocity(1.,1.,definition=2)
356.71131858379488